KLU Integration from scratch #4, changed files

This commit is contained in:
Francesco Lannutti 2014-04-21 14:11:28 +02:00 committed by rlar
parent f37e049ce3
commit 4677aa608c
20 changed files with 408 additions and 76 deletions

View File

@ -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

View File

@ -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 \

View File

@ -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",

View File

@ -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
};

View File

@ -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 */

View File

@ -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 */

View File

@ -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 <math.h>
#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 *);

View File

@ -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 */

View File

@ -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

View File

@ -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

View File

@ -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);
}

View File

@ -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)

View File

@ -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];

View File

@ -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);
}

View File

@ -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;
}

View File

@ -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;

View File

@ -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

View File

@ -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;

View File

@ -30,6 +30,20 @@ int nthreads;
return(E_NOMEM);\
}
#ifdef KLU
#include <stdlib.h>
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);
}

View File

@ -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);