diff --git a/src/include/ngspice/klu-binding.h b/src/include/ngspice/klu-binding.h
new file mode 100644
index 000000000..b3f5fa48d
--- /dev/null
+++ b/src/include/ngspice/klu-binding.h
@@ -0,0 +1,20 @@
+#ifndef _KLU_BINDING_H
+#define _KLU_BINDING_H
+
+#define CREATE_KLU_BINDING_TABLE(ptr, binding, a, b)                    \
+    if ((here->a != 0) && (here->b != 0)) {                             \
+        i = here->ptr ;                                                 \
+        matched = (BindElement *) bsearch (&i, BindStruct, nz, sizeof(BindElement), BindCompare) ; \
+        here->binding = matched ;                                       \
+        here->ptr = matched->CSC ;                                      \
+    }
+
+#define CONVERT_KLU_BINDING_TABLE_TO_COMPLEX(ptr, binding, a, b)        \
+    if ((here->a != 0) && (here->b != 0))                               \
+        here->ptr = here->binding->CSC_Complex ;
+
+#define CONVERT_KLU_BINDING_TABLE_TO_REAL(ptr, binding, a, b)   \
+    if ((here->a != 0) && (here->b != 0))                       \
+        here->ptr = here->binding->CSC ;
+
+#endif
diff --git a/src/maths/KLU/Makefile.am b/src/maths/KLU/Makefile.am
new file mode 100644
index 000000000..b972ae12c
--- /dev/null
+++ b/src/maths/KLU/Makefile.am
@@ -0,0 +1,97 @@
+## Process this file with automake to produce Makefile.in
+
+noinst_LTLIBRARIES = libKLU_real.la libKLU_complex.la libKLU.la
+
+libKLU_real_la_SOURCES = 	\
+	amd_1.c			\
+	amd_2.c			\
+	amd_aat.c		\
+	amd_control.c		\
+	amd_defaults.c		\
+	amd_dump.c		\
+	amd_global.c		\
+	amd_info.c		\
+	amd_order.c		\
+	amd_postorder.c		\
+	amd_post_tree.c		\
+	amd_preprocess.c	\
+	amd_valid.c		\
+	btf_maxtrans.c		\
+	btf_order.c		\
+	btf_strongcomp.c	\
+	colamd.c		\
+	colamd_global.c		\
+	klu.c			\
+	klu_analyze.c		\
+	klu_analyze_given.c	\
+	klu_defaults.c		\
+	klu_diagnostics.c	\
+	klu_dump.c		\
+	klu_extract.c		\
+	klu_factor.c		\
+	klu_free_numeric.c 	\
+	klu_free_symbolic.c	\
+	klu_kernel.c		\
+	klu_memory.c		\
+	klu_refactor.c		\
+	klu_scale.c		\
+	klu_solve.c		\
+	klu_sort.c		\
+	klu_tsolve.c
+
+
+libKLU_real_la_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include
+
+
+libKLU_complex_la_SOURCES = 	\
+	amd_1.c			\
+	amd_2.c			\
+	amd_aat.c		\
+	amd_control.c		\
+	amd_defaults.c		\
+	amd_dump.c		\
+	amd_global.c		\
+	amd_info.c		\
+	amd_order.c		\
+	amd_postorder.c		\
+	amd_post_tree.c		\
+	amd_preprocess.c	\
+	amd_valid.c		\
+	btf_maxtrans.c		\
+	btf_order.c		\
+	btf_strongcomp.c	\
+	colamd.c		\
+	colamd_global.c		\
+	klu.c			\
+	klu_analyze.c		\
+	klu_analyze_given.c	\
+	klu_defaults.c		\
+	klu_diagnostics.c	\
+	klu_dump.c		\
+	klu_extract.c		\
+	klu_factor.c		\
+	klu_free_numeric.c 	\
+	klu_free_symbolic.c	\
+	klu_kernel.c		\
+	klu_memory.c		\
+	klu_refactor.c		\
+	klu_scale.c		\
+	klu_solve.c		\
+	klu_sort.c		\
+	klu_tsolve.c
+
+
+libKLU_complex_la_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include -DCOMPLEX
+
+
+libKLU_la_SOURCES =		\
+	klusmp.c
+
+libKLU_la_LIBADD =		\
+	libKLU_real.la		\
+	libKLU_complex.la
+
+
+libKLU_la_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include
+
+MAINTAINERCLEANFILES = Makefile.in
diff --git a/src/maths/KLU/klusmp.c b/src/maths/KLU/klusmp.c
new file mode 100644
index 000000000..1d07fec03
--- /dev/null
+++ b/src/maths/KLU/klusmp.c
@@ -0,0 +1,748 @@
+/*
+ *  Spice3 COMPATIBILITY MODULE
+ *
+ *  Author:                     Advising professor:
+ *     Kenneth S. Kundert           Alberto Sangiovanni-Vincentelli
+ *     UC Berkeley
+ *
+ *  This module contains routines that make Sparse1.3 a direct
+ *  replacement for the SMP sparse matrix package in Spice3c1 or Spice3d1.
+ *  Sparse1.3 is in general a faster and more robust package than SMP.
+ *  These advantages become significant on large circuits.
+ *
+ *  >>> User accessible functions contained in this file:
+ *  SMPaddElt
+ *  SMPmakeElt
+ *  SMPcClear
+ *  SMPclear
+ *  SMPcLUfac
+ *  SMPluFac
+ *  SMPcReorder
+ *  SMPreorder
+ *  SMPcaSolve
+ *  SMPcSolve
+ *  SMPsolve
+ *  SMPmatSize
+ *  SMPnewMatrix
+ *  SMPdestroy
+ *  SMPpreOrder
+ *  SMPprint
+ *  SMPgetError
+ *  SMPcProdDiag
+ *  LoadGmin
+ *  SMPfindElt
+ *  SMPcombine
+ *  SMPcCombine
+ */
+
+/*
+ *  To replace SMP with Sparse, rename the file spSpice3.h to
+ *  spMatrix.h and place Sparse in a subdirectory of SPICE called
+ *  `sparse'.  Then on UNIX compile Sparse by executing `make spice'.
+ *  If not on UNIX, after compiling Sparse and creating the sparse.a
+ *  archive, compile this file (spSMP.c) and spSMP.o to the archive,
+ *  then copy sparse.a into the SPICE main directory and rename it
+ *  SMP.a.  Finally link SPICE.
+ *
+ *  To be compatible with SPICE, the following Sparse compiler options
+ *  (in spConfig.h) should be set as shown below:
+ *
+ *      EXPANDABLE                      YES
+ *      TRANSLATE                       NO
+ *      INITIALIZE                      NO or YES, YES for use with test prog.
+ *      DIAGONAL_PIVOTING               YES
+ *      MODIFIED_MARKOWITZ              NO
+ *      DELETE                          NO
+ *      STRIP                           NO
+ *      MODIFIED_NODAL                  YES
+ *      QUAD_ELEMENT                    NO
+ *      TRANSPOSE                       YES
+ *      SCALING                         NO
+ *      DOCUMENTATION                   YES
+ *      MULTIPLICATION                  NO
+ *      DETERMINANT                     YES
+ *      STABILITY                       NO
+ *      CONDITION                       NO
+ *      PSEUDOCONDITION                 NO
+ *      DEBUG                           YES
+ *
+ *      spREAL  double
+ */
+
+/*
+ *  Revision and copyright information.
+ *
+ *  Copyright (c) 1985,86,87,88,89,90
+ *  by Kenneth S. Kundert and the University of California.
+ *
+ *  Permission to use, copy, modify, and distribute this software and its
+ *  documentation for any purpose and without fee is hereby granted, provided
+ *  that the above copyright notice appear in all copies and supporting
+ *  documentation and that the authors and the University of California
+ *  are properly credited.  The authors and the University of California
+ *  make no representations as to the suitability of this software for
+ *  any purpose.  It is provided `as is', without express or implied warranty.
+ */
+
+/*
+ *  IMPORTS
+ *
+ *  >>> Import descriptions:
+ *  spMatrix.h
+ *     Sparse macros and declarations.
+ *  SMPdefs.h
+ *     Spice3's matrix macro definitions.
+ */
+
+#include "ngspice/config.h"
+#include <assert.h>
+#include <stdio.h>
+#include <math.h>
+#include "ngspice/spmatrix.h"
+#include "../sparse/spdefs.h"
+#include "ngspice/smpdefs.h"
+
+#if defined (_MSC_VER)
+extern double scalbn(double, int);
+#define logb _logb
+extern double logb(double);
+#endif
+
+static void LoadGmin_CSC (double **diag, int n, double Gmin) ;
+static void LoadGmin (SMPmatrix *eMatrix, double Gmin) ;
+
+void
+SMPmatrix_CSC (SMPmatrix *Matrix)
+{
+    spMatrix_CSC (Matrix->SPmatrix, Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx,
+                  Matrix->CKTkluAx_Complex, Matrix->CKTkluN, Matrix->CKTbindStruct, Matrix->CKTdiag_CSC) ;
+    return ;
+}
+
+void
+SMPnnz (SMPmatrix *Matrix)
+{
+    Matrix->CKTklunz = Matrix->SPmatrix->Elements ;
+
+    return ;
+}
+
+/*
+ * SMPaddElt()
+ */
+int
+SMPaddElt (SMPmatrix *Matrix, int Row, int Col, double Value)
+{
+    *spGetElement (Matrix->SPmatrix, Row, Col) = Value ;
+    return spError (Matrix->SPmatrix) ;
+}
+
+/*
+ * SMPmakeElt()
+ */
+double *
+SMPmakeElt (SMPmatrix *Matrix, int Row, int Col)
+{
+    return spGetElement (Matrix->SPmatrix, Row, Col) ;
+}
+
+/*
+ * SMPcClear()
+ */
+
+void
+SMPcClear (SMPmatrix *Matrix)
+{
+    int i ;
+    if (Matrix->CKTkluMODE)
+    {
+        spClear (Matrix->SPmatrix) ;
+        if (Matrix->CKTkluAx_Complex != NULL)
+        {
+            for (i = 0 ; i < 2 * Matrix->CKTklunz ; i++)
+                Matrix->CKTkluAx_Complex [i] = 0 ;
+        }
+    } else {
+        spClear (Matrix->SPmatrix) ;
+    }
+}
+
+/*
+ * SMPclear()
+ */
+
+void
+SMPclear (SMPmatrix *Matrix)
+{
+    int i ;
+    if (Matrix->CKTkluMODE)
+    {
+        spClear (Matrix->SPmatrix) ;
+        if (Matrix->CKTkluAx != NULL)
+        {
+            for (i = 0 ; i < Matrix->CKTklunz ; i++)
+                Matrix->CKTkluAx [i] = 0 ;
+        }
+    } else {
+        spClear (Matrix->SPmatrix) ;
+    }
+}
+
+#define NG_IGNORE(x)  (void)x
+
+/*
+ * SMPcLUfac()
+ */
+/*ARGSUSED*/
+
+int
+SMPcLUfac (SMPmatrix *Matrix, double PivTol)
+{
+    int ret ;
+
+    NG_IGNORE (PivTol) ;
+
+    if (Matrix->CKTkluMODE)
+    {
+        spSetComplex (Matrix->SPmatrix) ;
+        ret = klu_z_refactor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx_Complex,
+                              Matrix->CKTkluSymbolic, Matrix->CKTkluNumeric, Matrix->CKTkluCommon) ;
+        return (!ret) ;
+    } else {
+        spSetComplex (Matrix->SPmatrix) ;
+        return spFactor (Matrix->SPmatrix) ;
+    }
+}
+
+/*
+ * SMPluFac()
+ */
+/*ARGSUSED*/
+
+int
+SMPluFac (SMPmatrix *Matrix, double PivTol, double Gmin)
+{
+    int ret ;
+
+    NG_IGNORE (PivTol) ;
+
+    if (Matrix->CKTkluMODE)
+    {
+        spSetReal (Matrix->SPmatrix) ;
+        LoadGmin_CSC (Matrix->CKTdiag_CSC, Matrix->CKTkluN, Gmin) ;
+        ret = klu_refactor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx,
+                            Matrix->CKTkluSymbolic, Matrix->CKTkluNumeric, Matrix->CKTkluCommon) ;
+        return (!ret) ;
+
+//        if (ret == 1)
+//            return 0 ;
+//        else if (ret == 0)
+//            return (E_SINGULAR) ;
+//        else {
+//            fprintf (stderr, "KLU Error in re-factor!") ;
+//            return 1 ;
+//        }
+    } else {
+        spSetReal (Matrix->SPmatrix) ;
+        LoadGmin (Matrix, Gmin) ;
+        return spFactor (Matrix->SPmatrix) ;
+    }
+}
+
+/*
+ * SMPcReorder()
+ */
+
+int
+SMPcReorder (SMPmatrix *Matrix, double PivTol, double PivRel, int *NumSwaps)
+{
+    if (Matrix->CKTkluMODE)
+    {
+        *NumSwaps = 1 ;
+        spSetComplex (Matrix->SPmatrix) ;
+
+        if (Matrix->CKTkluNumeric != NULL)
+        {
+            klu_z_free_numeric (&(Matrix->CKTkluNumeric), Matrix->CKTkluCommon) ;
+            Matrix->CKTkluNumeric = klu_z_factor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx_Complex, Matrix->CKTkluSymbolic, Matrix->CKTkluCommon) ;
+        } else
+            Matrix->CKTkluNumeric = klu_z_factor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx_Complex, Matrix->CKTkluSymbolic, Matrix->CKTkluCommon) ;
+
+        if (Matrix->CKTkluNumeric == NULL)
+            return 1 ;
+        else
+            return 0 ;
+    } else {
+        *NumSwaps = 1 ;
+        spSetComplex (Matrix->SPmatrix) ;
+        return spOrderAndFactor (Matrix->SPmatrix, NULL, (spREAL)PivRel, (spREAL)PivTol, YES) ;
+    }
+}
+
+/*
+ * SMPreorder()
+ */
+
+int
+SMPreorder (SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin)
+{
+    if (Matrix->CKTkluMODE)
+    {
+        spSetReal (Matrix->SPmatrix) ;
+        LoadGmin_CSC (Matrix->CKTdiag_CSC, Matrix->CKTkluN, Gmin) ;
+
+        if (Matrix->CKTkluNumeric != NULL)
+        {
+            klu_free_numeric (&(Matrix->CKTkluNumeric), Matrix->CKTkluCommon) ;
+            Matrix->CKTkluNumeric = klu_factor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx, Matrix->CKTkluSymbolic, Matrix->CKTkluCommon) ;
+        } else
+            Matrix->CKTkluNumeric = klu_factor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx, Matrix->CKTkluSymbolic, Matrix->CKTkluCommon) ;
+
+        if (Matrix->CKTkluNumeric == NULL)
+            return 1 ;
+        else
+            return 0 ;
+    } else {
+        spSetReal (Matrix->SPmatrix) ;
+        LoadGmin (Matrix, Gmin) ;
+        return spOrderAndFactor (Matrix->SPmatrix, NULL, (spREAL)PivRel, (spREAL)PivTol, YES) ;
+    }
+}
+
+/*
+ * SMPcaSolve()
+ */
+void
+SMPcaSolve (SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[], double iSpare[])
+{
+    printf ("SMPcaSolve\n") ;
+    NG_IGNORE (iSpare) ;
+    NG_IGNORE (Spare) ;
+
+    spSolveTransposed (Matrix->SPmatrix, RHS, RHS, iRHS, iRHS) ;
+}
+
+/*
+ * SMPcSolve()
+ */
+
+void
+SMPcSolve (SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[], double iSpare[])
+{
+    int ret, i, *pExtOrder ;
+
+    NG_IGNORE (iSpare) ;
+    NG_IGNORE (Spare) ;
+
+    if (Matrix->CKTkluMODE)
+    {
+        pExtOrder = &Matrix->SPmatrix->IntToExtRowMap [Matrix->CKTkluN] ;
+        for (i = 2 * Matrix->CKTkluN - 1 ; i > 0 ; i -= 2)
+        {
+            Matrix->CKTkluIntermediate_Complex [i] = RHS [*(pExtOrder)] ;
+            Matrix->CKTkluIntermediate_Complex [i - 1] = iRHS [*(pExtOrder--)] ;
+        }
+
+        ret = klu_z_solve (Matrix->CKTkluSymbolic, Matrix->CKTkluNumeric, Matrix->CKTkluN, 1, Matrix->CKTkluIntermediate_Complex, Matrix->CKTkluCommon) ;
+
+        pExtOrder = &Matrix->SPmatrix->IntToExtColMap [Matrix->CKTkluN] ;
+        for (i = 2 * Matrix->CKTkluN - 1 ; i > 0 ; i -= 2)
+        {
+            RHS [*(pExtOrder)] = Matrix->CKTkluIntermediate_Complex [i] ;
+            iRHS [*(pExtOrder--)] = Matrix->CKTkluIntermediate_Complex [i - 1] ;
+        }
+
+    } else {
+
+        spSolve (Matrix->SPmatrix, RHS, RHS, iRHS, iRHS) ;
+    }
+}
+
+/*
+ * SMPsolve()
+ */
+
+void
+SMPsolve (SMPmatrix *Matrix, double RHS[], double Spare[])
+{
+    int ret, i, *pExtOrder ;
+
+    NG_IGNORE (Spare) ;
+
+    if (Matrix->CKTkluMODE) {
+
+        pExtOrder = &Matrix->SPmatrix->IntToExtRowMap [Matrix->CKTkluN] ;
+        for (i = Matrix->CKTkluN - 1 ; i >= 0 ; i--)
+            Matrix->CKTkluIntermediate [i] = RHS [*(pExtOrder--)] ;
+
+        ret = klu_solve (Matrix->CKTkluSymbolic, Matrix->CKTkluNumeric, Matrix->CKTkluN, 1, Matrix->CKTkluIntermediate, Matrix->CKTkluCommon) ;
+
+        pExtOrder = &Matrix->SPmatrix->IntToExtColMap [Matrix->CKTkluN] ;
+        for (i = Matrix->CKTkluN - 1 ; i >= 0 ; i--)
+            RHS [*(pExtOrder--)] = Matrix->CKTkluIntermediate [i] ;
+    } else {
+        spSolve (Matrix->SPmatrix, RHS, RHS, NULL, NULL) ;
+    }
+}
+
+/*
+ * SMPmatSize()
+ */
+int
+SMPmatSize (SMPmatrix *Matrix)
+{
+    return spGetSize (Matrix->SPmatrix, 1) ;
+}
+
+/*
+ * SMPnewMatrix()
+ */
+int
+SMPnewMatrix (SMPmatrix *Matrix, int size)
+{
+    int Error ;
+    Matrix->SPmatrix = spCreate (size, 1, &Error) ;
+    return Error ;
+}
+
+/*
+ * SMPdestroy()
+ */
+
+void
+SMPdestroy (SMPmatrix *Matrix)
+{
+    if (Matrix->CKTkluMODE)
+    {
+        spDestroy (Matrix->SPmatrix) ;
+        klu_free_numeric (&(Matrix->CKTkluNumeric), Matrix->CKTkluCommon) ;
+        klu_free_symbolic (&(Matrix->CKTkluSymbolic), Matrix->CKTkluCommon) ;
+        free (Matrix->CKTkluAp) ;
+        free (Matrix->CKTkluAi) ;
+        free (Matrix->CKTkluAx) ;
+        free (Matrix->CKTdiag_CSC) ;
+        free (Matrix->CKTkluIntermediate) ;
+        free (Matrix->CKTkluIntermediate_Complex) ;
+    } else {
+        spDestroy (Matrix->SPmatrix) ;
+    }
+}
+
+/*
+ * SMPpreOrder()
+ */
+
+int
+SMPpreOrder (SMPmatrix *Matrix)
+{
+    if (Matrix->CKTkluMODE)
+    {
+        Matrix->CKTkluSymbolic = klu_analyze (Matrix->CKTkluN, Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluCommon) ;
+
+        return 0 ;
+    } else {
+        spMNA_Preorder (Matrix->SPmatrix) ;
+        return spError (Matrix->SPmatrix) ;
+    }
+}
+
+/*
+ * SMPprintRHS()
+ */
+
+void
+SMPprintRHS (SMPmatrix *Matrix, char *Filename, RealVector RHS, RealVector iRHS)
+{
+    if (!Matrix->CKTkluMODE)
+        spFileVector (Matrix->SPmatrix, Filename, RHS, iRHS) ;
+}
+
+/*
+ * SMPprint()
+ */
+
+void
+SMPprint (SMPmatrix *Matrix, char *Filename)
+{
+    if (!Matrix->CKTkluMODE)
+    {
+        if (Filename)
+            spFileMatrix (Matrix->SPmatrix, Filename, "Circuit Matrix", 0, 1, 1) ;
+        else
+            spPrint (Matrix->SPmatrix, 0, 1, 1) ;
+    }
+}
+
+/*
+ * SMPgetError()
+ */
+void
+SMPgetError (SMPmatrix *Matrix, int *Col, int *Row)
+{
+    if (Matrix->CKTkluMODE)
+    {
+        *Row = Matrix->SPmatrix->IntToExtRowMap [Matrix->CKTkluCommon->singular_col] ;
+        *Col = Matrix->SPmatrix->IntToExtColMap [Matrix->CKTkluCommon->singular_col] ;
+    } else {
+        spWhereSingular (Matrix->SPmatrix, Row, Col) ;
+    }
+}
+
+/*
+ * SMPcProdDiag()
+ *    note: obsolete for Spice3d2 and later
+ */
+int
+SMPcProdDiag (SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
+{
+    spDeterminant (Matrix->SPmatrix, pExponent, &(pMantissa->real), &(pMantissa->imag)) ;
+    return spError (Matrix->SPmatrix) ;
+}
+
+/*
+ * SMPcDProd()
+ */
+int
+SMPcDProd (SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
+{
+    double	re, im, x, y, z;
+    int		p;
+
+    spDeterminant (Matrix->SPmatrix, &p, &re, &im) ;
+
+#ifndef M_LN2
+#define M_LN2   0.69314718055994530942
+#endif
+#ifndef M_LN10
+#define M_LN10  2.30258509299404568402
+#endif
+
+#ifdef debug_print
+    printf ("Determinant 10: (%20g,%20g)^%d\n", re, im, p) ;
+#endif
+
+    /* Convert base 10 numbers to base 2 numbers, for comparison */
+    y = p * M_LN10 / M_LN2;
+    x = (int) y;
+    y -= x;
+
+    /* ASSERT
+     *	x = integral part of exponent, y = fraction part of exponent
+     */
+
+    /* Fold in the fractional part */
+#ifdef debug_print
+    printf (" ** base10 -> base2 int =  %g, frac = %20g\n", x, y) ;
+#endif
+    z = pow (2.0, y) ;
+    re *= z ;
+    im *= z ;
+#ifdef debug_print
+    printf (" ** multiplier = %20g\n", z) ;
+#endif
+
+    /* Re-normalize (re or im may be > 2.0 or both < 1.0 */
+    if (re != 0.0)
+    {
+	y = logb (re) ;
+	if (im != 0.0)
+	    z = logb (im) ;
+	else
+	    z = 0 ;
+    } else if (im != 0.0) {
+	z = logb (im) ;
+	y = 0 ;
+    } else {
+	/* Singular */
+	/*printf("10 -> singular\n");*/
+	y = 0 ;
+	z = 0 ;
+    }
+
+#ifdef debug_print
+    printf (" ** renormalize changes = %g,%g\n", y, z) ;
+#endif
+    if (y < z)
+	y = z ;
+
+    *pExponent = (int)(x + y) ;
+    x = scalbn (re, (int) -y) ;
+    z = scalbn (im, (int) -y) ;
+#ifdef debug_print
+    printf (" ** values are: re %g, im %g, y %g, re' %g, im' %g\n", re, im, y, x, z) ;
+#endif
+    pMantissa->real = scalbn (re, (int) -y) ;
+    pMantissa->imag = scalbn (im, (int) -y) ;
+
+#ifdef debug_print
+    printf ("Determinant 10->2: (%20g,%20g)^%d\n", pMantissa->real, pMantissa->imag, *pExponent) ;
+#endif
+    return spError (Matrix->SPmatrix) ;
+}
+
+
+
+/*
+ *  The following routines need internal knowledge of the Sparse data
+ *  structures.
+ */
+
+/*
+ *  LOAD GMIN
+ *
+ *  This routine adds Gmin to each diagonal element.  Because Gmin is
+ *  added to the current diagonal, which may bear little relation to
+ *  what the outside world thinks is a diagonal, and because the
+ *  elements that are diagonals may change after calling spOrderAndFactor,
+ *  use of this routine is not recommended.  It is included here simply
+ *  for compatibility with Spice3.
+ */
+
+
+static void
+LoadGmin_CSC (double **diag, int n, double Gmin)
+{
+    int i ;
+
+    if (Gmin != 0.0)
+        for (i = 0 ; i < n ; i++)
+            if (diag [i] != NULL)
+                *(diag [i]) += Gmin ;
+}
+
+static void
+LoadGmin (SMPmatrix *eMatrix, double Gmin)
+{
+    MatrixPtr Matrix = eMatrix->SPmatrix ;
+    int I ;
+    ArrayOfElementPtrs Diag ;
+    ElementPtr diag ;
+
+    /* Begin `LoadGmin'. */
+    assert (IS_SPARSE (Matrix)) ;
+
+    if (Gmin != 0.0) {
+	Diag = Matrix->Diag ;
+	for (I = Matrix->Size ; I > 0 ; I--)
+        {
+	    if ((diag = Diag [I]) != NULL)
+		diag->Real += Gmin ;
+	}
+    }
+    return ;
+}
+
+
+
+
+/*
+ *  FIND ELEMENT
+ *
+ *  This routine finds an element in the matrix by row and column number.
+ *  If the element exists, a pointer to it is returned.  If not, then NULL
+ *  is returned unless the CreateIfMissing flag is TRUE, in which case a
+ *  pointer to the new element is returned.
+ */
+
+SMPelement *
+SMPfindElt (SMPmatrix *eMatrix, int Row, int Col, int CreateIfMissing)
+{
+    MatrixPtr Matrix = eMatrix->SPmatrix ;
+    ElementPtr Element ;
+
+    /* Begin `SMPfindElt'. */
+    assert (IS_SPARSE (Matrix)) ;
+    Row = Matrix->ExtToIntRowMap [Row] ;
+    Col = Matrix->ExtToIntColMap [Col] ;
+    Element = Matrix->FirstInCol [Col] ;
+    Element = spcFindElementInCol (Matrix, &Element, Row, Col, CreateIfMissing) ;
+    return (SMPelement *)Element ;
+}
+
+/* XXX The following should probably be implemented in spUtils */
+
+/*
+ * SMPcZeroCol()
+ */
+int
+SMPcZeroCol (SMPmatrix *eMatrix, int Col)
+{
+    MatrixPtr Matrix = eMatrix->SPmatrix ;
+    ElementPtr	Element ;
+
+    Col = Matrix->ExtToIntColMap [Col] ;
+
+    for (Element = Matrix->FirstInCol [Col] ; Element != NULL ; Element = Element->NextInCol)
+    {
+	Element->Real = 0.0 ;
+	Element->Imag = 0.0 ;
+    }
+
+    return spError (Matrix) ;
+}
+
+/*
+ * SMPcAddCol()
+ */
+int
+SMPcAddCol (SMPmatrix *eMatrix, int Accum_Col, int Addend_Col)
+{
+    MatrixPtr Matrix = eMatrix->SPmatrix ;
+    ElementPtr	Accum, Addend, *Prev ;
+
+    Accum_Col = Matrix->ExtToIntColMap [Accum_Col] ;
+    Addend_Col = Matrix->ExtToIntColMap [Addend_Col] ;
+
+    Addend = Matrix->FirstInCol [Addend_Col] ;
+    Prev = &Matrix->FirstInCol [Accum_Col] ;
+    Accum = *Prev;
+
+    while (Addend != NULL)
+    {
+	while (Accum && Accum->Row < Addend->Row)
+        {
+	    Prev = &Accum->NextInCol ;
+	    Accum = *Prev ;
+	}
+	if (!Accum || Accum->Row > Addend->Row)
+        {
+	    Accum = spcCreateElement (Matrix, Addend->Row, Accum_Col, Prev, 0) ;
+	}
+	Accum->Real += Addend->Real ;
+	Accum->Imag += Addend->Imag ;
+	Addend = Addend->NextInCol ;
+    }
+
+    return spError (Matrix) ;
+}
+
+/*
+ * SMPzeroRow()
+ */
+int
+SMPzeroRow (SMPmatrix *eMatrix, int Row)
+{
+    MatrixPtr Matrix = eMatrix->SPmatrix ;
+    ElementPtr	Element ;
+
+    Row = Matrix->ExtToIntColMap [Row] ;
+
+    if (Matrix->RowsLinked == NO)
+	spcLinkRows (Matrix) ;
+
+    if (Matrix->PreviousMatrixWasComplex || Matrix->Complex)
+    {
+	for (Element = Matrix->FirstInRow[Row] ; Element != NULL; Element = Element->NextInRow)
+	{
+	    Element->Real = 0.0 ;
+	    Element->Imag = 0.0 ;
+	}
+    } else {
+	for (Element = Matrix->FirstInRow [Row] ; Element != NULL ; Element = Element->NextInRow)
+	{
+	    Element->Real = 0.0 ;
+	}
+    }
+
+    return spError (Matrix) ;
+}
diff --git a/src/maths/sparse/spCSC.c b/src/maths/sparse/spCSC.c
new file mode 100644
index 000000000..2d7cd704e
--- /dev/null
+++ b/src/maths/sparse/spCSC.c
@@ -0,0 +1,145 @@
+/* Sparse Matrix to CSC Matrix Conversion Routines
+ * Including Dump Routines
+ *
+ * Author: Francesco Lannutti 2011-2012
+ *
+ * Instructions:
+ * spMatrix_CSC_dump and spRHS_CSC_dump are the dump routines;
+ * insert them in a point in your code after that the Sparse Matrix
+ * is filled in to dump the whole matrix in the CSC format.
+ * To solve correctly the resulting CSC linear system, it's crucial
+ * to perform another inversion of the Solution Vector following this code:
+ *
+ * pExtOrder = IntToExtColMap [n] ;
+ * for (i = n - 1 ; i >= 0 ; i--)
+ *     RHS [*(pExtOrder--)] = Intermediate [i] ;
+ */
+
+/* Includes */
+#include "ngspice/spmatrix.h"
+#include "spdefs.h"
+
+/* Body */
+int
+WriteCol_original (MatrixPtr Matrix, int Col, spREAL *CSC_Element, spREAL *CSC_Element_Complex, int *CSC_Row, BindElement *BindSparseCSC, spREAL **diag)
+{
+    int i ;
+    ElementPtr current ;
+
+    i = 0 ;
+    current = Matrix->FirstInCol [Col] ;
+
+    while (current != NULL) {
+        BindSparseCSC [i].Sparse = (double *)current ;
+        BindSparseCSC [i].CSC = &(CSC_Element [i]) ;
+        BindSparseCSC [i].CSC_Complex = &(CSC_Element_Complex [2 * i]) ;
+        CSC_Row [i] = (current->Row) - 1 ;
+        if (CSC_Row [i] == Col - 1)
+            diag [0] = &(CSC_Element [i]) ;
+        i++ ;
+        current = current->NextInCol ;
+    }
+
+    return i ;
+}
+
+int
+WriteCol_original_dump (MatrixPtr Matrix, int Col, spREAL *CSC_Element, int *CSC_Row)
+{
+    int i ;
+    ElementPtr current ;
+    i = 0 ;
+    current = Matrix->FirstInCol [Col] ;
+
+    while (current != NULL) {
+        CSC_Element [i] = current->Real ;
+        CSC_Row [i] = (current->Row) - 1 ;
+        i++ ;
+        current = current->NextInCol ;
+    }
+
+    return i ;
+}
+
+void
+spMatrix_CSC (MatrixPtr Matrix, int *Ap, int *Ai, double *Ax, double *Ax_Complex, int n, BindElement *BindSparseCSC, double **diag)
+{
+    int offset, i ;
+
+    offset = 0 ;
+    Ap[0] = offset ;
+    for (i = 1 ; i <= n ; i++) {
+        offset += WriteCol_original (Matrix, i, (spREAL *)(Ax + offset), (spREAL *)(Ax_Complex + 2 * offset),
+                                     (int *)(Ai + offset), BindSparseCSC + offset, (spREAL **)(diag + (i - 1))) ;
+
+        Ap[i] = offset ;
+    }
+}
+
+void
+spMatrix_CSC_dump (MatrixPtr Matrix, char *CSC_output)
+{
+    FILE *output ;
+    int offset, i, j, *Ap, *Ai, n, nz ;
+    double *Ax ;
+
+    n = spGetSize (Matrix, 1) ;
+    nz = Matrix->Elements ;
+    Ap = (int *) SP_MALLOC (int, n + 1) ;
+    Ai = (int *) SP_MALLOC (int, nz) ;
+    Ax = (double *) SP_MALLOC (double, nz) ;
+
+    offset = 0 ;
+    Ap[0] = offset ;
+    for (i = 1 ; i <= n ; i++) {
+        offset += WriteCol_original_dump (Matrix, i, (spREAL *)(Ax + offset), (int *)(Ai + offset)) ;
+        Ap[i] = offset ;
+    }
+
+    output = fopen (CSC_output, "w") ;
+    fprintf (output, "%%%%MatrixMarket matrix coordinate real general\n") ;
+    fprintf (output, "%%-------------------------------------------------------------------------------\n") ;
+    fprintf (output, "%% Transient Matrix Dump\n%% Family: ISCAS Circuit\n") ;
+    fprintf (output, "%%-------------------------------------------------------------------------------\n") ;
+    fprintf (output, "%d %d %d\n", n, n, offset) ;
+    for (i = 0 ; i < n ; i++)
+        for (j = Ap [i] ; j < Ap [i + 1] ; j++)
+            fprintf (output, "%d %d %-.9g\n", Ai [j] + 1, i + 1, Ax [j]) ;
+    fclose (output) ;
+
+    output = fopen ("IntToExtColMap.txt", "w") ;
+    for (i = 1 ; i <= n ; i++)
+        fprintf (output, "%d\n", Matrix->IntToExtColMap [i]) ;
+    fclose (output) ;
+
+    SP_FREE (Ap) ;
+    SP_FREE (Ai) ;
+    SP_FREE (Ax) ;
+}
+
+void
+spRHS_CSC_dump (RealNumber *RHS, char *CSC_output_b, MatrixPtr Matrix)
+{
+    FILE *output ;
+    int i, n, *pExtOrder ;
+    double *Intermediate ;
+
+    n = spGetSize (Matrix, 1) ;
+    Intermediate = (double *) SP_MALLOC (double, n) ;
+
+    pExtOrder = &Matrix->IntToExtRowMap [n] ;
+    for (i = n - 1 ; i >= 0 ; i--)
+        Intermediate [i] = RHS [*(pExtOrder--)] ;
+
+    output = fopen (CSC_output_b, "w") ;
+    fprintf (output, "%%%%MatrixMarket matrix array real general\n") ;
+    fprintf (output, "%%-------------------------------------------------------------------------------\n") ;
+    fprintf (output, "%% Transient RHS Vector Dump\n%% Family: ISCAS Circuit\n") ;
+    fprintf (output, "%%-------------------------------------------------------------------------------\n") ;
+    fprintf (output, "%d %d\n", n, 1) ;
+    for (i = 1 ; i < n + 1 ; i++)
+        fprintf (output, "%-.9g\n", Intermediate [i]) ;
+    fclose (output) ;
+
+    SP_FREE (Intermediate) ;
+}