1 /**********
2 Copyright 1991 Regents of the University of California. All rights reserved.
3 Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
4 Author: 1991 David A. Gates, U. C. Berkeley CAD Group
5 **********/
6
7 #include "../../maths/misc/norm.h"
8 #include "ngspice/bool.h"
9 #include "ngspice/cidersupt.h"
10 #include "ngspice/cpextern.h"
11 #include "ngspice/ifsim.h"
12 #include "ngspice/macros.h"
13 #include "ngspice/ngspice.h"
14 #include "ngspice/numenum.h"
15 #include "ngspice/numglobs.h"
16 #include "ngspice/spmatrix.h"
17 #include "ngspice/twodev.h"
18 #include "ngspice/twomesh.h"
19 #include "twoddefs.h"
20 #include "twodext.h"
21 extern IFfrontEnd *SPfrontEnd;
22
23
24 /* functions to calculate the 2D solutions */
25
26
27 /* the iteration driving loop and convergence check */
28 void
TWOdcSolve(TWOdevice * pDevice,int iterationLimit,BOOLEAN newSolver,BOOLEAN tranAnalysis,TWOtranInfo * info)29 TWOdcSolve(TWOdevice *pDevice, int iterationLimit, BOOLEAN newSolver,
30 BOOLEAN tranAnalysis, TWOtranInfo *info)
31 {
32 TWOnode *pNode;
33 TWOelem *pElem;
34 int size = pDevice->numEqns;
35 int index, eIndex, error;
36 int timesConverged = 0;
37 BOOLEAN quitLoop;
38 BOOLEAN debug;
39 BOOLEAN negConc = FALSE;
40 double *rhs = pDevice->rhs;
41 double *solution = pDevice->dcSolution;
42 double *delta = pDevice->dcDeltaSolution;
43 double poissNorm, contNorm;
44 double startTime, totalStartTime;
45 double totalTime, loadTime, factorTime, solveTime, updateTime, checkTime;
46 double orderTime = 0.0;
47
48 totalTime = loadTime = factorTime = solveTime = updateTime = checkTime = 0.0;
49 totalStartTime = SPfrontEnd->IFseconds();
50
51 quitLoop = FALSE;
52 debug = (!tranAnalysis && TWOdcDebug) || (tranAnalysis && TWOtranDebug);
53 pDevice->iterationNumber = 0;
54 pDevice->converged = FALSE;
55
56 if (debug) {
57 if (pDevice->poissonOnly) {
58 fprintf(stdout, "Equilibrium Solution:\n");
59 } else {
60 fprintf(stdout, "Bias Solution:\n");
61 }
62 fprintf(stdout, "Iteration RHS Norm\n");
63 }
64 while (!(pDevice->converged || (pDevice->iterationNumber > iterationLimit)
65 || quitLoop)) {
66 pDevice->iterationNumber++;
67
68 if ((!pDevice->poissonOnly) && (iterationLimit > 0)
69 &&(!tranAnalysis)) {
70 TWOjacCheck(pDevice, tranAnalysis, info);
71 }
72
73 /* LOAD */
74 startTime = SPfrontEnd->IFseconds();
75 if (pDevice->poissonOnly) {
76 TWOQsysLoad(pDevice);
77 } else if (!OneCarrier) {
78 TWO_sysLoad(pDevice, tranAnalysis, info);
79 } else if (OneCarrier == N_TYPE) {
80 TWONsysLoad(pDevice, tranAnalysis, info);
81 } else if (OneCarrier == P_TYPE) {
82 TWOPsysLoad(pDevice, tranAnalysis, info);
83 }
84 pDevice->rhsNorm = maxNorm(rhs, size);
85 loadTime += SPfrontEnd->IFseconds() - startTime;
86 if (debug) {
87 fprintf(stdout, "%7d %11.4e%s\n",
88 pDevice->iterationNumber - 1, pDevice->rhsNorm,
89 negConc ? " negative conc encountered" : "");
90 negConc = FALSE;
91 }
92
93 /* FACTOR */
94 startTime = SPfrontEnd->IFseconds();
95 error = spFactor(pDevice->matrix);
96 factorTime += SPfrontEnd->IFseconds() - startTime;
97 if (newSolver) {
98 if (pDevice->iterationNumber == 1) {
99 orderTime = factorTime;
100 }
101 else if (pDevice->iterationNumber == 2) {
102 orderTime -= factorTime - orderTime;
103 factorTime -= orderTime;
104 if (pDevice->poissonOnly) {
105 pDevice->pStats->orderTime[STAT_SETUP] += orderTime;
106 }
107 else {
108 pDevice->pStats->orderTime[STAT_DC] += orderTime;
109 }
110 /* After first two iterations, no special handling for a
111 * new solver */
112 newSolver = FALSE;
113 } /* end of case of iteratin 2 */
114 } /* end of special processing for a new solver */
115
116 if (foundError(error)) {
117 if (error == spSINGULAR) {
118 int badRow, badCol;
119 spWhereSingular(pDevice->matrix, &badRow, &badCol);
120 printf("***** singular at (%d,%d)\n", badRow, badCol);
121 }
122 pDevice->converged = FALSE;
123 quitLoop = TRUE;
124 continue;
125 }
126
127 /* SOLVE */
128 startTime = SPfrontEnd->IFseconds();
129 spSolve(pDevice->matrix, rhs, delta, NULL, NULL);
130 solveTime += SPfrontEnd->IFseconds() - startTime;
131
132 /* UPDATE */
133 startTime = SPfrontEnd->IFseconds();
134 /* Use norm reducing Newton method for DC bias solutions only. */
135 if ((!pDevice->poissonOnly) && (iterationLimit > 0)
136 && (!tranAnalysis) && (pDevice->rhsNorm > 1e-1)) {
137 error = TWOnewDelta(pDevice, tranAnalysis, info);
138 if (error) {
139 pDevice->converged = FALSE;
140 quitLoop = TRUE;
141 updateTime += SPfrontEnd->IFseconds() - startTime;
142 continue;
143 }
144 }
145 for (index = 1; index <= size; index++) {
146 solution[index] += delta[index];
147 }
148 updateTime += SPfrontEnd->IFseconds() - startTime;
149
150 /* CHECK CONVERGENCE */
151 startTime = SPfrontEnd->IFseconds();
152 /* Check if updates have gotten sufficiently small. */
153 if (pDevice->iterationNumber != 1) {
154 pDevice->converged = TWOdeltaConverged(pDevice);
155 }
156 /* Check if the residual is smaller than abstol. */
157 if (pDevice->converged && (!pDevice->poissonOnly)
158 && (!tranAnalysis)) {
159 if (!OneCarrier) {
160 TWO_rhsLoad(pDevice, tranAnalysis, info);
161 } else if (OneCarrier == N_TYPE) {
162 TWONrhsLoad(pDevice, tranAnalysis, info);
163 } else if (OneCarrier == P_TYPE) {
164 TWOPrhsLoad(pDevice, tranAnalysis, info);
165 }
166 pDevice->rhsNorm = maxNorm(rhs, size);
167 if (pDevice->rhsNorm > pDevice->abstol) {
168 pDevice->converged = FALSE;
169 }
170 if ((++timesConverged >= 2)
171 && (pDevice->rhsNorm < 1e3 * pDevice->abstol)) {
172 pDevice->converged = TRUE;
173 } else if (timesConverged >= 5) {
174 pDevice->converged = FALSE;
175 quitLoop = TRUE;
176 continue;
177 }
178 } else if (pDevice->converged && pDevice->poissonOnly) {
179 TWOQrhsLoad(pDevice);
180 pDevice->rhsNorm = maxNorm(rhs, size);
181 if (pDevice->rhsNorm > pDevice->abstol) {
182 pDevice->converged = FALSE;
183 }
184 if (++timesConverged >= 5) {
185 pDevice->converged = TRUE;
186 }
187 }
188 /* Check if the carrier concentrations are negative. */
189 if (pDevice->converged && (!pDevice->poissonOnly)) {
190 /* Clear garbage entry since carrier-free elements reference it. */
191 solution[0] = 0.0;
192 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
193 pElem = pDevice->elements[eIndex];
194 for (index = 0; index <= 3; index++) {
195 if (pElem->evalNodes[index]) {
196 pNode = pElem->pNodes[index];
197 if (solution[pNode->nEqn] < 0.0) {
198 pDevice->converged = FALSE;
199 negConc = TRUE;
200 if (tranAnalysis) {
201 quitLoop = TRUE;
202 } else {
203 solution[pNode->nEqn] = 0.0;
204 }
205 }
206 if (solution[pNode->pEqn] < 0.0) {
207 pDevice->converged = FALSE;
208 negConc = TRUE;
209 if (tranAnalysis) {
210 quitLoop = TRUE;
211 } else {
212 solution[pNode->pEqn] = 0.0;
213 }
214 }
215 }
216 }
217 }
218 if (!pDevice->converged) {
219 if (!OneCarrier) {
220 TWO_rhsLoad(pDevice, tranAnalysis, info);
221 } else if (OneCarrier == N_TYPE) {
222 TWONrhsLoad(pDevice, tranAnalysis, info);
223 } else if (OneCarrier == P_TYPE) {
224 TWOPrhsLoad(pDevice, tranAnalysis, info);
225 }
226 pDevice->rhsNorm = maxNorm(rhs, size);
227 }
228 }
229 checkTime += SPfrontEnd->IFseconds() - startTime;
230 }
231 totalTime += SPfrontEnd->IFseconds() - totalStartTime;
232
233 if (tranAnalysis) {
234 pDevice->pStats->loadTime[STAT_TRAN] += loadTime;
235 pDevice->pStats->factorTime[STAT_TRAN] += factorTime;
236 pDevice->pStats->solveTime[STAT_TRAN] += solveTime;
237 pDevice->pStats->updateTime[STAT_TRAN] += updateTime;
238 pDevice->pStats->checkTime[STAT_TRAN] += checkTime;
239 pDevice->pStats->numIters[STAT_TRAN] += pDevice->iterationNumber;
240 } else if (pDevice->poissonOnly) {
241 pDevice->pStats->loadTime[STAT_SETUP] += loadTime;
242 pDevice->pStats->factorTime[STAT_SETUP] += factorTime;
243 pDevice->pStats->solveTime[STAT_SETUP] += solveTime;
244 pDevice->pStats->updateTime[STAT_SETUP] += updateTime;
245 pDevice->pStats->checkTime[STAT_SETUP] += checkTime;
246 pDevice->pStats->numIters[STAT_SETUP] += pDevice->iterationNumber;
247 } else {
248 pDevice->pStats->loadTime[STAT_DC] += loadTime;
249 pDevice->pStats->factorTime[STAT_DC] += factorTime;
250 pDevice->pStats->solveTime[STAT_DC] += solveTime;
251 pDevice->pStats->updateTime[STAT_DC] += updateTime;
252 pDevice->pStats->checkTime[STAT_DC] += checkTime;
253 pDevice->pStats->numIters[STAT_DC] += pDevice->iterationNumber;
254 }
255
256 if (debug) {
257 if (!tranAnalysis) {
258 pDevice->rhsNorm = maxNorm(rhs, size);
259 fprintf(stdout, "%7d %11.4e%s\n",
260 pDevice->iterationNumber, pDevice->rhsNorm,
261 negConc ? " negative conc in solution" : "");
262 }
263 if (pDevice->converged) {
264 if (!pDevice->poissonOnly) {
265 poissNorm = contNorm = 0.0;
266 rhs[0] = 0.0; /* Make sure garbage entry is clear. */
267 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
268 pElem = pDevice->elements[eIndex];
269 for (index = 0; index <= 3; index++) {
270 if (pElem->evalNodes[index]) {
271 pNode = pElem->pNodes[index];
272 poissNorm = MAX(poissNorm,ABS(rhs[pNode->psiEqn]));
273 contNorm = MAX(contNorm,ABS(rhs[pNode->nEqn]));
274 contNorm = MAX(contNorm,ABS(rhs[pNode->pEqn]));
275 }
276 }
277 }
278 fprintf(stdout,
279 "Residual: %11.4e C/um poisson, %11.4e A/um continuity\n",
280 poissNorm * EpsNorm * VNorm * 1e-4,
281 contNorm * JNorm * LNorm * 1e-4);
282 } else {
283 fprintf(stdout, "Residual: %11.4e C/um poisson\n",
284 pDevice->rhsNorm * EpsNorm * VNorm * 1e-4);
285 }
286 }
287 }
288 }
289
290 BOOLEAN
TWOdeltaConverged(TWOdevice * pDevice)291 TWOdeltaConverged(TWOdevice *pDevice)
292 {
293 /* This function returns a TRUE if converged else a FALSE. */
294 int index;
295 BOOLEAN converged = TRUE;
296 double xOld, xNew, tol;
297
298 for (index = 1; index <= pDevice->numEqns; index++) {
299 xOld = pDevice->dcSolution[index];
300 xNew = xOld + pDevice->dcDeltaSolution[index];
301 tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew));
302 if (ABS(xOld - xNew) > tol) {
303 converged = FALSE;
304 break;
305 }
306 }
307 return (converged);
308 }
309
310 BOOLEAN
TWOdeviceConverged(TWOdevice * pDevice)311 TWOdeviceConverged(TWOdevice *pDevice)
312 {
313 int index, eIndex;
314 BOOLEAN converged = TRUE;
315 double *solution = pDevice->dcSolution;
316 TWOnode *pNode;
317 TWOelem *pElem;
318 double startTime;
319
320 /* If the update is sufficiently small, and the carrier concentrations
321 * are all positive, then return TRUE, else return FALSE.
322 */
323
324 /* CHECK CONVERGENCE */
325 startTime = SPfrontEnd->IFseconds();
326 if ((converged = TWOdeltaConverged(pDevice)) == TRUE) {
327 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
328 pElem = pDevice->elements[eIndex];
329 for (index = 0; index <= 3; index++) {
330 if (pElem->evalNodes[index]) {
331 pNode = pElem->pNodes[index];
332 if (pNode->nEqn != 0 && solution[pNode->nEqn] < 0.0) {
333 converged = FALSE;
334 solution[pNode->nEqn] = pNode->nConc = 0.0;
335 }
336 if (pNode->pEqn != 0 && solution[pNode->pEqn] < 0.0) {
337 converged = FALSE;
338 solution[pNode->pEqn] = pNode->pConc = 0.0;
339 }
340 }
341 }
342 }
343 }
344 pDevice->pStats->checkTime[STAT_TRAN] += SPfrontEnd->IFseconds() - startTime;
345
346 return (converged);
347 }
348
349 void
TWOresetJacobian(TWOdevice * pDevice)350 TWOresetJacobian(TWOdevice *pDevice)
351 {
352 int error;
353
354
355 if (!OneCarrier) {
356 TWO_jacLoad(pDevice);
357 } else if (OneCarrier == N_TYPE) {
358 TWONjacLoad(pDevice);
359 } else if (OneCarrier == P_TYPE) {
360 TWOPjacLoad(pDevice);
361 } else {
362 printf("TWOresetJacobian: unknown carrier type\n");
363 exit(-1);
364 }
365 error = spFactor(pDevice->matrix);
366 if (foundError(error)) {
367 exit(-1);
368 }
369 }
370
371 /*
372 * Compute the device state assuming charge neutrality exists everywhere in
373 * the device. In insulators, where there is no charge, assign the potential
374 * at half the insulator band gap (refPsi).
375 */
376 void
TWOstoreNeutralGuess(TWOdevice * pDevice)377 TWOstoreNeutralGuess(TWOdevice *pDevice)
378 {
379 int nIndex, eIndex;
380 TWOelem *pElem;
381 TWOnode *pNode;
382 double nie, conc, absConc, refPsi, psi, ni, pi, sign;
383
384 /* assign the initial guess for Poisson's equation */
385 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
386 pElem = pDevice->elements[eIndex];
387 refPsi = pElem->matlInfo->refPsi;
388 if (pElem->elemType == INSULATOR) {
389 for (nIndex = 0; nIndex <= 3; nIndex++) {
390 if (pElem->evalNodes[nIndex]) {
391 pNode = pElem->pNodes[nIndex];
392 if (pNode->nodeType == CONTACT) {
393 /*
394 * a metal contact to insulator domain so use work function
395 * difference as the value of psi
396 */
397 pNode->psi = RefPsi - pNode->eaff;
398 } else {
399 pNode->psi = refPsi;
400 }
401 }
402 }
403 }
404 if (pElem->elemType == SEMICON) {
405 for (nIndex = 0; nIndex <= 3; nIndex++) {
406 if (pElem->evalNodes[nIndex]) {
407 pNode = pElem->pNodes[nIndex];
408 /* silicon nodes */
409 nie = pNode->nie;
410 conc = pNode->netConc / pNode->nie;
411 psi = 0.0;
412 ni = nie;
413 pi = nie;
414 sign = SGN(conc);
415 absConc = ABS(conc);
416 if (conc != 0.0) {
417 psi = sign * log(0.5 * absConc
418 + sqrt(1.0 + 0.25 * absConc * absConc));
419 ni = nie * exp(psi);
420 pi = nie * exp(-psi);
421 }
422 pNode->psi = refPsi + psi;
423 pNode->nConc = ni;
424 pNode->pConc = pi;
425 /* store the initial guess in the dc solution vector */
426 if (pNode->nodeType != CONTACT) {
427 pDevice->dcSolution[pNode->poiEqn] = pNode->psi;
428 }
429 }
430 }
431 }
432 }
433 }
434
435 /* computing the equilibrium solution; solution of Poisson's eqn */
436 /* the solution is used as an initial starting point for bias conditions */
TWOequilSolve(TWOdevice * pDevice)437 int TWOequilSolve(TWOdevice *pDevice)
438 {
439 BOOLEAN newSolver = FALSE;
440 int error;
441 int nIndex, eIndex;
442 TWOelem *pElem;
443 TWOnode *pNode;
444 double startTime, setupTime, miscTime;
445
446 setupTime = miscTime = 0.0;
447
448 /* SETUP */
449 startTime = SPfrontEnd->IFseconds();
450
451 /* Set up pDevice to compute the equilibrium solution. If the solver
452 * is for bias, the arrays must be freed and allocated to the correct
453 * sizes for an equilibrium solution; if it is a new solver, they are
454 * only allocated; and if already an equilibrium solve, nothing
455 * needs to be done */
456 /* FALLTHROUGH added to suppress GCC warning due to
457 * -Wimplicit-fallthrough flag */
458 switch (pDevice->solverType) {
459 case SLV_SMSIG:
460 case SLV_BIAS:
461 /* Free memory allocated for the bias solution */
462 FREE(pDevice->dcSolution);
463 FREE(pDevice->dcDeltaSolution);
464 FREE(pDevice->copiedSolution);
465 FREE(pDevice->rhs);
466 FREE(pDevice->rhsImag);
467 spDestroy(pDevice->matrix);
468 /* FALLTHROUGH */
469 case SLV_NONE: {
470 /* Allocate memory needed for an equilibrium solution */
471 const int n_dim = pDevice->dimEquil;
472 const int n_eqn = n_dim - 1;
473 pDevice->poissonOnly = TRUE;
474 pDevice->numEqns = n_eqn;
475 XCALLOC(pDevice->dcSolution, double, n_dim);
476 XCALLOC(pDevice->dcDeltaSolution, double, n_dim);
477 XCALLOC(pDevice->copiedSolution, double, n_dim);
478 XCALLOC(pDevice->rhs, double, n_dim);
479 pDevice->matrix = spCreate(n_eqn, 0, &error);
480 if (error == spNO_MEMORY) {
481 (void) fprintf(cp_err, "TWOequilSolve: Out of Memory\n");
482 return E_NOMEM;
483 }
484 newSolver = TRUE;
485 spSetReal(pDevice->matrix); /* set to a real matrix */
486 TWOQjacBuild(pDevice);
487 pDevice->numOrigEquil = spElementCount(pDevice->matrix);
488 pDevice->numFillEquil = 0;
489 pDevice->solverType = SLV_EQUIL;
490 break;
491 }
492 case SLV_EQUIL: /* Nothing to do if already equilibrium solver */
493 break;
494 default: /* Invalid data */
495 fprintf(stderr, "Panic: Unknown solver type in equil solution.\n");
496 return E_PANIC;
497 } /* end of switch over solve type */
498 TWOstoreNeutralGuess(pDevice);
499 setupTime += SPfrontEnd->IFseconds() - startTime;
500
501 /* SOLVE */
502 TWOdcSolve(pDevice, MaxIterations, newSolver, FALSE, NULL);
503
504 /* MISCELLANEOUS */
505 startTime = SPfrontEnd->IFseconds();
506 if (newSolver) {
507 pDevice->numFillEquil = spFillinCount(pDevice->matrix);
508 }
509 if (pDevice->converged) {
510 TWOQcommonTerms(pDevice);
511
512 /* save equilibrium potential */
513 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
514 pElem = pDevice->elements[eIndex];
515 for (nIndex = 0; nIndex <= 3; nIndex++) {
516 if (pElem->evalNodes[nIndex]) {
517 pNode = pElem->pNodes[nIndex];
518 pNode->psi0 = pNode->psi;
519 }
520 }
521 }
522 } else {
523 printf("TWOequilSolve: No Convergence\n");
524 }
525 miscTime += SPfrontEnd->IFseconds() - startTime;
526 pDevice->pStats->setupTime[STAT_SETUP] += setupTime;
527 pDevice->pStats->miscTime[STAT_SETUP] += miscTime;
528
529 return 0;
530 }
531
532 /* compute the solution under an applied bias */
533 /* the equilibrium solution is taken as an initial guess */
534 void
TWObiasSolve(TWOdevice * pDevice,int iterationLimit,BOOLEAN tranAnalysis,TWOtranInfo * info)535 TWObiasSolve(TWOdevice *pDevice, int iterationLimit, BOOLEAN tranAnalysis,
536 TWOtranInfo *info)
537 {
538 BOOLEAN newSolver = FALSE;
539 int error;
540 int index, eIndex;
541 TWOelem *pElem;
542 TWOnode *pNode;
543 double refPsi;
544 double startTime, setupTime, miscTime;
545
546 setupTime = miscTime = 0.0;
547
548 /* SETUP */
549 startTime = SPfrontEnd->IFseconds();
550 /* Set up for solving for bias */
551 /* FALLTHROUGH added to suppress GCC warning due to
552 * -Wimplicit-fallthrough flag */
553 switch (pDevice->solverType) {
554 case SLV_EQUIL:
555 /* free up the vectors allocated in the equilibrium solution */
556 FREE(pDevice->dcSolution);
557 FREE(pDevice->dcDeltaSolution);
558 FREE(pDevice->copiedSolution);
559 FREE(pDevice->rhs);
560 spDestroy(pDevice->matrix);
561 /* FALLTHROUGH */
562 case SLV_NONE:
563 /* Set up for bias */
564 pDevice->poissonOnly = FALSE;
565 pDevice->numEqns = pDevice->dimBias - 1;
566 XCALLOC(pDevice->dcSolution, double, pDevice->dimBias);
567 XCALLOC(pDevice->dcDeltaSolution, double, pDevice->dimBias);
568 XCALLOC(pDevice->copiedSolution, double, pDevice->dimBias);
569 XCALLOC(pDevice->rhs, double, pDevice->dimBias);
570 XCALLOC(pDevice->rhsImag, double, pDevice->dimBias);
571 pDevice->matrix = spCreate(pDevice->numEqns, 1, &error);
572 if (error == spNO_MEMORY) {
573 printf("TWObiasSolve: Out of Memory\n");
574 exit(-1);
575 }
576 newSolver = TRUE;
577 if (!OneCarrier) {
578 TWO_jacBuild(pDevice);
579 } else if (OneCarrier == N_TYPE) {
580 TWONjacBuild(pDevice);
581 } else if (OneCarrier == P_TYPE) {
582 TWOPjacBuild(pDevice);
583 }
584 pDevice->numOrigBias = spElementCount(pDevice->matrix);
585 pDevice->numFillBias = 0;
586 TWOstoreInitialGuess(pDevice);
587 /* FALLTHROUGH */
588 case SLV_SMSIG:
589 spSetReal(pDevice->matrix);
590 pDevice->solverType = SLV_BIAS;
591 break;
592 case SLV_BIAS:
593 break;
594 default:
595 fprintf(stderr, "Panic: Unknown solver type in bias solution.\n");
596 exit(-1);
597 break;
598 }
599 setupTime += SPfrontEnd->IFseconds() - startTime;
600
601 /* SOLVE */
602 TWOdcSolve(pDevice, iterationLimit, newSolver, tranAnalysis, info);
603
604 /* MISCELLANEOUS */
605 startTime = SPfrontEnd->IFseconds();
606 if (newSolver) {
607 pDevice->numFillBias = spFillinCount(pDevice->matrix);
608 }
609 if ((!pDevice->converged) && iterationLimit > 1) {
610 printf("TWObiasSolve: No Convergence\n");
611 } else if (pDevice->converged) {
612 /* update the nodal quantities */
613 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
614 pElem = pDevice->elements[eIndex];
615 refPsi = pElem->matlInfo->refPsi;
616 for (index = 0; index <= 3; index++) {
617 if (pElem->evalNodes[index]) {
618 pNode = pElem->pNodes[index];
619 if (pNode->nodeType != CONTACT) {
620 pNode->psi = pDevice->dcSolution[pNode->psiEqn];
621 if (pElem->elemType == SEMICON) {
622 if (!OneCarrier) {
623 pNode->nConc = pDevice->dcSolution[pNode->nEqn];
624 pNode->pConc = pDevice->dcSolution[pNode->pEqn];
625 } else if (OneCarrier == N_TYPE) {
626 pNode->nConc = pDevice->dcSolution[pNode->nEqn];
627 pNode->pConc = pNode->nie * exp(-pNode->psi + refPsi);
628 } else if (OneCarrier == P_TYPE) {
629 pNode->pConc = pDevice->dcSolution[pNode->pEqn];
630 pNode->nConc = pNode->nie * exp(pNode->psi - refPsi);
631 }
632 }
633 }
634 }
635 }
636 }
637
638 /* update the current terms */
639 if (!OneCarrier) {
640 TWO_commonTerms(pDevice, FALSE, tranAnalysis, info);
641 } else if (OneCarrier == N_TYPE) {
642 TWONcommonTerms(pDevice, FALSE, tranAnalysis, info);
643 } else if (OneCarrier == P_TYPE) {
644 TWOPcommonTerms(pDevice, FALSE, tranAnalysis, info);
645 }
646 } else if (iterationLimit <= 1) {
647 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
648 pElem = pDevice->elements[eIndex];
649 refPsi = pElem->matlInfo->refPsi;
650 for (index = 0; index <= 3; index++) {
651 if (pElem->evalNodes[index]) {
652 pNode = pElem->pNodes[index];
653 if (pNode->nodeType != CONTACT) {
654 pNode->psi = pDevice->dcSolution[pNode->psiEqn];
655 pDevice->devState0 [pNode->nodePsi] = pNode->psi;
656 if (pElem->elemType == SEMICON) {
657 if (!OneCarrier) {
658 pNode->nConc = pDevice->dcSolution[pNode->nEqn];
659 pNode->pConc = pDevice->dcSolution[pNode->pEqn];
660 } else if (OneCarrier == N_TYPE) {
661 pNode->nConc = pDevice->dcSolution[pNode->nEqn];
662 pNode->pConc = pNode->nie * exp(-pNode->psi + refPsi);
663 } else if (OneCarrier == P_TYPE) {
664 pNode->pConc = pDevice->dcSolution[pNode->pEqn];
665 pNode->nConc = pNode->nie * exp(pNode->psi - refPsi);
666 }
667 pDevice->devState0 [pNode->nodeN] = pNode->nConc;
668 pDevice->devState0 [pNode->nodeP] = pNode->pConc;
669 }
670 }
671 }
672 }
673 }
674 }
675 miscTime += SPfrontEnd->IFseconds() - startTime;
676 if (tranAnalysis) {
677 pDevice->pStats->setupTime[STAT_TRAN] += setupTime;
678 pDevice->pStats->miscTime[STAT_TRAN] += miscTime;
679 } else {
680 pDevice->pStats->setupTime[STAT_DC] += setupTime;
681 pDevice->pStats->miscTime[STAT_DC] += miscTime;
682 }
683 }
684
685 /* Copy the device's equilibrium state into the solution vector. */
686 void
TWOstoreEquilibGuess(TWOdevice * pDevice)687 TWOstoreEquilibGuess(TWOdevice *pDevice)
688 {
689 int nIndex, eIndex;
690 double *solution = pDevice->dcSolution;
691 double refPsi;
692 TWOelem *pElem;
693 TWOnode *pNode;
694
695 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
696 pElem = pDevice->elements[eIndex];
697 refPsi = pElem->matlInfo->refPsi;
698 for (nIndex = 0; nIndex <= 3; nIndex++) {
699 if (pElem->evalNodes[nIndex]) {
700 pNode = pElem->pNodes[nIndex];
701 if (pNode->nodeType != CONTACT) {
702 pDevice->dcSolution[pNode->psiEqn] = pNode->psi0;
703 if (pElem->elemType == SEMICON) {
704 if (!OneCarrier) {
705 solution[pNode->nEqn] = pNode->nie * exp(pNode->psi0 - refPsi);
706 solution[pNode->pEqn] = pNode->nie * exp(-pNode->psi0 + refPsi);
707 } else if (OneCarrier == N_TYPE) {
708 solution[pNode->nEqn] = pNode->nie * exp(pNode->psi0 - refPsi);
709 } else if (OneCarrier == P_TYPE) {
710 solution[pNode->pEqn] = pNode->nie * exp(-pNode->psi0 + refPsi);
711 }
712 }
713 }
714 }
715 }
716 }
717 }
718
719 /* Copy the device's internal state into the solution vector. */
720 void
TWOstoreInitialGuess(TWOdevice * pDevice)721 TWOstoreInitialGuess(TWOdevice *pDevice)
722 {
723 int nIndex, eIndex;
724 double *solution = pDevice->dcSolution;
725 TWOelem *pElem;
726 TWOnode *pNode;
727
728 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
729 pElem = pDevice->elements[eIndex];
730 for (nIndex = 0; nIndex <= 3; nIndex++) {
731 if (pElem->evalNodes[nIndex]) {
732 pNode = pElem->pNodes[nIndex];
733 if (pNode->nodeType != CONTACT) {
734 solution[pNode->psiEqn] = pNode->psi;
735 if (pElem->elemType == SEMICON) {
736 if (!OneCarrier) {
737 solution[pNode->nEqn] = pNode->nConc;
738 solution[pNode->pEqn] = pNode->pConc;
739 } else if (OneCarrier == N_TYPE) {
740 solution[pNode->nEqn] = pNode->nConc;
741 } else if (OneCarrier == P_TYPE) {
742 solution[pNode->pEqn] = pNode->pConc;
743 }
744 }
745 }
746 }
747 }
748 }
749 }
750
751
752 /*
753 * function computeNewDelta computes an acceptable delta
754 */
755
756 void
oldTWOnewDelta(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)757 oldTWOnewDelta(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info)
758 {
759 int index;
760 double newNorm, fib, lambda, fibn, fibp;
761 BOOLEAN acceptable = FALSE;
762 lambda = 1.0;
763 fibn = 1.0;
764 fibp = 1.0;
765
766 /*
767 * copy the contents of dcSolution into copiedSolution and modify
768 * dcSolution by adding the deltaSolution
769 */
770
771 for (index = 1; index <= pDevice->numEqns; ++index) {
772 pDevice->copiedSolution[index] = pDevice->dcSolution[index];
773 pDevice->dcSolution[index] += pDevice->dcDeltaSolution[index];
774 }
775
776 /* compute L2 norm of the deltaX vector */
777 pDevice->rhsNorm = l2Norm(pDevice->dcDeltaSolution, pDevice->numEqns);
778
779 if (pDevice->poissonOnly) {
780 TWOQrhsLoad(pDevice);
781 } else if (!OneCarrier) {
782 TWO_rhsLoad(pDevice, tranAnalysis, info);
783 } else if (OneCarrier == N_TYPE) {
784 TWONrhsLoad(pDevice, tranAnalysis, info);
785 } else if (OneCarrier == P_TYPE) {
786 TWOPrhsLoad(pDevice, tranAnalysis, info);
787 }
788 newNorm = TWOnuNorm(pDevice);
789 if (newNorm <= pDevice->rhsNorm) {
790 acceptable = TRUE;
791 } else {
792 /* chop the step size so that deltax is acceptable */
793 for (; !acceptable;) {
794 fib = fibp;
795 fibp = fibn;
796 fibn += fib;
797 lambda *= (fibp / fibn);
798
799 for (index = 1; index <= pDevice->numEqns; index++) {
800 pDevice->dcSolution[index] = pDevice->copiedSolution[index]
801 + lambda * pDevice->dcDeltaSolution[index];
802 }
803 if (pDevice->poissonOnly) {
804 TWOQrhsLoad(pDevice);
805 } else if (!OneCarrier) {
806 TWO_rhsLoad(pDevice, tranAnalysis, info);
807 } else if (OneCarrier == N_TYPE) {
808 TWONrhsLoad(pDevice, tranAnalysis, info);
809 } else if (OneCarrier == P_TYPE) {
810 TWOPrhsLoad(pDevice, tranAnalysis, info);
811 }
812 newNorm = TWOnuNorm(pDevice);
813 if (newNorm <= pDevice->rhsNorm) {
814 acceptable = TRUE;
815 }
816 }
817 }
818 /* restore the previous dcSolution and the new deltaSolution */
819 pDevice->rhsNorm = newNorm;
820 for (index = 1; index <= pDevice->numEqns; index++) {
821 pDevice->dcSolution[index] = pDevice->copiedSolution[index];
822 pDevice->dcDeltaSolution[index] *= lambda;
823 }
824 }
825
826
827
828 int
TWOnewDelta(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)829 TWOnewDelta(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info)
830 {
831 int index, iterNum = 0;
832 double newNorm;
833 double fib, lambda, fibn, fibp;
834 BOOLEAN acceptable = FALSE, error = FALSE;
835
836 lambda = 1.0;
837 fibn = 1.0;
838 fibp = 1.0;
839
840 /*
841 * copy the contents of dcSolution into copiedSolution and modify
842 * dcSolution by adding the deltaSolution
843 */
844
845 for (index = 1; index <= pDevice->numEqns; ++index) {
846 pDevice->copiedSolution[index] = pDevice->dcSolution[index];
847 pDevice->dcSolution[index] += pDevice->dcDeltaSolution[index];
848 }
849
850 if (pDevice->poissonOnly) {
851 TWOQrhsLoad(pDevice);
852 } else if (!OneCarrier) {
853 TWO_rhsLoad(pDevice, tranAnalysis, info);
854 } else if (OneCarrier == N_TYPE) {
855 TWONrhsLoad(pDevice, tranAnalysis, info);
856 } else if (OneCarrier == P_TYPE) {
857 TWOPrhsLoad(pDevice, tranAnalysis, info);
858 }
859 newNorm = maxNorm(pDevice->rhs, pDevice->numEqns);
860 if (pDevice->rhsNorm <= pDevice->abstol) {
861 lambda = 0.0;
862 newNorm = pDevice->rhsNorm;
863 } else if (newNorm < pDevice->rhsNorm) {
864 acceptable = TRUE;
865 } else {
866 /* chop the step size so that deltax is acceptable */
867 if (TWOdcDebug) {
868 fprintf(stdout, " %11.4e %11.4e\n",
869 newNorm, lambda);
870 }
871 while (!acceptable) {
872 iterNum++;
873
874 if (iterNum > NORM_RED_MAXITERS) {
875 error = TRUE;
876 lambda = 0.0;
877 /* Don't break out until after we've reset the device. */
878 }
879 fib = fibp;
880 fibp = fibn;
881 fibn += fib;
882 lambda *= (fibp / fibn);
883
884 for (index = 1; index <= pDevice->numEqns; index++) {
885 pDevice->dcSolution[index] = pDevice->copiedSolution[index]
886 + lambda * pDevice->dcDeltaSolution[index];
887 }
888 if (pDevice->poissonOnly) {
889 TWOQrhsLoad(pDevice);
890 } else if (!OneCarrier) {
891 TWO_rhsLoad(pDevice, tranAnalysis, info);
892 } else if (OneCarrier == N_TYPE) {
893 TWONrhsLoad(pDevice, tranAnalysis, info);
894 } else if (OneCarrier == P_TYPE) {
895 TWOPrhsLoad(pDevice, tranAnalysis, info);
896 }
897 newNorm = maxNorm(pDevice->rhs, pDevice->numEqns);
898 if (error) {
899 break;
900 }
901 if (newNorm <= pDevice->rhsNorm) {
902 acceptable = TRUE;
903 }
904 if (TWOdcDebug) {
905 fprintf(stdout, " %11.4e %11.4e\n",
906 newNorm, lambda);
907 }
908 }
909 }
910 /* restore the previous dcSolution and the new deltaSolution */
911 pDevice->rhsNorm = newNorm;
912 for (index = 1; index <= pDevice->numEqns; index++) {
913 pDevice->dcSolution[index] = pDevice->copiedSolution[index];
914 pDevice->dcDeltaSolution[index] *= lambda;
915 }
916 return (error);
917 }
918
919
920 /* Predict the values of the internal variables at the next timepoint. */
921 /* Needed for Predictor-Corrector LTE estimation */
922 void
TWOpredict(TWOdevice * pDevice,TWOtranInfo * info)923 TWOpredict(TWOdevice *pDevice, TWOtranInfo *info)
924 {
925 int nIndex, eIndex;
926 TWOnode *pNode;
927 TWOelem *pElem;
928 double startTime, miscTime = 0.0;
929
930 /* TRANSIENT MISC */
931 startTime = SPfrontEnd->IFseconds();
932 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
933 pElem = pDevice->elements[eIndex];
934 for (nIndex = 0; nIndex <= 3; nIndex++) {
935 if (pElem->evalNodes[nIndex]) {
936 pNode = pElem->pNodes[nIndex];
937 pNode->psi = pDevice->devState1 [pNode->nodePsi];
938 if ((pElem->elemType == SEMICON) && (pNode->nodeType != CONTACT)) {
939 if (!OneCarrier) {
940 pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
941 pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
942 } else if (OneCarrier == N_TYPE) {
943 pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
944 pNode->pPred = pDevice->devState1 [pNode->nodeP];
945 } else if (OneCarrier == P_TYPE) {
946 pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
947 pNode->nPred = pDevice->devState1 [pNode->nodeN];
948 }
949 pNode->nConc = pNode->nPred;
950 pNode->pConc = pNode->pPred;
951 }
952 }
953 }
954 }
955 miscTime += SPfrontEnd->IFseconds() - startTime;
956 pDevice->pStats->miscTime[STAT_TRAN] += miscTime;
957 }
958
959
960 /* Estimate the device's overall truncation error. */
961 double
TWOtrunc(TWOdevice * pDevice,TWOtranInfo * info,double delta)962 TWOtrunc(TWOdevice *pDevice, TWOtranInfo *info, double delta)
963 {
964 int nIndex, eIndex;
965 TWOelem *pElem;
966 TWOnode *pNode;
967 double tolN, tolP, lte, relError, temp;
968 double lteCoeff = info->lteCoeff;
969 double mult = 10.0;
970 double reltol;
971 double startTime, lteTime = 0.0;
972
973 /* TRANSIENT LTE */
974 startTime = SPfrontEnd->IFseconds();
975
976 relError = 0.0;
977 reltol = pDevice->reltol * mult;
978
979 /* need to get the predictor for the current order */
980 /* the scheme currently used is very dumb. need to fix later */
981 computePredCoeff(info->method, info->order, info->predCoeff, info->delta);
982
983 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
984 pElem = pDevice->elements[eIndex];
985 for (nIndex = 0; nIndex <= 3; nIndex++) {
986 if (pElem->evalNodes[nIndex]) {
987 pNode = pElem->pNodes[nIndex];
988 if ((pElem->elemType == SEMICON) && (pNode->nodeType != CONTACT)) {
989 if (!OneCarrier) {
990 tolN = pDevice->abstol + reltol * ABS(pNode->nConc);
991 tolP = pDevice->abstol + reltol * ABS(pNode->pConc);
992 pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
993 pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
994 lte = lteCoeff * (pNode->nConc - pNode->nPred);
995 temp = lte / tolN;
996 relError += temp * temp;
997 lte = lteCoeff * (pNode->pConc - pNode->pPred);
998 temp = lte / tolP;
999 relError += temp * temp;
1000 } else if (OneCarrier == N_TYPE) {
1001 tolN = pDevice->abstol + reltol * ABS(pNode->nConc);
1002 pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
1003 lte = lteCoeff * (pNode->nConc - pNode->nPred);
1004 temp = lte / tolN;
1005 relError += temp * temp;
1006 } else if (OneCarrier == P_TYPE) {
1007 tolP = pDevice->abstol + reltol * ABS(pNode->pConc);
1008 pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
1009 lte = lteCoeff * (pNode->pConc - pNode->pPred);
1010 temp = lte / tolP;
1011 relError += temp * temp;
1012 }
1013 }
1014 }
1015 }
1016 }
1017
1018 relError = MAX(pDevice->abstol, relError); /* make sure it is non zero */
1019
1020 /* the total relative error has been calculated norm-2 squared */
1021 relError = sqrt(relError / pDevice->numEqns);
1022
1023 /* depending on the order of the integration method compute new delta */
1024 temp = delta / pow(relError, 1.0 / (info->order + 1));
1025
1026 lteTime += SPfrontEnd->IFseconds() - startTime;
1027 pDevice->pStats->lteTime += lteTime;
1028
1029 /* return the new delta (stored as temp) */
1030 return (temp);
1031 }
1032
1033 /* Save info from state table into the internal state. */
1034 void
TWOsaveState(TWOdevice * pDevice)1035 TWOsaveState(TWOdevice *pDevice)
1036 {
1037 int nIndex, eIndex;
1038 TWOnode *pNode;
1039 TWOelem *pElem;
1040
1041 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
1042 pElem = pDevice->elements[eIndex];
1043 for (nIndex = 0; nIndex <= 3; nIndex++) {
1044 if (pElem->evalNodes[nIndex]) {
1045 pNode = pElem->pNodes[nIndex];
1046 pNode->psi = pDevice->devState1 [pNode->nodePsi];
1047 if ((pElem->elemType == SEMICON) && (pNode->nodeType != CONTACT)) {
1048 pNode->nConc = pDevice->devState1 [pNode->nodeN];
1049 pNode->pConc = pDevice->devState1 [pNode->nodeP];
1050 }
1051 }
1052 }
1053 }
1054 }
1055
1056 /*
1057 * A function that checks convergence based on the convergence of the quasi
1058 * Fermi levels. This should be better since we are looking at potentials in
1059 * all (psi, phin, phip)
1060 */
1061 BOOLEAN
TWOpsiDeltaConverged(TWOdevice * pDevice)1062 TWOpsiDeltaConverged(TWOdevice *pDevice)
1063 {
1064 int index, nIndex, eIndex;
1065 TWOnode *pNode;
1066 TWOelem *pElem;
1067 double xOld, xNew, xDelta, tol;
1068 double psi, newPsi, nConc, pConc, newN, newP;
1069 double phiN, phiP, newPhiN, newPhiP;
1070 BOOLEAN converged = TRUE;
1071
1072 /* equilibrium solution */
1073 if (pDevice->poissonOnly) {
1074 for (index = 1; index <= pDevice->numEqns; index++) {
1075 xOld = pDevice->dcSolution[index];
1076 xDelta = pDevice->dcDeltaSolution[index];
1077 xNew = xOld + xDelta;
1078 tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew));
1079 if (ABS(xDelta) > tol) {
1080 converged = FALSE;
1081 goto done;
1082 }
1083 }
1084 return (converged);
1085 }
1086 /* bias solution. check convergence on psi, phin, phip */
1087 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
1088 pElem = pDevice->elements[eIndex];
1089 for (nIndex = 0; nIndex <= 3; nIndex++) {
1090 if (pElem->evalNodes[nIndex]) {
1091 pNode = pElem->pNodes[nIndex];
1092 if (pNode->nodeType != CONTACT) {
1093 /* check convergence on psi */
1094 xOld = pDevice->dcSolution[pNode->psiEqn];
1095 xDelta = pDevice->dcDeltaSolution[pNode->psiEqn];
1096 xNew = xOld + xDelta;
1097 tol = pDevice->abstol +
1098 pDevice->reltol * MAX(ABS(xOld), ABS(xNew));
1099 if (ABS(xDelta) > tol) {
1100 converged = FALSE;
1101 goto done;
1102 }
1103 }
1104 /* check convergence on phin and phip */
1105 if ((pElem->elemType == SEMICON) && (pNode->nodeType != CONTACT)) {
1106 psi = pDevice->dcSolution[pNode->psiEqn];
1107 nConc = pDevice->dcSolution[pNode->nEqn];
1108 pConc = pDevice->dcSolution[pNode->pEqn];
1109 newPsi = psi + pDevice->dcDeltaSolution[pNode->psiEqn];
1110 newN = nConc + pDevice->dcDeltaSolution[pNode->nEqn];
1111 newP = pConc + pDevice->dcDeltaSolution[pNode->pEqn];
1112 phiN = psi - log(nConc / pNode->nie);
1113 phiP = psi + log(pConc / pNode->nie);
1114 newPhiN = newPsi - log(newN / pNode->nie);
1115 newPhiP = newPsi + log(newP / pNode->nie);
1116 tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiN), ABS(newPhiN));
1117 if (ABS(newPhiN - phiN) > tol) {
1118 converged = FALSE;
1119 goto done;
1120 }
1121 tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiP), ABS(newPhiP));
1122 if (ABS(newPhiP - phiP) > tol) {
1123 converged = FALSE;
1124 goto done;
1125 }
1126 }
1127 }
1128 }
1129 }
1130 done:
1131 return (converged);
1132 }
1133
1134
1135 /* Function to compute Nu norm of the rhs vector. */
1136 /* nu-Norm calculation based upon work done at Stanford. */
1137 double
TWOnuNorm(TWOdevice * pDevice)1138 TWOnuNorm(TWOdevice *pDevice)
1139 {
1140 double norm = 0.0;
1141 double temp;
1142 int index;
1143
1144 /* the LU Decomposed matrix is available. use it to calculate x */
1145
1146 spSolve(pDevice->matrix, pDevice->rhs, pDevice->rhsImag,
1147 NULL, NULL);
1148
1149 /* the solution is in the rhsImag vector */
1150 /* compute L2-norm of the rhsImag vector */
1151
1152 for (index = 1; index <= pDevice->numEqns; index++) {
1153 temp = pDevice->rhsImag[index];
1154 norm += temp * temp;
1155 }
1156 norm = sqrt(norm);
1157 return (norm);
1158 }
1159
1160 /*
1161 * Check for numerical errors in the Jacobian. Useful debugging aid when new
1162 * models are being incorporated.
1163 */
1164 void
TWOjacCheck(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)1165 TWOjacCheck(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info)
1166 {
1167 int index, rIndex;
1168 double del, diff, tol, *dptr;
1169
1170 if (TWOjacDebug) {
1171 if (!OneCarrier) {
1172 TWO_sysLoad(pDevice, tranAnalysis, info);
1173 } else if (OneCarrier == N_TYPE) {
1174 TWONsysLoad(pDevice, tranAnalysis, info);
1175 } else if (OneCarrier == P_TYPE) {
1176 TWOPsysLoad(pDevice, tranAnalysis, info);
1177 }
1178 /*
1179 * spPrint( pDevice->matrix, 0, 1, 1 );
1180 */
1181 pDevice->rhsNorm = maxNorm(pDevice->rhs, pDevice->numEqns);
1182 for (index = 1; index <= pDevice->numEqns; index++) {
1183 if (1e3 * ABS(pDevice->rhs[index]) > pDevice->rhsNorm) {
1184 fprintf(stderr, "eqn %d: res %11.4e, norm %11.4e\n",
1185 index, pDevice->rhs[index], pDevice->rhsNorm);
1186 }
1187 }
1188 for (index = 1; index <= pDevice->numEqns; index++) {
1189 pDevice->rhsImag[index] = pDevice->rhs[index];
1190 }
1191 for (index = 1; index <= pDevice->numEqns; index++) {
1192 pDevice->copiedSolution[index] = pDevice->dcSolution[index];
1193 del = 1e-4 * pDevice->abstol + 1e-6 * ABS(pDevice->dcSolution[index]);
1194 pDevice->dcSolution[index] += del;
1195 if (!OneCarrier) {
1196 TWO_rhsLoad(pDevice, tranAnalysis, info);
1197 } else if (OneCarrier == N_TYPE) {
1198 TWONrhsLoad(pDevice, tranAnalysis, info);
1199 } else if (OneCarrier == P_TYPE) {
1200 TWOPrhsLoad(pDevice, tranAnalysis, info);
1201 }
1202 pDevice->dcSolution[index] = pDevice->copiedSolution[index];
1203 for (rIndex = 1; rIndex <= pDevice->numEqns; rIndex++) {
1204 diff = (pDevice->rhsImag[rIndex] - pDevice->rhs[rIndex]) / del;
1205 dptr = spFindElement(pDevice->matrix, rIndex, index);
1206 if (dptr != NULL) {
1207 tol = (1e-4 * pDevice->abstol) + (1e-2 * MAX(ABS(diff), ABS(*dptr)));
1208 if ((diff != 0.0) && (ABS(diff - *dptr) > tol)) {
1209 fprintf(stderr, "Mismatch[%d][%d]: FD = %11.4e, AJ = %11.4e\n\t FD-AJ = %11.4e vs. %11.4e\n",
1210 rIndex, index, diff, *dptr,
1211 ABS(diff - *dptr), tol);
1212 }
1213 } else {
1214 if (diff != 0.0) {
1215 fprintf(stderr, "Missing [%d][%d]: FD = %11.4e, AJ = 0.0\n",
1216 rIndex, index, diff);
1217 }
1218 }
1219 }
1220 }
1221 }
1222 }
1223