luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fixed Sensibility Analysis for KLU - First Trial

This commit is contained in:
Francesco Lannutti 2016-06-22 23:38:19 +02:00 committed by rlar
parent 4e42a93ba1
commit 74b6460326
7 changed files with 600 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

75
src/include/ngspice/klu.h
View File

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

8
src/maths/KLU/Makefile.am
View File

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

155
src/maths/KLU/klu_multiply.c Normal file
View File

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

132
src/maths/KLU/klu_utils.c Normal file
View File

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

6
src/maths/KLU/klu_version.h
View File

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

48
src/maths/KLU/klusmp.c
View File

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

193
src/spicelib/analysis/cktsens.c
View File

@ -45,6 +45,22 @@ static double inc_freq(double freq, int type, double step_size);
}
#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
/*
* 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;