diff --git a/configure.ac b/configure.ac index dbe9f2b0e..92ebfc65f 100644 --- a/configure.ac +++ b/configure.ac @@ -198,6 +198,17 @@ if test "x$with_editline" = xyes; then fi fi +# --enable-klu: Use KLU linear systems solver +AC_ARG_ENABLE([klu], + [AS_HELP_STRING([--enable-klu], [Use KLU linear systems solver])]) + +# Add KLU solver to ngspice +if test "x$enable_klu" = xyes; then + AC_DEFINE(KLU, [], [Define if we want KLU linear systems solver]) + AC_MSG_WARN([KLU solver enabled]) +fi +AM_CONDITIONAL([KLU_WANTED], [test "x$enable_klu" = xyes]) + # Enable maintainer commands only if requested AM_MAINTAINER_MODE([enable]) @@ -1152,6 +1163,7 @@ AC_CONFIG_FILES([Makefile src/maths/deriv/Makefile src/maths/poly/Makefile src/maths/sparse/Makefile + src/maths/KLU/Makefile src/misc/Makefile src/xspice/Makefile src/xspice/cm/Makefile diff --git a/src/Makefile.am b/src/Makefile.am index 1bdafe5e4..4b55a7aeb 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -177,7 +177,13 @@ ngspice_LDADD += \ maths/misc/libmathmisc.la \ maths/fft/libmathfft.la \ maths/poly/libpoly.la \ - maths/ni/libni.la \ + maths/ni/libni.la + +if KLU_WANTED +ngspice_LDADD += maths/KLU/libKLU.la +endif + +ngspice_LDADD += \ maths/sparse/libsparse.la \ misc/libmisc.la @@ -280,6 +286,10 @@ ngmultidec_LDADD = \ maths/sparse/libsparse.la \ misc/libmisc.la +if KLU_WANTED +ngmultidec_LDADD += maths/KLU/libKLU.la +endif + ## ngmakeidx: @@ -418,6 +428,10 @@ libspice_la_LIBADD += \ ciderlib/support/libcidersuprt.la endif +if KLU_WANTED +libspice_la_LIBADD += maths/KLU/libKLU.la +endif + libspice_la_LIBADD += \ maths/deriv/libderiv.la \ maths/cmaths/libcmaths.la \ diff --git a/src/frontend/misccoms.c b/src/frontend/misccoms.c index d41ff6774..40360a294 100644 --- a/src/frontend/misccoms.c +++ b/src/frontend/misccoms.c @@ -181,6 +181,14 @@ com_version(wordlist *wl) fprintf(cp_out, "******\n" "** %s-%s : %s\n" + +#ifdef KLU + "** Compiled with KLU Direct Linear Solver\n" + "** (ASRC, CPL, NDEV and URC models not supported yet)\n" +#else + "** Compiled with Sparse Direct Linear Solver\n" +#endif + "** The U. C. Berkeley CAD Group\n" "** Copyright 1985-1994, Regents of the University of California.\n" "** %s\n", @@ -213,6 +221,14 @@ com_version(wordlist *wl) fprintf(cp_out, "******\n" "** %s-%s : %s\n" + +#ifdef KLU + "** Compiled with KLU Direct Linear Solver\n" + "** (ASRC, CPL, NDEV and URC models not supported yet)\n" +#else + "** Compiled with Sparse Direct Linear Solver\n" +#endif + "** The U. C. Berkeley CAD Group\n" "** Copyright 1985-1994, Regents of the University of California.\n" "** %s\n", diff --git a/src/include/ngspice/cktdefs.h b/src/include/ngspice/cktdefs.h index f7777ae1d..5668d67f1 100644 --- a/src/include/ngspice/cktdefs.h +++ b/src/include/ngspice/cktdefs.h @@ -299,6 +299,10 @@ struct CKTcircuit { NGHASHPTR MODnameHash; GENinstance *noise_input; /* identify the input vsrc/isrc during noise analysis */ + +#ifdef KLU + unsigned int CKTkluMODE:1; +#endif }; diff --git a/src/include/ngspice/devdefs.h b/src/include/ngspice/devdefs.h index 0e338e348..be65de693 100644 --- a/src/include/ngspice/devdefs.h +++ b/src/include/ngspice/devdefs.h @@ -114,6 +114,15 @@ typedef struct SPICEdev { int *DEVinstSize; /* size of an instance */ int *DEVmodSize; /* size of a model */ +#ifdef KLU + int (*DEVbindCSC)(GENmodel *, CKTcircuit *); + /* routine to convert Sparse linked list to Real CSC array */ + int (*DEVbindCSCComplex)(GENmodel *, CKTcircuit *); + /* routine to convert Real CSC array to Complex CSC array */ + int (*DEVbindCSCComplexToReal)(GENmodel *, CKTcircuit *); + /* routine to convert Complex CSC array to Real CSC array */ +#endif + } SPICEdev; /* instance of structure for each possible type of device */ diff --git a/src/include/ngspice/optdefs.h b/src/include/ngspice/optdefs.h index bca139d06..06f127ebb 100644 --- a/src/include/ngspice/optdefs.h +++ b/src/include/ngspice/optdefs.h @@ -123,6 +123,9 @@ typedef struct { #define OPT_INDVERBOSITY 70 #define OPT_EPSMIN 71 +#ifdef KLU +#define OPT_SPARSE 72 +#endif #ifdef XSPICE /* gtri - begin - wbk - add new options */ diff --git a/src/include/ngspice/smpdefs.h b/src/include/ngspice/smpdefs.h index eda52847f..df001903b 100644 --- a/src/include/ngspice/smpdefs.h +++ b/src/include/ngspice/smpdefs.h @@ -1,8 +1,12 @@ #ifndef ngspice_SMPDEFS_H #define ngspice_SMPDEFS_H +/* Typedef removed by Francesco Lannutti (2012-02) to create the new SMPmatrix structure */ +/* typedef struct MatrixFrame SMPmatrix; -typedef struct MatrixElement SMPelement; +*/ +typedef struct MatrixFrame MatrixFrame; +typedef struct MatrixElement *SMPelement; /********** Copyright 1990 Regents of the University of California. All rights reserved. @@ -14,6 +18,44 @@ Modified: 2000 AlansFixes #include #include "ngspice/complex.h" +#ifdef KLU +#include "ngspice/klu.h" +#include "ngspice/spmatrix.h" +#endif + +/* SMPmatrix structure - Francesco Lannutti (2012-02) */ +typedef struct sSMPmatrix { + MatrixFrame *SPmatrix ; /* pointer to sparse matrix */ + +#ifdef KLU + klu_common *CKTkluCommon ; /* KLU common object */ + klu_symbolic *CKTkluSymbolic ; /* KLU symbolic object */ + klu_numeric *CKTkluNumeric ; /* KLU numeric object */ + int *CKTkluAp ; /* KLU column pointer */ + int *CKTkluAi ; /* KLU row pointer */ + double *CKTkluAx ; /* KLU Real Elements */ + double *CKTkluAx_Complex ; /* KLU Complex Elements */ + int CKTkluMatrixIsComplex ; /* KLU Matrix Is Complex Flag */ + #define CKTkluMatrixReal 0 /* KLU Matrix Real definition */ + #define CKTkluMatrixComplex 1 /* KLU Matrix Complex definition */ + double *CKTkluIntermediate ; /* KLU RHS Intermediate for Solve Real Step */ + double *CKTkluIntermediate_Complex ; /* KLU iRHS Intermediate for Solve Complex Step */ + BindElement *CKTbindStruct ; /* KLU - Sparse Binding Structure */ + double **CKTdiag_CSC ; /* KLU pointer to diagonal element to perform Gmin */ + int CKTkluN ; /* KLU N, copied */ + int CKTklunz ; /* KLU nz, copied for AC Analysis */ + 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 */ +#endif + +} SMPmatrix ; + + +#ifdef KLU +void SMPmatrix_CSC (SMPmatrix *) ; +void SMPnnz (SMPmatrix *) ; +#endif int SMPaddElt( SMPmatrix *, int , int , double ); double * SMPmakeElt( SMPmatrix * , int , int ); void SMPcClear( SMPmatrix *); @@ -27,7 +69,7 @@ void SMPcaSolve(SMPmatrix *Matrix, double RHS[], double iRHS[], void SMPcSolve( SMPmatrix *, double [], double [], double [], double []); void SMPsolve( SMPmatrix *, double [], double []); int SMPmatSize( SMPmatrix *); -int SMPnewMatrix( SMPmatrix **, int ); +int SMPnewMatrix( SMPmatrix *, int ); void SMPdestroy( SMPmatrix *); int SMPpreOrder( SMPmatrix *); void SMPprint( SMPmatrix * , char *); diff --git a/src/include/ngspice/spmatrix.h b/src/include/ngspice/spmatrix.h index 33780a56e..65c5488a8 100644 --- a/src/include/ngspice/spmatrix.h +++ b/src/include/ngspice/spmatrix.h @@ -35,6 +35,8 @@ */ +/* Francesco Lannutti 2012-09 - NGSPICE Configuration File Inclusion */ +#include "ngspice/config.h" #ifndef spOKAY @@ -294,4 +296,20 @@ extern void spMultTransposed(MatrixPtr,spREAL*,spREAL*,spREAL*,spREAL*); extern void spSolve( MatrixPtr, spREAL*, spREAL*, spREAL*, spREAL* ); extern void spSolveTransposed(MatrixPtr,spREAL*,spREAL*,spREAL*,spREAL*); +/* Francesco Lannutti - CSC Data Structure Conversion */ +#ifdef KLU +typedef struct sBindElement { + double *Sparse ; + double *CSC ; + double *CSC_Complex ; +} BindElement ; + +extern int WriteCol_original (MatrixPtr, int, spREAL *, spREAL *, int *, BindElement *, spREAL **) ; +extern int WriteCol_original_dump (MatrixPtr, int, spREAL *, int *) ; +extern void spMatrix_CSC (MatrixPtr, int *, int *, double *, double *, int, BindElement *, double **) ; +extern void spMatrix_CSC_dump (MatrixPtr, char *) ; +extern void spRHS_CSC_dump (spREAL *, char *, MatrixPtr) ; +#endif +/* ------------------------------------------------------ */ + #endif /* spOKAY */ diff --git a/src/include/ngspice/tskdefs.h b/src/include/ngspice/tskdefs.h index 1bba643ef..67744f7ed 100644 --- a/src/include/ngspice/tskdefs.h +++ b/src/include/ngspice/tskdefs.h @@ -69,6 +69,11 @@ struct TSKtask { double TSKrelDv; /* rel limit for iter-iter voltage change */ unsigned int TSKnoopac:1; /* flag for no OP calculation before AC */ double TSKepsmin; /* minimum value for log */ + +#ifdef KLU + unsigned int TSKkluMODE:1; +#endif + }; #endif diff --git a/src/maths/Makefile.am b/src/maths/Makefile.am index 9bc84a2cf..862326353 100644 --- a/src/maths/Makefile.am +++ b/src/maths/Makefile.am @@ -1,6 +1,10 @@ ## Process this file with automake SUBDIRS = cmaths ni sparse poly deriv misc fft -DIST_SUBDIRS = cmaths ni sparse poly deriv misc fft +DIST_SUBDIRS = cmaths ni sparse poly deriv misc fft KLU + +if KLU_WANTED +SUBDIRS += KLU +endif MAINTAINERCLEANFILES = Makefile.in diff --git a/src/maths/ni/niinit.c b/src/maths/ni/niinit.c index 1df37f66c..301269de5 100644 --- a/src/maths/ni/niinit.c +++ b/src/maths/ni/niinit.c @@ -16,6 +16,9 @@ Author: 1985 Thomas L. Quarles #include "ngspice/sperror.h" #include "ngspice/smpdefs.h" +#ifdef KLU +#include "ngspice/klu.h" +#endif int NIinit(CKTcircuit *ckt) @@ -24,6 +27,29 @@ NIinit(CKTcircuit *ckt) /* a concession to Ken Kundert's sparse matrix package - SMP doesn't need this*/ int Error; #endif /* SPARSE */ + +/* Allocation of the new SMPmatrix structure - Francesco Lannutti (2012-02) */ + ckt->CKTmatrix = TMALLOC (SMPmatrix, 1) ; + +#ifdef KLU + ckt->CKTmatrix->CKTkluCommon = TMALLOC (klu_common, 1) ; + ckt->CKTmatrix->CKTkluSymbolic = NULL ; + ckt->CKTmatrix->CKTkluNumeric = NULL ; + ckt->CKTmatrix->CKTkluAp = NULL ; + ckt->CKTmatrix->CKTkluAi = NULL ; + ckt->CKTmatrix->CKTkluAx = NULL ; + ckt->CKTmatrix->CKTkluMatrixIsComplex = CKTkluMatrixReal ; + ckt->CKTmatrix->CKTkluIntermediate = NULL ; + ckt->CKTmatrix->CKTkluIntermediate_Complex = NULL ; + ckt->CKTmatrix->CKTbindStruct = NULL ; + ckt->CKTmatrix->CKTdiag_CSC = NULL ; + ckt->CKTmatrix->CKTkluN = 0 ; + ckt->CKTmatrix->CKTklunz = 0 ; + ckt->CKTmatrix->CKTkluMODE = ckt->CKTkluMODE ; /* TO BE SUBSTITUTED WITH THE HEURISTICS */ + + klu_defaults (ckt->CKTmatrix->CKTkluCommon) ; +#endif + ckt->CKTniState = NIUNINITIALIZED; - return SMPnewMatrix(&(ckt->CKTmatrix), 0); + return SMPnewMatrix(ckt->CKTmatrix, 0); } diff --git a/src/maths/sparse/Makefile.am b/src/maths/sparse/Makefile.am index ccbdcf969..3ab29a6a4 100644 --- a/src/maths/sparse/Makefile.am +++ b/src/maths/sparse/Makefile.am @@ -2,7 +2,7 @@ noinst_LTLIBRARIES = libsparse.la -libsparse_la_SOURCES = \ +libsparse_la_SOURCES = \ spalloc.c \ spbuild.c \ spconfig.h \ @@ -15,7 +15,9 @@ libsparse_la_SOURCES = \ sputils.c - +if KLU_WANTED +libsparse_la_SOURCES += spCSC.c +endif AM_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include AM_CFLAGS = $(STATIC) diff --git a/src/maths/sparse/spsmp.c b/src/maths/sparse/spsmp.c index d0b476eed..a21df7b1a 100644 --- a/src/maths/sparse/spsmp.c +++ b/src/maths/sparse/spsmp.c @@ -106,7 +106,7 @@ extern double scalbn(double, int); extern double logb(double); #endif -static void LoadGmin(SMPmatrix *Matrix, double Gmin); +static void LoadGmin(SMPmatrix *eMatrix, double Gmin); /* @@ -115,8 +115,8 @@ static void LoadGmin(SMPmatrix *Matrix, double Gmin); int SMPaddElt(SMPmatrix *Matrix, int Row, int Col, double Value) { - *spGetElement( Matrix, Row, Col ) = Value; - return spError( Matrix ); + *spGetElement( Matrix->SPmatrix, Row, Col ) = Value; + return spError( Matrix->SPmatrix ); } /* @@ -125,7 +125,7 @@ SMPaddElt(SMPmatrix *Matrix, int Row, int Col, double Value) double * SMPmakeElt(SMPmatrix *Matrix, int Row, int Col) { - return spGetElement( Matrix, Row, Col ); + return spGetElement( Matrix->SPmatrix, Row, Col ); } /* @@ -134,7 +134,7 @@ SMPmakeElt(SMPmatrix *Matrix, int Row, int Col) void SMPcClear(SMPmatrix *Matrix) { - spClear( Matrix ); + spClear( Matrix->SPmatrix ); } /* @@ -143,7 +143,7 @@ SMPcClear(SMPmatrix *Matrix) void SMPclear(SMPmatrix *Matrix) { - spClear( Matrix ); + spClear( Matrix->SPmatrix ); } #define NG_IGNORE(x) (void)x @@ -157,8 +157,8 @@ SMPcLUfac(SMPmatrix *Matrix, double PivTol) { NG_IGNORE(PivTol); - spSetComplex( Matrix ); - return spFactor( Matrix ); + spSetComplex( Matrix->SPmatrix ); + return spFactor( Matrix->SPmatrix ); } /* @@ -169,9 +169,9 @@ int SMPluFac(SMPmatrix *Matrix, double PivTol, double Gmin) { NG_IGNORE(PivTol); - spSetReal( Matrix ); + spSetReal( Matrix->SPmatrix ); LoadGmin( Matrix, Gmin ); - return spFactor( Matrix ); + return spFactor( Matrix->SPmatrix ); } /* @@ -182,8 +182,8 @@ SMPcReorder(SMPmatrix *Matrix, double PivTol, double PivRel, int *NumSwaps) { *NumSwaps = 1; - spSetComplex( Matrix ); - return spOrderAndFactor( Matrix, NULL, + spSetComplex( Matrix->SPmatrix ); + return spOrderAndFactor( Matrix->SPmatrix, NULL, PivRel, PivTol, YES ); } @@ -193,9 +193,9 @@ SMPcReorder(SMPmatrix *Matrix, double PivTol, double PivRel, int SMPreorder(SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin) { - spSetReal( Matrix ); + spSetReal( Matrix->SPmatrix ); LoadGmin( Matrix, Gmin ); - return spOrderAndFactor( Matrix, NULL, + return spOrderAndFactor( Matrix->SPmatrix, NULL, PivRel, PivTol, YES ); } @@ -209,7 +209,7 @@ SMPcaSolve(SMPmatrix *Matrix, double RHS[], double iRHS[], NG_IGNORE(iSpare); NG_IGNORE(Spare); - spSolveTransposed( Matrix, RHS, RHS, iRHS, iRHS ); + spSolveTransposed( Matrix->SPmatrix, RHS, RHS, iRHS, iRHS ); } /* @@ -222,7 +222,7 @@ SMPcSolve(SMPmatrix *Matrix, double RHS[], double iRHS[], NG_IGNORE(iSpare); NG_IGNORE(Spare); - spSolve( Matrix, RHS, RHS, iRHS, iRHS ); + spSolve( Matrix->SPmatrix, RHS, RHS, iRHS, iRHS ); } /* @@ -233,7 +233,7 @@ SMPsolve(SMPmatrix *Matrix, double RHS[], double Spare[]) { NG_IGNORE(Spare); - spSolve( Matrix, RHS, RHS, NULL, NULL ); + spSolve( Matrix->SPmatrix, RHS, RHS, NULL, NULL ); } /* @@ -242,17 +242,17 @@ SMPsolve(SMPmatrix *Matrix, double RHS[], double Spare[]) int SMPmatSize(SMPmatrix *Matrix) { - return spGetSize( Matrix, 1 ); + return spGetSize( Matrix->SPmatrix, 1 ); } /* * SMPnewMatrix() */ int -SMPnewMatrix(SMPmatrix **pMatrix, int size) +SMPnewMatrix(SMPmatrix *Matrix, int size) { int Error; - *pMatrix = spCreate( size, 1, &Error ); + Matrix->SPmatrix = spCreate( size, 1, &Error ); return Error; } @@ -262,7 +262,7 @@ SMPnewMatrix(SMPmatrix **pMatrix, int size) void SMPdestroy(SMPmatrix *Matrix) { - spDestroy( Matrix ); + spDestroy( Matrix->SPmatrix ); } /* @@ -271,18 +271,18 @@ SMPdestroy(SMPmatrix *Matrix) int SMPpreOrder(SMPmatrix *Matrix) { - spMNA_Preorder( Matrix ); - return spError( Matrix ); + spMNA_Preorder( Matrix->SPmatrix ); + return spError( Matrix->SPmatrix ); } /* - * SMPprint() + * SMPprintRHS() */ void SMPprintRHS(SMPmatrix *Matrix, char *Filename, RealVector RHS, RealVector iRHS) { - spFileVector( Matrix, Filename, RHS, iRHS ); + spFileVector( Matrix->SPmatrix, Filename, RHS, iRHS ); } /* @@ -293,9 +293,9 @@ void SMPprint(SMPmatrix *Matrix, char *Filename) { if (Filename) - spFileMatrix(Matrix, Filename, "Circuit Matrix", 0, 1, 1 ); + spFileMatrix(Matrix->SPmatrix, Filename, "Circuit Matrix", 0, 1, 1 ); else - spPrint( Matrix, 0, 1, 1 ); + spPrint( Matrix->SPmatrix, 0, 1, 1 ); } /* @@ -304,7 +304,7 @@ SMPprint(SMPmatrix *Matrix, char *Filename) void SMPgetError(SMPmatrix *Matrix, int *Col, int *Row) { - spWhereSingular( Matrix, Row, Col ); + spWhereSingular( Matrix->SPmatrix, Row, Col ); } /* @@ -314,9 +314,9 @@ SMPgetError(SMPmatrix *Matrix, int *Col, int *Row) int SMPcProdDiag(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent) { - spDeterminant( Matrix, pExponent, &(pMantissa->real), + spDeterminant ( Matrix->SPmatrix, pExponent, &(pMantissa->real), &(pMantissa->imag) ); - return spError( Matrix ); + return spError( Matrix->SPmatrix ); } /* @@ -328,7 +328,7 @@ SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent) double re, im, x, y, z; int p; - spDeterminant( Matrix, &p, &re, &im); + spDeterminant( Matrix->SPmatrix, &p, &re, &im); #ifndef M_LN2 #define M_LN2 0.69314718055994530942 @@ -398,7 +398,7 @@ SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent) printf("Determinant 10->2: (%20g,%20g)^%d\n", pMantissa->real, pMantissa->imag, *pExponent); #endif - return spError( Matrix ); + return spError( Matrix->SPmatrix ); } @@ -420,8 +420,9 @@ SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent) */ static void -LoadGmin(SMPmatrix *Matrix, double Gmin) +LoadGmin(SMPmatrix *eMatrix, double Gmin) { + MatrixPtr Matrix = eMatrix->SPmatrix; int I; ArrayOfElementPtrs Diag; ElementPtr diag; @@ -452,8 +453,9 @@ LoadGmin(SMPmatrix *Matrix, double Gmin) */ SMPelement * -SMPfindElt(SMPmatrix *Matrix, int Row, int Col, int CreateIfMissing) +SMPfindElt(SMPmatrix *eMatrix, int Row, int Col, int CreateIfMissing) { + MatrixPtr Matrix = eMatrix->SPmatrix; ElementPtr Element; /* Begin `SMPfindElt'. */ @@ -462,7 +464,7 @@ SMPfindElt(SMPmatrix *Matrix, int Row, int Col, int CreateIfMissing) Col = Matrix->ExtToIntColMap[Col]; Element = Matrix->FirstInCol[Col]; Element = spcFindElementInCol(Matrix, &Element, Row, Col, CreateIfMissing); - return Element; + return (SMPelement *)Element; } /* XXX The following should probably be implemented in spUtils */ @@ -471,8 +473,9 @@ SMPfindElt(SMPmatrix *Matrix, int Row, int Col, int CreateIfMissing) * SMPcZeroCol() */ int -SMPcZeroCol(SMPmatrix *Matrix, int Col) +SMPcZeroCol(SMPmatrix *eMatrix, int Col) { + MatrixPtr Matrix = eMatrix->SPmatrix; ElementPtr Element; Col = Matrix->ExtToIntColMap[Col]; @@ -492,8 +495,9 @@ SMPcZeroCol(SMPmatrix *Matrix, int Col) * SMPcAddCol() */ int -SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col) +SMPcAddCol(SMPmatrix *eMatrix, int Accum_Col, int Addend_Col) { + MatrixPtr Matrix = eMatrix->SPmatrix; ElementPtr Accum, Addend, *Prev; Accum_Col = Matrix->ExtToIntColMap[Accum_Col]; @@ -523,8 +527,9 @@ SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col) * SMPzeroRow() */ int -SMPzeroRow(SMPmatrix *Matrix, int Row) +SMPzeroRow(SMPmatrix *eMatrix, int Row) { + MatrixPtr Matrix = eMatrix->SPmatrix; ElementPtr Element; Row = Matrix->ExtToIntColMap[Row]; diff --git a/src/spicelib/analysis/acan.c b/src/spicelib/analysis/acan.c index 1bef51747..87c8e39f9 100644 --- a/src/spicelib/analysis/acan.c +++ b/src/spicelib/analysis/acan.c @@ -251,6 +251,23 @@ ACan(CKTcircuit *ckt, int restart) ckt->CKTcurrentAnalysis = DOING_AC; +#ifdef KLU + int i ; + + if (ckt->CKTmatrix->CKTkluMODE) + { + /* Conversion from Real Matrix to Complex Matrix */ + if (!ckt->CKTmatrix->CKTkluMatrixIsComplex) + { + for (i = 0 ; i < DEVmaxnum ; i++) + if (DEVices [i] && DEVices [i]->DEVbindCSCComplex && ckt->CKThead [i]) + DEVices [i]->DEVbindCSCComplex (ckt->CKThead [i], ckt) ; + + ckt->CKTmatrix->CKTkluMatrixIsComplex = CKTkluMatrixComplex ; + } + } +#endif + /* main loop through all scheduled frequencies */ while (freq <= job->ACstopFreq + freqTol) { if(SPfrontEnd->IFpauseTest()) { @@ -262,6 +279,19 @@ ACan(CKTcircuit *ckt, int restart) /* Update opertating point, if variable 'hertz' is given */ if (ckt->CKTvarHertz) { + +#ifdef KLU + if (ckt->CKTmatrix->CKTkluMODE) + { + /* Conversion from Complex Matrix to Real Matrix */ + for (i = 0 ; i < DEVmaxnum ; i++) + if (DEVices [i] && DEVices [i]->DEVbindCSCComplexToReal && ckt->CKThead [i]) + DEVices [i]->DEVbindCSCComplexToReal (ckt->CKThead [i], ckt) ; + + ckt->CKTmatrix->CKTkluMatrixIsComplex = CKTkluMatrixReal ; + } +#endif + #ifdef XSPICE /* Call EVTop if event-driven instances exist */ @@ -290,6 +320,19 @@ ACan(CKTcircuit *ckt, int restart) ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG; error = CKTload(ckt); if(error) return(error); + +#ifdef KLU + if (ckt->CKTmatrix->CKTkluMODE) + { + /* Conversion from Real Matrix to Complex Matrix */ + for (i = 0 ; i < DEVmaxnum ; i++) + if (DEVices [i] && DEVices [i]->DEVbindCSCComplex && ckt->CKThead [i]) + DEVices [i]->DEVbindCSCComplex (ckt->CKThead [i], ckt) ; + + ckt->CKTmatrix->CKTkluMatrixIsComplex = CKTkluMatrixComplex ; + } +#endif + } ckt->CKTmode = (ckt->CKTmode&MODEUIC) | MODEAC; @@ -389,6 +432,19 @@ endsweep: SPfrontEnd->OUTendPlot (acPlot); acPlot = NULL; UPDATE_STATS(0); + +#ifdef KLU + if (ckt->CKTmatrix->CKTkluMODE) + { + /* Conversion from Complex Matrix to Real Matrix */ + for (i = 0 ; i < DEVmaxnum ; i++) + if (DEVices [i] && DEVices [i]->DEVbindCSCComplexToReal && ckt->CKThead [i]) + DEVices [i]->DEVbindCSCComplexToReal (ckt->CKThead [i], ckt) ; + + ckt->CKTmatrix->CKTkluMatrixIsComplex = CKTkluMatrixReal ; + } +#endif + return(0); } diff --git a/src/spicelib/analysis/cktacct.c b/src/spicelib/analysis/cktacct.c index 63b0a1063..161b68dc4 100644 --- a/src/spicelib/analysis/cktacct.c +++ b/src/spicelib/analysis/cktacct.c @@ -18,7 +18,14 @@ Author: 1985 Thomas L. Quarles #include "ngspice/cktdefs.h" #include "ngspice/spmatrix.h" +#ifdef KLU +#include "ngspice/klu.h" +#endif +/* Francesco Lannutti + * If ACCT is called before NIinit, the new SMPmatrix structure is not allocated + * so the control must be performed on the CKTmatrix pointer to the SMPmatrix structure +*/ /* ARGSUSED */ int @@ -33,21 +40,35 @@ CKTacct(CKTcircuit *ckt, JOB *anal, int which, IFvalue *val) break; case OPT_ORIGNZ: if ( ckt->CKTmatrix != NULL ) { - val->iValue = spOriginalCount(ckt->CKTmatrix); + val->iValue = spOriginalCount(ckt->CKTmatrix->SPmatrix); } else { val->iValue = 0; } break; case OPT_FILLNZ: if ( ckt->CKTmatrix != NULL ) { - val->iValue = spFillinCount(ckt->CKTmatrix); +#ifdef KLU + if (ckt->CKTmatrix->CKTkluMODE) + val->iValue = ckt->CKTmatrix->CKTkluNumeric->lnz + ckt->CKTmatrix->CKTkluNumeric->unz + ckt->CKTmatrix->CKTkluNumeric->nzoff - ckt->CKTmatrix->CKTklunz ; + else + val->iValue = spFillinCount(ckt->CKTmatrix->SPmatrix); +#else + val->iValue = spFillinCount(ckt->CKTmatrix->SPmatrix); +#endif } else { val->iValue = 0; } break; case OPT_TOTALNZ: if ( ckt->CKTmatrix != NULL ) { - val->iValue = spElementCount(ckt->CKTmatrix); +#ifdef KLU + if (ckt->CKTmatrix->CKTkluMODE) + val->iValue = ckt->CKTmatrix->CKTkluNumeric->lnz + ckt->CKTmatrix->CKTkluNumeric->unz + ckt->CKTmatrix->CKTkluNumeric->nzoff ; + else + val->iValue = spElementCount(ckt->CKTmatrix->SPmatrix); +#else + val->iValue = spElementCount(ckt->CKTmatrix->SPmatrix); +#endif } else { val->iValue = 0; } diff --git a/src/spicelib/analysis/cktdojob.c b/src/spicelib/analysis/cktdojob.c index fc52981bb..63f64786d 100644 --- a/src/spicelib/analysis/cktdojob.c +++ b/src/spicelib/analysis/cktdojob.c @@ -107,6 +107,11 @@ CKTdoJob(CKTcircuit *ckt, int reset, TSKtask *task) ckt->CKTtroubleElt = NULL; ckt->CKTnoopac = task->TSKnoopac && ckt->CKTisLinear; ckt->CKTepsmin = task->TSKepsmin; + +#ifdef KLU + ckt->CKTkluMODE = task->TSKkluMODE; +#endif + #ifdef NEWTRUNC ckt->CKTlteReltol = task->TSKlteReltol; ckt->CKTlteAbstol = task->TSKlteAbstol; diff --git a/src/spicelib/analysis/cktntask.c b/src/spicelib/analysis/cktntask.c index cb56373a3..491a80dad 100644 --- a/src/spicelib/analysis/cktntask.c +++ b/src/spicelib/analysis/cktntask.c @@ -75,6 +75,11 @@ CKTnewTask(CKTcircuit *ckt, TSKtask **taskPtr, IFuid taskName, TSKtask **defPtr) tsk->TSKrelDv = def->TSKrelDv; tsk->TSKnoopac = def->TSKnoopac; tsk->TSKepsmin = def->TSKepsmin; + +#ifdef KLU + tsk->TSKkluMODE = def->TSKkluMODE; +#endif + #ifdef NEWTRUNC tsk->TSKlteReltol = def->TSKlteReltol; tsk->TSKlteAbstol = def->TSKlteAbstol; @@ -132,6 +137,10 @@ CKTnewTask(CKTcircuit *ckt, TSKtask **taskPtr, IFuid taskName, TSKtask **defPtr) tsk->TSKrelDv = 2.0; tsk->TSKepsmin = 1e-28; +#ifdef KLU + tsk->TSKkluMODE = CKTkluON; +#endif + #if (1) /*CDHW*/ } #endif diff --git a/src/spicelib/analysis/cktsens.c b/src/spicelib/analysis/cktsens.c index e55c20d3b..55477cac4 100644 --- a/src/spicelib/analysis/cktsens.c +++ b/src/spicelib/analysis/cktsens.c @@ -71,7 +71,8 @@ int sens_sens(CKTcircuit *ckt, int restart) static double freq; static int nfreqs; static int i; - static SMPmatrix *delta_Y = NULL, *Y; + static SMPmatrix *delta_Y = NULL; + static MatrixFrame *Y; static double step_size; double *E, *iE; IFvalue value, nvalue; @@ -89,7 +90,9 @@ int sens_sens(CKTcircuit *ckt, int restart) int type; double *saved_rhs = NULL, *saved_irhs = NULL; - SMPmatrix *saved_matrix = NULL; + MatrixFrame *saved_matrix = NULL; + + delta_Y = TMALLOC(SMPmatrix, 1); #ifndef notdef #ifdef notdef @@ -139,10 +142,10 @@ int sens_sens(CKTcircuit *ckt, int restart) if (error) return error; - size = SMPmatSize(ckt->CKTmatrix); + size = SMPmatSize(ckt->CKTmatrix->SPmatrix); /* Create the perturbation matrix */ - error = SMPnewMatrix(&delta_Y, size); + error = SMPnewMatrix(delta_Y, size); if (error) return error; @@ -181,7 +184,9 @@ int sens_sens(CKTcircuit *ckt, int restart) sg->ptable[sg->param].keyword); } - SPfrontEnd->IFnewUid (ckt, output_names + k, NULL, namebuf, UID_OTHER, NULL); + SPfrontEnd->IFnewUid (ckt, + output_names + k, NULL, + namebuf, UID_OTHER, NULL); k += 1; } @@ -190,14 +195,16 @@ int sens_sens(CKTcircuit *ckt, int restart) freq_name = NULL; } else { type = IF_COMPLEX; - SPfrontEnd->IFnewUid (ckt, &freq_name, NULL, "frequency", UID_OTHER, NULL); + SPfrontEnd->IFnewUid (ckt, + &freq_name, NULL, + "frequency", UID_OTHER, NULL); } - error = SPfrontEnd->OUTpBeginPlot (ckt, ckt->CKTcurJob, - ckt->CKTcurJob->JOBname, - freq_name, IF_REAL, - num_vars, output_names, type, - &sen_data); + error = SPfrontEnd->OUTpBeginPlot ( + ckt, ckt->CKTcurJob, + ckt->CKTcurJob->JOBname, + freq_name, IF_REAL, + num_vars, output_names, type, &sen_data); if (error) return error; @@ -231,12 +238,12 @@ int sens_sens(CKTcircuit *ckt, int restart) bypass = ckt->CKTbypass; ckt->CKTbypass = 0; - /* CKTop solves into CKTrhs and CKTmatrix, + /* CKTop solves into CKTrhs and CKTmatrix->SPmatrix, * CKTirhs is hopefully zero (fresh allocated ?) */ E = ckt->CKTrhs; iE = ckt->CKTirhs; - Y = ckt->CKTmatrix; + Y = ckt->CKTmatrix->SPmatrix; #ifdef ASDEBUG DEBUG(1) { @@ -279,7 +286,7 @@ int sens_sens(CKTcircuit *ckt, int restart) if (error) return error; - /* XXX ckt->CKTmatrix = Y; */ + /* XXX ckt->CKTmatrix->SPmatrix = Y; */ error = CKTsetup(ckt); if (error) @@ -308,17 +315,17 @@ int sens_sens(CKTcircuit *ckt, int restart) } #endif - /* NIacIter solves into CKTrhsOld, CKTirhsOld and CKTmatrix */ + /* NIacIter solves into CKTrhsOld, CKTirhsOld and CKTmatrix->SPmatrix */ E = ckt->CKTrhsOld; iE = ckt->CKTirhsOld; - Y = ckt->CKTmatrix; + Y = ckt->CKTmatrix->SPmatrix; } /* Use a different vector & matrix */ save_context(ckt->CKTrhs, saved_rhs); save_context(ckt->CKTirhs, saved_irhs); - save_context(ckt->CKTmatrix, saved_matrix); + save_context(ckt->CKTmatrix->SPmatrix, saved_matrix); ckt->CKTrhs = delta_I; ckt->CKTirhs = delta_iI; @@ -351,7 +358,7 @@ int sens_sens(CKTcircuit *ckt, int restart) } #endif - SMPcClear(delta_Y); + SMPcClear(delta_Y->SPmatrix); for (j = 0; j < size; j++) { delta_I[j] = 0.0; @@ -399,7 +406,7 @@ int sens_sens(CKTcircuit *ckt, int restart) #ifdef ASDEBUG DEBUG(2) { printf("Effect of device:\n"); - SMPprint(delta_Y, NULL); + SMPprint(delta_Y->SPmatrix, NULL); printf("LHS:\n"); for (j = 0; j < size; j++) printf("%d: %g, %g\n", j, @@ -423,7 +430,7 @@ int sens_sens(CKTcircuit *ckt, int restart) if (error && error != E_BADPARM) return error; - SMPconstMult(delta_Y, -1.0); + SMPconstMult(delta_Y->SPmatrix, -1.0); for (j = 0; j < size; j++) { delta_I[j] *= -1.0; delta_iI[j] *= -1.0; @@ -432,7 +439,7 @@ int sens_sens(CKTcircuit *ckt, int restart) #ifdef ASDEBUG DEBUG(2) { printf("Effect of negating matrix:\n"); - SMPprint(delta_Y, NULL); + SMPprint(delta_Y->SPmatrix, NULL); for (j = 0; j < size; j++) printf("%d: %g, %g\n", j, delta_I[j], delta_iI[j]); @@ -458,7 +465,7 @@ int sens_sens(CKTcircuit *ckt, int restart) #ifdef ASDEBUG DEBUG(2) { printf("Effect of changing the parameter:\n"); - SMPprint(delta_Y, NULL); + SMPprint(delta_Y->SPmatrix, NULL); for (j = 0; j < size; j++) printf("%d: %g, %g\n", j, delta_I[j], delta_iI[j]); @@ -482,7 +489,7 @@ int sens_sens(CKTcircuit *ckt, int restart) #endif /* delta_Y E */ - SMPmultiply(delta_Y, delta_I_delta_Y, E, + SMPmultiply(delta_Y->SPmatrix, delta_I_delta_Y, E, delta_iI_delta_Y, iE); #ifdef ASDEBUG @@ -574,7 +581,7 @@ int sens_sens(CKTcircuit *ckt, int restart) release_context(ckt->CKTrhs, saved_rhs); release_context(ckt->CKTirhs, saved_irhs); - release_context(ckt->CKTmatrix, saved_matrix); + release_context(ckt->CKTmatrix->SPmatrix, saved_matrix); if (is_dc) nvalue.v.vec.rVec = output_values; @@ -599,9 +606,9 @@ int sens_sens(CKTcircuit *ckt, int restart) release_context(ckt->CKTrhs, saved_rhs); release_context(ckt->CKTirhs, saved_irhs); - release_context(ckt->CKTmatrix, saved_matrix); + release_context(ckt->CKTmatrix->SPmatrix, saved_matrix); - SMPdestroy(delta_Y); + SMPdestroy(delta_Y->SPmatrix); FREE(delta_I); FREE(delta_iI); @@ -625,7 +632,7 @@ int sens_sens(CKTcircuit *ckt, int restart) double inc_freq(double freq, int type, double step_size) { - if (type != SENS_LINEAR) + if (type != LINEAR) freq *= step_size; else freq += step_size; diff --git a/src/spicelib/analysis/cktsetup.c b/src/spicelib/analysis/cktsetup.c index 8804f24c4..3e50f1f05 100644 --- a/src/spicelib/analysis/cktsetup.c +++ b/src/spicelib/analysis/cktsetup.c @@ -30,6 +30,20 @@ int nthreads; return(E_NOMEM);\ } +#ifdef KLU +#include + +static +int +BindCompare (const void *a, const void *b) +{ + BindElement *A, *B ; + A = (BindElement *)a ; + B = (BindElement *)b ; + + return ((int)(A->Sparse - B->Sparse)) ; +} +#endif int CKTsetup(CKTcircuit *ckt) @@ -91,6 +105,53 @@ CKTsetup(CKTcircuit *ckt) if(error) return(error); } } + +#ifdef KLU + if (ckt->CKTmatrix->CKTkluMODE) + SMPnnz (ckt->CKTmatrix) ; +#endif + +#ifdef KLU + if (ckt->CKTmatrix->CKTkluMODE) + { + fprintf (stderr, "Using KLU as Direct Linear Solver\n") ; + + int i ; + int n = SMPmatSize (ckt->CKTmatrix) ; + ckt->CKTmatrix->CKTkluN = n ; + + int nz = ckt->CKTmatrix->CKTklunz ; + + ckt->CKTmatrix->CKTkluAp = TMALLOC (int, n + 1) ; + ckt->CKTmatrix->CKTkluAi = TMALLOC (int, nz) ; + ckt->CKTmatrix->CKTkluAx = TMALLOC (double, nz) ; + ckt->CKTmatrix->CKTkluIntermediate = TMALLOC (double, n) ; + + ckt->CKTmatrix->CKTbindStruct = TMALLOC (BindElement, nz) ; + + ckt->CKTmatrix->CKTdiag_CSC = TMALLOC (double *, n) ; + + /* Complex Stuff needed for AC Analysis */ + ckt->CKTmatrix->CKTkluAx_Complex = TMALLOC (double, 2 * nz) ; + ckt->CKTmatrix->CKTkluIntermediate_Complex = TMALLOC (double, 2 * n) ; + + /* Binding Table from Sparse to CSC Format Creation */ + SMPmatrix_CSC (ckt->CKTmatrix) ; + + /* Binding Table Sorting */ + qsort (ckt->CKTmatrix->CKTbindStruct, (size_t)nz, sizeof(BindElement), BindCompare) ; + + /* KLU Pointers Assignment */ + for (i = 0 ; i < DEVmaxnum ; i++) + if (DEVices [i] && DEVices [i]->DEVbindCSC && ckt->CKThead [i]) + DEVices [i]->DEVbindCSC (ckt->CKThead [i], ckt) ; + + ckt->CKTmatrix->CKTkluMatrixIsComplex = CKTkluMatrixReal ; + } else { + fprintf (stderr, "Using SPARSE 1.3 as Direct Linear Solver\n") ; + } +#endif + for(i=0;i<=MAX(2,ckt->CKTmaxOrder)+1;i++) { /* dctran needs 3 states as minimum */ CKALLOC(ckt->CKTstates[i],ckt->CKTnumStates,double); } diff --git a/src/spicelib/analysis/cktsopt.c b/src/spicelib/analysis/cktsopt.c index 88cb6fbc7..a802892be 100644 --- a/src/spicelib/analysis/cktsopt.c +++ b/src/spicelib/analysis/cktsopt.c @@ -172,6 +172,13 @@ CKTsetOpt(CKTcircuit *ckt, JOB *anal, int opt, IFvalue *val) case OPT_EPSMIN: task->TSKepsmin = val->rValue; break; + +#ifdef KLU + case OPT_SPARSE: + task->TSKkluMODE = (val->iValue == 0); + break; +#endif + /* gtri - begin - wbk - add new options */ #ifdef XSPICE case OPT_EVT_MAX_OP_ALTER: @@ -328,7 +335,13 @@ static IFparm OPTtbl[] = { { "noopac", OPT_NOOPAC, IF_SET|IF_FLAG, "No op calculation in ac if circuit is linear" }, { "epsmin", OPT_EPSMIN, IF_SET|IF_REAL, - "Minimum value for log" } + "Minimum value for log" }, + +#ifdef KLU + { "sparse", OPT_SPARSE, IF_SET|IF_FLAG, + "Set SPARSE 1.3 as Direct Linear Solver" }, +#endif + }; int OPTcount = NUMELEMS(OPTtbl);