From 74b6460326e46bc095c7b6febed62353dad34b87 Mon Sep 17 00:00:00 2001 From: Francesco Lannutti Date: Wed, 22 Jun 2016 23:38:19 +0200 Subject: [PATCH] Fixed Sensibility Analysis for KLU - First Trial --- src/include/ngspice/klu.h | 75 +++++++++++++ src/maths/KLU/Makefile.am | 8 +- src/maths/KLU/klu_multiply.c | 155 +++++++++++++++++++++++++ src/maths/KLU/klu_utils.c | 132 ++++++++++++++++++++++ src/maths/KLU/klu_version.h | 6 + src/maths/KLU/klusmp.c | 48 +++++++- src/spicelib/analysis/cktsens.c | 193 ++++++++++++++++++++++++++++++-- 7 files changed, 600 insertions(+), 17 deletions(-) create mode 100644 src/maths/KLU/klu_multiply.c create mode 100644 src/maths/KLU/klu_utils.c diff --git a/src/include/ngspice/klu.h b/src/include/ngspice/klu.h index 6e6238ace..aedacb518 100644 --- a/src/include/ngspice/klu.h +++ b/src/include/ngspice/klu.h @@ -861,6 +861,81 @@ int klu_z_print ) ; +int klu_constant_multiply +( + int *Ap, + double *Ax, + int n, + klu_common *Common, + double constant +) ; + +int klu_z_constant_multiply +( + int *Ap, + double *Ax, + int n, + klu_common *Common, + double constant +) ; + + +int klu_matrix_vector_multiply +( + int *Ap, /* CSR */ + int *Ai, /* CSR */ + double *Ax, /* CSR */ + double *RHS, + double *Solution, + int *IntToExtRowMap, + int *IntToExtColMap, + int n, + klu_common *Common +) ; + +int klu_z_matrix_vector_multiply +( + int *Ap, /* CSR */ + int *Ai, /* CSR */ + double *Ax, /* CSR */ + double *RHS, + double *Solution, + double *iRHS, + double *iSolution, + int *IntToExtRowMap, + int *IntToExtColMap, + int n, + klu_common *Common +) ; + + +int klu_convert_matrix_in_CSR +( + int *Ap_CSC, /* CSC */ + int *Ai_CSC, /* CSC */ + double *Ax_CSC, /* CSC */ + int *Ap_CSR, /* CSR */ + int *Ai_CSR, /* CSR */ + double *Ax_CSR, /* CSR */ + int n, + int nz, + klu_common *Common +) ; + +int klu_z_convert_matrix_in_CSR +( + int *Ap_CSC, /* CSC */ + int *Ai_CSC, /* CSC */ + double *Ax_CSC, /* CSC */ + int *Ap_CSR, /* CSR */ + int *Ai_CSR, /* CSR */ + double *Ax_CSR, /* CSR */ + int n, + int nz, + klu_common *Common +) ; + + /* ========================================================================== */ /* === KLU version ========================================================== */ /* ========================================================================== */ diff --git a/src/maths/KLU/Makefile.am b/src/maths/KLU/Makefile.am index aebc62484..66bedada8 100644 --- a/src/maths/KLU/Makefile.am +++ b/src/maths/KLU/Makefile.am @@ -10,11 +10,13 @@ libKLU_real_la_SOURCES = \ klu_factor.c \ klu_free_numeric.c \ klu_kernel.c \ + klu_multiply.c \ klu_refactor.c \ klu_scale.c \ klu_solve.c \ klu_sort.c \ - klu_tsolve.c + klu_tsolve.c \ + klu_utils.c libKLU_real_la_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include @@ -28,11 +30,13 @@ libKLU_complex_la_SOURCES = \ klu_factor.c \ klu_free_numeric.c \ klu_kernel.c \ + klu_multiply.c \ klu_refactor.c \ klu_scale.c \ klu_solve.c \ klu_sort.c \ - klu_tsolve.c + klu_tsolve.c \ + klu_utils.c libKLU_complex_la_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include -DCOMPLEX diff --git a/src/maths/KLU/klu_multiply.c b/src/maths/KLU/klu_multiply.c new file mode 100644 index 000000000..ca4d240c8 --- /dev/null +++ b/src/maths/KLU/klu_multiply.c @@ -0,0 +1,155 @@ +#include "klu_internal.h" + +Int KLU_constant_multiply /* return TRUE if successful, FALSE otherwise */ +( + Int *Ap, + double *Ax, + Int n, + KLU_common *Common, + double constant +) +{ + Entry *Az ; + int i, j ; + + /* ---------------------------------------------------------------------- */ + /* check inputs */ + /* ---------------------------------------------------------------------- */ + + if (Common == NULL) + { + return (FALSE) ; + } + + if (Ap == NULL || Ax == NULL) + { + Common->status = KLU_INVALID ; + return (FALSE) ; + } + + Common->status = KLU_OK ; + + Az = (Entry *)Ax ; + + for (i = 0 ; i < n ; i++) + { + for (j = Ap [i] ; j < Ap [i + 1] ; j++) + { + +#ifdef COMPLEX + Az [j].Real *= constant ; + Az [j].Imag *= constant ; +#else + Az [j] *= constant ; +#endif + + } + } + + return (TRUE) ; +} + +/* Macro function that multiplies two complex numbers and then adds them + * to another. to += mult_a * mult_b */ +#define CMPLX_MULT_ADD_ASSIGN(to,from_a,from_b) \ +{ (to).Real += (from_a).Real * (from_b).Real - \ + (from_a).Imag * (from_b).Imag; \ + (to).Imag += (from_a).Real * (from_b).Imag + \ + (from_a).Imag * (from_b).Real; \ +} + +Int KLU_matrix_vector_multiply /* return TRUE if successful, FALSE otherwise */ +( + Int *Ap, /* CSR */ + Int *Ai, /* CSR */ + double *Ax, /* CSR */ + double *RHS, + double *Solution, +#ifdef COMPLEX + double *iRHS, + double *iSolution, +#endif + Int *IntToExtRowMap, + Int *IntToExtColMap, + Int n, + KLU_common *Common +) +{ + Entry *Az, *Intermediate, Sum ; + Int i, j, *pExtOrder ; + + /* ---------------------------------------------------------------------- */ + /* check inputs */ + /* ---------------------------------------------------------------------- */ + + if (Common == NULL) + { + return (FALSE) ; + } + + if (Ap == NULL || Ai == NULL || Ax == NULL || RHS == NULL || Solution == NULL + +#ifdef COMPLEX + || iRHS == NULL || iSolution == NULL +#endif + + ) + { + Common->status = KLU_INVALID ; + return (FALSE) ; + } + + Common->status = KLU_OK ; + + + Intermediate = (Entry *) malloc ((size_t)n * sizeof (Entry)) ; + + Az = (Entry *)Ax ; + + pExtOrder = &IntToExtColMap [n] ; + for (i = n - 1 ; i >= 0 ; i--) + { + +#ifdef COMPLEX + Intermediate [i].Real = Solution [*(pExtOrder)] ; + Intermediate [i].Imag = iSolution [*(pExtOrder--)] ; +#else + Intermediate [i] = Solution [*(pExtOrder--)] ; +#endif + } + + pExtOrder = &IntToExtRowMap [n] ; + for (i = n - 1 ; i >= 0 ; i--) + { + +#ifdef COMPLEX + Sum.Real = 0.0 ; + Sum.Imag = 0.0 ; +#else + Sum = 0.0 ; +#endif + + for (j = Ap [i] ; j < Ap [i + 1] ; j++) + { + +#ifdef COMPLEX + CMPLX_MULT_ADD_ASSIGN (Sum, Az [j], Intermediate [Ai [j]]) ; +#else + Sum += Az [j] * Intermediate [Ai [j]] ; +#endif + + } + +#ifdef COMPLEX + RHS [*(pExtOrder)] = Sum.Real ; + iRHS [*(pExtOrder--)] = Sum.Imag ; +#else + RHS [*(pExtOrder--)] = Sum ; +#endif + + } + + free (Intermediate) ; + + return (TRUE) ; +} diff --git a/src/maths/KLU/klu_utils.c b/src/maths/KLU/klu_utils.c new file mode 100644 index 000000000..861792f20 --- /dev/null +++ b/src/maths/KLU/klu_utils.c @@ -0,0 +1,132 @@ +#include "klu_internal.h" + +typedef struct sElement { + Int row ; + Int col ; + Entry val ; +} 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 void +Compress (Int *Ai, Int *Bp, int num_rows, int n_COO) +{ + int i, j ; + + for (i = 0 ; i <= Ai [0] ; i++) + Bp [i] = 0 ; + + j = Ai [0] + 1 ; + for (i = 1 ; i < n_COO ; 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 <= num_rows ; j++) + Bp [j] = i ; +} + +Int +KLU_convert_matrix_in_CSR /* return TRUE if successful, FALSE otherwise */ +( + Int *Ap_CSC, /* CSC */ + Int *Ai_CSC, /* CSC */ + double *Ax_CSC, /* CSC */ + Int *Ap_CSR, /* CSR */ + Int *Ai_CSR, /* CSR */ + double *Ax_CSR, /* CSR */ + Int n, + Int nz, + KLU_common *Common +) +{ + Element *MatrixCOO ; + Entry *Az_CSC, *Az_CSR ; + Int *Ap_COO, i, j, k ; + + /* ---------------------------------------------------------------------- */ + /* check inputs */ + /* ---------------------------------------------------------------------- */ + + if (Common == NULL) + { + return (FALSE) ; + } + + if (Ap_CSC == NULL || Ai_CSC == NULL || Ax_CSC == NULL || Ap_CSR == NULL || Ai_CSR == NULL || Ax_CSR == NULL) + { + Common->status = KLU_INVALID ; + return (FALSE) ; + } + + Common->status = KLU_OK ; + + + MatrixCOO = (Element *) malloc ((size_t)nz * sizeof (Element)) ; + Ap_COO = (Int *) malloc ((size_t)nz * sizeof (Int)) ; + + Az_CSC = (Entry *)Ax_CSC ; + Az_CSR = (Entry *)Ax_CSR ; + + k = 0 ; + for (i = 0 ; i < n ; i++) + { + for (j = Ap_CSC [i] ; j < Ap_CSC [i + 1] ; j++) + { + MatrixCOO [k].row = Ai_CSC [j] ; + MatrixCOO [k].col = i ; + +#ifdef COMPLEX + MatrixCOO [k].val.Real = Az_CSC [j].Real ; + MatrixCOO [k].val.Imag = Az_CSC [j].Imag ; +#else + MatrixCOO [k].val = Az_CSC [j] ; +#endif + + k++ ; + } + } + + /* Order the MatrixCOO along the rows */ + qsort (MatrixCOO, (size_t)nz, sizeof(Element), CompareRow) ; + + for (i = 0 ; i < nz ; i++) + { + Ap_COO [i] = MatrixCOO [i].row ; + Ai_CSR [i] = MatrixCOO [i].col ; + +#ifdef COMPLEX + Az_CSR [i].Real = MatrixCOO [i].val.Real ; + Az_CSR [i].Imag = MatrixCOO [i].val.Imag ; +#else + Az_CSR [i] = MatrixCOO [i].val ; +#endif + + } + + Compress (Ap_COO, Ap_CSR, n, nz) ; + + free (MatrixCOO) ; + free (Ap_COO) ; + + return (TRUE) ; +} diff --git a/src/maths/KLU/klu_version.h b/src/maths/KLU/klu_version.h index 6ba7b5dfa..e64031c80 100644 --- a/src/maths/KLU/klu_version.h +++ b/src/maths/KLU/klu_version.h @@ -86,6 +86,9 @@ #define KLU_condest klu_z_condest #define KLU_flops klu_z_flops #define KLU_print klu_z_print +#define KLU_constant_multiply klu_z_constant_multiply +#define KLU_matrix_vector_multiply klu_z_matrix_vector_multiply +#define KLU_convert_matrix_in_CSR klu_z_convert_matrix_in_CSR #endif @@ -138,6 +141,9 @@ #define KLU_condest klu_condest #define KLU_flops klu_flops #define KLU_print klu_print +#define KLU_constant_multiply klu_constant_multiply +#define KLU_matrix_vector_multiply klu_matrix_vector_multiply +#define KLU_convert_matrix_in_CSR klu_convert_matrix_in_CSR #endif diff --git a/src/maths/KLU/klusmp.c b/src/maths/KLU/klusmp.c index 06e87f375..7816bfe06 100644 --- a/src/maths/KLU/klusmp.c +++ b/src/maths/KLU/klusmp.c @@ -523,7 +523,12 @@ SMPprint (SMPmatrix *Matrix, char *Filename) { if (Matrix->CKTkluMODE) { - klu_z_print (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx_Complex, Matrix->CKTkluN, Matrix->SPmatrix->IntToExtRowMap, Matrix->SPmatrix->IntToExtColMap) ; + if (Matrix->CKTkluMatrixIsComplex) + { + klu_z_print (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx_Complex, Matrix->CKTkluN, Matrix->SPmatrix->IntToExtRowMap, Matrix->SPmatrix->IntToExtColMap) ; + } else { + klu_print (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx, Matrix->CKTkluN, Matrix->SPmatrix->IntToExtRowMap, Matrix->SPmatrix->IntToExtColMap) ; + } } else { if (Filename) spFileMatrix (Matrix->SPmatrix, Filename, "Circuit Matrix", 0, 1, 1) ; @@ -1097,7 +1102,17 @@ SMPzeroRow (SMPmatrix *eMatrix, int Row) void SMPconstMult (SMPmatrix *Matrix, double constant) { - spConstMult (Matrix->SPmatrix, constant) ; + if (Matrix->CKTkluMODE) + { + if (Matrix->CKTkluMatrixIsComplex) + { + klu_z_constant_multiply (Matrix->CKTkluAp, Matrix->CKTkluAx_Complex, Matrix->CKTkluN, Matrix->CKTkluCommon, constant) ; + } else { + klu_constant_multiply (Matrix->CKTkluAp, Matrix->CKTkluAx, Matrix->CKTkluN, Matrix->CKTkluCommon, constant) ; + } + } else { + spConstMult (Matrix->SPmatrix, constant) ; + } } /* @@ -1106,5 +1121,32 @@ SMPconstMult (SMPmatrix *Matrix, double constant) void SMPmultiply (SMPmatrix *Matrix, double *RHS, double *Solution, double *iRHS, double *iSolution) { - spMultiply (Matrix->SPmatrix, RHS, Solution, iRHS, iSolution) ; + if (Matrix->CKTkluMODE) + { + int *Ap_CSR, *Ai_CSR ; + double *Ax_CSR ; + + Ap_CSR = (int *) malloc ((size_t)(Matrix->CKTkluN + 1) * sizeof (int)) ; + Ai_CSR = (int *) malloc ((size_t)Matrix->CKTklunz * sizeof (int)) ; + + if (Matrix->CKTkluMatrixIsComplex) + { + Ax_CSR = (double *) malloc ((size_t)(2 * Matrix->CKTklunz) * sizeof (double)) ; + klu_z_convert_matrix_in_CSR (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx_Complex, Ap_CSR, Ai_CSR, Ax_CSR, Matrix->CKTkluN, Matrix->CKTklunz, Matrix->CKTkluCommon) ; + klu_z_matrix_vector_multiply (Ap_CSR, Ai_CSR, Ax_CSR, RHS, Solution, iRHS, iSolution, Matrix->SPmatrix->IntToExtRowMap, + Matrix->SPmatrix->IntToExtColMap, Matrix->CKTkluN, Matrix->CKTkluCommon) ; + } else { + Ax_CSR = (double *) malloc ((size_t)Matrix->CKTklunz * sizeof (double)) ; + klu_convert_matrix_in_CSR (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx, Ap_CSR, Ai_CSR, Ax_CSR, Matrix->CKTkluN, Matrix->CKTklunz, Matrix->CKTkluCommon) ; + klu_matrix_vector_multiply (Ap_CSR, Ai_CSR, Ax_CSR, RHS, Solution, Matrix->SPmatrix->IntToExtRowMap, + Matrix->SPmatrix->IntToExtColMap, Matrix->CKTkluN, Matrix->CKTkluCommon) ; + iSolution = iRHS ; + } + + free (Ap_CSR) ; + free (Ai_CSR) ; + free (Ax_CSR) ; + } else { + spMultiply (Matrix->SPmatrix, RHS, Solution, iRHS, iSolution) ; + } } diff --git a/src/spicelib/analysis/cktsens.c b/src/spicelib/analysis/cktsens.c index da42b1ec3..dce4be602 100644 --- a/src/spicelib/analysis/cktsens.c +++ b/src/spicelib/analysis/cktsens.c @@ -45,6 +45,22 @@ static double inc_freq(double freq, int type, double step_size); } +#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 + + /* * Procedure: * @@ -62,7 +78,7 @@ static double inc_freq(double freq, int type, double step_size); static int error; int sens_sens(CKTcircuit *ckt, int restart) { - +/* #ifdef KLU if (ckt->CKTkluMODE) { @@ -70,7 +86,7 @@ int sens_sens(CKTcircuit *ckt, int restart) return OK ; } else { #endif - +*/ SENS_AN *job = (SENS_AN *) ckt->CKTcurJob; static int size; @@ -98,9 +114,11 @@ int sens_sens(CKTcircuit *ckt, int restart) int type; double *saved_rhs = NULL, *saved_irhs = NULL; - MatrixFrame *saved_matrix = NULL; + SMPmatrix *saved_matrix = NULL; - delta_Y = TMALLOC(SMPmatrix, 1); +#ifdef KLU + int nz, size_CSC; +#endif #ifndef notdef #ifdef notdef @@ -153,10 +171,38 @@ int sens_sens(CKTcircuit *ckt, int restart) size = SMPmatSize(ckt->CKTmatrix); /* Create the perturbation matrix */ + delta_Y = TMALLOC(SMPmatrix, 1); + +#ifdef KLU + delta_Y->CKTkluSymbolic = NULL; + delta_Y->CKTkluNumeric = NULL; + delta_Y->CKTkluAp = NULL; + delta_Y->CKTkluAi = NULL; + delta_Y->CKTkluAx = NULL; + delta_Y->CKTkluMatrixIsComplex = CKTkluMatrixReal; + delta_Y->CKTkluIntermediate = NULL; + delta_Y->CKTkluIntermediate_Complex = NULL; + delta_Y->CKTbindStruct = NULL; + delta_Y->CKTdiag_CSC = NULL; + delta_Y->CKTkluN = 0; + delta_Y->CKTklunz = 0; + delta_Y->CKTkluMODE = ckt->CKTkluMODE; + + if (ckt->CKTkluMODE) + { + delta_Y->CKTkluCommon = TMALLOC(klu_common, 1); + + klu_defaults(delta_Y->CKTkluCommon); + } else { + delta_Y->CKTkluCommon = NULL; + } +#endif + error = SMPnewMatrix(delta_Y, size); if (error) return error; + SMPprint(delta_Y, NULL); size += 1; /* Create an extra rhs */ @@ -253,6 +299,22 @@ int sens_sens(CKTcircuit *ckt, int restart) iE = ckt->CKTirhs; Y = ckt->CKTmatrix; +#ifdef KLU + if (ckt->CKTkluMODE) + { + /* Convert the KLU matrix to complex */ + for (i = 0 ; i < Y->CKTklunz ; i++) + { + Y->CKTkluAx_Complex [2 * i] = Y->CKTkluAx [i]; + Y->CKTkluAx_Complex [2 * i + 1] = 0.0; + } + + Y->CKTkluMatrixIsComplex = CKTkluMatrixComplex; + + SMPcReorder(Y, ckt->CKTpivotAbsTol, ckt->CKTpivotRelTol, &size_CSC); // size_CSC is just a placeholder here + } +#endif + #ifdef ASDEBUG DEBUG(1) { printf("Operating point:\n"); @@ -312,6 +374,22 @@ int sens_sens(CKTcircuit *ckt, int restart) error = CKTload(ckt); /* INITSMSIGS */ if (error) return error; + +//#ifdef KLU +// 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 + error = NIacIter(ckt); if (error) return error; @@ -333,7 +411,7 @@ int sens_sens(CKTcircuit *ckt, int restart) save_context(ckt->CKTrhs, saved_rhs); save_context(ckt->CKTirhs, saved_irhs); - save_context(ckt->CKTmatrix->SPmatrix, saved_matrix); + save_context(ckt->CKTmatrix, saved_matrix); ckt->CKTrhs = delta_I; ckt->CKTirhs = delta_iI; @@ -389,6 +467,50 @@ int sens_sens(CKTcircuit *ckt, int restart) } } +#ifdef KLU + if (ckt->CKTmatrix->CKTkluMODE) + { + size_CSC = SMPmatSize (ckt->CKTmatrix); + ckt->CKTmatrix->CKTkluN = size_CSC; + + SMPnnz (ckt->CKTmatrix); + nz = ckt->CKTmatrix->CKTklunz; + + ckt->CKTmatrix->CKTkluAp = TMALLOC (int, size_CSC + 1); + ckt->CKTmatrix->CKTkluAi = TMALLOC (int, nz); + ckt->CKTmatrix->CKTkluAx = TMALLOC (double, nz); + ckt->CKTmatrix->CKTkluIntermediate = TMALLOC (double, size_CSC); + + ckt->CKTmatrix->CKTbindStruct = TMALLOC (BindElement, nz); + + ckt->CKTmatrix->CKTdiag_CSC = TMALLOC (double *, size_CSC); + + /* Complex Stuff needed for AC Analysis */ + ckt->CKTmatrix->CKTkluAx_Complex = TMALLOC (double, 2 * nz); + ckt->CKTmatrix->CKTkluIntermediate_Complex = TMALLOC (double, 2 * size_CSC); + + /* 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 */ + if (DEVices [sg->dev]->DEVbindCSC) + DEVices [sg->dev]->DEVbindCSC (sg->model, ckt); + + ckt->CKTmatrix->CKTkluMatrixIsComplex = CKTkluMatrixReal; + + /* Clear KLU Vectors */ + for (i = 0; i < nz; i++) + { + ckt->CKTmatrix->CKTkluAx [i] = 0; + ckt->CKTmatrix->CKTkluAx_Complex [2 * i] = 0; + ckt->CKTmatrix->CKTkluAx_Complex [2 * i + 1] = 0; + } + } +#endif + /* ? CKTsetup would call NIreinit instead */ ckt->CKTniState = NISHOULDREORDER | NIACSHOULDREORDER; @@ -399,6 +521,18 @@ int sens_sens(CKTcircuit *ckt, int restart) /* XXX Leave original E until here!! so that temp reads * the right node voltages */ +#ifdef KLU + if (ckt->CKTkluMODE) + { + if (!is_dc) + { + if (DEVices [sg->dev]->DEVbindCSCComplex) + DEVices [sg->dev]->DEVbindCSCComplex (sg->model, ckt); + + ckt->CKTmatrix->CKTkluMatrixIsComplex = CKTkluMatrixComplex; + } + } +#endif if (sens_load(sg, ckt, is_dc)) { if (error && error != E_BADPARM) return error; /* XXX */ @@ -414,7 +548,7 @@ int sens_sens(CKTcircuit *ckt, int restart) #ifdef ASDEBUG DEBUG(2) { printf("Effect of device:\n"); - SMPprint(delta_Y->SPmatrix, NULL); + SMPprint(delta_Y, NULL); printf("LHS:\n"); for (j = 0; j < size; j++) printf("%d: %g, %g\n", j, @@ -439,6 +573,7 @@ int sens_sens(CKTcircuit *ckt, int restart) return error; SMPconstMult(delta_Y, -1.0); + for (j = 0; j < size; j++) { delta_I[j] *= -1.0; delta_iI[j] *= -1.0; @@ -497,8 +632,12 @@ int sens_sens(CKTcircuit *ckt, int restart) #endif /* delta_Y E */ +// fprintf(stderr, "\n\nPRIMA\n"); +// SMPprint(delta_Y, NULL); SMPmultiply(delta_Y, delta_I_delta_Y, E, delta_iI_delta_Y, iE); +// fprintf(stderr, "\n\nDOPO\n"); +// SMPprint(delta_Y, NULL); #ifdef ASDEBUG DEBUG(2) @@ -507,12 +646,24 @@ int sens_sens(CKTcircuit *ckt, int restart) j, delta_I_delta_Y[j]); #endif +// fprintf (stderr, "\n\nPRIMA 1\n"); +// for (j = 0; j < size; j++) +// { +// fprintf (stderr, "RHS [%d]: %-.9g j%-.9g\n", j, delta_I [j], delta_iI [j]); +// } + /* delta_I - delta_Y E */ for (j = 0; j < size; j++) { delta_I[j] -= delta_I_delta_Y[j]; delta_iI[j] -= delta_iI_delta_Y[j]; } +// fprintf (stderr, "\n\nDOPO 1\n"); +// for (j = 0; j < size; j++) +// { +// fprintf (stderr, "RHS [%d]: %-.9g j%-.9g\n", j, delta_I [j], delta_iI [j]); +// } + #ifdef ASDEBUG DEBUG(2) { printf(">>> Y:\n"); @@ -522,9 +673,22 @@ int sens_sens(CKTcircuit *ckt, int restart) delta_I[j], delta_iI[j]); } #endif + +// fprintf (stderr, "\n\nPRIMA\n"); +// for (j = 0; j < size; j++) +// { +// fprintf (stderr, "RHS [%d]: %-.14g j%-.14g\n", j, delta_I [j], delta_iI [j]); +// } + /* Solve; Y already factored */ SMPcSolve(Y, delta_I, delta_iI, NULL, NULL); +// fprintf (stderr, "\n\nDOPO\n"); +// for (j = 0; j < size; j++) +// { +// fprintf (stderr, "RHS [%d]: %-.14g j%-.14g\n", j, delta_I [j], delta_iI [j]); +// } + /* the special `0' node * the matrix indizes are [1..n] * yet the vector indizes are [0..n] @@ -557,14 +721,18 @@ int sens_sens(CKTcircuit *ckt, int restart) /* delta_I is now equal to delta_E */ if (is_dc) { - if (job->output_volt) + if (job->output_volt) { output_values[n] = delta_I [job->output_pos->number] - delta_I [job->output_neg->number]; +// fprintf (stderr, "Pos: %d\tNeg: %d\n", job->output_pos->number, job->output_neg->number); + } else { output_values[n] = delta_I[branch_eq]; } +// fprintf (stderr, "output_values real PRIMA: %-.14g\n", output_values [n]); output_values[n] /= delta_var; + fprintf(stderr, "output_values real DOPO: %-.14g - delta_var: %-.9g\n", output_values [n], delta_var); } else { if (job->output_volt) { output_cvalues[n].real = @@ -579,6 +747,7 @@ int sens_sens(CKTcircuit *ckt, int restart) output_cvalues[n].imag = delta_iI[branch_eq]; } +// fprintf (stderr, "output_values complex: %-.9g j%-.9g\n", output_cvalues [n].real, output_cvalues [n].imag); output_cvalues[n].real /= delta_var; output_cvalues[n].imag /= delta_var; } @@ -589,7 +758,7 @@ int sens_sens(CKTcircuit *ckt, int restart) release_context(ckt->CKTrhs, saved_rhs); release_context(ckt->CKTirhs, saved_irhs); - release_context(ckt->CKTmatrix->SPmatrix, saved_matrix); + release_context(ckt->CKTmatrix, saved_matrix); if (is_dc) nvalue.v.vec.rVec = output_values; @@ -614,7 +783,7 @@ int sens_sens(CKTcircuit *ckt, int restart) release_context(ckt->CKTrhs, saved_rhs); release_context(ckt->CKTirhs, saved_irhs); - release_context(ckt->CKTmatrix->SPmatrix, saved_matrix); + release_context(ckt->CKTmatrix, saved_matrix); SMPdestroy(delta_Y); FREE(delta_I); @@ -635,11 +804,11 @@ int sens_sens(CKTcircuit *ckt, int restart) #endif return OK; - +/* #ifdef KLU } #endif - +*/ } double @@ -724,7 +893,7 @@ count_steps(int type, double low, double high, int steps, double *stepsize) static int sens_load(sgen *sg, CKTcircuit *ckt, int is_dc) -{ +{//fprintf (stderr, "LOAD - is_dc: %d\n", is_dc); int (*fn) (GENmodel *, CKTcircuit *); error = 0;