epilectrik
/
pyNgSpice
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

First KLU support of CIDER ONED simulations

pre-master-46
Francesco Lannutti 6 years ago
committed by Holger Vogt
parent
4c351f51bd
commit
11fb209ee6
15 changed files with 1555 additions and 123 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 231
      src/ciderlib/oned/oneadmit.c
  2. 21
      src/ciderlib/oned/onecond.c
  3. 397
      src/ciderlib/oned/onecont.c
  4. 20
      src/ciderlib/oned/onedest.c
  5. 4
      src/ciderlib/oned/onedext.h
  6. 137
      src/ciderlib/oned/onepoiss.c
  7. 21
      src/ciderlib/oned/oneproj.c
  8. 203
      src/ciderlib/oned/onesolve.c
  9. 35
      src/include/ngspice/klu-binding.h
  10. 27
      src/include/ngspice/klu.h
  11. 45
      src/include/ngspice/onemesh.h
  12. 26
      src/include/ngspice/smpdefs.h
  13. 12
      src/maths/KLU/klu_extract.c
  14. 462
      src/maths/KLU/klusmp.c
  15. 37
      src/maths/sparse/spsmp.c

231
src/ciderlib/oned/oneadmit.c
View File

@ -116,30 +116,54 @@ NUMDadmittance(ONEdevice *pDevice, double omega, SPcomplex *yd)
pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
}
ONE_jacLoad(pDevice);
spSetComplex(pDevice->matrix);
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
if (pElem->elemType == SEMICON) {
for (i = 0; i <= 1; i++) {
pNode = pElem->pNodes[i];
if (pNode->nodeType != CONTACT) {
spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
}
}
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
} else {
#endif
spSetComplex (pDevice->matrix->SPmatrix) ;
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
if (pElem->elemType == SEMICON) {
for (i = 0; i <= 1; i++) {
pNode = pElem->pNodes[i];
if (pNode->nodeType != CONTACT) {
spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
}
}
}
}
#ifdef KLU
}
#endif
pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
/* FACTOR */
startTime = SPfrontEnd->IFseconds();
spFactor(pDevice->matrix);
#ifdef KLU
SMPluFacKLUforCIDER (pDevice->matrix) ;
#else
SMPcLUfac(pDevice->matrix, 0);
#endif
pDevice->pStats->factorTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
/* SOLVE */
startTime = SPfrontEnd->IFseconds();
spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
pDevice->rhsImag, solutionImag);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#else
SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#endif
pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
}
/* MISC */
@ -288,30 +312,54 @@ NBJTadmittance(ONEdevice *pDevice, double omega, SPcomplex *yIeVce,
pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
}
spSetComplex(pDevice->matrix);
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
if (pElem->elemType == SEMICON) {
for (i = 0; i <= 1; i++) {
pNode = pElem->pNodes[i];
if (pNode->nodeType != CONTACT) {
spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
}
}
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
} else {
#endif
spSetComplex (pDevice->matrix->SPmatrix) ;
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
if (pElem->elemType == SEMICON) {
for (i = 0; i <= 1; i++) {
pNode = pElem->pNodes[i];
if (pNode->nodeType != CONTACT) {
spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
}
}
}
}
#ifdef KLU
}
#endif
pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
/* FACTOR */
startTime = SPfrontEnd->IFseconds();
spFactor(pDevice->matrix);
#ifdef KLU
SMPluFacKLUforCIDER (pDevice->matrix) ;
#else
SMPcLUfac(pDevice->matrix, 0);
#endif
pDevice->pStats->factorTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
/* SOLVE */
startTime = SPfrontEnd->IFseconds();
spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
pDevice->rhsImag, solutionImag);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#else
SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#endif
pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
/* MISC */
@ -345,8 +393,13 @@ NBJTadmittance(ONEdevice *pDevice, double omega, SPcomplex *yIeVce,
/* SOLVE */
startTime = SPfrontEnd->IFseconds();
spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
pDevice->rhsImag, solutionImag);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#else
SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#endif
pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
}
/* MISC */
@ -418,7 +471,11 @@ ONEsorSolve(ONEdevice *pDevice, double *xReal, double *xImag, double omega)
rhsSOR[index] += pDevice->rhs[index];
}
/* compute xReal(k+1). solution stored in rhsSOR */
spSolve(pDevice->matrix, rhsSOR, rhsSOR, NULL, NULL);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, rhsSOR, rhsSOR, NULL, NULL) ;
#else
SMPsolve(pDevice->matrix, rhsSOR, rhsSOR);
#endif
/* modify solution when wRelax is not 1 */
if (wRelax != 1) {
@ -454,8 +511,12 @@ ONEsorSolve(ONEdevice *pDevice, double *xReal, double *xImag, double omega)
}
}
/* compute xImag(k+1) */
spSolve(pDevice->matrix, rhsSOR, rhsSOR,
NULL, NULL);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, rhsSOR, rhsSOR, NULL, NULL) ;
#else
SMPsolve(pDevice->matrix, rhsSOR, rhsSOR);
#endif
/* modify solution when wRelax is not 1 */
if (wRelax != 1) {
for (index = 1; index <= numEqns; index++) {
@ -516,25 +577,44 @@ NUMDys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yd)
pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
}
spSetComplex(pDevice->matrix);
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
if (pElem->elemType == SEMICON) {
for (i = 0; i <= 1; i++) {
pNode = pElem->pNodes[i];
if (pNode->nodeType != CONTACT) {
CMPLX_MULT_SCALAR(temp, cOmega, 0.5 * pElem->dx);
spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
}
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
} else {
#endif
spSetComplex (pDevice->matrix->SPmatrix) ;
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
if (pElem->elemType == SEMICON) {
for (i = 0; i <= 1; i++) {
pNode = pElem->pNodes[i];
if (pNode->nodeType != CONTACT) {
CMPLX_MULT_SCALAR(temp, cOmega, 0.5 * pElem->dx);
spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
}
}
}
}
#ifdef KLU
}
#endif
#ifdef KLU
SMPluFacKLUforCIDER (pDevice->matrix) ;
#else
SMPcLUfac (pDevice->matrix, 0) ;
#endif
spFactor(pDevice->matrix);
spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
pDevice->rhsImag, solutionImag);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#else
SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#endif
pElem = pDevice->elemArray[1];
pNode = pElem->pLeftNode;
@ -584,24 +664,45 @@ NBJTys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yIeVce,
pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
}
spSetComplex(pDevice->matrix);
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
if (pElem->elemType == SEMICON) {
for (i = 0; i <= 1; i++) {
pNode = pElem->pNodes[i];
if (pNode->nodeType != CONTACT) {
CMPLX_MULT_SCALAR(temp, cOmega, 0.5 * pElem->dx);
spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
}
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
} else {
#endif
spSetComplex (pDevice->matrix->SPmatrix) ;
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
if (pElem->elemType == SEMICON) {
for (i = 0; i <= 1; i++) {
pNode = pElem->pNodes[i];
if (pNode->nodeType != CONTACT) {
CMPLX_MULT_SCALAR(temp, cOmega, 0.5 * pElem->dx);
spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
}
}
}
}
#ifdef KLU
}
#endif
#ifdef KLU
SMPluFacKLUforCIDER (pDevice->matrix) ;
#else
SMPcLUfac (pDevice->matrix, 0) ;
#endif
#ifdef KLU
SMPsolveKLUforCIDER(pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag);
#else
SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#endif
spFactor(pDevice->matrix);
spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
pDevice->rhsImag, solutionImag);
pElem = pDevice->elemArray[1];
pNode = pElem->pLeftNode;
y = computeAdmittance(pNode, FALSE, solutionReal, solutionImag, &cOmega);
@ -624,8 +725,12 @@ NBJTys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yIeVce,
}
/* don't need to LU factor the jacobian since it exists */
spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
pDevice->rhsImag, solutionImag);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#else
SMPcSolveForCIDER (pDevice->matrix, pDevice->rhs, solutionReal, pDevice->rhsImag, solutionImag) ;
#endif
pElem = pDevice->elemArray[1];
pNode = pElem->pLeftNode;
y = computeAdmittance(pNode, FALSE, solutionReal, solutionImag, &cOmega);

21
src/ciderlib/oned/onecond.c
View File

@ -42,7 +42,12 @@ NUMDconductance(ONEdevice *pDevice, BOOLEAN tranAnalysis,
pDevice->rhs[pNode->pEqn] = -pEdge->dJpDpsiP1;
}
incVpn = pDevice->dcDeltaSolution;
spSolve(pDevice->matrix, pDevice->rhs, incVpn, NULL, NULL);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVpn, NULL, NULL) ;
#else
SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVpn) ;
#endif
pElem = pDevice->elemArray[1];
pNode = pElem->pRightNode;
@ -95,7 +100,12 @@ NBJTconductance(ONEdevice *pDevice, BOOLEAN tranAnalysis, double *intCoeff,
pDevice->rhs[pNode->pEqn] = -pEdge->dJpDpsiP1;
}
incVce = pDevice->dcDeltaSolution;
spSolve(pDevice->matrix, pDevice->rhs, incVce, NULL, NULL);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVce, NULL, NULL) ;
#else
SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVce) ;
#endif
/* zero the rhs before loading in the new rhs base contribution */
for (index = 1; index <= pDevice->numEqns; index++) {
@ -114,7 +124,12 @@ NBJTconductance(ONEdevice *pDevice, BOOLEAN tranAnalysis, double *intCoeff,
}
incVbe = pDevice->copiedSolution;
spSolve(pDevice->matrix, pDevice->rhs, incVbe, NULL, NULL);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVbe, NULL, NULL) ;
#else
SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVbe) ;
#endif
pElem = pDevice->elemArray[1];/* first element */
pEdge = pElem->pEdge;

397
src/ciderlib/oned/onecont.c
View File

@ -15,6 +15,10 @@ Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
#include "ngspice/cidersupt.h"
#include "../../maths/misc/bernoull.h"
#ifdef KLU
#include "ngspice/klu-binding.h"
#endif
/* functions to setup and solve the continuity equations */
/* Both continuity equations are solved */
@ -39,21 +43,76 @@ ONE_jacBuild(ONEdevice *pDevice)
pNode = pElem->pNodes[index];
/* get poisson pointer */
psiEqn = pNode->psiEqn;
pNode->fPsiPsi = spGetElement(matrix, psiEqn, psiEqn);
#ifdef KLU
pNode->fPsiPsi = SMPmakeEltKLUforCIDER (matrix, psiEqn, psiEqn) ;
pNode->fPsiPsiBinding = NULL ;
#else
pNode->fPsiPsi = SMPmakeElt(matrix, psiEqn, psiEqn);
#endif
if (pElem->elemType == SEMICON) {
/* get continuity-coupling terms */
nEqn = pNode->nEqn;
pEqn = pNode->pEqn;
/* pointers for additional terms */
pNode->fPsiN = spGetElement(matrix, psiEqn, nEqn);
pNode->fPsiP = spGetElement(matrix, psiEqn, pEqn);
pNode->fNPsi = spGetElement(matrix, nEqn, psiEqn);
pNode->fNN = spGetElement(matrix, nEqn, nEqn);
pNode->fNP = spGetElement(matrix, nEqn, pEqn);
pNode->fPPsi = spGetElement(matrix, pEqn, psiEqn);
pNode->fPP = spGetElement(matrix, pEqn, pEqn);
pNode->fPN = spGetElement(matrix, pEqn, nEqn);
#ifdef KLU
pNode->fPsiN = SMPmakeEltKLUforCIDER (matrix, psiEqn, nEqn) ;
pNode->fPsiNBinding = NULL ;
#else
pNode->fPsiN = SMPmakeElt(matrix, psiEqn, nEqn);
#endif
#ifdef KLU
pNode->fPsiP = SMPmakeEltKLUforCIDER (matrix, psiEqn, pEqn) ;
pNode->fPsiPBinding = NULL ;
#else
pNode->fPsiP = SMPmakeElt(matrix, psiEqn, pEqn);
#endif
#ifdef KLU
pNode->fNPsi = SMPmakeEltKLUforCIDER (matrix, nEqn, psiEqn) ;
pNode->fNPsiBinding = NULL ;
#else
pNode->fNPsi = SMPmakeElt(matrix, nEqn, psiEqn);
#endif
#ifdef KLU
pNode->fNN = SMPmakeEltKLUforCIDER (matrix, nEqn, nEqn) ;
pNode->fNNBinding = NULL ;
#else
pNode->fNN = SMPmakeElt(matrix, nEqn, nEqn);
#endif
#ifdef KLU
pNode->fNP = SMPmakeEltKLUforCIDER (matrix, nEqn, pEqn) ;
pNode->fNPBinding = NULL ;
#else
pNode->fNP = SMPmakeElt(matrix, nEqn, pEqn);
#endif
#ifdef KLU
pNode->fPPsi = SMPmakeEltKLUforCIDER (matrix, pEqn, psiEqn) ;
pNode->fPPsiBinding = NULL ;
#else
pNode->fPPsi = SMPmakeElt(matrix, pEqn, psiEqn);
#endif
#ifdef KLU
pNode->fPP = SMPmakeEltKLUforCIDER (matrix, pEqn, pEqn) ;
pNode->fPPBinding = NULL ;
#else
pNode->fPP = SMPmakeElt(matrix, pEqn, pEqn);
#endif
#ifdef KLU
pNode->fPN = SMPmakeEltKLUforCIDER (matrix, pEqn, nEqn) ;
pNode->fPNBinding = NULL ;
#else
pNode->fPN = SMPmakeElt(matrix, pEqn, nEqn);
#endif
} else {
nEqn = 0;
pEqn = 0;
@ -73,34 +132,302 @@ ONE_jacBuild(ONEdevice *pDevice)
/* now terms to couple to adjacent nodes */
pNode = pElem->pLeftNode;
pNode->fPsiPsiiP1 = spGetElement(matrix, psiEqnL, psiEqnR);
#ifdef KLU
pNode->fPsiPsiiP1 = SMPmakeEltKLUforCIDER (matrix, psiEqnL, psiEqnR) ;
pNode->fPsiPsiiP1Binding = NULL ;
#else
pNode->fPsiPsiiP1 = SMPmakeElt(matrix, psiEqnL, psiEqnR);
#endif
if (pElem->elemType == SEMICON) {
/* pointers for additional terms */
pNode->fNPsiiP1 = spGetElement(matrix, nEqnL, psiEqnR);
pNode->fNNiP1 = spGetElement(matrix, nEqnL, nEqnR);
pNode->fPPsiiP1 = spGetElement(matrix, pEqnL, psiEqnR);
pNode->fPPiP1 = spGetElement(matrix, pEqnL, pEqnR);
#ifdef KLU
pNode->fNPsiiP1 = SMPmakeEltKLUforCIDER (matrix, nEqnL, psiEqnR) ;
pNode->fNPsiiP1Binding = NULL ;
#else
pNode->fNPsiiP1 = SMPmakeElt(matrix, nEqnL, psiEqnR);
#endif
#ifdef KLU
pNode->fNNiP1 = SMPmakeEltKLUforCIDER (matrix, nEqnL, nEqnR) ;
pNode->fNNiP1Binding = NULL ;
#else
pNode->fNNiP1 = SMPmakeElt(matrix, nEqnL, nEqnR);
#endif
#ifdef KLU
pNode->fPPsiiP1 = SMPmakeEltKLUforCIDER (matrix, pEqnL, psiEqnR) ;
pNode->fPPsiiP1Binding = NULL ;
#else
pNode->fPPsiiP1 = SMPmakeElt(matrix, pEqnL, psiEqnR);
#endif
#ifdef KLU
pNode->fPPiP1 = SMPmakeEltKLUforCIDER (matrix, pEqnL, pEqnR) ;
pNode->fPPiP1Binding = NULL ;
#else
pNode->fPPiP1 = SMPmakeElt(matrix, pEqnL, pEqnR);
#endif
if (AvalancheGen) {
pNode->fNPiP1 = spGetElement(matrix, nEqnL, pEqnR);
pNode->fPNiP1 = spGetElement(matrix, pEqnL, nEqnR);
#ifdef KLU
pNode->fNPiP1 = SMPmakeEltKLUforCIDER (matrix, nEqnL, pEqnR) ;
pNode->fNPiP1Binding = NULL ;
#else
pNode->fNPiP1 = SMPmakeElt(matrix, nEqnL, pEqnR);
#endif
#ifdef KLU
pNode->fPNiP1 = SMPmakeEltKLUforCIDER (matrix, pEqnL, nEqnR) ;
pNode->fPNiP1Binding = NULL ;
#else
pNode->fPNiP1 = SMPmakeElt(matrix, pEqnL, nEqnR);
#endif
}
}
pNode = pElem->pRightNode;
pNode->fPsiPsiiM1 = spGetElement(matrix, psiEqnR, psiEqnL);
#ifdef KLU
pNode->fPsiPsiiM1 = SMPmakeEltKLUforCIDER (matrix, psiEqnR, psiEqnL) ;
pNode->fPsiPsiiM1Binding = NULL ;
#else
pNode->fPsiPsiiM1 = SMPmakeElt(matrix, psiEqnR, psiEqnL);
#endif
if (pElem->elemType == SEMICON) {
/* pointers for additional terms */
pNode->fNPsiiM1 = spGetElement(matrix, nEqnR, psiEqnL);
pNode->fNNiM1 = spGetElement(matrix, nEqnR, nEqnL);
pNode->fPPsiiM1 = spGetElement(matrix, pEqnR, psiEqnL);
pNode->fPPiM1 = spGetElement(matrix, pEqnR, pEqnL);
#ifdef KLU
pNode->fNPsiiM1 = SMPmakeEltKLUforCIDER (matrix, nEqnR, psiEqnL) ;
pNode->fNPsiiM1Binding = NULL ;
#else
pNode->fNPsiiM1 = SMPmakeElt(matrix, nEqnR, psiEqnL);
#endif
#ifdef KLU
pNode->fNNiM1 = SMPmakeEltKLUforCIDER (matrix, nEqnR, nEqnL) ;
pNode->fNNiM1Binding = NULL ;
#else
pNode->fNNiM1 = SMPmakeElt(matrix, nEqnR, nEqnL);
#endif
#ifdef KLU
pNode->fPPsiiM1 = SMPmakeEltKLUforCIDER (matrix, pEqnR, psiEqnL) ;
pNode->fPPsiiM1Binding = NULL ;
#else
pNode->fPPsiiM1 = SMPmakeElt(matrix, pEqnR, psiEqnL);
#endif
#ifdef KLU
pNode->fPPiM1 = SMPmakeEltKLUforCIDER (matrix, pEqnR, pEqnL) ;
pNode->fPPiM1Binding = NULL ;
#else
pNode->fPPiM1 = SMPmakeElt(matrix, pEqnR, pEqnL);
#endif
if (AvalancheGen) {
pNode->fNPiM1 = spGetElement(matrix, nEqnR, pEqnL);
pNode->fPNiM1 = spGetElement(matrix, pEqnR, nEqnL);
#ifdef KLU
pNode->fNPiM1 = SMPmakeEltKLUforCIDER (matrix, nEqnR, pEqnL) ;
pNode->fNPiM1Binding = NULL ;
#else
pNode->fNPiM1 = SMPmakeElt(matrix, nEqnR, pEqnL);
#endif
#ifdef KLU
pNode->fPNiM1 = SMPmakeEltKLUforCIDER (matrix, pEqnR, nEqnL) ;
pNode->fPNiM1Binding = NULL ;
#else
pNode->fPNiM1 = SMPmakeElt(matrix, pEqnR, nEqnL);
#endif
}
}
}
}
#ifdef KLU
/*
#define CREATE_KLU_BINDING_TABLE_CIDER(ptr, binding, a, b) \
printf ("Swapping Pointer %s: (%d,%d)\n", #ptr, a, b) ; \
if ((a > 0) && (b > 0)) { \
if (pNode->binding != NULL) { \
if (pNode->binding->CSC_Complex != NULL) { \
printf (" Looking for the Pointer: %p\n", pNode->binding->CSC_Complex) ; \
qsort (BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
i.COO = NULL ; \
i.CSC_Complex = pNode->binding->CSC_Complex ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
if (matched != NULL) { \
printf (" Found the Old Pointer\n") ; \
pNode->ptr = pNode->binding->CSC_Complex ; \
} else { \
i.COO = pNode->ptr ; \
i.CSC_Complex = NULL ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
if (matched != NULL) { \
printf (" Looking for the Pointer 1\n") ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} else { \
printf (" Leaving the Pointer as is\n") ; \
} \
} \
} else { \
printf (" Looking for the Pointer 2\n") ; \
i.COO = pNode->ptr ; \
i.CSC_Complex = NULL ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} \
} else { \
printf (" Looking for the Pointer 3\n") ; \
i.COO = pNode->ptr ; \
i.CSC_Complex = NULL ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} \
}
*/
/*
#define CREATE_KLU_BINDING_TABLE_CIDER_TO_REAL(ptr, binding, a, b) \
if ((a > 0) && (b > 0)) { \
printf ("Macro\n") ; \
if (pNode->binding) { \
printf ("IF: %p\n", pNode->binding) ; \
printf ("COO: %p\n", pNode->binding->COO) ; \
printf ("CSC: %p\n", pNode->binding->CSC) ; \
if (pNode->binding->CSC_Complex) { \
printf ("CSC_Complex: %p\n", pNode->binding->CSC_Complex) ; \
pNode->ptr = pNode->binding->CSC_Complex ; \
} else { \
i = pNode->ptr ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} \
} else { \
i = pNode->ptr ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} \
}
*/
void
ONEbindCSC (ONEdevice *pDevice)
{
ONEelem *pElem;
ONEnode *pNode;
int index, eIndex;
int psiEqn, nEqn, pEqn; /* scratch for deref'd eqn numbers */
int psiEqnL=0, nEqnL=0, pEqnL=0;
int psiEqnR=0, nEqnR=0, pEqnR=0;
BindKluElementCOO i, *matched, *BindStruct, *BindStructCSC ;
size_t nz ;
BindStruct = pDevice->matrix->SMPkluMatrix->KLUmatrixBindStructCOO ;
nz = pDevice->matrix->SMPkluMatrix->KLUmatrixNZ ;
BindStructCSC = (BindKluElementCOO *) malloc (nz * sizeof(BindKluElementCOO)) ;
for (index = 0 ; index < (int)nz ; index++) {
BindStructCSC [index] = BindStruct [index] ;
}
for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
pElem = pDevice->elemArray[eIndex];
/* first the self terms */
for (index = 0; index <= 1; index++) {
pNode = pElem->pNodes[index];
/* get poisson pointer */
psiEqn = pNode->psiEqn;
CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsi, fPsiPsiBinding, psiEqn, psiEqn) ;
if (pElem->elemType == SEMICON) {
/* get continuity-coupling terms */
nEqn = pNode->nEqn;
pEqn = pNode->pEqn;
/* pointers for additional terms */
CREATE_KLU_BINDING_TABLE_CIDER(fPsiN, fPsiNBinding, psiEqn, nEqn) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPsiP, fPsiPBinding, psiEqn, pEqn) ;
CREATE_KLU_BINDING_TABLE_CIDER(fNPsi, fNPsiBinding, nEqn, psiEqn) ;
CREATE_KLU_BINDING_TABLE_CIDER(fNN, fNNBinding, nEqn, nEqn) ;
CREATE_KLU_BINDING_TABLE_CIDER(fNP, fNPBinding, nEqn, pEqn) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPPsi, fPPsiBinding, pEqn, psiEqn) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPP, fPPBinding, pEqn, pEqn) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPN, fPNBinding, pEqn, nEqn) ;
} else {
nEqn = 0;
pEqn = 0;
}
/* save indices */
if (index == 0) { /* left node */
psiEqnL = psiEqn;
nEqnL = nEqn;
pEqnL = pEqn;
} else {
psiEqnR = psiEqn;
nEqnR = nEqn;
pEqnR = pEqn;
}
}
/* now terms to couple to adjacent nodes */
pNode = pElem->pLeftNode;
CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsiiP1, fPsiPsiiP1Binding, psiEqnL, psiEqnR) ;
if (pElem->elemType == SEMICON) {
/* pointers for additional terms */
CREATE_KLU_BINDING_TABLE_CIDER(fNPsiiP1, fNPsiiP1Binding, nEqnL, psiEqnR) ;
CREATE_KLU_BINDING_TABLE_CIDER(fNNiP1, fNNiP1Binding, nEqnL, nEqnR) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPPsiiP1, fPPsiiP1Binding, pEqnL, psiEqnR) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPPiP1, fPPiP1Binding, pEqnL, pEqnR) ;
if (AvalancheGen) {
CREATE_KLU_BINDING_TABLE_CIDER(fNPiP1, fNPiP1Binding, nEqnL, pEqnR) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPNiP1, fPNiP1Binding, pEqnL, nEqnR) ;
}
}
pNode = pElem->pRightNode;
CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsiiM1, fPsiPsiiM1Binding, psiEqnR, psiEqnL) ;
if (pElem->elemType == SEMICON) {
/* pointers for additional terms */
CREATE_KLU_BINDING_TABLE_CIDER(fNPsiiM1, fNPsiiM1Binding, nEqnR, psiEqnL) ;
CREATE_KLU_BINDING_TABLE_CIDER(fNNiM1, fNNiM1Binding, nEqnR, nEqnL) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPPsiiM1, fPPsiiM1Binding, pEqnR, psiEqnL) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPPiM1, fPPiM1Binding, pEqnR, pEqnL) ;
if (AvalancheGen) {
CREATE_KLU_BINDING_TABLE_CIDER(fNPiM1, fNPiM1Binding, nEqnR, pEqnL) ;
CREATE_KLU_BINDING_TABLE_CIDER(fPNiM1, fPNiM1Binding, pEqnR, nEqnL) ;
}
}
}
free (BindStructCSC) ;
}
#endif
void
ONE_sysLoad(ONEdevice *pDevice, BOOLEAN tranAnalysis,
@ -131,7 +458,17 @@ ONE_sysLoad(ONEdevice *pDevice, BOOLEAN tranAnalysis,
}
/* zero the matrix */
spClear(pDevice->matrix);
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
SMPclearKLUforCIDER (pDevice->matrix) ;
} else {
#endif
SMPclear(pDevice->matrix);
#ifdef KLU
}
#endif
for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
pElem = pDevice->elemArray[eIndex];
@ -275,7 +612,17 @@ ONE_jacLoad(ONEdevice *pDevice)
ONE_commonTerms(pDevice, FALSE, FALSE, NULL);
/* zero the matrix */
spClear(pDevice->matrix);
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
SMPclearKLUforCIDER (pDevice->matrix) ;
} else {
#endif
SMPclear(pDevice->matrix);
#ifdef KLU
}
#endif
for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
pElem = pDevice->elemArray[eIndex];

20
src/ciderlib/oned/onedest.c
View File

@ -35,7 +35,15 @@ ONEdestroy(ONEdevice *pDevice)
FREE(pDevice->copiedSolution);
FREE(pDevice->rhs);
FREE(pDevice->rhsImag);
spDestroy(pDevice->matrix);
#ifdef KLU
SMPdestroyKLUforCIDER (pDevice->matrix) ;
#else
SMPdestroy (pDevice->matrix) ;
#endif
FREE (pDevice->matrix) ;
break;
case SLV_EQUIL:
/* free up the vectors allocated in the equilibrium solution */
@ -43,7 +51,15 @@ ONEdestroy(ONEdevice *pDevice)
FREE(pDevice->dcDeltaSolution);
FREE(pDevice->copiedSolution);
FREE(pDevice->rhs);
spDestroy(pDevice->matrix);
#ifdef KLU
SMPdestroyKLUforCIDER (pDevice->matrix) ;
#else
SMPdestroy (pDevice->matrix) ;
#endif
FREE (pDevice->matrix) ;
break;
case SLV_NONE:
break;

4
src/ciderlib/oned/onedext.h
View File

@ -115,5 +115,9 @@ extern void ONEjacCheck(ONEdevice *, BOOLEAN, ONEtranInfo *);
extern void ONEpredict(ONEdevice *, ONEtranInfo *);
extern BOOLEAN ONEdeviceConverged(ONEdevice *);
#ifdef KLU
void ONEbindCSC (ONEdevice *) ;
void ONEQbindCSC (ONEdevice *) ;
#endif
#endif

137
src/ciderlib/oned/onepoiss.c
View File

@ -12,6 +12,10 @@ Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
#include "oneddefs.h"
#include "ngspice/spmatrix.h"
#ifdef KLU
#include "ngspice/klu-binding.h"
#endif
/* Functions to setup and solve the 1D poisson equation. */
@ -27,16 +31,129 @@ ONEQjacBuild(ONEdevice *pDevice)
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
pNode = pElem->pLeftNode;
pNode->fPsiPsi = spGetElement(matrix, pNode->poiEqn, pNode->poiEqn);
#ifdef KLU
pNode->fPsiPsi = SMPmakeEltKLUforCIDER (matrix, pNode->poiEqn, pNode->poiEqn) ;
pNode->fPsiPsiBinding = NULL ;
#else
pNode->fPsiPsi = SMPmakeElt(matrix, pNode->poiEqn, pNode->poiEqn);
#endif
pNode1 = pElem->pRightNode;
pNode->fPsiPsiiP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn);
#ifdef KLU
pNode->fPsiPsiiP1 = SMPmakeEltKLUforCIDER (matrix, pNode->poiEqn, pNode1->poiEqn) ;
pNode->fPsiPsiiP1Binding = NULL ;
#else
pNode->fPsiPsiiP1 = SMPmakeElt(matrix, pNode->poiEqn, pNode1->poiEqn);
#endif
pNode = pElem->pRightNode;
pNode->fPsiPsi = spGetElement(matrix, pNode->poiEqn, pNode->poiEqn);
#ifdef KLU
pNode->fPsiPsi = SMPmakeEltKLUforCIDER (matrix, pNode->poiEqn, pNode->poiEqn) ;
pNode->fPsiPsiBinding = NULL ;
#else
pNode->fPsiPsi = SMPmakeElt(matrix, pNode->poiEqn, pNode->poiEqn);
#endif
pNode1 = pElem->pLeftNode;
pNode->fPsiPsiiM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn);
#ifdef KLU
pNode->fPsiPsiiM1 = SMPmakeEltKLUforCIDER (matrix, pNode->poiEqn, pNode1->poiEqn) ;
pNode->fPsiPsiiM1Binding = NULL ;
#else
pNode->fPsiPsiiM1 = SMPmakeElt(matrix, pNode->poiEqn, pNode1->poiEqn);
#endif
}
}
#ifdef KLU
/*
#define CREATE_KLU_BINDING_TABLE_CIDER(ptr, binding, a, b) \
printf ("Swapping Pointer %s: (%d,%d)\n", #ptr, a, b) ; \
if ((a > 0) && (b > 0)) { \
if (pNode->binding != NULL) { \
if (pNode->binding->CSC_Complex != NULL) { \
printf (" Looking for the Pointer: %p\n", pNode->binding->CSC_Complex) ; \
qsort (BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
i.COO = NULL ; \
i.CSC_Complex = pNode->binding->CSC_Complex ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
if (matched != NULL) { \
printf (" Found the Old Pointer\n") ; \
pNode->ptr = pNode->binding->CSC_Complex ; \
} else { \
i.COO = pNode->ptr ; \
i.CSC_Complex = NULL ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
if (matched != NULL) { \
printf (" Looking for the Pointer 1\n") ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} else { \
printf (" Leaving the Pointer as is\n") ; \
} \
} \
} else { \
printf (" Looking for the Pointer 2\n") ; \
i.COO = pNode->ptr ; \
i.CSC_Complex = NULL ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} \
} else { \
printf (" Looking for the Pointer 3\n") ; \
i.COO = pNode->ptr ; \
i.CSC_Complex = NULL ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} \
}
*/
void
ONEQbindCSC (ONEdevice *pDevice)
{
ONEelem *pElem ;
ONEnode *pNode, *pNode1 ;
int index ;
BindKluElementCOO i, *matched, *BindStruct, *BindStructCSC ;
size_t nz ;
BindStruct = pDevice->matrix->SMPkluMatrix->KLUmatrixBindStructCOO ;
nz = pDevice->matrix->SMPkluMatrix->KLUmatrixNZ ;
BindStructCSC = (BindKluElementCOO *) malloc (nz * sizeof(BindKluElementCOO)) ;
for (index = 0 ; index < (int)nz ; index++) {
BindStructCSC [index] = BindStruct [index] ;
}
for (index = 1 ; index < pDevice->numNodes ; index++) {
pElem = pDevice->elemArray [index] ;
pNode = pElem->pLeftNode ;
CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsi, fPsiPsiBinding, pNode->poiEqn, pNode->poiEqn) ;
pNode1 = pElem->pRightNode ;
CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsiiP1, fPsiPsiiP1Binding, pNode->poiEqn, pNode1->poiEqn) ;
pNode = pElem->pRightNode ;
CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsi, fPsiPsiBinding, pNode->poiEqn, pNode->poiEqn) ;
pNode1 = pElem->pLeftNode ;
CREATE_KLU_BINDING_TABLE_CIDER(fPsiPsiiM1, fPsiPsiiM1Binding, pNode->poiEqn, pNode1->poiEqn) ;
}
free (BindStructCSC) ;
}
#endif
void
ONEQsysLoad(ONEdevice *pDevice)
@ -58,7 +175,17 @@ ONEQsysLoad(ONEdevice *pDevice)
}
/* zero the matrix */
spClear(pDevice->matrix);
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
SMPclearKLUforCIDER (pDevice->matrix) ;
} else {
#endif
SMPclear(pDevice->matrix);
#ifdef KLU
}
#endif
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];

21
src/ciderlib/oned/oneproj.c
View File

@ -52,7 +52,12 @@ NUMDproject(ONEdevice *pDevice, double delV)
pDevice->rhs[pNode->pEqn] = -pEdge->dJpDpsiP1;
}
incVpn = pDevice->dcDeltaSolution;
spSolve(pDevice->matrix, pDevice->rhs, incVpn, NULL, NULL);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVpn, NULL, NULL) ;
#else
SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVpn) ;
#endif
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
@ -128,7 +133,12 @@ NBJTproject(ONEdevice *pDevice, double delVce, double delVbe,
pDevice->rhs[pNode->pEqn] = -pEdge->dJpDpsiP1;
}
incVce = pDevice->dcDeltaSolution;
spSolve(pDevice->matrix, pDevice->rhs, incVce, NULL, NULL);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVce, NULL, NULL) ;
#else
SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVce) ;
#endif
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];
@ -179,7 +189,12 @@ NBJTproject(ONEdevice *pDevice, double delVce, double delVbe,
}
incVbe = pDevice->copiedSolution;
spSolve(pDevice->matrix, pDevice->rhs, incVbe, NULL, NULL);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, incVbe, NULL, NULL) ;
#else
SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, incVbe) ;
#endif
for (index = 1; index < pDevice->numNodes; index++) {
pElem = pDevice->elemArray[index];

203
src/ciderlib/oned/onesolve.c
View File

@ -91,7 +91,13 @@ ONEdcSolve(ONEdevice *pDevice, int iterationLimit, BOOLEAN newSolver,
}
/* FACTOR */
startTime = SPfrontEnd->IFseconds();
error = spFactor(pDevice->matrix);
#ifdef KLU
error = SMPreorderKLUforCIDER (pDevice->matrix) ;
#else
error = SMPluFacForCIDER (pDevice->matrix) ;
#endif
factorTime += SPfrontEnd->IFseconds() - startTime;
if (newSolver) {
if (pDevice->iterationNumber == 1) {
@ -110,7 +116,7 @@ ONEdcSolve(ONEdevice *pDevice, int iterationLimit, BOOLEAN newSolver,
if (foundError(error)) {
if (error == spSINGULAR) {
int badRow, badCol;
spWhereSingular(pDevice->matrix, &badRow, &badCol);
SMPgetError(pDevice->matrix, &badCol, &badRow);
printf("***** singular at (%d,%d)\n", badRow, badCol);
}
/* Should probably try to recover now, but we'll punt instead. */
@ -118,7 +124,13 @@ ONEdcSolve(ONEdevice *pDevice, int iterationLimit, BOOLEAN newSolver,
}
/* SOLVE */
startTime = SPfrontEnd->IFseconds();
spSolve(pDevice->matrix, rhs, delta, NULL, NULL);
#ifdef KLU
SMPsolveKLUforCIDER (pDevice->matrix, rhs, delta, NULL, NULL) ;
#else
SMPsolveForCIDER (pDevice->matrix, rhs, delta) ;
#endif
solveTime += SPfrontEnd->IFseconds() - startTime;
/* UPDATE */
@ -441,7 +453,13 @@ ONEresetJacobian(ONEdevice *pDevice)
ONE_jacLoad(pDevice);
error = spFactor(pDevice->matrix);
#ifdef KLU
error = SMPreorderKLUforCIDER (pDevice->matrix) ;
#else
error = SMPluFacForCIDER (pDevice->matrix) ;
#endif
if (foundError(error)) {
exit(-1);
}
@ -539,7 +557,6 @@ ONEequilSolve(ONEdevice *pDevice)
ONEnode *pNode;
double startTime, setupTime, miscTime;
setupTime = miscTime = 0.0;
/* SETUP */
@ -553,7 +570,15 @@ ONEequilSolve(ONEdevice *pDevice)
FREE(pDevice->copiedSolution);
FREE(pDevice->rhs);
FREE(pDevice->rhsImag);
spDestroy(pDevice->matrix);
#ifdef KLU
SMPdestroyKLUforCIDER (pDevice->matrix) ;
#else
SMPdestroy (pDevice->matrix) ;
#endif
FREE (pDevice->matrix) ;
/* FALLTHRU */
case SLV_NONE:
pDevice->poissonOnly = TRUE;
@ -562,15 +587,63 @@ ONEequilSolve(ONEdevice *pDevice)
XCALLOC(pDevice->dcDeltaSolution, double, pDevice->dimEquil);
XCALLOC(pDevice->copiedSolution, double, pDevice->dimEquil);
XCALLOC(pDevice->rhs, double, pDevice->dimEquil);
pDevice->matrix = spCreate(pDevice->numEqns, 0, &error);
pDevice->matrix = TMALLOC (SMPmatrix, 1) ;
#ifdef KLU
pDevice->matrix->CKTkluMODE = CKTkluON ; /* Francesco Lannutti - To be sustitued with a value coming from the uplevel */
error = SMPnewMatrixKLUforCIDER (pDevice->matrix, pDevice->numEqns, CKTkluMatrixReal) ;
#else
error = SMPnewMatrixForCIDER (pDevice->matrix, pDevice->numEqns, 0) ;
#endif
if (error == spNO_MEMORY) {
printf("ONEequilSolve: Out of Memory\n");
exit(-1);
}
newSolver = TRUE;
spSetReal(pDevice->matrix);
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixReal ;
} else {
#endif
spSetReal (pDevice->matrix->SPmatrix) ;
#ifdef KLU
}
#endif
ONEQjacBuild(pDevice);
pDevice->numOrigEquil = spElementCount(pDevice->matrix);
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
/* Convert the COO Storage to CSC for KLU and Fill the Binding Table */
SMPconvertCOOtoCSCKLUforCIDER (pDevice->matrix) ;
/* Bind the KLU Pointers */
ONEQbindCSC (pDevice) ;
/* Perform KLU Matrix Analysis */
pDevice->matrix->SMPkluMatrix->KLUmatrixSymbolic = klu_analyze ((int)pDevice->matrix->SMPkluMatrix->KLUmatrixN, pDevice->matrix->SMPkluMatrix->KLUmatrixAp,
pDevice->matrix->SMPkluMatrix->KLUmatrixAi, pDevice->matrix->SMPkluMatrix->KLUmatrixCommon) ;
if (pDevice->matrix->SMPkluMatrix->KLUmatrixSymbolic == NULL) {
printf ("CIDER: KLU Failed\n") ;
if (pDevice->matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) {
return ; // Francesco Lannutti - Fix KLU return values
}
}
printf ("CIDER: KLU to be fixed: spElementCount\n") ;
pDevice->numOrigEquil = 0 ; // Francesco Lannutti - Fix for KLU
} else {
pDevice->numOrigEquil = spElementCount (pDevice->matrix->SPmatrix) ;
}
#else
pDevice->numOrigEquil = spElementCount (pDevice->matrix->SPmatrix) ;
#endif
pDevice->numFillEquil = 0;
/* FALLTHRU */
case SLV_EQUIL:
@ -590,7 +663,21 @@ ONEequilSolve(ONEdevice *pDevice)
/* MISCELLANEOUS */
startTime = SPfrontEnd->IFseconds();
if (newSolver) {
pDevice->numFillEquil = spFillinCount(pDevice->matrix);
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
// Francesco Lannutti - Fix for KLU
printf ("CIDER: KLU to be fixed: spFillinCount\n") ;
pDevice->numFillEquil = 0 ;
} else {
#endif
pDevice->numFillEquil = spFillinCount(pDevice->matrix->SPmatrix);
#ifdef KLU
}
#endif
}
if (pDevice->converged) {
ONEQcommonTerms(pDevice);
@ -642,7 +729,15 @@ ONEbiasSolve(ONEdevice *pDevice, int iterationLimit,
FREE(pDevice->dcDeltaSolution);
FREE(pDevice->copiedSolution);
FREE(pDevice->rhs);
spDestroy(pDevice->matrix);
#ifdef KLU
SMPdestroyKLUforCIDER (pDevice->matrix) ;
#else
SMPdestroy (pDevice->matrix) ;
#endif
FREE (pDevice->matrix) ;
/* FALLTHRU */
case SLV_NONE:
pDevice->poissonOnly = FALSE;
@ -652,18 +747,65 @@ ONEbiasSolve(ONEdevice *pDevice, int iterationLimit,
XCALLOC(pDevice->copiedSolution, double, pDevice->dimBias);
XCALLOC(pDevice->rhs, double, pDevice->dimBias);
XCALLOC(pDevice->rhsImag, double, pDevice->dimBias);
pDevice->matrix = spCreate(pDevice->numEqns, 1, &error);
pDevice->matrix = TMALLOC (SMPmatrix, 1) ;
#ifdef KLU
pDevice->matrix->CKTkluMODE = CKTkluON ; /* Francesco Lannutti - To be sustitued with a value coming from the uplevel */
error = SMPnewMatrixKLUforCIDER (pDevice->matrix, pDevice->numEqns, CKTkluMatrixComplex) ;
#else
error = SMPnewMatrixForCIDER (pDevice->matrix, pDevice->numEqns, 1) ;
#endif
if (error == spNO_MEMORY) {
exit(-1);
}
newSolver = TRUE;
ONE_jacBuild(pDevice);
pDevice->numOrigBias = spElementCount(pDevice->matrix);
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
/* Convert the COO Storage to CSC for KLU and Fill the Binding Table */
SMPconvertCOOtoCSCKLUforCIDER (pDevice->matrix) ;
/* Bind the KLU Pointers */
ONEbindCSC (pDevice) ;
/* Perform KLU Matrix Analysis */
pDevice->matrix->SMPkluMatrix->KLUmatrixSymbolic = klu_analyze ((int)pDevice->matrix->SMPkluMatrix->KLUmatrixN, pDevice->matrix->SMPkluMatrix->KLUmatrixAp,
pDevice->matrix->SMPkluMatrix->KLUmatrixAi, pDevice->matrix->SMPkluMatrix->KLUmatrixCommon) ;
if (pDevice->matrix->SMPkluMatrix->KLUmatrixSymbolic == NULL) {
if (pDevice->matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) {
printf ("CIDER: KLU failed\n") ;
return ; // Francesco Lannutti - Fix KLU return values
}
}
pDevice->numOrigBias = 0 ; // Francesco Lannutti - Fix for KLU
} else {
pDevice->numOrigBias = spElementCount(pDevice->matrix->SPmatrix);
}
#else
pDevice->numOrigBias = spElementCount (pDevice->matrix->SPmatrix) ;
#endif
pDevice->numFillBias = 0;
ONEstoreInitialGuess(pDevice);
/* FALLTHRU */
case SLV_SMSIG:
spSetReal(pDevice->matrix);
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixReal ;
} else {
#endif
spSetReal (pDevice->matrix->SPmatrix) ;
#ifdef KLU
}
#endif
/* FALLTHRU */
case SLV_BIAS:
pDevice->solverType = SLV_BIAS;
@ -681,7 +823,21 @@ ONEbiasSolve(ONEdevice *pDevice, int iterationLimit,
/* MISCELLANEOUS */
startTime = SPfrontEnd->IFseconds();
if (newSolver) {
pDevice->numFillBias = spFillinCount(pDevice->matrix);
#ifdef KLU
if (pDevice->matrix->CKTkluMODE) {
// Francesco Lannutti - Fix for KLU
printf ("CIDER: KLU to be fixed: spFillinCount\n") ;
pDevice->numFillBias = 0 ;
} else {
#endif
pDevice->numFillBias = spFillinCount (pDevice->matrix->SPmatrix) ; // Francesco Lannutti - Fix for KLU
#ifdef KLU
}
#endif
}
solution = pDevice->dcSolution;
if ((!pDevice->converged) && iterationLimit > 1) {
@ -999,8 +1155,12 @@ double
ONEnuNorm(ONEdevice *pDevice)
{
/* The LU Decomposed matrix is available. Use it to calculate x. */
spSolve(pDevice->matrix, pDevice->rhs, pDevice->rhsImag,
NULL, NULL);
#ifdef KLU
printf ("CIDER: KLU to be fixed ONEnuNorm\n") ;
SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, pDevice->rhsImag, NULL, NULL) ;
#else
SMPsolveForCIDER (pDevice->matrix, pDevice->rhs, pDevice->rhsImag) ;
#endif
/* Compute L2-norm of the solution vector (stored in rhsImag) */
return (l2Norm(pDevice->rhsImag, pDevice->numEqns));
@ -1041,7 +1201,14 @@ ONEjacCheck(ONEdevice *pDevice, BOOLEAN tranAnalysis, ONEtranInfo *info)
pDevice->dcSolution[index] = pDevice->copiedSolution[index];
for (rIndex = 1; rIndex <= pDevice->numEqns; rIndex++) {
diff = (pDevice->rhsImag[rIndex] - pDevice->rhs[rIndex]) / del;
dptr = spFindElement(pDevice->matrix, rIndex, index);
#ifdef KLU
printf ("CIDER: KLU to be fixed: spFindElement") ;
dptr = NULL ;
#else
dptr = spFindElement(pDevice->matrix->SPmatrix, rIndex, index); // Francesco Lannutti - Fix for KLU
#endif
/*
* if ( dptr ISNOT NULL ) { fprintf( stderr, "[%d][%d]: FD =
* %11.4e, AJ = %11.4e\n", rIndex, index, diff, *dptr ); } else {

35
src/include/ngspice/klu-binding.h
View File

@ -99,4 +99,39 @@
smp_data_out->input[k].port[l].ptr = smp_data_out->input[k].port[l].h.binding->CSC ;
#endif
#ifdef CIDER
#define CREATE_KLU_BINDING_TABLE_CIDER(ptr, binding, a, b) \
if ((a > 0) && (b > 0)) { \
if (pNode->binding != NULL) { \
if (pNode->binding->CSC_Complex != NULL) { \
qsort (BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
i.COO = NULL ; \
i.CSC_Complex = pNode->binding->CSC_Complex ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStructCSC, nz, sizeof(BindKluElementCOO), BindKluCompareCSC) ; \
if (matched == NULL) { \
i.COO = pNode->ptr ; \
i.CSC_Complex = NULL ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
if (matched != NULL) { \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} \
} \
} else { \
i.COO = pNode->ptr ; \
i.CSC_Complex = NULL ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} \
} else { \
i.COO = pNode->ptr ; \
i.CSC_Complex = NULL ; \
matched = (BindKluElementCOO *) bsearch (&i, BindStruct, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ; \
pNode->binding = matched ; \
pNode->ptr = matched->CSC_Complex ; \
} \
}
#endif
#endif

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

@ -935,6 +935,33 @@ int klu_z_convert_matrix_in_CSR
klu_common *Common
) ;
#ifdef CIDER
typedef struct sBindKluElementCOO {
double *COO ;
double *CSC_Complex ;
} BindKluElementCOO ;
int BindKluCompareCOO (const void *a, const void *b) ;
int BindKluCompareCSC (const void *a, const void *b) ;
typedef struct sKLUmatrix {
klu_common *KLUmatrixCommon ; /* KLU common object */
klu_symbolic *KLUmatrixSymbolic ; /* KLU symbolic object */
klu_numeric *KLUmatrixNumeric ; /* KLU numeric object */
int *KLUmatrixAp ; /* KLU column pointer */
int *KLUmatrixAi ; /* KLU row pointer */
double *KLUmatrixAxComplex ; /* KLU Complex Elements */
unsigned int KLUmatrixIsComplex:1 ; /* KLU Matrix Is Complex Flag */
double *KLUmatrixIntermediateComplex ; /* KLU iRHS Intermediate for Solve Complex Step */
unsigned int KLUmatrixN ; /* KLU N, copied */
unsigned int KLUmatrixNZ ; /* KLU nz, copied for AC Analysis */
int *KLUmatrixColCOO ; /* KLU Col Index for COO storage */
int *KLUmatrixRowCOO ; /* KLU Row Index for COO storage */
double *KLUmatrixValueComplexCOO ; /* KLU Complex Elements for COO storage */
BindKluElementCOO *KLUmatrixBindStructCOO ; /* KLU COO Binding Structure */
double *KLUmatrixTrashCOO ; /* KLU COO Trash Pointer for Ground Node not Stored in the Matrix */
} KLUmatrix ;
#endif
/* ========================================================================== */
/* === KLU version ========================================================== */

45
src/include/ngspice/onemesh.h
View File

@ -16,6 +16,10 @@ Authors: 1987 Karti Mayaram, 1991 David Gates
#include "ngspice/material.h"
#ifdef KLU
#include "ngspice/klu.h"
#endif
typedef struct sONEelem {
struct sONEelem *pElems[2]; /* array to store neighbor elements */
struct sONEnode *pNodes[2]; /* array to store the element nodes */
@ -27,6 +31,20 @@ typedef struct sONEelem {
ONEmaterial *matlInfo; /* material information */
double epsRel; /* epsilon of material rel to Semicon */
int evalNodes[2]; /* nodes to be evaluated in elem */
#ifdef KLU
double *KLUleftLeftNode ;
double *KLUleftRightNode ;
double *KLUrightLeftNode ;
double *KLUrightRightNode ;
BindKluElementCOO *KLUleftLeftNodeBinding ;
BindKluElementCOO *KLUleftRightNodeBinding ;
BindKluElementCOO *KLUrightLeftNodeBinding ;
BindKluElementCOO *KLUrightRightNodeBinding ;
#endif
} ONEelem;
#define pLeftElem pElems[0]
@ -112,6 +130,33 @@ typedef struct sONEnode {
double *fPNiM1;
double *fPN;
double *fPNiP1;
#ifdef KLU
BindKluElementCOO *fPsiPsiiM1Binding ;
BindKluElementCOO *fPsiPsiBinding ;
BindKluElementCOO *fPsiPsiiP1Binding ;
BindKluElementCOO *fPsiNBinding ;
BindKluElementCOO *fPsiPBinding ;
BindKluElementCOO *fNPsiiM1Binding ;
BindKluElementCOO *fNPsiBinding ;
BindKluElementCOO *fNPsiiP1Binding ;
BindKluElementCOO *fNNiM1Binding ;
BindKluElementCOO *fNNBinding ;
BindKluElementCOO *fNNiP1Binding ;
BindKluElementCOO *fNPiM1Binding ;
BindKluElementCOO *fNPBinding ;
BindKluElementCOO *fNPiP1Binding ;
BindKluElementCOO *fPPsiiM1Binding ;
BindKluElementCOO *fPPsiBinding ;
BindKluElementCOO *fPPsiiP1Binding ;
BindKluElementCOO *fPPiM1Binding ;
BindKluElementCOO *fPPBinding ;
BindKluElementCOO *fPPiP1Binding ;
BindKluElementCOO *fPNiM1Binding ;
BindKluElementCOO *fPNBinding ;
BindKluElementCOO *fPNiP1Binding ;
#endif
} ONEnode;

26
src/include/ngspice/smpdefs.h
View File

@ -47,6 +47,11 @@ typedef struct sSMPmatrix {
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 */
#ifdef CIDER
KLUmatrix *SMPkluMatrix ; /* KLU Pointer to the KLU Matrix Data Structure (only for CIDER, for the moment) */
#endif
#endif
} SMPmatrix ;
@ -56,7 +61,21 @@ typedef struct sSMPmatrix {
void SMPmatrix_CSC (SMPmatrix *) ;
void SMPnnz (SMPmatrix *) ;
void spDeterminant_KLU (SMPmatrix *, int *, double *, double *) ;
#ifdef CIDER
void SMPsolveKLUforCIDER (SMPmatrix *, double [], double [], double [], double []) ;
int SMPreorderKLUforCIDER (SMPmatrix *) ;
double *SMPmakeEltKLUforCIDER (SMPmatrix *, int, int) ;
void SMPclearKLUforCIDER (SMPmatrix *) ;
void SMPconvertCOOtoCSCKLUforCIDER (SMPmatrix *) ;
void SMPdestroyKLUforCIDER (SMPmatrix *) ;
int SMPnewMatrixKLUforCIDER (SMPmatrix *, int, unsigned int) ;
int SMPluFacKLUforCIDER (SMPmatrix *) ;
void SMPprintKLUforCIDER (SMPmatrix *, char *) ;
#endif
#endif
int SMPaddElt( SMPmatrix *, int , int , double );
double * SMPmakeElt( SMPmatrix * , int , int );
void SMPcClear( SMPmatrix *);
@ -85,4 +104,11 @@ int SMPzeroRow(SMPmatrix *Matrix, int Row);
void SMPconstMult(SMPmatrix *, double);
void SMPmultiply(SMPmatrix *, double *, double *, double *, double *);
#ifdef CIDER
void SMPcSolveForCIDER (SMPmatrix *, double [], double [], double [], double []) ;
int SMPluFacForCIDER (SMPmatrix *) ;
int SMPnewMatrixForCIDER (SMPmatrix *, int, int) ;
void SMPsolveForCIDER (SMPmatrix *, double [], double []) ;
#endif
#endif

12
src/maths/KLU/klu_extract.c
View File

@ -448,9 +448,17 @@ Int KLU_print
{
#ifdef COMPLEX
fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j].Real, Az [j].Imag) ;
if (IntToExtRowMap && IntToExtColMap) {
fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j].Real, Az [j].Imag) ;
} else {
fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", Ai [j] + 1, i + 1, Az [j].Real, Az [j].Imag) ;
}
#else
fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j]) ;
if (IntToExtRowMap && IntToExtColMap) {
fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j]) ;
} else {
fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", Ai [j] + 1, i + 1, Az [j]) ;
}
#endif
}

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

@ -130,6 +130,170 @@ SMPnnz (SMPmatrix *Matrix)
return ;
}
#ifdef CIDER
typedef struct sElement {
unsigned int row ;
unsigned int col ;
double *pointer ;
} 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 int
CompareColumn (const void *a, const void *b)
{
Element *A = (Element *) a ;
Element *B = (Element *) b ;
return
(A->col > B->col) ? 1 :
(A->col < B->col) ? -1 :
0 ;
}
static void
Compress (unsigned int *Ai, unsigned int *Bp, unsigned int n, unsigned int nz)
{
unsigned int i, j ;
for (i = 0 ; i <= Ai [0] ; i++)
Bp [i] = 0 ;
j = Ai [0] + 1 ;
for (i = 1 ; i < nz ; 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 <= n ; j++)
Bp [j] = i ;
}
int
BindKluCompareCOO (const void *a, const void *b)
{
BindKluElementCOO *A = (BindKluElementCOO *) a ;
BindKluElementCOO *B = (BindKluElementCOO *) b ;
return
(A->COO > B->COO) ? 1 :
(A->COO < B->COO) ? -1 :
0 ;
}
int
BindKluCompareCSC (const void *a, const void *b)
{
BindKluElementCOO *A = (BindKluElementCOO *) a ;
BindKluElementCOO *B = (BindKluElementCOO *) b ;
return
(A->CSC_Complex > B->CSC_Complex) ? 1 :
(A->CSC_Complex < B->CSC_Complex) ? -1 :
0 ;
}
void SMPconvertCOOtoCSCKLUforCIDER (SMPmatrix *Matrix)
{
Element *MatrixCOO ;
unsigned int *Ap_COO, i, j, nz ;
/* Count the non-zero elements and store it */
nz = 0 ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
if ((Matrix->SMPkluMatrix->KLUmatrixRowCOO [i] != -1) && (Matrix->SMPkluMatrix->KLUmatrixColCOO [i] != -1)) {
nz++ ;
}
}
Matrix->SMPkluMatrix->KLUmatrixNZ = nz ;
/* Allocate the compressed COO elements */
MatrixCOO = (Element *) malloc (nz * sizeof(Element)) ;
/* Allocate the temporary COO Column Index */
Ap_COO = (unsigned int *) malloc (nz * sizeof(unsigned int)) ;
/* Allocate the needed KLU data structures */
Matrix->SMPkluMatrix->KLUmatrixAp = (int *) malloc ((Matrix->SMPkluMatrix->KLUmatrixN + 1) * sizeof(int)) ;
Matrix->SMPkluMatrix->KLUmatrixAi = (int *) malloc (nz * sizeof(int)) ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO = (BindKluElementCOO *) malloc (nz * sizeof(BindKluElementCOO)) ;
Matrix->SMPkluMatrix->KLUmatrixAxComplex = (double *) malloc (2 * nz * sizeof(double)) ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = (double *) malloc (2 * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(double)) ;
/* Populate the compressed COO elements and COO value of Binding Table */
j = 0 ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
if ((Matrix->SMPkluMatrix->KLUmatrixRowCOO [i] != -1) && (Matrix->SMPkluMatrix->KLUmatrixColCOO [i] != -1)) {
MatrixCOO [j].row = (unsigned int)Matrix->SMPkluMatrix->KLUmatrixRowCOO [i] ;
MatrixCOO [j].col = (unsigned int)Matrix->SMPkluMatrix->KLUmatrixColCOO [i] ;
MatrixCOO [j].pointer = &(Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO [2 * i]) ;
j++ ;
}
}
/* Order the MatrixCOO along the columns */
qsort (MatrixCOO, nz, sizeof(Element), CompareColumn) ;
/* Order the MatrixCOO along the rows */
i = 0 ;
while (i < nz)
{
/* Look for the next column */
for (j = i + 1 ; j < nz ; j++)
{
if (MatrixCOO [j].col != MatrixCOO [i].col)
{
break ;
}
}
qsort (MatrixCOO + i, j - i, sizeof(Element), CompareRow) ;
i = j ;
}
/* Copy back the Matrix in partial CSC */
for (i = 0 ; i < nz ; i++)
{
Ap_COO [i] = MatrixCOO [i].col ;
Matrix->SMPkluMatrix->KLUmatrixAi [i] = (int)MatrixCOO [i].row ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO [i].COO = MatrixCOO [i].pointer ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO [i].CSC_Complex = &(Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i]) ;
}
/* Compress the COO Column Index to CSC Column Index */
Compress (Ap_COO, (unsigned int *)Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixN, nz) ;
/* Free the temporary stuff */
free (Ap_COO) ;
free (MatrixCOO) ;
/* Sort the Binding Table */
qsort (Matrix->SMPkluMatrix->KLUmatrixBindStructCOO, nz, sizeof(BindKluElementCOO), BindKluCompareCOO) ;
return ;
}
#endif
/*
* SMPaddElt()
*/
@ -149,6 +313,26 @@ SMPmakeElt (SMPmatrix *Matrix, int Row, int Col)
return spGetElement (Matrix->SPmatrix, Row, Col) ;
}
#ifdef CIDER
double *
SMPmakeEltKLUforCIDER (SMPmatrix *Matrix, int Row, int Col)
{
if (Matrix->CKTkluMODE) {
if ((Row > 0) && (Col > 0)) {
Row = Row - 1 ;
Col = Col - 1 ;
Matrix->SMPkluMatrix->KLUmatrixRowCOO [Row * (int)Matrix->SMPkluMatrix->KLUmatrixN + Col] = Row ;
Matrix->SMPkluMatrix->KLUmatrixColCOO [Row * (int)Matrix->SMPkluMatrix->KLUmatrixN + Col] = Col ;
return &(Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO [2 * (Row * (int)Matrix->SMPkluMatrix->KLUmatrixN + Col)]) ;
} else {
return Matrix->SMPkluMatrix->KLUmatrixTrashCOO ;
}
} else {
return spGetElement (Matrix->SPmatrix, Row, Col) ;
}
}
#endif
/*
* SMPcClear()
*/
@ -191,6 +375,18 @@ SMPclear (SMPmatrix *Matrix)
}
}
#ifdef CIDER
void
SMPclearKLUforCIDER (SMPmatrix *Matrix)
{
unsigned int i ;
for (i = 0 ; i < 2 * Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
Matrix->SMPkluMatrix->KLUmatrixAxComplex [i] = 0 ;
}
}
#endif
#define NG_IGNORE(x) (void)x
/*
@ -262,6 +458,46 @@ SMPluFac (SMPmatrix *Matrix, double PivTol, double Gmin)
}
}
#ifdef CIDER
int
SMPluFacKLUforCIDER (SMPmatrix *Matrix)
{
unsigned int i ;
double *KLUmatrixAx ;
int ret ;
if (Matrix->CKTkluMODE)
{
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex) {
ret = klu_z_refactor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAxComplex,
Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
} else {
/* Allocate the Real Matrix */
KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof(double)) ;
/* Copy the Complex Matrix into the Real Matrix */
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
KLUmatrixAx [i] = Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] ;
}
/* Re-Factor the Real Matrix */
ret = klu_refactor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, KLUmatrixAx,
Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
/* Free the Real Matrix Storage */
free (KLUmatrixAx) ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) {
printf ("CIDER: KLU Empty Matrix\n") ;
return 0 ;
}
return (!ret) ;
} else {
return spFactor (Matrix->SPmatrix) ;
}
}
#endif
/*
* SMPcReorder()
*/
@ -336,6 +572,58 @@ SMPreorder (SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin)
}
}
#ifdef CIDER
int
SMPreorderKLUforCIDER (SMPmatrix *Matrix)
{
unsigned int i ;
double *KLUmatrixAx ;
if (Matrix->CKTkluMODE)
{
// Matrix->CKTkluCommon->tol = PivTol ;
if (Matrix->SMPkluMatrix->KLUmatrixNumeric != NULL) {
klu_free_numeric (&(Matrix->SMPkluMatrix->KLUmatrixNumeric), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
}
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex) {
Matrix->SMPkluMatrix->KLUmatrixNumeric = klu_z_factor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi,
Matrix->SMPkluMatrix->KLUmatrixAxComplex, Matrix->SMPkluMatrix->KLUmatrixSymbolic,
Matrix->SMPkluMatrix->KLUmatrixCommon) ;
} else {
/* Allocate the Real Matrix */
KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof(double)) ;
/* Copy the Complex Matrix into the Real Matrix */
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
KLUmatrixAx [i] = Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] ;
}
/* Factor the Real Matrix */
Matrix->SMPkluMatrix->KLUmatrixNumeric = klu_factor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi,
KLUmatrixAx, Matrix->SMPkluMatrix->KLUmatrixSymbolic,
Matrix->SMPkluMatrix->KLUmatrixCommon) ;
/* Free the Real Matrix Storage */
free (KLUmatrixAx) ;
}
if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL)
{
printf ("CIDER: KLU Factorization Error\n") ;
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX)
{
return 0 ;
}
return 1 ;
}
else
return 0 ;
} else {
return spFactor (Matrix->SPmatrix) ;
}
}
#endif
/*
* SMPcaSolve()
*/
@ -405,6 +693,57 @@ SMPcSolve (SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[], doubl
}
}
#ifdef CIDER
void
SMPsolveKLUforCIDER (SMPmatrix *Matrix, double RHS[], double RHSsolution[], double iRHS[], double iRHSsolution[])
{
int ret ;
unsigned int i ;
double *KLUmatrixIntermediate ;
if (Matrix->CKTkluMODE)
{
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex) {
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
{
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] = RHS [i + 1] ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] = iRHS [i + 1] ;
}
ret = klu_z_solve (Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, (int)Matrix->SMPkluMatrix->KLUmatrixN, 1,
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
{
RHSsolution [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] ;
iRHSsolution [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] ;
}
} else {
/* Allocate the Intermediate Vector */
KLUmatrixIntermediate = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * sizeof(double)) ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
KLUmatrixIntermediate [i] = RHS [i + 1] ;
}
ret = klu_solve (Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, (int)Matrix->SMPkluMatrix->KLUmatrixN, 1,
KLUmatrixIntermediate, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
RHSsolution [i + 1] = KLUmatrixIntermediate [i] ;
}
/* Free the Intermediate Vector */
free (KLUmatrixIntermediate) ;
}
} else {
spSolve (Matrix->SPmatrix, RHS, RHSsolution, iRHS, iRHSsolution) ;
}
}
#endif
/*
* SMPsolve()
*/
@ -452,6 +791,58 @@ SMPnewMatrix (SMPmatrix *Matrix, int size)
return Error ;
}
#ifdef CIDER
int
SMPnewMatrixKLUforCIDER (SMPmatrix *Matrix, int size, unsigned int KLUmatrixIsComplex)
{
int Error ;
unsigned int i ;
if (Matrix->CKTkluMODE) {
/* Allocate the KLU Matrix Data Structure */
Matrix->SMPkluMatrix = (KLUmatrix *) malloc (sizeof (KLUmatrix)) ;
/* Initialize the KLU Matrix Internal Pointers */
Matrix->SMPkluMatrix->KLUmatrixCommon = (klu_common *) malloc (sizeof (klu_common)) ; ;
Matrix->SMPkluMatrix->KLUmatrixSymbolic = NULL ;
Matrix->SMPkluMatrix->KLUmatrixNumeric = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAp = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAi = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAxComplex = NULL ;
if (KLUmatrixIsComplex) {
Matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixComplex ;
} else {
Matrix->SMPkluMatrix->KLUmatrixIsComplex = CKTkluMatrixReal ;
}
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixNZ = 0 ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO = NULL ;
Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO = NULL ;
/* Initialize the KLU Common Data Structure */
klu_defaults (Matrix->SMPkluMatrix->KLUmatrixCommon) ;
/* Allocate KLU data structures */
Matrix->SMPkluMatrix->KLUmatrixN = (unsigned int)size ;
Matrix->SMPkluMatrix->KLUmatrixColCOO = (int *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(int)) ;
Matrix->SMPkluMatrix->KLUmatrixRowCOO = (int *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(int)) ;
Matrix->SMPkluMatrix->KLUmatrixTrashCOO = (double *) malloc (sizeof(double)) ;
Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO = (double *) malloc (2 * Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(double)) ;
/* Pre-set the values of Row and Col */
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
Matrix->SMPkluMatrix->KLUmatrixRowCOO [i] = -1 ;
Matrix->SMPkluMatrix->KLUmatrixColCOO [i] = -1 ;
}
return spOKAY ;
} else {
Matrix->SPmatrix = spCreate (size, (int)KLUmatrixIsComplex, &Error) ;
return Error ;
}
}
#endif
/*
* SMPdestroy()
*/
@ -486,6 +877,38 @@ SMPdestroy (SMPmatrix *Matrix)
}
}
#ifdef CIDER
void
SMPdestroyKLUforCIDER (SMPmatrix *Matrix)
{
if (Matrix->CKTkluMODE)
{
klu_free_numeric (&(Matrix->SMPkluMatrix->KLUmatrixNumeric), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
klu_free_symbolic (&(Matrix->SMPkluMatrix->KLUmatrixSymbolic), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
free (Matrix->SMPkluMatrix->KLUmatrixAp) ;
free (Matrix->SMPkluMatrix->KLUmatrixAi) ;
free (Matrix->SMPkluMatrix->KLUmatrixAxComplex) ;
free (Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex) ;
free (Matrix->SMPkluMatrix->KLUmatrixBindStructCOO) ;
free (Matrix->SMPkluMatrix->KLUmatrixColCOO) ;
free (Matrix->SMPkluMatrix->KLUmatrixRowCOO) ;
free (Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO) ;
free (Matrix->SMPkluMatrix->KLUmatrixTrashCOO) ;
Matrix->SMPkluMatrix->KLUmatrixAp = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAi = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAxComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO = NULL ;
Matrix->SMPkluMatrix->KLUmatrixColCOO = NULL ;
Matrix->SMPkluMatrix->KLUmatrixRowCOO = NULL ;
Matrix->SMPkluMatrix->KLUmatrixValueComplexCOO = NULL ;
Matrix->SMPkluMatrix->KLUmatrixTrashCOO = NULL ;
} else {
spDestroy (Matrix->SPmatrix) ;
}
}
#endif
/*
* SMPpreOrder()
*/
@ -547,6 +970,43 @@ SMPprint (SMPmatrix *Matrix, char *Filename)
}
}
#ifdef CIDER
void
SMPprintKLUforCIDER (SMPmatrix *Matrix, char *Filename)
{
unsigned int i ;
double *KLUmatrixAx ;
if (Matrix->CKTkluMODE)
{
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex)
{
klu_z_print (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAxComplex,
(int)Matrix->SMPkluMatrix->KLUmatrixN, NULL, NULL) ;
} else {
/* Allocate the Real Matrix */
KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof(double)) ;
/* Copy the Complex Matrix into the Real Matrix */
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
KLUmatrixAx [i] = Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] ;
}
/* Print the Real Matrix */
klu_print (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, KLUmatrixAx, (int)Matrix->SMPkluMatrix->KLUmatrixN, NULL, NULL) ;
/* Free the Real Matrix Storage */
free (KLUmatrixAx) ;
}
} else {
if (Filename)
spFileMatrix (Matrix->SPmatrix, Filename, "Circuit Matrix", 0, 1, 1) ;
else
spPrint (Matrix->SPmatrix, 0, 1, 1) ;
}
}
#endif
/*
* SMPgetError()
*/
@ -562,7 +1022,6 @@ SMPgetError (SMPmatrix *Matrix, int *Col, int *Row)
}
}
#ifdef KLU
void
spDeterminant_KLU (SMPmatrix *Matrix, int *pExponent, RealNumber *pDeterminant, RealNumber *piDeterminant)
{
@ -819,7 +1278,6 @@ spDeterminant_KLU (SMPmatrix *Matrix, int *pExponent, RealNumber *pDeterminant,
free (Ux) ;
free (Rs) ;
}
#endif
/*
* SMPcProdDiag()

37
src/maths/sparse/spsmp.c
View File

@ -174,6 +174,14 @@ SMPluFac(SMPmatrix *Matrix, double PivTol, double Gmin)
return spFactor( Matrix->SPmatrix );
}
#ifdef CIDER
int
SMPluFacForCIDER (SMPmatrix *Matrix)
{
return spFactor (Matrix->SPmatrix) ;
}
#endif
/*
* SMPcReorder()
*/
@ -225,6 +233,14 @@ SMPcSolve(SMPmatrix *Matrix, double RHS[], double iRHS[],
spSolve( Matrix->SPmatrix, RHS, RHS, iRHS, iRHS );
}
#ifdef CIDER
void
SMPcSolveForCIDER (SMPmatrix *Matrix, double RHS[], double RHSsolution[], double iRHS[], double iRHSsolution[])
{
spSolve (Matrix->SPmatrix, RHS, RHSsolution, iRHS, iRHSsolution) ;
}
#endif
/*
* SMPsolve()
*/
@ -236,6 +252,14 @@ SMPsolve(SMPmatrix *Matrix, double RHS[], double Spare[])
spSolve( Matrix->SPmatrix, RHS, RHS, NULL, NULL );
}
#ifdef CIDER
void
SMPsolveForCIDER (SMPmatrix *Matrix, double RHS[], double RHSsolution[])
{
spSolve (Matrix->SPmatrix, RHS, RHSsolution, NULL, NULL) ;
}
#endif
/*
* SMPmatSize()
*/
@ -256,6 +280,19 @@ SMPnewMatrix(SMPmatrix *Matrix, int size)
return Error;
}
#ifdef CIDER
int
SMPnewMatrixForCIDER (SMPmatrix *Matrix, int size, int complex)
{
int error ;
Matrix->SPmatrix = spCreate (size, complex, &error) ;
return error ;
}
#endif
/*
* SMPdestroy()
*/

Write Preview
Loading…
Cancel
Save