From 11fb209ee6ebdcad909e0bc04cd9bc49639faf26 Mon Sep 17 00:00:00 2001
From: Francesco Lannutti <nicolati@users.sf.net>
Date: Sat, 28 Mar 2020 16:41:04 +0100
Subject: [PATCH] First KLU support of CIDER ONED simulations

---
 src/ciderlib/oned/oneadmit.c      | 231 +++++++++++----
 src/ciderlib/oned/onecond.c       |  21 +-
 src/ciderlib/oned/onecont.c       | 397 +++++++++++++++++++++++--
 src/ciderlib/oned/onedest.c       |  20 +-
 src/ciderlib/oned/onedext.h       |   4 +
 src/ciderlib/oned/onepoiss.c      | 137 ++++++++-
 src/ciderlib/oned/oneproj.c       |  21 +-
 src/ciderlib/oned/onesolve.c      | 203 +++++++++++--
 src/include/ngspice/klu-binding.h |  35 +++
 src/include/ngspice/klu.h         |  27 ++
 src/include/ngspice/onemesh.h     |  45 +++
 src/include/ngspice/smpdefs.h     |  26 ++
 src/maths/KLU/klu_extract.c       |  12 +-
 src/maths/KLU/klusmp.c            | 462 +++++++++++++++++++++++++++++-
 src/maths/sparse/spsmp.c          |  37 +++
 15 files changed, 1555 insertions(+), 123 deletions(-)

diff --git a/src/ciderlib/oned/oneadmit.c b/src/ciderlib/oned/oneadmit.c
index 13df2f602..2e0b607fa 100644
--- a/src/ciderlib/oned/oneadmit.c
+++ b/src/ciderlib/oned/oneadmit.c
@@ -116,30 +116,54 @@ NUMDadmittance(ONEdevice *pDevice, double omega, SPcomplex *yd)
       pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
     }
     ONE_jacLoad(pDevice);
-    spSetComplex(pDevice->matrix);
-    for (index = 1; index < pDevice->numNodes; index++) {
-      pElem = pDevice->elemArray[index];
-      if (pElem->elemType == SEMICON) {
-	for (i = 0; i <= 1; i++) {
-	  pNode = pElem->pNodes[i];
-	  if (pNode->nodeType != CONTACT) {
-	    spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
-	    spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
-	  }
-	}
+
+#ifdef KLU
+    if (pDevice->matrix->CKTkluMODE) {
+      pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
+    } else {
+#endif
+
+      spSetComplex (pDevice->matrix->SPmatrix) ;
+
+      for (index = 1; index < pDevice->numNodes; index++) {
+        pElem = pDevice->elemArray[index];
+        if (pElem->elemType == SEMICON) {
+          for (i = 0; i <= 1; i++) {
+            pNode = pElem->pNodes[i];
+	    if (pNode->nodeType != CONTACT) {
+	      spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
+	      spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
+            }
+          }
+        }
       }
+
+#ifdef KLU
     }
+#endif
+
     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
 
     /* FACTOR */
     startTime = SPfrontEnd->IFseconds();
-    spFactor(pDevice->matrix);
+
+#ifdef KLU
+    SMPluFacKLUforCIDER (pDevice->matrix) ;
+#else
+    SMPcLUfac(pDevice->matrix, 0);
+#endif
+
     pDevice->pStats->factorTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
 
     /* SOLVE */
     startTime = SPfrontEnd->IFseconds();
-    spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
-	pDevice->rhsImag, solutionImag);
+
+#ifdef KLU
+    SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#else
+    SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#endif
+
     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
   }
   /* MISC */
@@ -288,30 +312,54 @@ NBJTadmittance(ONEdevice *pDevice, double omega, SPcomplex *yIeVce,
       pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
       pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
     }
-    spSetComplex(pDevice->matrix);
-    for (index = 1; index < pDevice->numNodes; index++) {
-      pElem = pDevice->elemArray[index];
-      if (pElem->elemType == SEMICON) {
-	for (i = 0; i <= 1; i++) {
-	  pNode = pElem->pNodes[i];
-	  if (pNode->nodeType != CONTACT) {
-	    spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
-	    spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
-	  }
-	}
+
+#ifdef KLU
+    if (pDevice->matrix->CKTkluMODE) {
+      pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
+    } else {
+#endif
+
+      spSetComplex (pDevice->matrix->SPmatrix) ;
+
+      for (index = 1; index < pDevice->numNodes; index++) {
+        pElem = pDevice->elemArray[index];
+        if (pElem->elemType == SEMICON) {
+          for (i = 0; i <= 1; i++) {
+            pNode = pElem->pNodes[i];
+            if (pNode->nodeType != CONTACT) {
+              spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
+              spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
+            }
+          }
+        }
       }
+
+#ifdef KLU
     }
+#endif
+
     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
 
     /* FACTOR */
     startTime = SPfrontEnd->IFseconds();
-    spFactor(pDevice->matrix);
+
+#ifdef KLU
+    SMPluFacKLUforCIDER (pDevice->matrix) ;
+#else
+    SMPcLUfac(pDevice->matrix, 0);
+#endif
+
     pDevice->pStats->factorTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
 
     /* SOLVE */
     startTime = SPfrontEnd->IFseconds();
-    spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
-	pDevice->rhsImag, solutionImag);
+
+#ifdef KLU
+    SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#else
+    SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#endif
+
     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
 
     /* MISC */
@@ -345,8 +393,13 @@ NBJTadmittance(ONEdevice *pDevice, double omega, SPcomplex *yIeVce,
 
     /* SOLVE */
     startTime = SPfrontEnd->IFseconds();
-    spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
-	pDevice->rhsImag, solutionImag);
+
+#ifdef KLU
+    SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#else
+    SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#endif
+
     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
   }
   /* MISC */
@@ -418,7 +471,11 @@ ONEsorSolve(ONEdevice *pDevice, double *xReal, double *xImag, double omega)
       rhsSOR[index] += pDevice->rhs[index];
     }
     /* compute xReal(k+1). solution stored in rhsSOR */
-    spSolve(pDevice->matrix, rhsSOR, rhsSOR, NULL, NULL);
+#ifdef KLU
+    SMPsolveKLUforCIDER (pDevice->matrix, rhsSOR, rhsSOR, NULL, NULL) ;
+#else
+    SMPsolve(pDevice->matrix, rhsSOR, rhsSOR);
+#endif
 
     /* modify solution when wRelax is not 1 */
     if (wRelax != 1) {
@@ -454,8 +511,12 @@ ONEsorSolve(ONEdevice *pDevice, double *xReal, double *xImag, double omega)
       }
     }
     /* compute xImag(k+1) */
-    spSolve(pDevice->matrix, rhsSOR, rhsSOR,
-	NULL, NULL);
+#ifdef KLU
+    SMPsolveKLUforCIDER (pDevice->matrix, rhsSOR, rhsSOR, NULL, NULL) ;
+#else
+    SMPsolve(pDevice->matrix, rhsSOR, rhsSOR);
+#endif
+
     /* modify solution when wRelax is not 1 */
     if (wRelax != 1) {
       for (index = 1; index <= numEqns; index++) {
@@ -516,25 +577,44 @@ NUMDys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yd)
     pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
     pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
   }
-  spSetComplex(pDevice->matrix);
-  for (index = 1; index < pDevice->numNodes; index++) {
-    pElem = pDevice->elemArray[index];
-    if (pElem->elemType == SEMICON) {
-      for (i = 0; i <= 1; i++) {
-	pNode = pElem->pNodes[i];
-	if (pNode->nodeType != CONTACT) {
-	  CMPLX_MULT_SCALAR(temp, cOmega, 0.5 * pElem->dx);
-	  spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
-	  spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
-	}
+
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
+  } else {
+#endif
+
+    spSetComplex (pDevice->matrix->SPmatrix) ;
+
+    for (index = 1; index < pDevice->numNodes; index++) {
+      pElem = pDevice->elemArray[index];
+      if (pElem->elemType == SEMICON) {
+        for (i = 0; i <= 1; i++) {
+          pNode = pElem->pNodes[i];
+          if (pNode->nodeType != CONTACT) {
+            CMPLX_MULT_SCALAR(temp, cOmega, 0.5 * pElem->dx);
+            spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
+            spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
+          }
+        }
       }
     }
+
+#ifdef KLU
   }
+#endif
 
+#ifdef KLU
+  SMPluFacKLUforCIDER (pDevice->matrix) ;
+#else
+  SMPcLUfac (pDevice->matrix, 0) ;
+#endif
 
-  spFactor(pDevice->matrix);
-  spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
-      pDevice->rhsImag, solutionImag);
+#ifdef KLU
+  SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#else
+  SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#endif
 
   pElem = pDevice->elemArray[1];
   pNode = pElem->pLeftNode;
@@ -584,24 +664,45 @@ NBJTys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yIeVce,
     pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
     pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
   }
-  spSetComplex(pDevice->matrix);
-  for (index = 1; index < pDevice->numNodes; index++) {
-    pElem = pDevice->elemArray[index];
-    if (pElem->elemType == SEMICON) {
-      for (i = 0; i <= 1; i++) {
-	pNode = pElem->pNodes[i];
-	if (pNode->nodeType != CONTACT) {
-	  CMPLX_MULT_SCALAR(temp, cOmega, 0.5 * pElem->dx);
-	  spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
-	  spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
-	}
+
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
+  } else {
+#endif
+
+    spSetComplex (pDevice->matrix->SPmatrix) ;
+
+    for (index = 1; index < pDevice->numNodes; index++) {
+      pElem = pDevice->elemArray[index];
+      if (pElem->elemType == SEMICON) {
+        for (i = 0; i <= 1; i++) {
+          pNode = pElem->pNodes[i];
+          if (pNode->nodeType != CONTACT) {
+            CMPLX_MULT_SCALAR(temp, cOmega, 0.5 * pElem->dx);
+            spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
+            spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
+          }
+        }
       }
     }
+
+#ifdef KLU
   }
+#endif
+
+#ifdef KLU
+  SMPluFacKLUforCIDER (pDevice->matrix) ;
+#else
+  SMPcLUfac (pDevice->matrix, 0) ;
+#endif
+
+#ifdef KLU
+  SMPsolveKLUforCIDER(pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag);
+#else
+  SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#endif
 
-  spFactor(pDevice->matrix);
-  spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
-      pDevice->rhsImag, solutionImag);
   pElem = pDevice->elemArray[1];
   pNode = pElem->pLeftNode;
   y = computeAdmittance(pNode, FALSE, solutionReal, solutionImag, &cOmega);
@@ -624,8 +725,12 @@ NBJTys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yIeVce,
   }
 
   /* don't need to LU factor the jacobian since it exists */
-  spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
-      pDevice->rhsImag, solutionImag);
+#ifdef KLU
+    SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#else
+    SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
+#endif
+
   pElem = pDevice->elemArray[1];
   pNode = pElem->pLeftNode;
   y = computeAdmittance(pNode, FALSE, solutionReal, solutionImag, &cOmega);
diff --git a/src/ciderlib/oned/onecond.c b/src/ciderlib/oned/onecond.c
index e70078c50..bdf984f4a 100644
--- a/src/ciderlib/oned/onecond.c
+++ b/src/ciderlib/oned/onecond.c
@@ -42,7 +42,12 @@ NUMDconductance(ONEdevice *pDevice, BOOLEAN tranAnalysis,
     pDevice->rhs[pNode->pEqn] = -pEdge->dJpDpsiP1;
   }
   incVpn = pDevice->dcDeltaSolution;
-  spSolve(pDevice->matrix, pDevice->rhs, incVpn, NULL, NULL);
+
+#ifdef KLU
+  SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVpn, NULL, NULL) ;
+#else
+  SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVpn) ;
+#endif
 
   pElem = pDevice->elemArray[1];
   pNode = pElem->pRightNode;
@@ -95,7 +100,12 @@ NBJTconductance(ONEdevice *pDevice, BOOLEAN tranAnalysis, double *intCoeff,
     pDevice->rhs[pNode->pEqn] = -pEdge->dJpDpsiP1;
   }
   incVce = pDevice->dcDeltaSolution;
-  spSolve(pDevice->matrix, pDevice->rhs, incVce, NULL, NULL);
+
+#ifdef KLU
+  SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVce, NULL, NULL) ;
+#else
+  SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVce) ;
+#endif
 
   /* zero the rhs before loading in the new rhs base contribution */
   for (index = 1; index <= pDevice->numEqns; index++) {
@@ -114,7 +124,12 @@ NBJTconductance(ONEdevice *pDevice, BOOLEAN tranAnalysis, double *intCoeff,
   }
 
   incVbe = pDevice->copiedSolution;
-  spSolve(pDevice->matrix, pDevice->rhs, incVbe, NULL, NULL);
+
+#ifdef KLU
+  SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVbe, NULL, NULL) ;
+#else
+  SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVbe) ;
+#endif
 
   pElem = pDevice->elemArray[1];/* first element */
   pEdge = pElem->pEdge;
diff --git a/src/ciderlib/oned/onecont.c b/src/ciderlib/oned/onecont.c
index 4bcffdc33..705dec16e 100644
--- a/src/ciderlib/oned/onecont.c
+++ b/src/ciderlib/oned/onecont.c
@@ -15,6 +15,10 @@ Author:	1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
 #include "ngspice/cidersupt.h"
 #include "../../maths/misc/bernoull.h"
 
+#ifdef KLU
+#include "ngspice/klu-binding.h"
+#endif
+
 
 /* functions to setup and solve the continuity equations */
 /* Both continuity equations are solved */
@@ -39,21 +43,76 @@ ONE_jacBuild(ONEdevice *pDevice)
       pNode = pElem->pNodes[index];
       /* get poisson pointer */
       psiEqn = pNode->psiEqn;
-      pNode->fPsiPsi = spGetElement(matrix, psiEqn, psiEqn);
+
+#ifdef KLU
+      pNode->fPsiPsi = SMPmakeEltKLUforCIDER (matrix, psiEqn, psiEqn) ;
+      pNode->fPsiPsiBinding = NULL ;
+#else
+      pNode->fPsiPsi = SMPmakeElt(matrix, psiEqn, psiEqn);
+#endif
 
       if (pElem->elemType == SEMICON) {
 	/* get continuity-coupling terms */
 	nEqn = pNode->nEqn;
 	pEqn = pNode->pEqn;
 	/* pointers for additional terms */
-	pNode->fPsiN = spGetElement(matrix, psiEqn, nEqn);
-	pNode->fPsiP = spGetElement(matrix, psiEqn, pEqn);
-	pNode->fNPsi = spGetElement(matrix, nEqn, psiEqn);
-	pNode->fNN = spGetElement(matrix, nEqn, nEqn);
-	pNode->fNP = spGetElement(matrix, nEqn, pEqn);
-	pNode->fPPsi = spGetElement(matrix, pEqn, psiEqn);
-	pNode->fPP = spGetElement(matrix, pEqn, pEqn);
-	pNode->fPN = spGetElement(matrix, pEqn, nEqn);
+
+#ifdef KLU
+        pNode->fPsiN = SMPmakeEltKLUforCIDER (matrix, psiEqn, nEqn) ;
+        pNode->fPsiNBinding = NULL ;
+#else
+	pNode->fPsiN = SMPmakeElt(matrix, psiEqn, nEqn);
+#endif
+
+#ifdef KLU
+        pNode->fPsiP = SMPmakeEltKLUforCIDER (matrix, psiEqn, pEqn) ;
+        pNode->fPsiPBinding = NULL ;
+#else
+	pNode->fPsiP = SMPmakeElt(matrix, psiEqn, pEqn);
+#endif
+
+#ifdef KLU
+        pNode->fNPsi = SMPmakeEltKLUforCIDER (matrix, nEqn, psiEqn) ;
+        pNode->fNPsiBinding = NULL ;
+#else
+	pNode->fNPsi = SMPmakeElt(matrix, nEqn, psiEqn);
+#endif
+
+#ifdef KLU
+        pNode->fNN = SMPmakeEltKLUforCIDER (matrix, nEqn, nEqn) ;
+        pNode->fNNBinding = NULL ;
+#else
+	pNode->fNN = SMPmakeElt(matrix, nEqn, nEqn);
+#endif
+
+#ifdef KLU
+        pNode->fNP = SMPmakeEltKLUforCIDER (matrix, nEqn, pEqn) ;
+        pNode->fNPBinding = NULL ;
+#else
+	pNode->fNP = SMPmakeElt(matrix, nEqn, pEqn);
+#endif
+
+#ifdef KLU
+        pNode->fPPsi = SMPmakeEltKLUforCIDER (matrix, pEqn, psiEqn) ;
+        pNode->fPPsiBinding = NULL ;
+#else
+	pNode->fPPsi = SMPmakeElt(matrix, pEqn, psiEqn);
+#endif
+
+#ifdef KLU
+        pNode->fPP = SMPmakeEltKLUforCIDER (matrix, pEqn, pEqn) ;
+        pNode->fPPBinding = NULL ;
+#else
+	pNode->fPP = SMPmakeElt(matrix, pEqn, pEqn);
+#endif
+
+#ifdef KLU
+        pNode->fPN = SMPmakeEltKLUforCIDER (matrix, pEqn, nEqn) ;
+        pNode->fPNBinding = NULL ;
+#else
+	pNode->fPN = SMPmakeElt(matrix, pEqn, nEqn);
+#endif
+
       } else {
 	nEqn = 0;
 	pEqn = 0;
@@ -73,34 +132,302 @@ ONE_jacBuild(ONEdevice *pDevice)
 
     /* now terms to couple to adjacent nodes */
     pNode = pElem->pLeftNode;
-    pNode->fPsiPsiiP1 = spGetElement(matrix, psiEqnL, psiEqnR);
+
+#ifdef KLU
+    pNode->fPsiPsiiP1 = SMPmakeEltKLUforCIDER (matrix, psiEqnL, psiEqnR) ;
+    pNode->fPsiPsiiP1Binding = NULL ;
+#else
+    pNode->fPsiPsiiP1 = SMPmakeElt(matrix, psiEqnL, psiEqnR);
+#endif
+
     if (pElem->elemType == SEMICON) {
       /* pointers for additional terms */
-      pNode->fNPsiiP1 = spGetElement(matrix, nEqnL, psiEqnR);
-      pNode->fNNiP1 = spGetElement(matrix, nEqnL, nEqnR);
-      pNode->fPPsiiP1 = spGetElement(matrix, pEqnL, psiEqnR);
-      pNode->fPPiP1 = spGetElement(matrix, pEqnL, pEqnR);
+
+#ifdef KLU
+      pNode->fNPsiiP1 = SMPmakeEltKLUforCIDER (matrix, nEqnL, psiEqnR) ;
+      pNode->fNPsiiP1Binding = NULL ;
+#else
+      pNode->fNPsiiP1 = SMPmakeElt(matrix, nEqnL, psiEqnR);
+#endif
+
+#ifdef KLU
+      pNode->fNNiP1 = SMPmakeEltKLUforCIDER (matrix, nEqnL, nEqnR) ;
+      pNode->fNNiP1Binding = NULL ;
+#else
+      pNode->fNNiP1 = SMPmakeElt(matrix, nEqnL, nEqnR);
+#endif
+
+#ifdef KLU
+      pNode->fPPsiiP1 = SMPmakeEltKLUforCIDER (matrix, pEqnL, psiEqnR) ;
+      pNode->fPPsiiP1Binding = NULL ;
+#else
+      pNode->fPPsiiP1 = SMPmakeElt(matrix, pEqnL, psiEqnR);
+#endif
+
+#ifdef KLU
+      pNode->fPPiP1 = SMPmakeEltKLUforCIDER (matrix, pEqnL, pEqnR) ;
+      pNode->fPPiP1Binding = NULL ;
+#else
+      pNode->fPPiP1 = SMPmakeElt(matrix, pEqnL, pEqnR);
+#endif
+
       if (AvalancheGen) {
-	pNode->fNPiP1 = spGetElement(matrix, nEqnL, pEqnR);
-	pNode->fPNiP1 = spGetElement(matrix, pEqnL, nEqnR);
+
+#ifdef KLU
+        pNode->fNPiP1 = SMPmakeEltKLUforCIDER (matrix, nEqnL, pEqnR) ;
+        pNode->fNPiP1Binding = NULL ;
+#else
+	pNode->fNPiP1 = SMPmakeElt(matrix, nEqnL, pEqnR);
+#endif
+
+#ifdef KLU
+        pNode->fPNiP1 = SMPmakeEltKLUforCIDER (matrix, pEqnL, nEqnR) ;
+        pNode->fPNiP1Binding = NULL ;
+#else
+	pNode->fPNiP1 = SMPmakeElt(matrix, pEqnL, nEqnR);
+#endif
+
       }
     }
     pNode = pElem->pRightNode;
-    pNode->fPsiPsiiM1 = spGetElement(matrix, psiEqnR, psiEqnL);
+
+#ifdef KLU
+    pNode->fPsiPsiiM1 = SMPmakeEltKLUforCIDER (matrix, psiEqnR, psiEqnL) ;
+    pNode->fPsiPsiiM1Binding = NULL ;
+#else
+    pNode->fPsiPsiiM1 = SMPmakeElt(matrix, psiEqnR, psiEqnL);
+#endif
+
     if (pElem->elemType == SEMICON) {
       /* pointers for additional terms */
-      pNode->fNPsiiM1 = spGetElement(matrix, nEqnR, psiEqnL);
-      pNode->fNNiM1 = spGetElement(matrix, nEqnR, nEqnL);
-      pNode->fPPsiiM1 = spGetElement(matrix, pEqnR, psiEqnL);
-      pNode->fPPiM1 = spGetElement(matrix, pEqnR, pEqnL);
+
+#ifdef KLU
+      pNode->fNPsiiM1 = SMPmakeEltKLUforCIDER (matrix, nEqnR, psiEqnL) ;
+      pNode->fNPsiiM1Binding = NULL ;
+#else
+      pNode->fNPsiiM1 = SMPmakeElt(matrix, nEqnR, psiEqnL);
+#endif
+
+#ifdef KLU
+      pNode->fNNiM1 = SMPmakeEltKLUforCIDER (matrix, nEqnR, nEqnL) ;
+      pNode->fNNiM1Binding = NULL ;
+#else
+      pNode->fNNiM1 = SMPmakeElt(matrix, nEqnR, nEqnL);
+#endif
+
+#ifdef KLU
+      pNode->fPPsiiM1 = SMPmakeEltKLUforCIDER (matrix, pEqnR, psiEqnL) ;
+      pNode->fPPsiiM1Binding = NULL ;
+#else
+      pNode->fPPsiiM1 = SMPmakeElt(matrix, pEqnR, psiEqnL);
+#endif
+
+#ifdef KLU
+      pNode->fPPiM1 = SMPmakeEltKLUforCIDER (matrix, pEqnR, pEqnL) ;
+      pNode->fPPiM1Binding = NULL ;
+#else
+      pNode->fPPiM1 = SMPmakeElt(matrix, pEqnR, pEqnL);
+#endif
+
       if (AvalancheGen) {
-	pNode->fNPiM1 = spGetElement(matrix, nEqnR, pEqnL);
-	pNode->fPNiM1 = spGetElement(matrix, pEqnR, nEqnL);
+
+#ifdef KLU
+        pNode->fNPiM1 = SMPmakeEltKLUforCIDER (matrix, nEqnR, pEqnL) ;
+        pNode->fNPiM1Binding = NULL ;
+#else
+	pNode->fNPiM1 = SMPmakeElt(matrix, nEqnR, pEqnL);
+#endif
+
+#ifdef KLU
+        pNode->fPNiM1 = SMPmakeEltKLUforCIDER (matrix, pEqnR, nEqnL) ;
+        pNode->fPNiM1Binding = NULL ;
+#else
+	pNode->fPNiM1 = SMPmakeElt(matrix, pEqnR, nEqnL);
+#endif
+
       }
     }
   }
 }
 
+#ifdef KLU
+/*
+#define CREATE_KLU_BINDING_TABLE_CIDER(ptr, binding, a, b)                             \
+    printf ("Swapping Pointer %s: (%d,%d)\n", #ptr, a, b) ; \
+    if ((a > 0) && (b > 0)) {                                      \
+        if (pNode->binding != NULL) { \
+            if (pNode->binding->CSC_Complex != NULL) { \
+                printf ("  Looking for the Pointer: %p\n", pNode->binding->CSC_Complex) ; \
+                qsort (BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
+                i.COO = NULL ; \
+                i.CSC_Complex = pNode->binding->CSC_Complex ; \
+                matched = (BindKluElementCOO *) bsearch (&i, BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
+                if (matched != NULL) { \
+                    printf ("  Found the Old Pointer\n") ; \
+                    pNode->ptr = pNode->binding->CSC_Complex ; \
+                } else { \
+                    i.COO = pNode->ptr ;                                                          \
+                    i.CSC_Complex = NULL ; \
+                    matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+                    if (matched != NULL) { \
+                        printf ("  Looking for the Pointer 1\n") ; \
+                        pNode->binding = matched ;                                                \
+                        pNode->ptr = matched->CSC_Complex ;                                               \
+                    } else { \
+                        printf ("  Leaving the Pointer as is\n") ; \
+                    } \
+                } \
+            } else { \
+                printf ("  Looking for the Pointer 2\n") ; \
+                i.COO = pNode->ptr ;                                                          \
+                i.CSC_Complex = NULL ; \
+                matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+                pNode->binding = matched ;                                                \
+                pNode->ptr = matched->CSC_Complex ;                                               \
+            } \
+        } else { \
+            printf ("  Looking for the Pointer 3\n") ; \
+            i.COO = pNode->ptr ;                                                          \
+            i.CSC_Complex = NULL ; \
+            matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+            pNode->binding = matched ;                                                \
+            pNode->ptr = matched->CSC_Complex ;                                               \
+        } \
+    }
+*/
+
+/*
+#define CREATE_KLU_BINDING_TABLE_CIDER_TO_REAL(ptr, binding, a, b)                             \
+    if ((a > 0) && (b > 0)) {                                      \
+        printf ("Macro\n") ; \
+        if (pNode->binding) { \
+            printf ("IF: %p\n", pNode->binding) ; \
+            printf ("COO: %p\n", pNode->binding->COO) ; \
+            printf ("CSC: %p\n", pNode->binding->CSC) ; \
+            if (pNode->binding->CSC_Complex) { \
+                printf ("CSC_Complex: %p\n", pNode->binding->CSC_Complex) ; \
+                pNode->ptr = pNode->binding->CSC_Complex ; \
+            } else { \
+                i = pNode->ptr ;                                                          \
+                matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+                pNode->binding = matched ;                                                \
+                pNode->ptr = matched->CSC_Complex ;                                               \
+            } \
+        } else { \
+            i = pNode->ptr ;                                                          \
+            matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+            pNode->binding = matched ;                                                \
+            pNode->ptr = matched->CSC_Complex ;                                               \
+        } \
+    }
+*/
+
+void
+ONEbindCSC (ONEdevice *pDevice)
+{
+  ONEelem *pElem;
+  ONEnode *pNode;
+  int index, eIndex;
+  int psiEqn, nEqn, pEqn;	/* scratch for deref'd eqn numbers */
+  int psiEqnL=0, nEqnL=0, pEqnL=0;
+  int psiEqnR=0, nEqnR=0, pEqnR=0;
+  BindKluElementCOO i, *matched, *BindStruct, *BindStructCSC ;
+  size_t nz ;
+
+  BindStruct = pDevice->matrix->SMPkluMatrix->KLUmatrixBindStructCOO ;
+  nz = pDevice->matrix->SMPkluMatrix->KLUmatrixNZ ;
+
+  BindStructCSC = (BindKluElementCOO *) malloc (nz * sizeof(BindKluElementCOO)) ;
+  for (index = 0 ; index < (int)nz ; index++) {
+    BindStructCSC [index] = BindStruct [index] ;
+  }
+
+  for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
+    pElem = pDevice->elemArray[eIndex];
+    /* first the self terms */
+    for (index = 0; index <= 1; index++) {
+      pNode = pElem->pNodes[index];
+      /* get poisson pointer */
+      psiEqn = pNode->psiEqn;
+
+      CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsi, fPsiPsiBinding, psiEqn, psiEqn) ;
+
+      if (pElem->elemType == SEMICON) {
+	/* get continuity-coupling terms */
+	nEqn = pNode->nEqn;
+	pEqn = pNode->pEqn;
+	/* pointers for additional terms */
+
+        CREATE_KLU_BINDING_TABLE_CIDER(fPsiN, fPsiNBinding, psiEqn, nEqn) ;
+        CREATE_KLU_BINDING_TABLE_CIDER(fPsiP, fPsiPBinding, psiEqn, pEqn) ;
+        CREATE_KLU_BINDING_TABLE_CIDER(fNPsi, fNPsiBinding, nEqn, psiEqn) ;
+        CREATE_KLU_BINDING_TABLE_CIDER(fNN, fNNBinding, nEqn, nEqn) ;
+        CREATE_KLU_BINDING_TABLE_CIDER(fNP, fNPBinding, nEqn, pEqn) ;
+        CREATE_KLU_BINDING_TABLE_CIDER(fPPsi, fPPsiBinding, pEqn, psiEqn) ;
+        CREATE_KLU_BINDING_TABLE_CIDER(fPP, fPPBinding, pEqn, pEqn) ;
+        CREATE_KLU_BINDING_TABLE_CIDER(fPN, fPNBinding, pEqn, nEqn) ;
+
+      } else {
+	nEqn = 0;
+	pEqn = 0;
+      }
+
+      /* save indices */
+      if (index == 0) {		/* left node */
+	psiEqnL = psiEqn;
+	nEqnL = nEqn;
+	pEqnL = pEqn;
+      } else {
+	psiEqnR = psiEqn;
+	nEqnR = nEqn;
+	pEqnR = pEqn;
+      }
+    }
+
+    /* now terms to couple to adjacent nodes */
+    pNode = pElem->pLeftNode;
+
+    CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsiiP1, fPsiPsiiP1Binding, psiEqnL, psiEqnR) ;
+
+    if (pElem->elemType == SEMICON) {
+      /* pointers for additional terms */
+
+      CREATE_KLU_BINDING_TABLE_CIDER(fNPsiiP1, fNPsiiP1Binding, nEqnL, psiEqnR) ;
+      CREATE_KLU_BINDING_TABLE_CIDER(fNNiP1, fNNiP1Binding, nEqnL, nEqnR) ;
+      CREATE_KLU_BINDING_TABLE_CIDER(fPPsiiP1, fPPsiiP1Binding, pEqnL, psiEqnR) ;
+      CREATE_KLU_BINDING_TABLE_CIDER(fPPiP1, fPPiP1Binding, pEqnL, pEqnR) ;
+
+      if (AvalancheGen) {
+
+        CREATE_KLU_BINDING_TABLE_CIDER(fNPiP1, fNPiP1Binding, nEqnL, pEqnR) ;
+        CREATE_KLU_BINDING_TABLE_CIDER(fPNiP1, fPNiP1Binding, pEqnL, nEqnR) ;
+
+      }
+    }
+    pNode = pElem->pRightNode;
+
+    CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsiiM1, fPsiPsiiM1Binding, psiEqnR, psiEqnL) ;
+
+    if (pElem->elemType == SEMICON) {
+      /* pointers for additional terms */
+
+      CREATE_KLU_BINDING_TABLE_CIDER(fNPsiiM1, fNPsiiM1Binding, nEqnR, psiEqnL) ;
+      CREATE_KLU_BINDING_TABLE_CIDER(fNNiM1, fNNiM1Binding, nEqnR, nEqnL) ;
+      CREATE_KLU_BINDING_TABLE_CIDER(fPPsiiM1, fPPsiiM1Binding, pEqnR, psiEqnL) ;
+      CREATE_KLU_BINDING_TABLE_CIDER(fPPiM1, fPPiM1Binding, pEqnR, pEqnL) ;
+
+      if (AvalancheGen) {
+
+        CREATE_KLU_BINDING_TABLE_CIDER(fNPiM1, fNPiM1Binding, nEqnR, pEqnL) ;
+        CREATE_KLU_BINDING_TABLE_CIDER(fPNiM1, fPNiM1Binding, pEqnR, nEqnL) ;
+
+      }
+    }
+  }
+
+  free (BindStructCSC) ;
+}
+#endif
 
 void 
 ONE_sysLoad(ONEdevice *pDevice, BOOLEAN tranAnalysis, 
@@ -131,7 +458,17 @@ ONE_sysLoad(ONEdevice *pDevice, BOOLEAN tranAnalysis,
   }
 
   /* zero the matrix */
-  spClear(pDevice->matrix);
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    SMPclearKLUforCIDER (pDevice->matrix) ;
+  } else {
+#endif
+
+    SMPclear(pDevice->matrix);
+
+#ifdef KLU
+  }
+#endif
 
   for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
     pElem = pDevice->elemArray[eIndex];
@@ -275,7 +612,17 @@ ONE_jacLoad(ONEdevice *pDevice)
   ONE_commonTerms(pDevice, FALSE, FALSE, NULL);
 
   /* zero the matrix */
-  spClear(pDevice->matrix);
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    SMPclearKLUforCIDER (pDevice->matrix) ;
+  } else {
+#endif
+
+    SMPclear(pDevice->matrix);
+
+#ifdef KLU
+  }
+#endif
 
   for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
     pElem = pDevice->elemArray[eIndex];
diff --git a/src/ciderlib/oned/onedest.c b/src/ciderlib/oned/onedest.c
index 56a544ed2..a21d8b824 100644
--- a/src/ciderlib/oned/onedest.c
+++ b/src/ciderlib/oned/onedest.c
@@ -35,7 +35,15 @@ ONEdestroy(ONEdevice *pDevice)
     FREE(pDevice->copiedSolution);
     FREE(pDevice->rhs);
     FREE(pDevice->rhsImag);
-    spDestroy(pDevice->matrix);
+
+#ifdef KLU
+    SMPdestroyKLUforCIDER (pDevice->matrix) ;
+#else
+    SMPdestroy (pDevice->matrix) ;
+#endif
+
+    FREE (pDevice->matrix) ;
+
     break;
   case SLV_EQUIL:
     /* free up the vectors allocated in the equilibrium solution */
@@ -43,7 +51,15 @@ ONEdestroy(ONEdevice *pDevice)
     FREE(pDevice->dcDeltaSolution);
     FREE(pDevice->copiedSolution);
     FREE(pDevice->rhs);
-    spDestroy(pDevice->matrix);
+
+#ifdef KLU
+    SMPdestroyKLUforCIDER (pDevice->matrix) ;
+#else
+    SMPdestroy (pDevice->matrix) ;
+#endif
+
+    FREE (pDevice->matrix) ;
+
     break;
   case SLV_NONE:
     break;
diff --git a/src/ciderlib/oned/onedext.h b/src/ciderlib/oned/onedext.h
index 8548e7b36..2edbf0e82 100644
--- a/src/ciderlib/oned/onedext.h
+++ b/src/ciderlib/oned/onedext.h
@@ -115,5 +115,9 @@ extern void ONEjacCheck(ONEdevice *, BOOLEAN, ONEtranInfo *);
 extern void ONEpredict(ONEdevice *, ONEtranInfo *);
 extern BOOLEAN ONEdeviceConverged(ONEdevice *);
 
+#ifdef KLU
+void ONEbindCSC (ONEdevice *) ;
+void ONEQbindCSC (ONEdevice *) ;
+#endif
 
 #endif
diff --git a/src/ciderlib/oned/onepoiss.c b/src/ciderlib/oned/onepoiss.c
index 0d993e006..30552e1d1 100644
--- a/src/ciderlib/oned/onepoiss.c
+++ b/src/ciderlib/oned/onepoiss.c
@@ -12,6 +12,10 @@ Author:	1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
 #include "oneddefs.h"
 #include "ngspice/spmatrix.h"
 
+#ifdef KLU
+#include "ngspice/klu-binding.h"
+#endif
+
 /* Functions to setup and solve the 1D poisson equation. */
 
 
@@ -27,16 +31,129 @@ ONEQjacBuild(ONEdevice *pDevice)
   for (index = 1; index < pDevice->numNodes; index++) {
     pElem = pDevice->elemArray[index];
     pNode = pElem->pLeftNode;
-    pNode->fPsiPsi = spGetElement(matrix, pNode->poiEqn, pNode->poiEqn);
+
+#ifdef KLU
+    pNode->fPsiPsi = SMPmakeEltKLUforCIDER (matrix, pNode->poiEqn, pNode->poiEqn) ;
+    pNode->fPsiPsiBinding = NULL ;
+#else
+    pNode->fPsiPsi = SMPmakeElt(matrix, pNode->poiEqn, pNode->poiEqn);
+#endif
+
     pNode1 = pElem->pRightNode;
-    pNode->fPsiPsiiP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn);
+
+#ifdef KLU
+    pNode->fPsiPsiiP1 = SMPmakeEltKLUforCIDER (matrix, pNode->poiEqn, pNode1->poiEqn) ;
+    pNode->fPsiPsiiP1Binding = NULL ;
+#else
+    pNode->fPsiPsiiP1 = SMPmakeElt(matrix, pNode->poiEqn, pNode1->poiEqn);
+#endif
+
     pNode = pElem->pRightNode;
-    pNode->fPsiPsi = spGetElement(matrix, pNode->poiEqn, pNode->poiEqn);
+
+#ifdef KLU
+    pNode->fPsiPsi = SMPmakeEltKLUforCIDER (matrix, pNode->poiEqn, pNode->poiEqn) ;
+    pNode->fPsiPsiBinding = NULL ;
+#else
+    pNode->fPsiPsi = SMPmakeElt(matrix, pNode->poiEqn, pNode->poiEqn);
+#endif
+
     pNode1 = pElem->pLeftNode;
-    pNode->fPsiPsiiM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn);
+
+#ifdef KLU
+    pNode->fPsiPsiiM1 = SMPmakeEltKLUforCIDER (matrix, pNode->poiEqn, pNode1->poiEqn) ;
+    pNode->fPsiPsiiM1Binding = NULL ;
+#else
+    pNode->fPsiPsiiM1 = SMPmakeElt(matrix, pNode->poiEqn, pNode1->poiEqn);
+#endif
+
   }
 }
 
+#ifdef KLU
+/*
+#define CREATE_KLU_BINDING_TABLE_CIDER(ptr, binding, a, b)                             \
+    printf ("Swapping Pointer %s: (%d,%d)\n", #ptr, a, b) ; \
+    if ((a > 0) && (b > 0)) {                                      \
+        if (pNode->binding != NULL) { \
+            if (pNode->binding->CSC_Complex != NULL) { \
+                printf ("  Looking for the Pointer: %p\n", pNode->binding->CSC_Complex) ; \
+                qsort (BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
+                i.COO = NULL ; \
+                i.CSC_Complex = pNode->binding->CSC_Complex ; \
+                matched = (BindKluElementCOO *) bsearch (&i, BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
+                if (matched != NULL) { \
+                    printf ("  Found the Old Pointer\n") ; \
+                    pNode->ptr = pNode->binding->CSC_Complex ; \
+                } else { \
+                    i.COO = pNode->ptr ;                                                          \
+                    i.CSC_Complex = NULL ; \
+                    matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+                    if (matched != NULL) { \
+                        printf ("  Looking for the Pointer 1\n") ; \
+                        pNode->binding = matched ;                                                \
+                        pNode->ptr = matched->CSC_Complex ;                                               \
+                    } else { \
+                        printf ("  Leaving the Pointer as is\n") ; \
+                    } \
+                } \
+            } else { \
+                printf ("  Looking for the Pointer 2\n") ; \
+                i.COO = pNode->ptr ;                                                          \
+                i.CSC_Complex = NULL ; \
+                matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+                pNode->binding = matched ;                                                \
+                pNode->ptr = matched->CSC_Complex ;                                               \
+            } \
+        } else { \
+            printf ("  Looking for the Pointer 3\n") ; \
+            i.COO = pNode->ptr ;                                                          \
+            i.CSC_Complex = NULL ; \
+            matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+            pNode->binding = matched ;                                                \
+            pNode->ptr = matched->CSC_Complex ;                                               \
+        } \
+    }
+*/
+
+void
+ONEQbindCSC (ONEdevice *pDevice)
+{
+  ONEelem *pElem ;
+  ONEnode *pNode, *pNode1 ;
+  int index ;
+  BindKluElementCOO i, *matched, *BindStruct, *BindStructCSC ;
+  size_t nz ;
+
+  BindStruct = pDevice->matrix->SMPkluMatrix->KLUmatrixBindStructCOO ;
+  nz = pDevice->matrix->SMPkluMatrix->KLUmatrixNZ ;
+
+  BindStructCSC = (BindKluElementCOO *) malloc (nz * sizeof(BindKluElementCOO)) ;
+  for (index = 0 ; index < (int)nz ; index++) {
+    BindStructCSC [index] = BindStruct [index] ;
+  }
+
+  for (index = 1 ; index < pDevice->numNodes ; index++) {
+    pElem = pDevice->elemArray [index] ;
+    pNode = pElem->pLeftNode ;
+
+    CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsi, fPsiPsiBinding, pNode->poiEqn, pNode->poiEqn) ;
+
+    pNode1 = pElem->pRightNode ;
+
+    CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsiiP1, fPsiPsiiP1Binding, pNode->poiEqn, pNode1->poiEqn) ;
+
+    pNode = pElem->pRightNode ;
+
+    CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsi, fPsiPsiBinding, pNode->poiEqn, pNode->poiEqn) ;
+
+    pNode1 = pElem->pLeftNode ;
+
+    CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsiiM1, fPsiPsiiM1Binding, pNode->poiEqn, pNode1->poiEqn) ;
+  }
+
+  free (BindStructCSC) ;
+}
+#endif
 
 void
 ONEQsysLoad(ONEdevice *pDevice)
@@ -58,7 +175,17 @@ ONEQsysLoad(ONEdevice *pDevice)
   }
 
   /* zero the matrix */
-  spClear(pDevice->matrix);
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    SMPclearKLUforCIDER (pDevice->matrix) ;
+  } else {
+#endif
+
+    SMPclear(pDevice->matrix);
+
+#ifdef KLU
+  }
+#endif
 
   for (index = 1; index < pDevice->numNodes; index++) {
     pElem = pDevice->elemArray[index];
diff --git a/src/ciderlib/oned/oneproj.c b/src/ciderlib/oned/oneproj.c
index c7c88502f..7e21d7bf1 100644
--- a/src/ciderlib/oned/oneproj.c
+++ b/src/ciderlib/oned/oneproj.c
@@ -52,7 +52,12 @@ NUMDproject(ONEdevice *pDevice, double delV)
     pDevice->rhs[pNode->pEqn] = -pEdge->dJpDpsiP1;
   }
   incVpn = pDevice->dcDeltaSolution;
-  spSolve(pDevice->matrix, pDevice->rhs, incVpn, NULL, NULL);
+
+#ifdef KLU
+  SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVpn, NULL, NULL) ;
+#else
+  SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVpn) ;
+#endif
 
   for (index = 1; index < pDevice->numNodes; index++) {
     pElem = pDevice->elemArray[index];
@@ -128,7 +133,12 @@ NBJTproject(ONEdevice *pDevice, double delVce, double delVbe,
       pDevice->rhs[pNode->pEqn] = -pEdge->dJpDpsiP1;
     }
     incVce = pDevice->dcDeltaSolution;
-    spSolve(pDevice->matrix, pDevice->rhs, incVce, NULL, NULL);
+
+#ifdef KLU
+    SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVce, NULL, NULL) ;
+#else
+    SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVce) ;
+#endif
 
     for (index = 1; index < pDevice->numNodes; index++) {
       pElem = pDevice->elemArray[index];
@@ -179,7 +189,12 @@ NBJTproject(ONEdevice *pDevice, double delVce, double delVbe,
     }
 
     incVbe = pDevice->copiedSolution;
-    spSolve(pDevice->matrix, pDevice->rhs, incVbe, NULL, NULL);
+
+#ifdef KLU
+    SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVbe, NULL, NULL) ;
+#else
+    SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVbe) ;
+#endif
 
     for (index = 1; index < pDevice->numNodes; index++) {
       pElem = pDevice->elemArray[index];
diff --git a/src/ciderlib/oned/onesolve.c b/src/ciderlib/oned/onesolve.c
index 7586100c2..e360d1044 100644
--- a/src/ciderlib/oned/onesolve.c
+++ b/src/ciderlib/oned/onesolve.c
@@ -91,7 +91,13 @@ ONEdcSolve(ONEdevice *pDevice, int iterationLimit, BOOLEAN newSolver,
     }
     /* FACTOR */
     startTime = SPfrontEnd->IFseconds();
-    error = spFactor(pDevice->matrix);
+
+#ifdef KLU
+    error = SMPreorderKLUforCIDER (pDevice->matrix) ;
+#else
+    error = SMPluFacForCIDER (pDevice->matrix) ;
+#endif
+
     factorTime += SPfrontEnd->IFseconds() - startTime;
     if (newSolver) {
       if (pDevice->iterationNumber == 1) {
@@ -110,7 +116,7 @@ ONEdcSolve(ONEdevice *pDevice, int iterationLimit, BOOLEAN newSolver,
     if (foundError(error)) {
       if (error == spSINGULAR) {
 	int badRow, badCol;
-	spWhereSingular(pDevice->matrix, &badRow, &badCol);
+	SMPgetError(pDevice->matrix, &badCol, &badRow);
 	printf("*****  singular at (%d,%d)\n", badRow, badCol);
       }
       /* Should probably try to recover now, but we'll punt instead. */
@@ -118,7 +124,13 @@ ONEdcSolve(ONEdevice *pDevice, int iterationLimit, BOOLEAN newSolver,
     }
     /* SOLVE */
     startTime = SPfrontEnd->IFseconds();
-    spSolve(pDevice->matrix, rhs, delta, NULL, NULL);
+
+#ifdef KLU
+    SMPsolveKLUforCIDER (pDevice->matrix, rhs, delta, NULL, NULL) ;
+#else
+    SMPsolveForCIDER (pDevice->matrix, rhs, delta) ;
+#endif
+
     solveTime += SPfrontEnd->IFseconds() - startTime;
 
     /* UPDATE */
@@ -441,7 +453,13 @@ ONEresetJacobian(ONEdevice *pDevice)
   
   
   ONE_jacLoad(pDevice);
-  error = spFactor(pDevice->matrix);
+
+#ifdef KLU
+  error = SMPreorderKLUforCIDER (pDevice->matrix) ;
+#else
+  error = SMPluFacForCIDER (pDevice->matrix) ;
+#endif
+
   if (foundError(error)) {
     exit(-1);
   }
@@ -539,7 +557,6 @@ ONEequilSolve(ONEdevice *pDevice)
   ONEnode *pNode;
   double startTime, setupTime, miscTime;
 
-
   setupTime = miscTime = 0.0;
 
   /* SETUP */
@@ -553,7 +570,15 @@ ONEequilSolve(ONEdevice *pDevice)
     FREE(pDevice->copiedSolution);
     FREE(pDevice->rhs);
     FREE(pDevice->rhsImag);
-    spDestroy(pDevice->matrix);
+
+#ifdef KLU
+    SMPdestroyKLUforCIDER (pDevice->matrix) ;
+#else
+    SMPdestroy (pDevice->matrix) ;
+#endif
+
+    FREE (pDevice->matrix) ;
+
     /* FALLTHRU */
   case SLV_NONE:
     pDevice->poissonOnly = TRUE;
@@ -562,15 +587,63 @@ ONEequilSolve(ONEdevice *pDevice)
     XCALLOC(pDevice->dcDeltaSolution, double, pDevice->dimEquil);
     XCALLOC(pDevice->copiedSolution, double, pDevice->dimEquil);
     XCALLOC(pDevice->rhs, double, pDevice->dimEquil);
-    pDevice->matrix = spCreate(pDevice->numEqns, 0, &error);
+
+    pDevice->matrix = TMALLOC (SMPmatrix, 1) ;
+
+#ifdef KLU
+    pDevice->matrix->CKTkluMODE = CKTkluON ; /* Francesco Lannutti - To be sustitued with a value coming from the uplevel */
+    error = SMPnewMatrixKLUforCIDER (pDevice->matrix, pDevice->numEqns, CKTkluMatrixReal) ;
+#else
+    error = SMPnewMatrixForCIDER (pDevice->matrix, pDevice->numEqns, 0) ;
+#endif
+
     if (error == spNO_MEMORY) {
       printf("ONEequilSolve: Out of Memory\n");
       exit(-1);
     }
     newSolver = TRUE;
-    spSetReal(pDevice->matrix);
+
+#ifdef KLU
+    if (pDevice->matrix->CKTkluMODE) {
+      pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixReal ;
+    } else {
+#endif
+
+      spSetReal (pDevice->matrix->SPmatrix) ;
+
+#ifdef KLU
+    }
+#endif
+
     ONEQjacBuild(pDevice);
-    pDevice->numOrigEquil = spElementCount(pDevice->matrix);
+
+#ifdef KLU
+    if (pDevice->matrix->CKTkluMODE) {
+
+      /* Convert the COO Storage to CSC for KLU and Fill the Binding Table */
+      SMPconvertCOOtoCSCKLUforCIDER (pDevice->matrix) ;
+
+      /* Bind the KLU Pointers */
+      ONEQbindCSC (pDevice) ;
+
+      /* Perform KLU Matrix Analysis */
+      pDevice->matrix->SMPkluMatrix->KLUmatrixSymbolic = klu_analyze ((int)pDevice->matrix->SMPkluMatrix->KLUmatrixN, pDevice->matrix->SMPkluMatrix->KLUmatrixAp,
+                                                                      pDevice->matrix->SMPkluMatrix->KLUmatrixAi, pDevice->matrix->SMPkluMatrix->KLUmatrixCommon) ;
+      if (pDevice->matrix->SMPkluMatrix->KLUmatrixSymbolic == NULL) {
+        printf ("CIDER: KLU Failed\n") ;
+        if (pDevice->matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) {
+          return ; // Francesco Lannutti - Fix KLU return values
+        }
+      }
+      printf ("CIDER: KLU to be fixed: spElementCount\n") ;
+      pDevice->numOrigEquil = 0 ; // Francesco Lannutti - Fix for KLU
+    } else {
+      pDevice->numOrigEquil = spElementCount (pDevice->matrix->SPmatrix) ;
+    }
+#else
+    pDevice->numOrigEquil = spElementCount (pDevice->matrix->SPmatrix) ;
+#endif
+
     pDevice->numFillEquil = 0;
     /* FALLTHRU */
   case SLV_EQUIL:
@@ -590,7 +663,21 @@ ONEequilSolve(ONEdevice *pDevice)
   /* MISCELLANEOUS */
   startTime = SPfrontEnd->IFseconds();
   if (newSolver) {
-    pDevice->numFillEquil = spFillinCount(pDevice->matrix);
+
+#ifdef KLU
+    if (pDevice->matrix->CKTkluMODE) {
+      // Francesco Lannutti - Fix for KLU
+      printf ("CIDER: KLU to be fixed: spFillinCount\n") ;
+      pDevice->numFillEquil = 0 ;
+    } else {
+#endif
+
+      pDevice->numFillEquil = spFillinCount(pDevice->matrix->SPmatrix);
+
+#ifdef KLU
+    }
+#endif
+
   }
   if (pDevice->converged) {
     ONEQcommonTerms(pDevice);
@@ -642,7 +729,15 @@ ONEbiasSolve(ONEdevice *pDevice, int iterationLimit,
     FREE(pDevice->dcDeltaSolution);
     FREE(pDevice->copiedSolution);
     FREE(pDevice->rhs);
-    spDestroy(pDevice->matrix);
+
+#ifdef KLU
+    SMPdestroyKLUforCIDER (pDevice->matrix) ;
+#else
+    SMPdestroy (pDevice->matrix) ;
+#endif
+
+    FREE (pDevice->matrix) ;
+
     /* FALLTHRU */
   case SLV_NONE:
     pDevice->poissonOnly = FALSE;
@@ -652,18 +747,65 @@ ONEbiasSolve(ONEdevice *pDevice, int iterationLimit,
     XCALLOC(pDevice->copiedSolution, double, pDevice->dimBias);
     XCALLOC(pDevice->rhs, double, pDevice->dimBias);
     XCALLOC(pDevice->rhsImag, double, pDevice->dimBias);
-    pDevice->matrix = spCreate(pDevice->numEqns, 1, &error);
+
+    pDevice->matrix = TMALLOC (SMPmatrix, 1) ;
+
+#ifdef KLU
+    pDevice->matrix->CKTkluMODE = CKTkluON ; /* Francesco Lannutti - To be sustitued with a value coming from the uplevel */
+    error = SMPnewMatrixKLUforCIDER (pDevice->matrix, pDevice->numEqns, CKTkluMatrixComplex) ;
+#else
+    error = SMPnewMatrixForCIDER (pDevice->matrix, pDevice->numEqns, 1) ;
+#endif
+
     if (error == spNO_MEMORY) {
       exit(-1);
     }
     newSolver = TRUE;
     ONE_jacBuild(pDevice);
-    pDevice->numOrigBias = spElementCount(pDevice->matrix);
+
+#ifdef KLU
+    if (pDevice->matrix->CKTkluMODE) {
+
+      /* Convert the COO Storage to CSC for KLU and Fill the Binding Table */
+      SMPconvertCOOtoCSCKLUforCIDER (pDevice->matrix) ;
+
+      /* Bind the KLU Pointers */
+      ONEbindCSC (pDevice) ;
+
+      /* Perform KLU Matrix Analysis */
+      pDevice->matrix->SMPkluMatrix->KLUmatrixSymbolic = klu_analyze ((int)pDevice->matrix->SMPkluMatrix->KLUmatrixN, pDevice->matrix->SMPkluMatrix->KLUmatrixAp,
+                                                                      pDevice->matrix->SMPkluMatrix->KLUmatrixAi, pDevice->matrix->SMPkluMatrix->KLUmatrixCommon) ;
+      if (pDevice->matrix->SMPkluMatrix->KLUmatrixSymbolic == NULL) {
+        if (pDevice->matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) {
+          printf ("CIDER: KLU failed\n") ;
+          return ; // Francesco Lannutti - Fix KLU return values
+        }
+      }
+      pDevice->numOrigBias = 0 ; // Francesco Lannutti - Fix for KLU
+    } else {
+      pDevice->numOrigBias = spElementCount(pDevice->matrix->SPmatrix);
+    }
+#else
+    pDevice->numOrigBias = spElementCount (pDevice->matrix->SPmatrix) ;
+#endif
+
     pDevice->numFillBias = 0;
     ONEstoreInitialGuess(pDevice);
     /* FALLTHRU */
   case SLV_SMSIG:
-    spSetReal(pDevice->matrix);
+
+#ifdef KLU
+    if (pDevice->matrix->CKTkluMODE) {
+      pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixReal ;
+    } else {
+#endif
+
+      spSetReal (pDevice->matrix->SPmatrix) ;
+
+#ifdef KLU
+    }
+#endif
+
     /* FALLTHRU */
   case SLV_BIAS:
     pDevice->solverType = SLV_BIAS;
@@ -681,7 +823,21 @@ ONEbiasSolve(ONEdevice *pDevice, int iterationLimit,
   /* MISCELLANEOUS */
   startTime = SPfrontEnd->IFseconds();
   if (newSolver) {
-    pDevice->numFillBias = spFillinCount(pDevice->matrix);
+
+#ifdef KLU
+    if (pDevice->matrix->CKTkluMODE) {
+      // Francesco Lannutti - Fix for KLU
+      printf ("CIDER: KLU to be fixed: spFillinCount\n") ;
+      pDevice->numFillBias = 0 ;
+    } else {
+#endif
+
+      pDevice->numFillBias = spFillinCount (pDevice->matrix->SPmatrix) ; // Francesco Lannutti - Fix for KLU
+
+#ifdef KLU
+    }
+#endif
+
   }
   solution = pDevice->dcSolution;
   if ((!pDevice->converged) && iterationLimit > 1) {
@@ -999,8 +1155,12 @@ double
 ONEnuNorm(ONEdevice *pDevice)
 {
   /* The LU Decomposed matrix is available.  Use it to calculate x. */
-  spSolve(pDevice->matrix, pDevice->rhs, pDevice->rhsImag,
-      NULL, NULL);
+#ifdef KLU
+  printf ("CIDER: KLU to be fixed ONEnuNorm\n") ;
+  SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, pDevice->rhsImag, NULL, NULL) ;
+#else
+  SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, pDevice->rhsImag) ;
+#endif
 
   /* Compute L2-norm of the solution vector (stored in rhsImag) */
   return (l2Norm(pDevice->rhsImag, pDevice->numEqns));
@@ -1041,7 +1201,14 @@ ONEjacCheck(ONEdevice *pDevice, BOOLEAN tranAnalysis, ONEtranInfo *info)
       pDevice->dcSolution[index] = pDevice->copiedSolution[index];
       for (rIndex = 1; rIndex <= pDevice->numEqns; rIndex++) {
 	diff = (pDevice->rhsImag[rIndex] - pDevice->rhs[rIndex]) / del;
-	dptr = spFindElement(pDevice->matrix, rIndex, index);
+
+#ifdef KLU
+        printf ("CIDER: KLU to be fixed: spFindElement") ;
+        dptr = NULL ;
+#else
+	dptr = spFindElement(pDevice->matrix->SPmatrix, rIndex, index); // Francesco Lannutti - Fix for KLU
+#endif
+
 	/*
 	 * if ( dptr ISNOT NULL ) { fprintf( stderr, "[%d][%d]: FD =
 	 * %11.4e, AJ = %11.4e\n", rIndex, index, diff, *dptr ); } else {
diff --git a/src/include/ngspice/klu-binding.h b/src/include/ngspice/klu-binding.h
index b9c53d330..5166399bc 100644
--- a/src/include/ngspice/klu-binding.h
+++ b/src/include/ngspice/klu-binding.h
@@ -99,4 +99,39 @@
         smp_data_out->input[k].port[l].ptr = smp_data_out->input[k].port[l].h.binding->CSC ;
 #endif
 
+#ifdef CIDER
+#define CREATE_KLU_BINDING_TABLE_CIDER(ptr, binding, a, b)                                \
+    if ((a > 0) && (b > 0)) {                                                             \
+        if (pNode->binding != NULL) {                                                     \
+            if (pNode->binding->CSC_Complex != NULL) {                                    \
+                qsort (BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
+                i.COO = NULL ; \
+                i.CSC_Complex = pNode->binding->CSC_Complex ; \
+                matched = (BindKluElementCOO *) bsearch (&i, BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
+                if (matched == NULL) { \
+                    i.COO = pNode->ptr ;                                                          \
+                    i.CSC_Complex = NULL ; \
+                    matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+                    if (matched != NULL) { \
+                        pNode->binding = matched ;                                                \
+                        pNode->ptr = matched->CSC_Complex ;                                               \
+                    } \
+                } \
+            } else { \
+                i.COO = pNode->ptr ;                                                          \
+                i.CSC_Complex = NULL ; \
+                matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+                pNode->binding = matched ;                                                \
+                pNode->ptr = matched->CSC_Complex ;                                               \
+            } \
+        } else { \
+            i.COO = pNode->ptr ;                                                          \
+            i.CSC_Complex = NULL ; \
+            matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
+            pNode->binding = matched ;                                                \
+            pNode->ptr = matched->CSC_Complex ;                                               \
+        } \
+    }
+#endif
+
 #endif
diff --git a/src/include/ngspice/klu.h b/src/include/ngspice/klu.h
index aedacb518..0a54106e2 100644
--- a/src/include/ngspice/klu.h
+++ b/src/include/ngspice/klu.h
@@ -935,6 +935,33 @@ int klu_z_convert_matrix_in_CSR
     klu_common *Common
 ) ;
 
+#ifdef CIDER
+typedef struct sBindKluElementCOO {
+    double *COO ;
+    double *CSC_Complex ;
+} BindKluElementCOO ;
+
+int BindKluCompareCOO (const void *a, const void *b) ;
+int BindKluCompareCSC (const void *a, const void *b) ;
+
+typedef struct sKLUmatrix {
+    klu_common *KLUmatrixCommon ;               /* KLU common object */
+    klu_symbolic *KLUmatrixSymbolic ;           /* KLU symbolic object */
+    klu_numeric *KLUmatrixNumeric ;             /* KLU numeric object */
+    int *KLUmatrixAp ;                          /* KLU column pointer */
+    int *KLUmatrixAi ;                          /* KLU row pointer */
+    double *KLUmatrixAxComplex ;                /* KLU Complex Elements */
+    unsigned int KLUmatrixIsComplex:1 ;         /* KLU Matrix Is Complex Flag */
+    double *KLUmatrixIntermediateComplex ;      /* KLU iRHS Intermediate for Solve Complex Step */
+    unsigned int KLUmatrixN ;                   /* KLU N, copied */
+    unsigned int KLUmatrixNZ ;                  /* KLU nz, copied for AC Analysis */
+    int *KLUmatrixColCOO ;                      /* KLU Col Index for COO storage */
+    int *KLUmatrixRowCOO ;                      /* KLU Row Index for COO storage */
+    double *KLUmatrixValueComplexCOO ;          /* KLU Complex Elements for COO storage */
+    BindKluElementCOO *KLUmatrixBindStructCOO ; /* KLU COO Binding Structure */
+    double *KLUmatrixTrashCOO ;                 /* KLU COO Trash Pointer for Ground Node not Stored in the Matrix */
+} KLUmatrix ;
+#endif
 
 /* ========================================================================== */
 /* === KLU version ========================================================== */
diff --git a/src/include/ngspice/onemesh.h b/src/include/ngspice/onemesh.h
index 41238f71e..6228ab825 100644
--- a/src/include/ngspice/onemesh.h
+++ b/src/include/ngspice/onemesh.h
@@ -16,6 +16,10 @@ Authors: 1987 Karti Mayaram, 1991 David Gates
 
 #include "ngspice/material.h"
 
+#ifdef KLU
+#include "ngspice/klu.h"
+#endif
+
 typedef struct sONEelem {
     struct sONEelem *pElems[2];     /* array to store neighbor elements */
     struct sONEnode *pNodes[2];     /* array to store the element nodes */
@@ -27,6 +31,20 @@ typedef struct sONEelem {
     ONEmaterial *matlInfo;          /* material information */
     double epsRel;                  /* epsilon of material rel to Semicon */
     int evalNodes[2];               /* nodes to be evaluated in elem */
+
+#ifdef KLU
+    double *KLUleftLeftNode ;
+    double *KLUleftRightNode ;
+    double *KLUrightLeftNode ;
+    double *KLUrightRightNode ;
+
+    BindKluElementCOO *KLUleftLeftNodeBinding ;
+    BindKluElementCOO *KLUleftRightNodeBinding ;
+    BindKluElementCOO *KLUrightLeftNodeBinding ;
+    BindKluElementCOO *KLUrightRightNodeBinding ;
+
+#endif
+
 } ONEelem;
 
 #define pLeftElem pElems[0]
@@ -112,6 +130,33 @@ typedef struct sONEnode {
     double *fPNiM1;
     double *fPN;
     double *fPNiP1;
+
+#ifdef KLU
+    BindKluElementCOO *fPsiPsiiM1Binding ;
+    BindKluElementCOO *fPsiPsiBinding ;
+    BindKluElementCOO *fPsiPsiiP1Binding ;
+    BindKluElementCOO *fPsiNBinding ;
+    BindKluElementCOO *fPsiPBinding ;
+    BindKluElementCOO *fNPsiiM1Binding ;
+    BindKluElementCOO *fNPsiBinding ;
+    BindKluElementCOO *fNPsiiP1Binding ;
+    BindKluElementCOO *fNNiM1Binding ;
+    BindKluElementCOO *fNNBinding ;
+    BindKluElementCOO *fNNiP1Binding ;
+    BindKluElementCOO *fNPiM1Binding ;
+    BindKluElementCOO *fNPBinding ;
+    BindKluElementCOO *fNPiP1Binding ;
+    BindKluElementCOO *fPPsiiM1Binding ;
+    BindKluElementCOO *fPPsiBinding ;
+    BindKluElementCOO *fPPsiiP1Binding ;
+    BindKluElementCOO *fPPiM1Binding ;
+    BindKluElementCOO *fPPBinding ;
+    BindKluElementCOO *fPPiP1Binding ;
+    BindKluElementCOO *fPNiM1Binding ;
+    BindKluElementCOO *fPNBinding ;
+    BindKluElementCOO *fPNiP1Binding ;
+#endif
+
 } ONEnode;
 
 
diff --git a/src/include/ngspice/smpdefs.h b/src/include/ngspice/smpdefs.h
index a4fd24745..2b8a40ff4 100644
--- a/src/include/ngspice/smpdefs.h
+++ b/src/include/ngspice/smpdefs.h
@@ -47,6 +47,11 @@ typedef struct sSMPmatrix {
     unsigned int CKTkluMODE:1 ;            /* KLU MODE parameter to enable KLU or not from the heuristic */
     #define CKTkluON 1                     /* KLU MODE ON definition */
     #define CKTkluOFF 0                    /* KLU MODE OFF definition */
+
+#ifdef CIDER
+    KLUmatrix *SMPkluMatrix ;              /* KLU Pointer to the KLU Matrix Data Structure (only for CIDER, for the moment) */
+#endif
+
 #endif
 
 } SMPmatrix ;
@@ -56,7 +61,21 @@ typedef struct sSMPmatrix {
 void SMPmatrix_CSC (SMPmatrix *) ;
 void SMPnnz (SMPmatrix *) ;
 void spDeterminant_KLU (SMPmatrix *, int *, double *, double *) ;
+
+#ifdef CIDER
+void SMPsolveKLUforCIDER (SMPmatrix *, double [], double [], double [], double []) ;
+int SMPreorderKLUforCIDER (SMPmatrix *) ;
+double *SMPmakeEltKLUforCIDER (SMPmatrix *, int, int) ;
+void SMPclearKLUforCIDER (SMPmatrix *) ;
+void SMPconvertCOOtoCSCKLUforCIDER (SMPmatrix *) ;
+void SMPdestroyKLUforCIDER (SMPmatrix *) ;
+int SMPnewMatrixKLUforCIDER (SMPmatrix *, int, unsigned int) ;
+int SMPluFacKLUforCIDER (SMPmatrix *) ;
+void SMPprintKLUforCIDER (SMPmatrix *, char *) ;
 #endif
+
+#endif
+
 int SMPaddElt( SMPmatrix *, int , int , double );
 double * SMPmakeElt( SMPmatrix * , int , int );
 void SMPcClear( SMPmatrix *);
@@ -85,4 +104,11 @@ int SMPzeroRow(SMPmatrix *Matrix, int Row);
 void SMPconstMult(SMPmatrix *, double);
 void SMPmultiply(SMPmatrix *, double *, double *, double *, double *);
 
+#ifdef CIDER
+void SMPcSolveForCIDER (SMPmatrix *, double [], double [], double [], double []) ;
+int SMPluFacForCIDER (SMPmatrix *) ;
+int SMPnewMatrixForCIDER (SMPmatrix *, int, int) ;
+void SMPsolveForCIDER (SMPmatrix *, double [], double []) ;
+#endif
+
 #endif
diff --git a/src/maths/KLU/klu_extract.c b/src/maths/KLU/klu_extract.c
index f51041b8b..3fdf50712 100644
--- a/src/maths/KLU/klu_extract.c
+++ b/src/maths/KLU/klu_extract.c
@@ -448,9 +448,17 @@ Int KLU_print
         {
 
 #ifdef COMPLEX
-            fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j].Real, Az [j].Imag) ;
+            if (IntToExtRowMap && IntToExtColMap) {
+                fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j].Real, Az [j].Imag) ;
+            } else {
+                fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", Ai [j] + 1, i + 1, Az [j].Real, Az [j].Imag) ;
+            }
 #else
-            fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j]) ;
+            if (IntToExtRowMap && IntToExtColMap) {
+                fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j]) ;
+            } else {
+                fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", Ai [j] + 1, i + 1, Az [j]) ;
+            }
 #endif
 
         }
diff --git a/src/maths/KLU/klusmp.c b/src/maths/KLU/klusmp.c
index a22c4669d..ff375393e 100644
--- a/src/maths/KLU/klusmp.c
+++ b/src/maths/KLU/klusmp.c
@@ -130,6 +130,170 @@ SMPnnz (SMPmatrix *Matrix)
     return ;
 }
 
+#ifdef CIDER
+typedef struct sElement {
+  unsigned int row ;
+  unsigned int col ;
+  double *pointer ;
+} Element ;
+
+static int
+CompareRow (const void *a, const void *b)
+{
+    Element *A = (Element *) a ;
+    Element *B = (Element *) b ;
+
+    return
+        (A->row > B->row) ?  1 :
+        (A->row < B->row) ? -1 :
+        0 ;
+}
+
+static int
+CompareColumn (const void *a, const void *b)
+{
+    Element *A = (Element *) a ;
+    Element *B = (Element *) b ;
+
+    return
+        (A->col > B->col) ?  1 :
+        (A->col < B->col) ? -1 :
+        0 ;
+}
+
+static void
+Compress (unsigned int *Ai, unsigned int *Bp, unsigned int n, unsigned int nz)
+{
+    unsigned int i, j ;
+
+    for (i = 0 ; i <= Ai [0] ; i++)
+        Bp [i] = 0 ;
+
+    j = Ai [0] + 1 ;
+    for (i = 1 ; i < nz ; i++)
+    {
+        if (Ai [i] == Ai [i - 1] + 1)
+        {
+            Bp [j] = i ;
+            j++ ;
+        }
+        else if (Ai [i] > Ai [i - 1] + 1)
+        {
+            for ( ; j <= Ai [i] ; j++)
+                Bp [j] = i ;
+        }
+    }
+
+    for ( ; j <= n ; j++)
+        Bp [j] = i ;
+}
+
+int
+BindKluCompareCOO (const void *a, const void *b)
+{
+    BindKluElementCOO *A = (BindKluElementCOO *) a ;
+    BindKluElementCOO *B = (BindKluElementCOO *) b ;
+
+    return
+        (A->COO > B->COO) ?  1 :
+        (A->COO < B->COO) ? -1 :
+        0 ;
+}
+
+int
+BindKluCompareCSC (const void *a, const void *b)
+{
+    BindKluElementCOO *A = (BindKluElementCOO *) a ;
+    BindKluElementCOO *B = (BindKluElementCOO *) b ;
+
+    return
+        (A->CSC_Complex > B->CSC_Complex) ?  1 :
+        (A->CSC_Complex < B->CSC_Complex) ? -1 :
+        0 ;
+}
+
+void SMPconvertCOOtoCSCKLUforCIDER (SMPmatrix *Matrix)
+{
+    Element *MatrixCOO ;
+    unsigned int *Ap_COO, i, j, nz ;
+
+    /* Count the non-zero elements and store it */
+    nz = 0 ;
+    for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
+        if ((Matrix->SMPkluMatrix->KLUmatrixRowCOO [i] != -1) && (Matrix->SMPkluMatrix->KLUmatrixColCOO [i] != -1)) {
+            nz++ ;
+        }
+    }
+    Matrix->SMPkluMatrix->KLUmatrixNZ = nz ;
+
+    /* Allocate the compressed COO elements */
+    MatrixCOO = (Element *) malloc (nz * sizeof(Element)) ;
+
+    /* Allocate the temporary COO Column Index */
+    Ap_COO = (unsigned int *) malloc (nz * sizeof(unsigned int)) ;
+
+    /* Allocate the needed KLU data structures */
+    Matrix->SMPkluMatrix->KLUmatrixAp = (int *) malloc ((Matrix->SMPkluMatrix->KLUmatrixN + 1) * sizeof(int)) ;
+    Matrix->SMPkluMatrix->KLUmatrixAi = (int *) malloc (nz * sizeof(int)) ;
+    Matrix->SMPkluMatrix->KLUmatrixBindStructCOO = (BindKluElementCOO *) malloc (nz * sizeof(BindKluElementCOO)) ;
+    Matrix->SMPkluMatrix->KLUmatrixAxComplex = (double *) malloc (2 * nz * sizeof(double)) ;
+    Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = (double *) malloc (2 * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(double)) ;
+
+    /* Populate the compressed COO elements and COO value of Binding Table */
+    j = 0 ;
+    for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
+        if ((Matrix->SMPkluMatrix->KLUmatrixRowCOO [i] != -1) && (Matrix->SMPkluMatrix->KLUmatrixColCOO [i] != -1)) {
+            MatrixCOO [j].row = (unsigned int)Matrix->SMPkluMatrix->KLUmatrixRowCOO [i] ;
+            MatrixCOO [j].col = (unsigned int)Matrix->SMPkluMatrix->KLUmatrixColCOO [i] ;
+            MatrixCOO [j].pointer = &(Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO [2 * i]) ;
+            j++ ;
+        }
+    }
+
+    /* Order the MatrixCOO along the columns */
+    qsort (MatrixCOO, nz, sizeof(Element), CompareColumn) ;
+
+    /* Order the MatrixCOO along the rows */
+    i = 0 ;
+    while (i < nz)
+    {
+        /* Look for the next column */
+        for (j = i + 1 ; j < nz ; j++)
+        {
+            if (MatrixCOO [j].col != MatrixCOO [i].col)
+            {
+                break ;
+            }
+        }
+
+        qsort (MatrixCOO + i, j - i, sizeof(Element), CompareRow) ;
+
+        i = j ;
+    }
+
+    /* Copy back the Matrix in partial CSC */
+    for (i = 0 ; i < nz ; i++)
+    {
+        Ap_COO [i] = MatrixCOO [i].col ;
+        Matrix->SMPkluMatrix->KLUmatrixAi [i] = (int)MatrixCOO [i].row ;
+        Matrix->SMPkluMatrix->KLUmatrixBindStructCOO [i].COO = MatrixCOO [i].pointer ;
+        Matrix->SMPkluMatrix->KLUmatrixBindStructCOO [i].CSC_Complex = &(Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i]) ;
+    }
+
+    /* Compress the COO Column Index to CSC Column Index */
+    Compress (Ap_COO, (unsigned int *)Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixN, nz) ;
+
+    /* Free the temporary stuff */
+    free (Ap_COO) ;
+    free (MatrixCOO) ;
+
+    /* Sort the Binding Table */
+    qsort (Matrix->SMPkluMatrix->KLUmatrixBindStructCOO, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ;
+
+    return ;
+}
+#endif
+
 /*
  * SMPaddElt()
  */
@@ -149,6 +313,26 @@ SMPmakeElt (SMPmatrix *Matrix, int Row, int Col)
     return spGetElement (Matrix->SPmatrix, Row, Col) ;
 }
 
+#ifdef CIDER
+double *
+SMPmakeEltKLUforCIDER (SMPmatrix *Matrix, int Row, int Col)
+{
+    if (Matrix->CKTkluMODE) {
+        if ((Row > 0) && (Col > 0)) {
+            Row = Row - 1 ;
+            Col = Col - 1 ;
+            Matrix->SMPkluMatrix->KLUmatrixRowCOO [Row * (int)Matrix->SMPkluMatrix->KLUmatrixN + Col] = Row ;
+            Matrix->SMPkluMatrix->KLUmatrixColCOO [Row * (int)Matrix->SMPkluMatrix->KLUmatrixN + Col] = Col ;
+            return &(Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO [2 * (Row * (int)Matrix->SMPkluMatrix->KLUmatrixN + Col)]) ;
+        } else {
+            return Matrix->SMPkluMatrix->KLUmatrixTrashCOO ;
+        }
+    } else {
+        return spGetElement (Matrix->SPmatrix, Row, Col) ;
+    }
+}
+#endif
+
 /*
  * SMPcClear()
  */
@@ -191,6 +375,18 @@ SMPclear (SMPmatrix *Matrix)
     }
 }
 
+#ifdef CIDER
+void
+SMPclearKLUforCIDER (SMPmatrix *Matrix)
+{
+    unsigned int i ;
+
+    for (i = 0 ; i < 2 * Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
+        Matrix->SMPkluMatrix->KLUmatrixAxComplex [i] = 0 ;
+    }
+}
+#endif
+
 #define NG_IGNORE(x)  (void)x
 
 /*
@@ -262,6 +458,46 @@ SMPluFac (SMPmatrix *Matrix, double PivTol, double Gmin)
     }
 }
 
+#ifdef CIDER
+int
+SMPluFacKLUforCIDER (SMPmatrix *Matrix)
+{
+    unsigned int i ;
+    double *KLUmatrixAx ;
+    int ret ;
+
+    if (Matrix->CKTkluMODE)
+    {
+        if (Matrix->SMPkluMatrix->KLUmatrixIsComplex) {
+            ret = klu_z_refactor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAxComplex,
+                                  Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+        } else {
+            /* Allocate the Real Matrix */
+            KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof(double)) ;
+
+            /* Copy the Complex Matrix into the Real Matrix */
+            for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
+                KLUmatrixAx [i] = Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] ;
+            }
+
+            /* Re-Factor the Real Matrix */
+            ret = klu_refactor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, KLUmatrixAx,
+                                Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+
+            /* Free the Real Matrix Storage */
+            free (KLUmatrixAx) ;
+        }
+        if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) {
+            printf ("CIDER: KLU Empty Matrix\n") ;
+            return 0 ;
+        }
+        return (!ret) ;
+    } else {
+        return spFactor (Matrix->SPmatrix) ;
+    }
+}
+#endif
+
 /*
  * SMPcReorder()
  */
@@ -336,6 +572,58 @@ SMPreorder (SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin)
     }
 }
 
+#ifdef CIDER
+int
+SMPreorderKLUforCIDER (SMPmatrix *Matrix)
+{
+    unsigned int i ;
+    double *KLUmatrixAx ;
+
+    if (Matrix->CKTkluMODE)
+    {
+//        Matrix->CKTkluCommon->tol = PivTol ;
+
+        if (Matrix->SMPkluMatrix->KLUmatrixNumeric != NULL) {
+            klu_free_numeric (&(Matrix->SMPkluMatrix->KLUmatrixNumeric), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+        }
+        if (Matrix->SMPkluMatrix->KLUmatrixIsComplex) {
+            Matrix->SMPkluMatrix->KLUmatrixNumeric = klu_z_factor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi,
+                                                                   Matrix->SMPkluMatrix->KLUmatrixAxComplex, Matrix->SMPkluMatrix->KLUmatrixSymbolic,
+                                                                   Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+        } else {
+            /* Allocate the Real Matrix */
+            KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof(double)) ;
+
+            /* Copy the Complex Matrix into the Real Matrix */
+            for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
+                KLUmatrixAx [i] = Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] ;
+            }
+
+            /* Factor the Real Matrix */
+            Matrix->SMPkluMatrix->KLUmatrixNumeric = klu_factor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi,
+                                                                 KLUmatrixAx, Matrix->SMPkluMatrix->KLUmatrixSymbolic,
+                                                                 Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+
+            /* Free the Real Matrix Storage */
+            free (KLUmatrixAx) ;
+        }
+        if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL)
+        {
+            printf ("CIDER: KLU Factorization Error\n") ;
+            if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX)
+            {
+                return 0 ;
+            }
+            return 1 ;
+        }
+        else
+            return 0 ;
+    } else {
+        return spFactor (Matrix->SPmatrix) ;
+    }
+}
+#endif
+
 /*
  * SMPcaSolve()
  */
@@ -405,6 +693,57 @@ SMPcSolve (SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[], doubl
     }
 }
 
+#ifdef CIDER
+void
+SMPsolveKLUforCIDER (SMPmatrix *Matrix, double RHS[], double RHSsolution[], double iRHS[], double iRHSsolution[])
+{
+    int ret ;
+    unsigned int i ;
+    double *KLUmatrixIntermediate ;
+
+    if (Matrix->CKTkluMODE)
+    {
+        if (Matrix->SMPkluMatrix->KLUmatrixIsComplex) {
+            for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
+            {
+                Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] = RHS [i + 1] ;
+                Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] = iRHS [i + 1] ;
+            }
+
+            ret = klu_z_solve (Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, (int)Matrix->SMPkluMatrix->KLUmatrixN, 1,
+                               Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+
+            for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
+            {
+                RHSsolution [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] ;
+                iRHSsolution [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] ;
+            }
+        } else {
+            /* Allocate the Intermediate Vector */
+            KLUmatrixIntermediate = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * sizeof(double)) ;
+
+            for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
+                KLUmatrixIntermediate [i] = RHS [i + 1] ;
+            }
+
+            ret = klu_solve (Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, (int)Matrix->SMPkluMatrix->KLUmatrixN, 1,
+                             KLUmatrixIntermediate, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+
+            for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
+                RHSsolution [i + 1] = KLUmatrixIntermediate [i] ;
+            }
+
+            /* Free the Intermediate Vector */
+            free (KLUmatrixIntermediate) ;
+        }
+
+    } else {
+
+        spSolve (Matrix->SPmatrix, RHS, RHSsolution, iRHS, iRHSsolution) ;
+    }
+}
+#endif
+
 /*
  * SMPsolve()
  */
@@ -452,6 +791,58 @@ SMPnewMatrix (SMPmatrix *Matrix, int size)
     return Error ;
 }
 
+#ifdef CIDER
+int
+SMPnewMatrixKLUforCIDER (SMPmatrix *Matrix, int size, unsigned int KLUmatrixIsComplex)
+{
+    int Error ;
+    unsigned int i ;
+
+    if (Matrix->CKTkluMODE) {
+        /* Allocate the KLU Matrix Data Structure */
+        Matrix->SMPkluMatrix = (KLUmatrix *) malloc (sizeof (KLUmatrix)) ;
+
+        /* Initialize the KLU Matrix Internal Pointers */
+        Matrix->SMPkluMatrix->KLUmatrixCommon = (klu_common *) malloc (sizeof (klu_common)) ; ;
+        Matrix->SMPkluMatrix->KLUmatrixSymbolic = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixNumeric = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixAp = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixAi = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixAxComplex = NULL ;
+        if (KLUmatrixIsComplex) {
+            Matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
+        } else {
+            Matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixReal ;
+        }
+        Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixNZ = 0 ;
+        Matrix->SMPkluMatrix->KLUmatrixBindStructCOO = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO = NULL ;
+
+        /* Initialize the KLU Common Data Structure */
+        klu_defaults (Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+
+        /* Allocate KLU data structures */
+        Matrix->SMPkluMatrix->KLUmatrixN = (unsigned int)size ;
+        Matrix->SMPkluMatrix->KLUmatrixColCOO = (int *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(int)) ;
+        Matrix->SMPkluMatrix->KLUmatrixRowCOO = (int *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(int)) ;
+        Matrix->SMPkluMatrix->KLUmatrixTrashCOO = (double *) malloc (sizeof(double)) ;
+        Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO = (double *) malloc (2 * Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(double)) ;
+
+        /* Pre-set the values of Row and Col */
+        for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
+            Matrix->SMPkluMatrix->KLUmatrixRowCOO [i] = -1 ;
+            Matrix->SMPkluMatrix->KLUmatrixColCOO [i] = -1 ;
+        }
+
+        return spOKAY ;
+    } else {
+        Matrix->SPmatrix = spCreate (size, (int)KLUmatrixIsComplex, &Error) ;
+        return Error ;
+    }
+}
+#endif
+
 /*
  * SMPdestroy()
  */
@@ -486,6 +877,38 @@ SMPdestroy (SMPmatrix *Matrix)
     }
 }
 
+#ifdef CIDER
+void
+SMPdestroyKLUforCIDER (SMPmatrix *Matrix)
+{
+    if (Matrix->CKTkluMODE)
+    {
+        klu_free_numeric (&(Matrix->SMPkluMatrix->KLUmatrixNumeric), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+        klu_free_symbolic (&(Matrix->SMPkluMatrix->KLUmatrixSymbolic), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
+        free (Matrix->SMPkluMatrix->KLUmatrixAp) ;
+        free (Matrix->SMPkluMatrix->KLUmatrixAi) ;
+        free (Matrix->SMPkluMatrix->KLUmatrixAxComplex) ;
+        free (Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex) ;
+        free (Matrix->SMPkluMatrix->KLUmatrixBindStructCOO) ;
+        free (Matrix->SMPkluMatrix->KLUmatrixColCOO) ;
+        free (Matrix->SMPkluMatrix->KLUmatrixRowCOO) ;
+        free (Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO) ;
+        free (Matrix->SMPkluMatrix->KLUmatrixTrashCOO) ;
+        Matrix->SMPkluMatrix->KLUmatrixAp = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixAi = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixAxComplex = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixBindStructCOO = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixColCOO = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixRowCOO = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO = NULL ;
+        Matrix->SMPkluMatrix->KLUmatrixTrashCOO = NULL ;
+    } else {
+        spDestroy (Matrix->SPmatrix) ;
+    }
+}
+#endif
+
 /*
  * SMPpreOrder()
  */
@@ -547,6 +970,43 @@ SMPprint (SMPmatrix *Matrix, char *Filename)
     }
 }
 
+#ifdef CIDER
+void
+SMPprintKLUforCIDER (SMPmatrix *Matrix, char *Filename)
+{
+    unsigned int i ;
+    double *KLUmatrixAx ;
+
+    if (Matrix->CKTkluMODE)
+    {
+        if (Matrix->SMPkluMatrix->KLUmatrixIsComplex)
+        {
+            klu_z_print (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAxComplex,
+                         (int)Matrix->SMPkluMatrix->KLUmatrixN, NULL, NULL) ;
+        } else {
+            /* Allocate the Real Matrix */
+            KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof(double)) ;
+
+            /* Copy the Complex Matrix into the Real Matrix */
+            for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
+                KLUmatrixAx [i] = Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] ;
+            }
+
+            /* Print the Real Matrix */
+            klu_print (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, KLUmatrixAx, (int)Matrix->SMPkluMatrix->KLUmatrixN, NULL, NULL) ;
+
+            /* Free the Real Matrix Storage */
+            free (KLUmatrixAx) ;
+        }
+    } else {
+        if (Filename)
+            spFileMatrix (Matrix->SPmatrix, Filename, "Circuit Matrix", 0, 1, 1) ;
+        else
+            spPrint (Matrix->SPmatrix, 0, 1, 1) ;
+    }
+}
+#endif
+
 /*
  * SMPgetError()
  */
@@ -562,7 +1022,6 @@ SMPgetError (SMPmatrix *Matrix, int *Col, int *Row)
     }
 }
 
-#ifdef KLU
 void
 spDeterminant_KLU (SMPmatrix *Matrix, int *pExponent, RealNumber *pDeterminant, RealNumber *piDeterminant)
 {
@@ -819,7 +1278,6 @@ spDeterminant_KLU (SMPmatrix *Matrix, int *pExponent, RealNumber *pDeterminant,
     free (Ux) ;
     free (Rs) ;
 }
-#endif
 
 /*
  * SMPcProdDiag()
diff --git a/src/maths/sparse/spsmp.c b/src/maths/sparse/spsmp.c
index a7d64d17c..7a230dfc7 100644
--- a/src/maths/sparse/spsmp.c
+++ b/src/maths/sparse/spsmp.c
@@ -174,6 +174,14 @@ SMPluFac(SMPmatrix *Matrix, double PivTol, double Gmin)
     return spFactor( Matrix->SPmatrix );
 }
 
+#ifdef CIDER
+int
+SMPluFacForCIDER (SMPmatrix *Matrix)
+{
+    return spFactor (Matrix->SPmatrix) ;
+}
+#endif
+
 /*
  * SMPcReorder()
  */
@@ -225,6 +233,14 @@ SMPcSolve(SMPmatrix *Matrix, double RHS[], double iRHS[],
     spSolve( Matrix->SPmatrix, RHS, RHS, iRHS, iRHS );
 }
 
+#ifdef CIDER
+void
+SMPcSolveForCIDER (SMPmatrix *Matrix, double RHS[], double RHSsolution[], double iRHS[], double iRHSsolution[])
+{
+    spSolve (Matrix->SPmatrix, RHS, RHSsolution, iRHS, iRHSsolution) ;
+}
+#endif
+
 /*
  * SMPsolve()
  */
@@ -236,6 +252,14 @@ SMPsolve(SMPmatrix *Matrix, double RHS[], double Spare[])
     spSolve( Matrix->SPmatrix, RHS, RHS, NULL, NULL );
 }
 
+#ifdef CIDER
+void
+SMPsolveForCIDER (SMPmatrix *Matrix, double RHS[], double RHSsolution[])
+{
+    spSolve (Matrix->SPmatrix, RHS, RHSsolution, NULL, NULL) ;
+}
+#endif
+
 /*
  * SMPmatSize()
  */
@@ -256,6 +280,19 @@ SMPnewMatrix(SMPmatrix *Matrix, int size)
     return Error;
 }
 
+#ifdef CIDER
+int
+SMPnewMatrixForCIDER (SMPmatrix *Matrix, int size, int complex)
+{
+    int error ;
+
+    Matrix->SPmatrix = spCreate (size, complex, &error) ;
+
+    return error ;
+}
+
+#endif
+
 /*
  * SMPdestroy()
  */