diff --git a/src/ciderlib/oned/onesolve.c b/src/ciderlib/oned/onesolve.c index 7825d946a..128cd06fc 100644 --- a/src/ciderlib/oned/onesolve.c +++ b/src/ciderlib/oned/onesolve.c @@ -1,7 +1,7 @@ /********** Copyright 1991 Regents of the University of California. All rights reserved. -Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group -Author: 1991 David A. Gates, U. C. Berkeley CAD Group +Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group +Author: 1991 David A. Gates, U. C. Berkeley CAD Group **********/ /* @@ -25,269 +25,280 @@ Author: 1991 David A. Gates, U. C. Berkeley CAD Group #include "ngspice/cktdefs.h" #include "ngspice/ftedefs.h" -extern IFfrontEnd *SPfrontEnd; +extern IFfrontEnd* SPfrontEnd; /* The iteration driving loop and convergence check */ void -ONEdcSolve(ONEdevice *pDevice, int iterationLimit, bool newSolver, - bool tranAnalysis, ONEtranInfo *info) +ONEdcSolve(ONEdevice* pDevice, int iterationLimit, bool newSolver, + bool tranAnalysis, ONEtranInfo* info) { - ONEnode *pNode; - ONEelem *pElem; - int index, eIndex, error; - int timesConverged = 0, negConc = FALSE; - int size = pDevice->numEqns; - bool quitLoop; - bool debug = FALSE; - double *rhs = pDevice->rhs; -/* double *intermediate = pDevice->copiedSolution; */ - double *solution = pDevice->dcSolution; - double *delta = pDevice->dcDeltaSolution; - double poissNorm, contNorm; - double startTime, totalStartTime; - double totalTime, loadTime, factorTime, solveTime, updateTime, checkTime; - double orderTime = 0.0; - - quitLoop = FALSE; - debug = (!tranAnalysis && ONEdcDebug) - ||(tranAnalysis && ONEtranDebug); - pDevice->iterationNumber = 0; - pDevice->converged = FALSE; - - - totalTime = loadTime = factorTime = solveTime = updateTime = checkTime = 0.0; - totalStartTime = SPfrontEnd->IFseconds(); - - if (debug) { - if (pDevice->poissonOnly) { - fprintf(stdout, "Equilibrium Solution:\n"); - } else { - fprintf(stdout, "Bias Solution:\n"); - } - fprintf(stdout, "Iteration RHS Norm\n"); - } - while (! (pDevice->converged - || pDevice->iterationNumber > iterationLimit - || quitLoop)) { - pDevice->iterationNumber++; - - if ((!pDevice->poissonOnly) && (iterationLimit > 0) - &&(!tranAnalysis)) { - ONEjacCheck(pDevice, tranAnalysis, info); - } - /* LOAD */ - startTime = SPfrontEnd->IFseconds(); - if (pDevice->poissonOnly) { - ONEQsysLoad(pDevice); - } else { - ONE_sysLoad(pDevice, tranAnalysis, info); - } - pDevice->rhsNorm = maxNorm(rhs, size); - loadTime += SPfrontEnd->IFseconds() - startTime; + ONEnode* pNode; + ONEelem* pElem; + int index, eIndex, error; + int timesConverged = 0, negConc = FALSE; + int size = pDevice->numEqns; + bool quitLoop; + bool debug = FALSE; + double* rhs = pDevice->rhs; + /* double *intermediate = pDevice->copiedSolution; */ + double* solution = pDevice->dcSolution; + double* delta = pDevice->dcDeltaSolution; + double poissNorm, contNorm; + double startTime, totalStartTime; + double totalTime, loadTime, factorTime, solveTime, updateTime, checkTime; + double orderTime = 0.0; + + quitLoop = FALSE; + debug = (!tranAnalysis && ONEdcDebug) + || (tranAnalysis && ONEtranDebug); + pDevice->iterationNumber = 0; + pDevice->converged = FALSE; + + + totalTime = loadTime = factorTime = solveTime = updateTime = checkTime = 0.0; + totalStartTime = SPfrontEnd->IFseconds(); + if (debug) { - fprintf(stdout, "%7d %11.4e%s\n", - pDevice->iterationNumber - 1, pDevice->rhsNorm, - negConc ? " negative conc encountered" : ""); - negConc = FALSE; + if (pDevice->poissonOnly) { + fprintf(stdout, "Equilibrium Solution:\n"); + } + else { + fprintf(stdout, "Bias Solution:\n"); + } + fprintf(stdout, "Iteration RHS Norm\n"); } - /* FACTOR */ - startTime = SPfrontEnd->IFseconds(); + while (!(pDevice->converged + || pDevice->iterationNumber > iterationLimit + || quitLoop)) { + pDevice->iterationNumber++; + + if ((!pDevice->poissonOnly) && (iterationLimit > 0) + && (!tranAnalysis)) { + ONEjacCheck(pDevice, tranAnalysis, info); + } + /* LOAD */ + startTime = SPfrontEnd->IFseconds(); + if (pDevice->poissonOnly) { + ONEQsysLoad(pDevice); + } + else { + ONE_sysLoad(pDevice, tranAnalysis, info); + } + pDevice->rhsNorm = maxNorm(rhs, size); + loadTime += SPfrontEnd->IFseconds() - startTime; + if (debug) { + fprintf(stdout, "%7d %11.4e%s\n", + pDevice->iterationNumber - 1, pDevice->rhsNorm, + negConc ? " negative conc encountered" : ""); + negConc = FALSE; + } + /* FACTOR */ + startTime = SPfrontEnd->IFseconds(); #ifdef KLU - error = SMPreorderKLUforCIDER (pDevice->matrix) ; + error = SMPreorderKLUforCIDER(pDevice->matrix); #else - error = SMPluFacForCIDER (pDevice->matrix) ; + error = SMPluFacForCIDER(pDevice->matrix); #endif - factorTime += SPfrontEnd->IFseconds() - startTime; - if (newSolver) { - if (pDevice->iterationNumber == 1) { - orderTime = factorTime; - } else if (pDevice->iterationNumber == 2) { - orderTime -= factorTime - orderTime; - factorTime -= orderTime; - if (pDevice->poissonOnly) { - pDevice->pStats->orderTime[STAT_SETUP] += orderTime; - } else { - pDevice->pStats->orderTime[STAT_DC] += orderTime; - } - newSolver = FALSE; - } - } - if (foundError(error)) { - if (error == spSINGULAR) { - int 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. */ - exit(-1); - } - /* SOLVE */ - startTime = SPfrontEnd->IFseconds(); + factorTime += SPfrontEnd->IFseconds() - startTime; + if (newSolver) { + if (pDevice->iterationNumber == 1) { + orderTime = factorTime; + } + else if (pDevice->iterationNumber == 2) { + orderTime -= factorTime - orderTime; + factorTime -= orderTime; + if (pDevice->poissonOnly) { + pDevice->pStats->orderTime[STAT_SETUP] += orderTime; + } + else { + pDevice->pStats->orderTime[STAT_DC] += orderTime; + } + newSolver = FALSE; + } + } + if (foundError(error)) { + if (error == spSINGULAR) { + int 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. */ + exit(-1); + } + /* SOLVE */ + startTime = SPfrontEnd->IFseconds(); #ifdef KLU - SMPsolveKLUforCIDER (pDevice->matrix, rhs, delta, NULL, NULL) ; + SMPsolveKLUforCIDER(pDevice->matrix, rhs, delta, NULL, NULL); #else - SMPsolveForCIDER (pDevice->matrix, rhs, delta) ; + SMPsolveForCIDER(pDevice->matrix, rhs, delta); #endif - solveTime += SPfrontEnd->IFseconds() - startTime; - - /* UPDATE */ - startTime = SPfrontEnd->IFseconds(); - /* - * Use norm reducing Newton method only for DC bias solutions. Since norm - * reducing can get trapped by numerical errors, turn it off when we are - * somewhat close to the solution. - */ - if ((!pDevice->poissonOnly) && (iterationLimit > 0) - && (!tranAnalysis) && (pDevice->rhsNorm > 1e-6)) { - error = ONEnewDelta(pDevice, tranAnalysis, info); - if (error) { - pDevice->converged = FALSE; - quitLoop = TRUE; - updateTime += SPfrontEnd->IFseconds() - startTime; - continue; - } - } - for (index = 1; index <= size; index++) { - solution[index] += delta[index]; - } - updateTime += SPfrontEnd->IFseconds() - startTime; - - /* CHECK CONVERGENCE */ - startTime = SPfrontEnd->IFseconds(); - /* Check if updates have gotten sufficiently small. */ - if (pDevice->iterationNumber != 1) { - /* - * pDevice->converged = ONEdeltaConverged(pDevice); - */ - pDevice->converged = ONEpsiDeltaConverged(pDevice, &negConc); + solveTime += SPfrontEnd->IFseconds() - startTime; + + /* UPDATE */ + startTime = SPfrontEnd->IFseconds(); + /* + * Use norm reducing Newton method only for DC bias solutions. Since norm + * reducing can get trapped by numerical errors, turn it off when we are + * somewhat close to the solution. + */ + if ((!pDevice->poissonOnly) && (iterationLimit > 0) + && (!tranAnalysis) && (pDevice->rhsNorm > 1e-6)) { + error = ONEnewDelta(pDevice, tranAnalysis, info); + if (error) { + pDevice->converged = FALSE; + quitLoop = TRUE; + updateTime += SPfrontEnd->IFseconds() - startTime; + continue; + } + } + for (index = 1; index <= size; index++) { + solution[index] += delta[index]; + } + updateTime += SPfrontEnd->IFseconds() - startTime; + + /* CHECK CONVERGENCE */ + startTime = SPfrontEnd->IFseconds(); + /* Check if updates have gotten sufficiently small. */ + if (pDevice->iterationNumber != 1) { + /* + * pDevice->converged = ONEdeltaConverged(pDevice); + */ + pDevice->converged = ONEpsiDeltaConverged(pDevice, &negConc); + } + /* Check if the rhs residual is smaller than abstol. */ + if (pDevice->converged && (!pDevice->poissonOnly) + && (!tranAnalysis)) { + ONE_rhsLoad(pDevice, tranAnalysis, info); + pDevice->rhsNorm = maxNorm(rhs, size); + if (pDevice->rhsNorm > pDevice->abstol) { + pDevice->converged = FALSE; + } + if ((++timesConverged >= 2) + && (pDevice->rhsNorm < 1e3 * pDevice->abstol)) { + pDevice->converged = TRUE; + } + else if (timesConverged >= 5) { + pDevice->converged = FALSE; + quitLoop = TRUE; + } + } + else if (pDevice->converged && pDevice->poissonOnly) { + ONEQrhsLoad(pDevice); + pDevice->rhsNorm = maxNorm(rhs, size); + if (pDevice->rhsNorm > pDevice->abstol) { + pDevice->converged = FALSE; + } + if (++timesConverged >= 5) { + pDevice->converged = TRUE; + } + } + /* Check if any of the carrier concentrations are negative. */ + if (pDevice->converged && (!pDevice->poissonOnly)) { + /* Clear garbage entry since carrier-free elements reference it */ + solution[0] = 0.0; + + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (index = 0; index <= 1; index++) { + if (pElem->evalNodes[index]) { + pNode = pElem->pNodes[index]; + if (solution[pNode->nEqn] < 0.0) { + pDevice->converged = FALSE; + negConc = TRUE; + if (tranAnalysis) { + quitLoop = TRUE; + } + else { + solution[pNode->nEqn] = 0.0; + } + } + if (solution[pNode->pEqn] < 0.0) { + pDevice->converged = FALSE; + negConc = TRUE; + if (tranAnalysis) { + quitLoop = TRUE; + } + else { + solution[pNode->pEqn] = 0.0; + } + } + } + } + } + /* Set to a consistent state if negative conc was encountered. */ + if (!pDevice->converged) { + ONE_rhsLoad(pDevice, tranAnalysis, info); + pDevice->rhsNorm = maxNorm(rhs, size); + } + } + checkTime += SPfrontEnd->IFseconds() - startTime; } - /* Check if the rhs residual is smaller than abstol. */ - if (pDevice->converged &&(!pDevice->poissonOnly) - &&(!tranAnalysis)) { - ONE_rhsLoad(pDevice, tranAnalysis, info); - pDevice->rhsNorm = maxNorm(rhs, size); - if (pDevice->rhsNorm > pDevice->abstol) { - pDevice->converged = FALSE; - } - if ((++timesConverged >= 2) - &&(pDevice->rhsNorm < 1e3 * pDevice->abstol)) { - pDevice->converged = TRUE; - } else if (timesConverged >= 5) { - pDevice->converged = FALSE; - quitLoop = TRUE; - } - } else if (pDevice->converged && pDevice->poissonOnly) { - ONEQrhsLoad(pDevice); - pDevice->rhsNorm = maxNorm(rhs, size); - if (pDevice->rhsNorm > pDevice->abstol) { - pDevice->converged = FALSE; - } - if (++timesConverged >= 5) { - pDevice->converged = TRUE; - } + totalTime += SPfrontEnd->IFseconds() - totalStartTime; + + if (tranAnalysis) { + pDevice->pStats->loadTime[STAT_TRAN] += loadTime; + pDevice->pStats->factorTime[STAT_TRAN] += factorTime; + pDevice->pStats->solveTime[STAT_TRAN] += solveTime; + pDevice->pStats->updateTime[STAT_TRAN] += updateTime; + pDevice->pStats->checkTime[STAT_TRAN] += checkTime; + pDevice->pStats->numIters[STAT_TRAN] += pDevice->iterationNumber; } - /* Check if any of the carrier concentrations are negative. */ - if (pDevice->converged &&(!pDevice->poissonOnly)) { - /* Clear garbage entry since carrier-free elements reference it */ - solution[0] = 0.0; - - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (index = 0; index <= 1; index++) { - if (pElem->evalNodes[index]) { - pNode = pElem->pNodes[index]; - if (solution[pNode->nEqn] < 0.0) { - pDevice->converged = FALSE; - negConc = TRUE; - if (tranAnalysis) { - quitLoop = TRUE; - } else { - solution[pNode->nEqn] = 0.0; - } - } - if (solution[pNode->pEqn] < 0.0) { - pDevice->converged = FALSE; - negConc = TRUE; - if (tranAnalysis) { - quitLoop = TRUE; - } else { - solution[pNode->pEqn] = 0.0; - } - } - } - } - } - /* Set to a consistent state if negative conc was encountered. */ - if (!pDevice->converged) { - ONE_rhsLoad(pDevice, tranAnalysis, info); - pDevice->rhsNorm = maxNorm(rhs, size); - } + else if (pDevice->poissonOnly) { + pDevice->pStats->loadTime[STAT_SETUP] += loadTime; + pDevice->pStats->factorTime[STAT_SETUP] += factorTime; + pDevice->pStats->solveTime[STAT_SETUP] += solveTime; + pDevice->pStats->updateTime[STAT_SETUP] += updateTime; + pDevice->pStats->checkTime[STAT_SETUP] += checkTime; + pDevice->pStats->numIters[STAT_SETUP] += pDevice->iterationNumber; } - checkTime += SPfrontEnd->IFseconds() - startTime; - } - totalTime += SPfrontEnd->IFseconds() - totalStartTime; - - if (tranAnalysis) { - pDevice->pStats->loadTime[STAT_TRAN] += loadTime; - pDevice->pStats->factorTime[STAT_TRAN] += factorTime; - pDevice->pStats->solveTime[STAT_TRAN] += solveTime; - pDevice->pStats->updateTime[STAT_TRAN] += updateTime; - pDevice->pStats->checkTime[STAT_TRAN] += checkTime; - pDevice->pStats->numIters[STAT_TRAN] += pDevice->iterationNumber; - } else if (pDevice->poissonOnly) { - pDevice->pStats->loadTime[STAT_SETUP] += loadTime; - pDevice->pStats->factorTime[STAT_SETUP] += factorTime; - pDevice->pStats->solveTime[STAT_SETUP] += solveTime; - pDevice->pStats->updateTime[STAT_SETUP] += updateTime; - pDevice->pStats->checkTime[STAT_SETUP] += checkTime; - pDevice->pStats->numIters[STAT_SETUP] += pDevice->iterationNumber; - } else { - pDevice->pStats->loadTime[STAT_DC] += loadTime; - pDevice->pStats->factorTime[STAT_DC] += factorTime; - pDevice->pStats->solveTime[STAT_DC] += solveTime; - pDevice->pStats->updateTime[STAT_DC] += updateTime; - pDevice->pStats->checkTime[STAT_DC] += checkTime; - pDevice->pStats->numIters[STAT_DC] += pDevice->iterationNumber; - } - - if (debug) { - if (!tranAnalysis) { - pDevice->rhsNorm = maxNorm(rhs, size); - fprintf(stdout, "%7d %11.4e%s\n", - pDevice->iterationNumber, pDevice->rhsNorm, - negConc ? " negative conc in solution" : ""); + else { + pDevice->pStats->loadTime[STAT_DC] += loadTime; + pDevice->pStats->factorTime[STAT_DC] += factorTime; + pDevice->pStats->solveTime[STAT_DC] += solveTime; + pDevice->pStats->updateTime[STAT_DC] += updateTime; + pDevice->pStats->checkTime[STAT_DC] += checkTime; + pDevice->pStats->numIters[STAT_DC] += pDevice->iterationNumber; } - if (pDevice->converged) { - if (!pDevice->poissonOnly) { - rhs[0] = 0.0; - poissNorm = contNorm = 0.0; - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (index = 0; index <= 1; index++) { - if (pElem->evalNodes[index]) { - pNode = pElem->pNodes[index]; - poissNorm = MAX(poissNorm, ABS(rhs[pNode->psiEqn])); - contNorm = MAX(contNorm, ABS(rhs[pNode->nEqn])); - contNorm = MAX(contNorm, ABS(rhs[pNode->pEqn])); - } - } - } - fprintf(stdout, - "Residual: %11.4e C/um^2 poisson, %11.4e A/um^2 continuity\n", - poissNorm * EpsNorm * ENorm * 1e-8, - contNorm * JNorm * 1e-8); - } else { - fprintf(stdout, "Residual: %11.4e C/um^2 poisson\n", - pDevice->rhsNorm * EpsNorm * ENorm * 1e-8); - } + + if (debug) { + if (!tranAnalysis) { + pDevice->rhsNorm = maxNorm(rhs, size); + fprintf(stdout, "%7d %11.4e%s\n", + pDevice->iterationNumber, pDevice->rhsNorm, + negConc ? " negative conc in solution" : ""); + } + if (pDevice->converged) { + if (!pDevice->poissonOnly) { + rhs[0] = 0.0; + poissNorm = contNorm = 0.0; + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (index = 0; index <= 1; index++) { + if (pElem->evalNodes[index]) { + pNode = pElem->pNodes[index]; + poissNorm = MAX(poissNorm, ABS(rhs[pNode->psiEqn])); + contNorm = MAX(contNorm, ABS(rhs[pNode->nEqn])); + contNorm = MAX(contNorm, ABS(rhs[pNode->pEqn])); + } + } + } + fprintf(stdout, + "Residual: %11.4e C/um^2 poisson, %11.4e A/um^2 continuity\n", + poissNorm * EpsNorm * ENorm * 1e-8, + contNorm * JNorm * 1e-8); + } + else { + fprintf(stdout, "Residual: %11.4e C/um^2 poisson\n", + pDevice->rhsNorm * EpsNorm * ENorm * 1e-8); + } + } } - } } /* @@ -296,108 +307,108 @@ ONEdcSolve(ONEdevice *pDevice, int iterationLimit, bool newSolver, * being used since we are always looking at potentials: (psi, phin, phip). */ bool -ONEpsiDeltaConverged(ONEdevice *pDevice, int *pNegConc) +ONEpsiDeltaConverged(ONEdevice* pDevice, int* pNegConc) { - int index, nIndex, eIndex; - ONEnode *pNode; - ONEelem *pElem; - double xOld, xNew, xDelta, tol; - double psi, newPsi, nConc, pConc, newN, newP; - double phiN, phiP, newPhiN, newPhiP; - bool converged = TRUE; - - /* Equilibrium solution. */ - if (pDevice->poissonOnly) { - for (index = 1; index <= pDevice->numEqns; index++) { - xOld = pDevice->dcSolution[index]; - xDelta = pDevice->dcDeltaSolution[index]; - xNew = xOld + xDelta; - tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew)); - if (ABS(xDelta) > tol) { - converged = FALSE; - goto done; - } + int index, nIndex, eIndex; + ONEnode* pNode; + ONEelem* pElem; + double xOld, xNew, xDelta, tol; + double psi, newPsi, nConc, pConc, newN, newP; + double phiN, phiP, newPhiN, newPhiP; + bool converged = TRUE; + + /* Equilibrium solution. */ + if (pDevice->poissonOnly) { + for (index = 1; index <= pDevice->numEqns; index++) { + xOld = pDevice->dcSolution[index]; + xDelta = pDevice->dcDeltaSolution[index]; + xNew = xOld + xDelta; + tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew)); + if (ABS(xDelta) > tol) { + converged = FALSE; + goto done; + } + } + return (converged); + } - return (converged); - - } - /* Bias solution. Check convergence on psi, phin, phip. */ - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (nIndex = 0; nIndex <= 1; nIndex++) { - if (pElem->evalNodes[nIndex]) { - pNode = pElem->pNodes[nIndex]; - if (pNode->nodeType != CONTACT) { - /* check convergence on psi */ - xOld = pDevice->dcSolution[pNode->psiEqn]; - xDelta = pDevice->dcDeltaSolution[pNode->psiEqn]; - xNew = xOld + xDelta; - tol = pDevice->abstol + - pDevice->reltol * MAX(ABS(xOld), ABS(xNew)); - if (ABS(xDelta) > tol) { - converged = FALSE; - goto done; - } - } - /* check convergence on phin and phip */ - if (pElem->elemType == SEMICON && pNode->nodeType != CONTACT) { - psi = pDevice->dcSolution[pNode->psiEqn]; - nConc = pDevice->dcSolution[pNode->nEqn]; - pConc = pDevice->dcSolution[pNode->pEqn]; - newPsi = psi + pDevice->dcDeltaSolution[pNode->psiEqn]; - newN = nConc + pDevice->dcDeltaSolution[pNode->nEqn]; - newP = pConc + pDevice->dcDeltaSolution[pNode->pEqn]; - if (newN <= 0.0 || newP <= 0.0) { - *pNegConc = TRUE; - converged = FALSE; - goto done; - } - phiN = psi - log(nConc / pNode->nie); - phiP = psi + log(pConc / pNode->nie); - newPhiN = newPsi - log(newN / pNode->nie); - newPhiP = newPsi + log(newP / pNode->nie); - tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiN), ABS(newPhiN)); - if (ABS(newPhiN - phiN) > tol) { - converged = FALSE; - goto done; - } - tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiP), ABS(newPhiP)); - if (ABS(newPhiP - phiP) > tol) { - converged = FALSE; - goto done; - } - } - } + /* Bias solution. Check convergence on psi, phin, phip. */ + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (nIndex = 0; nIndex <= 1; nIndex++) { + if (pElem->evalNodes[nIndex]) { + pNode = pElem->pNodes[nIndex]; + if (pNode->nodeType != CONTACT) { + /* check convergence on psi */ + xOld = pDevice->dcSolution[pNode->psiEqn]; + xDelta = pDevice->dcDeltaSolution[pNode->psiEqn]; + xNew = xOld + xDelta; + tol = pDevice->abstol + + pDevice->reltol * MAX(ABS(xOld), ABS(xNew)); + if (ABS(xDelta) > tol) { + converged = FALSE; + goto done; + } + } + /* check convergence on phin and phip */ + if (pElem->elemType == SEMICON && pNode->nodeType != CONTACT) { + psi = pDevice->dcSolution[pNode->psiEqn]; + nConc = pDevice->dcSolution[pNode->nEqn]; + pConc = pDevice->dcSolution[pNode->pEqn]; + newPsi = psi + pDevice->dcDeltaSolution[pNode->psiEqn]; + newN = nConc + pDevice->dcDeltaSolution[pNode->nEqn]; + newP = pConc + pDevice->dcDeltaSolution[pNode->pEqn]; + if (newN <= 0.0 || newP <= 0.0) { + *pNegConc = TRUE; + converged = FALSE; + goto done; + } + phiN = psi - log(nConc / pNode->nie); + phiP = psi + log(pConc / pNode->nie); + newPhiN = newPsi - log(newN / pNode->nie); + newPhiP = newPsi + log(newP / pNode->nie); + tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiN), ABS(newPhiN)); + if (ABS(newPhiN - phiN) > tol) { + converged = FALSE; + goto done; + } + tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiP), ABS(newPhiP)); + if (ABS(newPhiP - phiP) > tol) { + converged = FALSE; + goto done; + } + } + } + } } - } done: - return (converged); + return (converged); } /* * See if the update to the solution is small enough. Returns TRUE if it is. */ bool -ONEdeltaConverged(ONEdevice *pDevice) +ONEdeltaConverged(ONEdevice* pDevice) { - int index; - bool converged = TRUE; - double *solution = pDevice->dcSolution; - double *delta = pDevice->dcDeltaSolution; - double xOld, xNew, tol; - - - for (index = 1; index <= pDevice->numEqns; index++) { - xOld = solution[index]; - xNew = xOld + delta[index]; - tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew)); - if (ABS(xOld - xNew) > tol) { - converged = FALSE; - break; + int index; + bool converged = TRUE; + double* solution = pDevice->dcSolution; + double* delta = pDevice->dcDeltaSolution; + double xOld, xNew, tol; + + + for (index = 1; index <= pDevice->numEqns; index++) { + xOld = solution[index]; + xNew = xOld + delta[index]; + tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew)); + if (ABS(xOld - xNew) > tol) { + converged = FALSE; + break; + } } - } - return (converged); + return (converged); } /* @@ -405,44 +416,44 @@ ONEdeltaConverged(ONEdevice *pDevice) * physically reasonable. Returns TRUE if it is. */ bool -ONEdeviceConverged(ONEdevice *pDevice) +ONEdeviceConverged(ONEdevice* pDevice) { - int index, eIndex; - bool converged = TRUE; - double *solution = pDevice->dcSolution; - ONEnode *pNode; - ONEelem *pElem; - double startTime; - - /* - * If the update is sufficently small, and the carrier concentrations are - * all positive, then return TRUE, else return FALSE. - */ - - - /* CHECK CONVERGENCE */ - startTime = SPfrontEnd->IFseconds(); - if ((converged = ONEdeltaConverged(pDevice)) == TRUE) { - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (index = 0; index <= 1; index++) { - if (pElem->evalNodes[index]) { - pNode = pElem->pNodes[index]; - if (pNode->nEqn != 0 && solution[pNode->nEqn] < 0.0) { - converged = FALSE; - solution[pNode->nEqn] = pNode->nConc = 0.0; - } - if (pNode->pEqn != 0 && solution[pNode->pEqn] < 0.0) { - converged = FALSE; - solution[pNode->pEqn] = pNode->pConc = 0.0; - } - } - } + int index, eIndex; + bool converged = TRUE; + double* solution = pDevice->dcSolution; + ONEnode* pNode; + ONEelem* pElem; + double startTime; + + /* + * If the update is sufficently small, and the carrier concentrations are + * all positive, then return TRUE, else return FALSE. + */ + + + /* CHECK CONVERGENCE */ + startTime = SPfrontEnd->IFseconds(); + if ((converged = ONEdeltaConverged(pDevice)) == TRUE) { + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (index = 0; index <= 1; index++) { + if (pElem->evalNodes[index]) { + pNode = pElem->pNodes[index]; + if (pNode->nEqn != 0 && solution[pNode->nEqn] < 0.0) { + converged = FALSE; + solution[pNode->nEqn] = pNode->nConc = 0.0; + } + if (pNode->pEqn != 0 && solution[pNode->pEqn] < 0.0) { + converged = FALSE; + solution[pNode->pEqn] = pNode->pConc = 0.0; + } + } + } + } } - } - pDevice->pStats->checkTime[STAT_TRAN] += SPfrontEnd->IFseconds() - startTime; + pDevice->pStats->checkTime[STAT_TRAN] += SPfrontEnd->IFseconds() - startTime; - return (converged); + return (converged); } /* @@ -450,22 +461,22 @@ ONEdeviceConverged(ONEdevice *pDevice) * solution. */ void -ONEresetJacobian(ONEdevice *pDevice) +ONEresetJacobian(ONEdevice* pDevice) { - int error; - - - ONE_jacLoad(pDevice); + int error; + + + ONE_jacLoad(pDevice); #ifdef KLU - error = SMPreorderKLUforCIDER (pDevice->matrix) ; + error = SMPreorderKLUforCIDER(pDevice->matrix); #else - error = SMPluFacForCIDER (pDevice->matrix) ; + error = SMPluFacForCIDER(pDevice->matrix); #endif - if (foundError(error)) { - exit(-1); - } + if (foundError(error)) { + exit(-1); + } } /* @@ -474,75 +485,76 @@ ONEresetJacobian(ONEdevice *pDevice) * at half the insulator band gap (refPsi). */ void -ONEstoreNeutralGuess(ONEdevice *pDevice) +ONEstoreNeutralGuess(ONEdevice* pDevice) { - int nIndex, eIndex; - ONEelem *pElem; - ONEnode *pNode; - double nie, conc, absConc, refPsi, psi, ni, pi, sign; - - - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - refPsi = pElem->matlInfo->refPsi; - if (pElem->elemType == INSULATOR) { - for (nIndex = 0; nIndex <= 1; nIndex++) { - if (pElem->evalNodes[nIndex]) { - pNode = pElem->pNodes[nIndex]; - if (pNode->nodeType == CONTACT) { - /* - * A metal contact to insulating domain, so use work function - * difference as the value of psi. - */ - pNode->psi = RefPsi - pNode->eaff; - } else { - pNode->psi = refPsi; - } - } - } - } - if (pElem->elemType == SEMICON) { - for (nIndex = 0; nIndex <= 1; nIndex++) { - if (pElem->evalNodes[nIndex]) { - pNode = pElem->pNodes[nIndex]; - nie = pNode->nie; - conc = pNode->netConc / nie; - psi = 0.0; - ni = nie; - pi = nie; - sign = SGN(conc); - absConc = ABS(conc); - if (conc != 0.0) { - psi = sign * log(0.5 * absConc + sqrt(1.0 + 0.25 * absConc * absConc)); - ni = nie * exp(psi); - pi = nie * exp(-psi); - if (FreezeOut) { - /* Use Newton's Method to solve for psi. */ - int ctr, maxiter = 10; - double rhs, deriv, fNa, fNd, fdNa, fdNd; - for (ctr = 0; ctr < maxiter; ctr++) { - pNode->nConc = ni; - pNode->pConc = pi; - ONEQfreezeOut(pNode, &fNd, &fNa, &fdNd, &fdNa); - rhs = pi - ni + pNode->nd * fNd - pNode->na * fNa; - deriv = pi + ni - pNode->nd * fdNd + pNode->na * fdNa; - psi += rhs / deriv; - ni = nie * exp(psi); - pi = nie * exp(-psi); - } - } - } - pNode->psi = refPsi + psi; - pNode->psi0 = pNode->psi; - pNode->vbe = refPsi; - pNode->nConc = ni; - pNode->pConc = pi; - /* Now store the initial guess in the dc solution vector. */ - pDevice->dcSolution[pNode->poiEqn] = pNode->psi; - } - } + int nIndex, eIndex; + ONEelem* pElem; + ONEnode* pNode; + double nie, conc, absConc, refPsi, psi, ni, pi, sign; + + + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + refPsi = pElem->matlInfo->refPsi; + if (pElem->elemType == INSULATOR) { + for (nIndex = 0; nIndex <= 1; nIndex++) { + if (pElem->evalNodes[nIndex]) { + pNode = pElem->pNodes[nIndex]; + if (pNode->nodeType == CONTACT) { + /* + * A metal contact to insulating domain, so use work function + * difference as the value of psi. + */ + pNode->psi = RefPsi - pNode->eaff; + } + else { + pNode->psi = refPsi; + } + } + } + } + if (pElem->elemType == SEMICON) { + for (nIndex = 0; nIndex <= 1; nIndex++) { + if (pElem->evalNodes[nIndex]) { + pNode = pElem->pNodes[nIndex]; + nie = pNode->nie; + conc = pNode->netConc / nie; + psi = 0.0; + ni = nie; + pi = nie; + sign = SGN(conc); + absConc = ABS(conc); + if (conc != 0.0) { + psi = sign * log(0.5 * absConc + sqrt(1.0 + 0.25 * absConc * absConc)); + ni = nie * exp(psi); + pi = nie * exp(-psi); + if (FreezeOut) { + /* Use Newton's Method to solve for psi. */ + int ctr, maxiter = 10; + double rhs, deriv, fNa, fNd, fdNa, fdNd; + for (ctr = 0; ctr < maxiter; ctr++) { + pNode->nConc = ni; + pNode->pConc = pi; + ONEQfreezeOut(pNode, &fNd, &fNa, &fdNd, &fdNa); + rhs = pi - ni + pNode->nd * fNd - pNode->na * fNa; + deriv = pi + ni - pNode->nd * fdNd + pNode->na * fdNa; + psi += rhs / deriv; + ni = nie * exp(psi); + pi = nie * exp(-psi); + } + } + } + pNode->psi = refPsi + psi; + pNode->psi0 = pNode->psi; + pNode->vbe = refPsi; + pNode->nConc = ni; + pNode->pConc = pi; + /* Now store the initial guess in the dc solution vector. */ + pDevice->dcSolution[pNode->poiEqn] = pNode->psi; + } + } + } } - } } /* @@ -551,154 +563,158 @@ ONEstoreNeutralGuess(ONEdevice *pDevice) * taken as an initial guess. */ void -ONEequilSolve(ONEdevice *pDevice) +ONEequilSolve(ONEdevice* pDevice) { - bool newSolver = FALSE; - int error; - int nIndex, eIndex; - ONEelem *pElem; - ONEnode *pNode; - double startTime, setupTime, miscTime; - - setupTime = miscTime = 0.0; - - /* SETUP */ - startTime = SPfrontEnd->IFseconds(); - switch (pDevice->solverType) { - case SLV_SMSIG: - case SLV_BIAS: - /* free up memory allocated for the bias solution */ - FREE(pDevice->dcSolution); - FREE(pDevice->dcDeltaSolution); - FREE(pDevice->copiedSolution); - FREE(pDevice->rhs); - FREE(pDevice->rhsImag); + bool newSolver = FALSE; + int error; + int nIndex, eIndex; + ONEelem* pElem; + ONEnode* pNode; + double startTime, setupTime, miscTime; + + setupTime = miscTime = 0.0; + + /* SETUP */ + startTime = SPfrontEnd->IFseconds(); + switch (pDevice->solverType) { + case SLV_SMSIG: + case SLV_BIAS: + /* free up memory allocated for the bias solution */ + FREE(pDevice->dcSolution); + FREE(pDevice->dcDeltaSolution); + FREE(pDevice->copiedSolution); + FREE(pDevice->rhs); + FREE(pDevice->rhsImag); #ifdef KLU - SMPdestroyKLUforCIDER (pDevice->matrix) ; + SMPdestroyKLUforCIDER(pDevice->matrix); #else - SMPdestroy (pDevice->matrix) ; + SMPdestroy(pDevice->matrix); #endif - FREE (pDevice->matrix) ; + FREE(pDevice->matrix); - /* FALLTHRU */ - case SLV_NONE: - pDevice->poissonOnly = TRUE; - pDevice->numEqns = pDevice->dimEquil - 1; - XCALLOC(pDevice->dcSolution, double, pDevice->dimEquil); - XCALLOC(pDevice->dcDeltaSolution, double, pDevice->dimEquil); - XCALLOC(pDevice->copiedSolution, double, pDevice->dimEquil); - XCALLOC(pDevice->rhs, double, pDevice->dimEquil); + /* FALLTHRU */ + case SLV_NONE: + pDevice->poissonOnly = TRUE; + pDevice->numEqns = pDevice->dimEquil - 1; + XCALLOC(pDevice->dcSolution, double, pDevice->dimEquil); + XCALLOC(pDevice->dcDeltaSolution, double, pDevice->dimEquil); + XCALLOC(pDevice->copiedSolution, double, pDevice->dimEquil); + XCALLOC(pDevice->rhs, double, pDevice->dimEquil); - pDevice->matrix = TMALLOC (SMPmatrix, 1) ; + pDevice->matrix = TMALLOC(SMPmatrix, 1); #ifdef KLU - pDevice->matrix->CKTkluMODE = ft_curckt->ci_ckt->CKTkluMODE ; - error = SMPnewMatrixKLUforCIDER (pDevice->matrix, pDevice->numEqns, KLUmatrixReal) ; + pDevice->matrix->CKTkluMODE = ft_curckt->ci_ckt->CKTkluMODE; + error = SMPnewMatrixKLUforCIDER(pDevice->matrix, pDevice->numEqns, KLUmatrixReal); #else - error = SMPnewMatrixForCIDER (pDevice->matrix, pDevice->numEqns, 0) ; + error = SMPnewMatrixForCIDER(pDevice->matrix, pDevice->numEqns, 0); #endif - if (error == spNO_MEMORY) { - fprintf(stderr, "Error: ONEequilSolve: Out of Memory\n"); - exit(-1); - } - newSolver = TRUE; + if (error == spNO_MEMORY) { + fprintf(stderr, "Error: ONEequilSolve: Out of Memory\n"); + exit(-1); + } + newSolver = TRUE; #ifdef KLU - if (pDevice->matrix->CKTkluMODE) { - pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUmatrixReal ; - } else { + if (pDevice->matrix->CKTkluMODE) { + pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUmatrixReal; + } + else { #endif - spSetReal (pDevice->matrix->SPmatrix) ; + spSetReal(pDevice->matrix->SPmatrix); #ifdef KLU - } + } #endif - ONEQjacBuild(pDevice); + ONEQjacBuild(pDevice); #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) { - fprintf (stderr, "Error: CIDER: KLU Failed\n") ; - if (pDevice->matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) { - return ; // Francesco Lannutti - Fix KLU return values + 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) { + fprintf(stderr, "Error: CIDER: KLU Failed\n"); + if (pDevice->matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) { + return; // Francesco Lannutti - Fix KLU return values + } + } + pDevice->numOrigEquil = (int)pDevice->matrix->SMPkluMatrix->KLUmatrixNZ; + } + else { + pDevice->numOrigEquil = spElementCount(pDevice->matrix->SPmatrix); } - } - pDevice->numOrigEquil = (int)pDevice->matrix->SMPkluMatrix->KLUmatrixNZ ; - } else { - pDevice->numOrigEquil = spElementCount (pDevice->matrix->SPmatrix) ; - } #else - pDevice->numOrigEquil = spElementCount (pDevice->matrix->SPmatrix) ; + pDevice->numOrigEquil = spElementCount(pDevice->matrix->SPmatrix); #endif - pDevice->numFillEquil = 0; - /* FALLTHRU */ - case SLV_EQUIL: - pDevice->solverType = SLV_EQUIL; - break; - default: - fprintf(stderr, "Error: Panic: Unknown solver type in equil solution.\n"); - exit(-1); - break; - } - ONEstoreNeutralGuess(pDevice); - setupTime += SPfrontEnd->IFseconds() - startTime; - - /* SOLVE */ - ONEdcSolve(pDevice, MaxIterations, newSolver, FALSE, NULL); - - /* MISCELLANEOUS */ - startTime = SPfrontEnd->IFseconds(); - if (newSolver) { + pDevice->numFillEquil = 0; + /* FALLTHRU */ + case SLV_EQUIL: + pDevice->solverType = SLV_EQUIL; + break; + default: + fprintf(stderr, "Error: Panic: Unknown solver type in equil solution.\n"); + exit(-1); + break; + } + ONEstoreNeutralGuess(pDevice); + setupTime += SPfrontEnd->IFseconds() - startTime; + + /* SOLVE */ + ONEdcSolve(pDevice, MaxIterations, newSolver, FALSE, NULL); + + /* MISCELLANEOUS */ + startTime = SPfrontEnd->IFseconds(); + if (newSolver) { #ifdef KLU - if (pDevice->matrix->CKTkluMODE) { - pDevice->numFillEquil = pDevice->matrix->SMPkluMatrix->KLUmatrixNumeric->lnz + pDevice->matrix->SMPkluMatrix->KLUmatrixNumeric->unz - - (int)pDevice->matrix->SMPkluMatrix->KLUmatrixNZ ; - } else { + if (pDevice->matrix->CKTkluMODE) { + pDevice->numFillEquil = pDevice->matrix->SMPkluMatrix->KLUmatrixNumeric->lnz + pDevice->matrix->SMPkluMatrix->KLUmatrixNumeric->unz + - (int)pDevice->matrix->SMPkluMatrix->KLUmatrixNZ; + } + else { #endif - pDevice->numFillEquil = spFillinCount(pDevice->matrix->SPmatrix); + pDevice->numFillEquil = spFillinCount(pDevice->matrix->SPmatrix); #ifdef KLU - } + } #endif - } - if (pDevice->converged) { - ONEQcommonTerms(pDevice); - - /* Save equilibrium potential. */ - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (nIndex = 0; nIndex <= 1; nIndex++) { - if (pElem->evalNodes[nIndex]) { - pNode = pElem->pNodes[nIndex]; - pNode->psi0 = pNode->psi; - } - } } - } else { - fprintf(stderr, "Error: ONEequilSolve: No Convergence\n"); - } - miscTime += SPfrontEnd->IFseconds() - startTime; - pDevice->pStats->setupTime[STAT_SETUP] += setupTime; - pDevice->pStats->miscTime[STAT_SETUP] += miscTime; + if (pDevice->converged) { + ONEQcommonTerms(pDevice); + + /* Save equilibrium potential. */ + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (nIndex = 0; nIndex <= 1; nIndex++) { + if (pElem->evalNodes[nIndex]) { + pNode = pElem->pNodes[nIndex]; + pNode->psi0 = pNode->psi; + } + } + } + } + else { + fprintf(stderr, "Error: ONEequilSolve: No Convergence\n"); + } + miscTime += SPfrontEnd->IFseconds() - startTime; + pDevice->pStats->setupTime[STAT_SETUP] += setupTime; + pDevice->pStats->miscTime[STAT_SETUP] += miscTime; } /* @@ -706,464 +722,474 @@ ONEequilSolve(ONEdevice *pDevice) * is taken as an initial guess the first time this is called. */ void -ONEbiasSolve(ONEdevice *pDevice, int iterationLimit, - bool tranAnalysis, ONEtranInfo *info) +ONEbiasSolve(ONEdevice* pDevice, int iterationLimit, + bool tranAnalysis, ONEtranInfo* info) { - bool newSolver = FALSE; - int error; - int index, eIndex; - double *solution; - ONEelem *pElem; - ONEnode *pNode; - double startTime, setupTime, miscTime; + bool newSolver = FALSE; + int error; + int index, eIndex; + double* solution; + ONEelem* pElem; + ONEnode* pNode; + double startTime, setupTime, miscTime; - setupTime = miscTime = 0.0; + setupTime = miscTime = 0.0; - /* SETUP */ - startTime = SPfrontEnd->IFseconds(); - switch (pDevice->solverType) { - case SLV_EQUIL: - /* Free up the vectors allocated in the equilibrium solution. */ - FREE(pDevice->dcSolution); - FREE(pDevice->dcDeltaSolution); - FREE(pDevice->copiedSolution); - FREE(pDevice->rhs); + /* SETUP */ + startTime = SPfrontEnd->IFseconds(); + switch (pDevice->solverType) { + case SLV_EQUIL: + /* Free up the vectors allocated in the equilibrium solution. */ + FREE(pDevice->dcSolution); + FREE(pDevice->dcDeltaSolution); + FREE(pDevice->copiedSolution); + FREE(pDevice->rhs); #ifdef KLU - SMPdestroyKLUforCIDER (pDevice->matrix) ; + SMPdestroyKLUforCIDER(pDevice->matrix); #else - SMPdestroy (pDevice->matrix) ; + SMPdestroy(pDevice->matrix); #endif - FREE (pDevice->matrix) ; + FREE(pDevice->matrix); - /* FALLTHRU */ - case SLV_NONE: - pDevice->poissonOnly = FALSE; - pDevice->numEqns = pDevice->dimBias - 1; - XCALLOC(pDevice->dcSolution, double, pDevice->dimBias); - XCALLOC(pDevice->dcDeltaSolution, double, pDevice->dimBias); - XCALLOC(pDevice->copiedSolution, double, pDevice->dimBias); - XCALLOC(pDevice->rhs, double, pDevice->dimBias); - XCALLOC(pDevice->rhsImag, double, pDevice->dimBias); + /* FALLTHRU */ + case SLV_NONE: + pDevice->poissonOnly = FALSE; + pDevice->numEqns = pDevice->dimBias - 1; + XCALLOC(pDevice->dcSolution, double, pDevice->dimBias); + XCALLOC(pDevice->dcDeltaSolution, double, pDevice->dimBias); + XCALLOC(pDevice->copiedSolution, double, pDevice->dimBias); + XCALLOC(pDevice->rhs, double, pDevice->dimBias); + XCALLOC(pDevice->rhsImag, double, pDevice->dimBias); - pDevice->matrix = TMALLOC (SMPmatrix, 1) ; + pDevice->matrix = TMALLOC(SMPmatrix, 1); #ifdef KLU - pDevice->matrix->CKTkluMODE = ft_curckt->ci_ckt->CKTkluMODE ; - error = SMPnewMatrixKLUforCIDER (pDevice->matrix, pDevice->numEqns, KLUMatrixComplex) ; + pDevice->matrix->CKTkluMODE = ft_curckt->ci_ckt->CKTkluMODE; + error = SMPnewMatrixKLUforCIDER(pDevice->matrix, pDevice->numEqns, KLUMatrixComplex); #else - error = SMPnewMatrixForCIDER (pDevice->matrix, pDevice->numEqns, 1) ; + error = SMPnewMatrixForCIDER(pDevice->matrix, pDevice->numEqns, 1); #endif - if (error == spNO_MEMORY) { - exit(-1); - } - newSolver = TRUE; - ONE_jacBuild(pDevice); + if (error == spNO_MEMORY) { + exit(-1); + } + newSolver = TRUE; + ONE_jacBuild(pDevice); #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) { - fprintf (stderr, "Error: CIDER: KLU failed\n") ; - return ; // Francesco Lannutti - Fix KLU return values + 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) { + fprintf(stderr, "Error: CIDER: KLU failed\n"); + return; // Francesco Lannutti - Fix KLU return values + } + } + pDevice->numOrigBias = (int)pDevice->matrix->SMPkluMatrix->KLUmatrixNZ; + } + else { + pDevice->numOrigBias = spElementCount(pDevice->matrix->SPmatrix); } - } - pDevice->numOrigBias = (int)pDevice->matrix->SMPkluMatrix->KLUmatrixNZ ; - } else { - pDevice->numOrigBias = spElementCount(pDevice->matrix->SPmatrix); - } #else - pDevice->numOrigBias = spElementCount (pDevice->matrix->SPmatrix) ; + pDevice->numOrigBias = spElementCount(pDevice->matrix->SPmatrix); #endif - pDevice->numFillBias = 0; - ONEstoreInitialGuess(pDevice); - /* FALLTHRU */ - case SLV_SMSIG: + pDevice->numFillBias = 0; + ONEstoreInitialGuess(pDevice); + /* FALLTHRU */ + case SLV_SMSIG: #ifdef KLU - if (pDevice->matrix->CKTkluMODE) { - pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUmatrixReal ; - } else { + if (pDevice->matrix->CKTkluMODE) { + pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUmatrixReal; + } + else { #endif - spSetReal (pDevice->matrix->SPmatrix) ; + spSetReal(pDevice->matrix->SPmatrix); #ifdef KLU - } + } #endif - /* FALLTHRU */ - case SLV_BIAS: - pDevice->solverType = SLV_BIAS; - break; - default: - fprintf(stderr, "Error: Panic: Unknown solver type in bias solution.\n"); - exit(-1); - break; - } - setupTime += SPfrontEnd->IFseconds() - startTime; + /* FALLTHRU */ + case SLV_BIAS: + pDevice->solverType = SLV_BIAS; + break; + default: + fprintf(stderr, "Error: Panic: Unknown solver type in bias solution.\n"); + exit(-1); + break; + } + setupTime += SPfrontEnd->IFseconds() - startTime; - /* SOLVE */ - ONEdcSolve(pDevice, iterationLimit, newSolver, tranAnalysis, info); + /* SOLVE */ + ONEdcSolve(pDevice, iterationLimit, newSolver, tranAnalysis, info); - /* MISCELLANEOUS */ - startTime = SPfrontEnd->IFseconds(); - if (newSolver) { + /* MISCELLANEOUS */ + startTime = SPfrontEnd->IFseconds(); + if (newSolver) { #ifdef KLU - if (pDevice->matrix->CKTkluMODE) { - pDevice->numFillBias = pDevice->matrix->SMPkluMatrix->KLUmatrixNumeric->lnz + pDevice->matrix->SMPkluMatrix->KLUmatrixNumeric->unz - - (int)pDevice->matrix->SMPkluMatrix->KLUmatrixNZ ; - } else { + if (pDevice->matrix->CKTkluMODE) { + pDevice->numFillBias = pDevice->matrix->SMPkluMatrix->KLUmatrixNumeric->lnz + pDevice->matrix->SMPkluMatrix->KLUmatrixNumeric->unz + - (int)pDevice->matrix->SMPkluMatrix->KLUmatrixNZ; + } + else { #endif - pDevice->numFillBias = spFillinCount (pDevice->matrix->SPmatrix) ; // Francesco Lannutti - Fix for KLU + pDevice->numFillBias = spFillinCount(pDevice->matrix->SPmatrix); // Francesco Lannutti - Fix for KLU #ifdef KLU - } + } #endif - } - solution = pDevice->dcSolution; - if ((!pDevice->converged) && iterationLimit > 1) { - } else if (pDevice->converged) { - /* Update the nodal quantities. */ - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (index = 0; index <= 1; index++) { - if (pElem->evalNodes[index]) { - pNode = pElem->pNodes[index]; - if (pNode->psiEqn != 0) { - pNode->psi = solution[pNode->psiEqn]; - } - if (pNode->nEqn != 0) { - pNode->nConc = solution[pNode->nEqn]; - } - if (pNode->pEqn != 0) { - pNode->pConc = solution[pNode->pEqn]; - } - } - } } - /* Update the current terms. */ - ONE_commonTerms(pDevice, FALSE, tranAnalysis, info); - } else if (iterationLimit <= 1) { - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (index = 0; index <= 1; index++) { - if (pElem->evalNodes[index]) { - pNode = pElem->pNodes[index]; - if (pNode->nodeType != CONTACT) { - pNode->psi = solution[pNode->psiEqn]; - pDevice->devState0 [pNode->nodePsi] = pNode->psi; - if (pElem->elemType == SEMICON) { - pNode->nConc = solution[pNode->nEqn]; - pNode->pConc = solution[pNode->pEqn]; - pDevice->devState0 [pNode->nodeN] = pNode->nConc; - pDevice->devState0 [pNode->nodeP] = pNode->pConc; - } - } - } - } + solution = pDevice->dcSolution; + if ((!pDevice->converged) && iterationLimit > 1) { + } + else if (pDevice->converged) { + /* Update the nodal quantities. */ + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (index = 0; index <= 1; index++) { + if (pElem->evalNodes[index]) { + pNode = pElem->pNodes[index]; + if (pNode->psiEqn != 0) { + pNode->psi = solution[pNode->psiEqn]; + } + if (pNode->nEqn != 0) { + pNode->nConc = solution[pNode->nEqn]; + } + if (pNode->pEqn != 0) { + pNode->pConc = solution[pNode->pEqn]; + } + } + } + } + /* Update the current terms. */ + ONE_commonTerms(pDevice, FALSE, tranAnalysis, info); + } + else if (iterationLimit <= 1) { + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (index = 0; index <= 1; index++) { + if (pElem->evalNodes[index]) { + pNode = pElem->pNodes[index]; + if (pNode->nodeType != CONTACT) { + pNode->psi = solution[pNode->psiEqn]; + pDevice->devState0[pNode->nodePsi] = pNode->psi; + if (pElem->elemType == SEMICON) { + pNode->nConc = solution[pNode->nEqn]; + pNode->pConc = solution[pNode->pEqn]; + pDevice->devState0[pNode->nodeN] = pNode->nConc; + pDevice->devState0[pNode->nodeP] = pNode->pConc; + } + } + } + } + } + } + miscTime += SPfrontEnd->IFseconds() - startTime; + if (tranAnalysis) { + pDevice->pStats->setupTime[STAT_TRAN] += setupTime; + pDevice->pStats->miscTime[STAT_TRAN] += miscTime; + } + else { + pDevice->pStats->setupTime[STAT_DC] += setupTime; + pDevice->pStats->miscTime[STAT_DC] += miscTime; } - } - miscTime += SPfrontEnd->IFseconds() - startTime; - if (tranAnalysis) { - pDevice->pStats->setupTime[STAT_TRAN] += setupTime; - pDevice->pStats->miscTime[STAT_TRAN] += miscTime; - } else { - pDevice->pStats->setupTime[STAT_DC] += setupTime; - pDevice->pStats->miscTime[STAT_DC] += miscTime; - } } /* Copy the device's equilibrium state into the solution vector. */ void -ONEstoreEquilibGuess(ONEdevice *pDevice) +ONEstoreEquilibGuess(ONEdevice* pDevice) { - int nIndex, eIndex; - double *solution = pDevice->dcSolution; - double refPsi; - ONEelem *pElem; - ONEnode *pNode; - - - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - refPsi = pElem->matlInfo->refPsi; - for (nIndex = 0; nIndex <= 1; nIndex++) { - if (pElem->evalNodes[nIndex]) { - pNode = pElem->pNodes[nIndex]; - if (pNode->nodeType != CONTACT) { - solution[pNode->psiEqn] = pNode->psi0; - if (pElem->elemType == SEMICON) { - solution[pNode->nEqn] = pNode->nie * exp(pNode->psi0 - refPsi); - solution[pNode->pEqn] = pNode->nie * exp(-pNode->psi0 + refPsi); - } - } - } + int nIndex, eIndex; + double* solution = pDevice->dcSolution; + double refPsi; + ONEelem* pElem; + ONEnode* pNode; + + + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + refPsi = pElem->matlInfo->refPsi; + for (nIndex = 0; nIndex <= 1; nIndex++) { + if (pElem->evalNodes[nIndex]) { + pNode = pElem->pNodes[nIndex]; + if (pNode->nodeType != CONTACT) { + solution[pNode->psiEqn] = pNode->psi0; + if (pElem->elemType == SEMICON) { + solution[pNode->nEqn] = pNode->nie * exp(pNode->psi0 - refPsi); + solution[pNode->pEqn] = pNode->nie * exp(-pNode->psi0 + refPsi); + } + } + } + } } - } } /* Copy the device's internal state into the solution vector. */ void -ONEstoreInitialGuess(ONEdevice *pDevice) +ONEstoreInitialGuess(ONEdevice* pDevice) { - int nIndex, eIndex; - double *solution = pDevice->dcSolution; - ONEelem *pElem; - ONEnode *pNode; - - - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (nIndex = 0; nIndex <= 1; nIndex++) { - if (pElem->evalNodes[nIndex]) { - pNode = pElem->pNodes[nIndex]; - if (pNode->nodeType != CONTACT) { - solution[pNode->psiEqn] = pNode->psi; - if (pElem->elemType == SEMICON) { - solution[pNode->nEqn] = pNode->nConc; - solution[pNode->pEqn] = pNode->pConc; - } - } - } + int nIndex, eIndex; + double* solution = pDevice->dcSolution; + ONEelem* pElem; + ONEnode* pNode; + + + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (nIndex = 0; nIndex <= 1; nIndex++) { + if (pElem->evalNodes[nIndex]) { + pNode = pElem->pNodes[nIndex]; + if (pNode->nodeType != CONTACT) { + solution[pNode->psiEqn] = pNode->psi; + if (pElem->elemType == SEMICON) { + solution[pNode->nEqn] = pNode->nConc; + solution[pNode->pEqn] = pNode->pConc; + } + } + } + } } - } } int -ONEnewDelta(ONEdevice *pDevice, bool tranAnalysis, ONEtranInfo *info) +ONEnewDelta(ONEdevice* pDevice, bool tranAnalysis, ONEtranInfo* info) { - int index, iterNum; - double newNorm, fib, lambda, fibn, fibp; - bool acceptable = FALSE, error = FALSE; - iterNum = 0; - lambda = 1.0; - fibn = 1.0; - fibp = 1.0; - - - /* - * Copy the contents of dcSolution into copiedSolution and modify - * dcSolution by adding the deltaSolution. - */ - for (index = 1; index <= pDevice->numEqns; index++) { - pDevice->copiedSolution[index] = pDevice->dcSolution[index]; - pDevice->dcSolution[index] += pDevice->dcDeltaSolution[index]; - } - - if (pDevice->poissonOnly) { - ONEQrhsLoad(pDevice); - } else { - ONE_rhsLoad(pDevice, tranAnalysis, info); - } - newNorm = maxNorm(pDevice->rhs, pDevice->numEqns); - if (pDevice->rhsNorm <= pDevice->abstol) { - lambda = 0.0; - newNorm = pDevice->rhsNorm; - } else if (newNorm < pDevice->rhsNorm) { - acceptable = TRUE; - } else { - /* chop the step size so that deltax is acceptable */ - if (ONEdcDebug) { - fprintf(stdout, " %11.4e %11.4e\n", - newNorm, lambda); + int index, iterNum; + double newNorm, fib, lambda, fibn, fibp; + bool acceptable = FALSE, error = FALSE; + iterNum = 0; + lambda = 1.0; + fibn = 1.0; + fibp = 1.0; + + + /* + * Copy the contents of dcSolution into copiedSolution and modify + * dcSolution by adding the deltaSolution. + */ + for (index = 1; index <= pDevice->numEqns; index++) { + pDevice->copiedSolution[index] = pDevice->dcSolution[index]; + pDevice->dcSolution[index] += pDevice->dcDeltaSolution[index]; + } + + if (pDevice->poissonOnly) { + ONEQrhsLoad(pDevice); + } + else { + ONE_rhsLoad(pDevice, tranAnalysis, info); + } + newNorm = maxNorm(pDevice->rhs, pDevice->numEqns); + if (pDevice->rhsNorm <= pDevice->abstol) { + lambda = 0.0; + newNorm = pDevice->rhsNorm; + } + else if (newNorm < pDevice->rhsNorm) { + acceptable = TRUE; } - while (!acceptable) { - iterNum++; - - if (iterNum > NORM_RED_MAXITERS) { - error = TRUE; - lambda = 0.0; - /* Don't break out until after we've reset the device. */ - } - fib = fibp; - fibp = fibn; - fibn += fib; - lambda *= (fibp / fibn); - - for (index = 1; index <= pDevice->numEqns; index++) { - pDevice->dcSolution[index] = pDevice->copiedSolution[index] - + lambda * pDevice->dcDeltaSolution[index]; - } - if (pDevice->poissonOnly) { - ONEQrhsLoad(pDevice); - } else { - ONE_rhsLoad(pDevice, tranAnalysis, info); - } - newNorm = maxNorm(pDevice->rhs, pDevice->numEqns); - if (error) { - break; - } - if (newNorm <= pDevice->rhsNorm) { - acceptable = TRUE; - } - if (ONEdcDebug) { - fprintf(stdout, " %11.4e %11.4e\n", - newNorm, lambda); - } + else { + /* chop the step size so that deltax is acceptable */ + if (ONEdcDebug) { + fprintf(stdout, " %11.4e %11.4e\n", + newNorm, lambda); + } + while (!acceptable) { + iterNum++; + + if (iterNum > NORM_RED_MAXITERS) { + error = TRUE; + lambda = 0.0; + /* Don't break out until after we've reset the device. */ + } + fib = fibp; + fibp = fibn; + fibn += fib; + lambda *= (fibp / fibn); + + for (index = 1; index <= pDevice->numEqns; index++) { + pDevice->dcSolution[index] = pDevice->copiedSolution[index] + + lambda * pDevice->dcDeltaSolution[index]; + } + if (pDevice->poissonOnly) { + ONEQrhsLoad(pDevice); + } + else { + ONE_rhsLoad(pDevice, tranAnalysis, info); + } + newNorm = maxNorm(pDevice->rhs, pDevice->numEqns); + if (error) { + break; + } + if (newNorm <= pDevice->rhsNorm) { + acceptable = TRUE; + } + if (ONEdcDebug) { + fprintf(stdout, " %11.4e %11.4e\n", + newNorm, lambda); + } + } + } + /* Restore the previous dcSolution and store the new deltaSolution. */ + pDevice->rhsNorm = newNorm; + for (index = 1; index <= pDevice->numEqns; index++) { + pDevice->dcSolution[index] = pDevice->copiedSolution[index]; + pDevice->dcDeltaSolution[index] *= lambda; } - } - /* Restore the previous dcSolution and store the new deltaSolution. */ - pDevice->rhsNorm = newNorm; - for (index = 1; index <= pDevice->numEqns; index++) { - pDevice->dcSolution[index] = pDevice->copiedSolution[index]; - pDevice->dcDeltaSolution[index] *= lambda; - } - return(error); + return(error); } /* Predict the values of the internal variables at the next timepoint. */ /* Needed for Predictor-Corrector LTE estimation */ void -ONEpredict(ONEdevice *pDevice, ONEtranInfo *info) +ONEpredict(ONEdevice* pDevice, ONEtranInfo* info) { - int nIndex, eIndex; - ONEnode *pNode; - ONEelem *pElem; - double startTime, miscTime = 0.0; - - - /* TRANSIENT MISC */ - startTime = SPfrontEnd->IFseconds(); - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (nIndex = 0; nIndex <= 1; nIndex++) { - if (pElem->evalNodes[nIndex]) { - pNode = pElem->pNodes[nIndex]; - pNode->psi = pDevice->devState1 [pNode->nodePsi]; - if (pElem->elemType == SEMICON && pNode->nodeType != CONTACT) { - pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN); - pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP); - pNode->nConc = pNode->nPred; - pNode->pConc = pNode->pPred; - } - } + int nIndex, eIndex; + ONEnode* pNode; + ONEelem* pElem; + double startTime, miscTime = 0.0; + + + /* TRANSIENT MISC */ + startTime = SPfrontEnd->IFseconds(); + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (nIndex = 0; nIndex <= 1; nIndex++) { + if (pElem->evalNodes[nIndex]) { + pNode = pElem->pNodes[nIndex]; + pNode->psi = pDevice->devState1[pNode->nodePsi]; + if (pElem->elemType == SEMICON && pNode->nodeType != CONTACT) { + pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN); + pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP); + pNode->nConc = pNode->nPred; + pNode->pConc = pNode->pPred; + } + } + } } - } - miscTime += SPfrontEnd->IFseconds() - startTime; - pDevice->pStats->miscTime[STAT_TRAN] += miscTime; + miscTime += SPfrontEnd->IFseconds() - startTime; + pDevice->pStats->miscTime[STAT_TRAN] += miscTime; } /* Estimate the device's overall truncation error. */ double -ONEtrunc(ONEdevice *pDevice, ONEtranInfo *info, double delta) +ONEtrunc(ONEdevice* pDevice, ONEtranInfo* info, double delta) { - int nIndex, eIndex; - ONEelem *pElem; - ONEnode *pNode; - double tolN, tolP, lte, relError, temp; - double lteCoeff = info->lteCoeff; - double mult = 10.0; - double reltol; - double startTime, lteTime = 0.0; - - - /* TRANSIENT LTE */ - startTime = SPfrontEnd->IFseconds(); - - relError = 0.0; - reltol = pDevice->reltol * mult; - - /* Need to get the predictor for the current order. */ - /* The scheme currently used is very dumb. Need to fix later. */ - /* XXX Why is the scheme dumb? Never understood this. */ - computePredCoeff(info->method, info->order, info->predCoeff, info->delta); - - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (nIndex = 0; nIndex <= 1; nIndex++) { - if (pElem->evalNodes[nIndex]) { - pNode = pElem->pNodes[nIndex]; - if (pElem->elemType == SEMICON && pNode->nodeType != CONTACT) { - tolN = pDevice->abstol + reltol * ABS(pNode->nConc); - tolP = pDevice->abstol + reltol * ABS(pNode->pConc); - pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN); - pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP); - lte = lteCoeff * (pNode->nConc - pNode->nPred); - temp = lte / tolN; - relError += temp * temp; - lte = lteCoeff * (pNode->pConc - pNode->pPred); - temp = lte / tolP; - relError += temp * temp; - } - } + int nIndex, eIndex; + ONEelem* pElem; + ONEnode* pNode; + double tolN, tolP, lte, relError, temp; + double lteCoeff = info->lteCoeff; + double mult = 10.0; + double reltol; + double startTime, lteTime = 0.0; + + + /* TRANSIENT LTE */ + startTime = SPfrontEnd->IFseconds(); + + relError = 0.0; + reltol = pDevice->reltol * mult; + + /* Need to get the predictor for the current order. */ + /* The scheme currently used is very dumb. Need to fix later. */ + /* XXX Why is the scheme dumb? Never understood this. */ + computePredCoeff(info->method, info->order, info->predCoeff, info->delta); + + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (nIndex = 0; nIndex <= 1; nIndex++) { + if (pElem->evalNodes[nIndex]) { + pNode = pElem->pNodes[nIndex]; + if (pElem->elemType == SEMICON && pNode->nodeType != CONTACT) { + tolN = pDevice->abstol + reltol * ABS(pNode->nConc); + tolP = pDevice->abstol + reltol * ABS(pNode->pConc); + pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN); + pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP); + lte = lteCoeff * (pNode->nConc - pNode->nPred); + temp = lte / tolN; + relError += temp * temp; + lte = lteCoeff * (pNode->pConc - pNode->pPred); + temp = lte / tolP; + relError += temp * temp; + } + } + } } - } - /* Make sure error is non-zero. */ - relError = MAX(pDevice->abstol, relError); + /* Make sure error is non-zero. */ + relError = MAX(pDevice->abstol, relError); - /* The total relative error has been calculated norm-2 squared. */ - relError = sqrt(relError / pDevice->numEqns); + /* The total relative error has been calculated norm-2 squared. */ + relError = sqrt(relError / pDevice->numEqns); - /* Use the order of the integration method to compute new delta. */ - temp = delta / pow(relError, 1.0 / (info->order + 1)); + /* Use the order of the integration method to compute new delta. */ + temp = delta / pow(relError, 1.0 / (info->order + 1)); - lteTime += SPfrontEnd->IFseconds() - startTime; - pDevice->pStats->lteTime += lteTime; + lteTime += SPfrontEnd->IFseconds() - startTime; + pDevice->pStats->lteTime += lteTime; - return (temp); + return (temp); } /* Save info from state table into the internal state. */ void -ONEsaveState(ONEdevice *pDevice) +ONEsaveState(ONEdevice* pDevice) { - int nIndex, eIndex; - ONEnode *pNode; - ONEelem *pElem; - - for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { - pElem = pDevice->elemArray[eIndex]; - for (nIndex = 0; nIndex <= 1; nIndex++) { - if (pElem->evalNodes[nIndex]) { - pNode = pElem->pNodes[nIndex]; - pNode->psi = pDevice->devState1 [pNode->nodePsi]; - if (pElem->elemType == SEMICON && pNode->nodeType != CONTACT) { - pNode->nConc = pDevice->devState1 [pNode->nodeN]; - pNode->pConc = pDevice->devState1 [pNode->nodeP]; - } - } + int nIndex, eIndex; + ONEnode* pNode; + ONEelem* pElem; + + for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) { + pElem = pDevice->elemArray[eIndex]; + for (nIndex = 0; nIndex <= 1; nIndex++) { + if (pElem->evalNodes[nIndex]) { + pNode = pElem->pNodes[nIndex]; + pNode->psi = pDevice->devState1[pNode->nodePsi]; + if (pElem->elemType == SEMICON && pNode->nodeType != CONTACT) { + pNode->nConc = pDevice->devState1[pNode->nodeN]; + pNode->pConc = pDevice->devState1[pNode->nodeP]; + } + } + } } - } } /* Function to compute Nu norm of the rhs vector. */ /* Nu-norm calculation based upon work done at Stanford. */ double -ONEnuNorm(ONEdevice *pDevice) +ONEnuNorm(ONEdevice* pDevice) { - /* The LU Decomposed matrix is available. Use it to calculate x. */ + /* The LU Decomposed matrix is available. Use it to calculate x. */ #ifdef KLU - printf ("CIDER: KLU to be fixed ONEnuNorm\n") ; - SMPsolveKLUforCIDER (pDevice->matrix, pDevice->rhs, pDevice->rhsImag, NULL, NULL) ; + 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) ; + SMPsolveForCIDER(pDevice->matrix, pDevice->rhs, pDevice->rhsImag); #endif - /* Compute L2-norm of the solution vector (stored in rhsImag) */ - return (l2Norm(pDevice->rhsImag, pDevice->numEqns)); + /* Compute L2-norm of the solution vector (stored in rhsImag) */ + return (l2Norm(pDevice->rhsImag, pDevice->numEqns)); } @@ -1172,63 +1198,64 @@ ONEnuNorm(ONEdevice *pDevice) * models are being incorporated. */ void -ONEjacCheck(ONEdevice *pDevice, bool tranAnalysis, ONEtranInfo *info) +ONEjacCheck(ONEdevice* pDevice, bool tranAnalysis, ONEtranInfo* info) { - int index, rIndex; - double del, diff, tol, *dptr; - - - if (ONEjacDebug) { - ONE_sysLoad(pDevice, tranAnalysis, info); - /* - * spPrint( pDevice->matrix, 0, 1, 1 ); - */ - pDevice->rhsNorm = maxNorm(pDevice->rhs, pDevice->numEqns); - for (index = 1; index <= pDevice->numEqns; index++) { - if (1e3 * ABS(pDevice->rhs[index]) > pDevice->rhsNorm) { - fprintf(stderr, "eqn %d: res %11.4e, norm %11.4e\n", - index, pDevice->rhs[index], pDevice->rhsNorm); - } - } - for (index = 1; index <= pDevice->numEqns; index++) { - pDevice->rhsImag[index] = pDevice->rhs[index]; - } - for (index = 1; index <= pDevice->numEqns; index++) { - pDevice->copiedSolution[index] = pDevice->dcSolution[index]; - del = 1e-4 * pDevice->abstol + 1e-6 * ABS(pDevice->dcSolution[index]); - pDevice->dcSolution[index] += del; - ONE_rhsLoad(pDevice, tranAnalysis, info); - pDevice->dcSolution[index] = pDevice->copiedSolution[index]; - for (rIndex = 1; rIndex <= pDevice->numEqns; rIndex++) { - diff = (pDevice->rhsImag[rIndex] - pDevice->rhs[rIndex]) / del; + int index, rIndex; + double del, diff, tol, * dptr; + + + if (ONEjacDebug) { + ONE_sysLoad(pDevice, tranAnalysis, info); + /* + * spPrint( pDevice->matrix, 0, 1, 1 ); + */ + pDevice->rhsNorm = maxNorm(pDevice->rhs, pDevice->numEqns); + for (index = 1; index <= pDevice->numEqns; index++) { + if (1e3 * ABS(pDevice->rhs[index]) > pDevice->rhsNorm) { + fprintf(stderr, "eqn %d: res %11.4e, norm %11.4e\n", + index, pDevice->rhs[index], pDevice->rhsNorm); + } + } + for (index = 1; index <= pDevice->numEqns; index++) { + pDevice->rhsImag[index] = pDevice->rhs[index]; + } + for (index = 1; index <= pDevice->numEqns; index++) { + pDevice->copiedSolution[index] = pDevice->dcSolution[index]; + del = 1e-4 * pDevice->abstol + 1e-6 * ABS(pDevice->dcSolution[index]); + pDevice->dcSolution[index] += del; + ONE_rhsLoad(pDevice, tranAnalysis, info); + pDevice->dcSolution[index] = pDevice->copiedSolution[index]; + for (rIndex = 1; rIndex <= pDevice->numEqns; rIndex++) { + diff = (pDevice->rhsImag[rIndex] - pDevice->rhs[rIndex]) / del; #ifdef KLU - printf ("CIDER: KLU to be fixed: spFindElement") ; - dptr = NULL ; + printf("CIDER: KLU to be fixed: spFindElement"); + dptr = NULL; #else - dptr = spFindElement(pDevice->matrix->SPmatrix, rIndex, index); // Francesco Lannutti - Fix for KLU + 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 { - * fprintf( stderr, "[%d][%d]: FD = %11.4e, AJ = %11.4e\n", rIndex, - * index, diff, 0.0 ); } - */ - if (dptr != NULL) { - tol = (1e-4 * pDevice->abstol) + (1e-2 * MAX(ABS(diff), ABS(*dptr))); - if ((diff != 0.0) && (ABS(diff - *dptr) > tol)) { - fprintf(stderr, "Mismatch[%d][%d]: FD = %11.4e, AJ = %11.4e\n\t FD-AJ = %11.4e vs. %11.4e\n", - rIndex, index, diff, *dptr, - ABS(diff - *dptr), tol); - } - } else { - if (diff != 0.0) { - fprintf(stderr, "Missing [%d][%d]: FD = %11.4e, AJ = 0.0\n", - rIndex, index, diff); - } - } - } + /* + * if ( dptr ISNOT NULL ) { fprintf( stderr, "[%d][%d]: FD = + * %11.4e, AJ = %11.4e\n", rIndex, index, diff, *dptr ); } else { + * fprintf( stderr, "[%d][%d]: FD = %11.4e, AJ = %11.4e\n", rIndex, + * index, diff, 0.0 ); } + */ + if (dptr != NULL) { + tol = (1e-4 * pDevice->abstol) + (1e-2 * MAX(ABS(diff), ABS(*dptr))); + if ((diff != 0.0) && (ABS(diff - *dptr) > tol)) { + fprintf(stderr, "Mismatch[%d][%d]: FD = %11.4e, AJ = %11.4e\n\t FD-AJ = %11.4e vs. %11.4e\n", + rIndex, index, diff, *dptr, + ABS(diff - *dptr), tol); + } + } + else { + if (diff != 0.0) { + fprintf(stderr, "Missing [%d][%d]: FD = %11.4e, AJ = 0.0\n", + rIndex, index, diff); + } + } + } + } } - } }