1 /**********
2 Copyright 1991 Regents of the University of California.  All rights reserved.
3 Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
4 **********/
5 
6 /* Functions to compute the ac admittances of a device. */
7 
8 #include "ngspice/ngspice.h"
9 #include "ngspice/numglobs.h"
10 #include "ngspice/numenum.h"
11 #include "ngspice/numconst.h"
12 #include "ngspice/twodev.h"
13 #include "ngspice/twomesh.h"
14 #include "ngspice/complex.h"
15 #include "ngspice/spmatrix.h"
16 #include "ngspice/bool.h"
17 #include "ngspice/macros.h"
18 #include "ngspice/ifsim.h"
19 #include "twoddefs.h"
20 #include "twodext.h"
21 #include "ngspice/cidersupt.h"
22 
23 extern IFfrontEnd *SPfrontEnd;
24 
25 /*
26  * mmhhh this may cause troubles
27  * Paolo Nenzi 2002
28  */
29  SPcomplex yTotal;
30 
31 int
NUMD2admittance(TWOdevice * pDevice,double omega,SPcomplex * yd)32 NUMD2admittance(TWOdevice *pDevice, double omega, SPcomplex *yd)
33 {
34   TWOnode *pNode;
35   TWOelem *pElem;
36   int index, eIndex;
37   double dxdy;
38   double *solnReal, *solnImag;
39   double *rhsReal, *rhsImag;
40   SPcomplex yAc, cOmega, *y;
41   BOOLEAN deltaVContact = FALSE;
42   BOOLEAN SORFailed;
43   double startTime;
44 
45   /* Each time we call this counts as one AC iteration. */
46   pDevice->pStats->numIters[STAT_AC] += 1;
47 
48   /*
49    * change context names of solution vectors for ac analysis dcDeltaSolution
50    * stores the real part and copiedSolution stores the imaginary part of the
51    * ac solution vector
52    */
53   pDevice->solverType = SLV_SMSIG;
54   rhsReal = pDevice->rhs;
55   rhsImag = pDevice->rhsImag;
56   solnReal = pDevice->dcDeltaSolution;
57   solnImag = pDevice->copiedSolution;
58 
59   /* use a normalized radian frequency */
60   omega *= TNorm;
61   CMPLX_ASSIGN_VALUE(cOmega, 0.0, omega);
62 
63   if ((AcAnalysisMethod == SOR) || (AcAnalysisMethod == SOR_ONLY)) {
64     /* LOAD */
65     startTime = SPfrontEnd->IFseconds();
66     /* zero the rhsImag */
67     for (index = 1; index <= pDevice->numEqns; index++) {
68       rhsImag[index] = 0.0;
69     }
70     /* store the new rhs vector */
71     storeNewRhs(pDevice, pDevice->pLastContact);
72     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
73 
74     /* SOLVE */
75     startTime = SPfrontEnd->IFseconds();
76     SORFailed = TWOsorSolve(pDevice, solnReal, solnImag, omega);
77     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
78     if (SORFailed && AcAnalysisMethod == SOR) {
79       AcAnalysisMethod = DIRECT;
80       printf("SOR failed at %g Hz, switching to direct-method ac analysis.\n",
81 	  omega / (2 * M_PI * TNorm) );
82     } else if (SORFailed) {	/* Told to only do SOR, so give up. */
83       printf("SOR failed at %g Hz, returning null admittance.\n",
84 	  omega / (2 * M_PI * TNorm) );
85       CMPLX_ASSIGN_VALUE(*yd, 0.0, 0.0);
86       return (AcAnalysisMethod);
87     }
88   }
89   if (AcAnalysisMethod == DIRECT) {
90     /* LOAD */
91     startTime = SPfrontEnd->IFseconds();
92     for (index = 1; index <= pDevice->numEqns; index++) {
93       rhsImag[index] = 0.0;
94     }
95     /* solve the system of equations directly */
96     if (!OneCarrier) {
97       TWO_jacLoad(pDevice);
98     } else if (OneCarrier == N_TYPE) {
99       TWONjacLoad(pDevice);
100     } else if (OneCarrier == P_TYPE) {
101       TWOPjacLoad(pDevice);
102     }
103     storeNewRhs(pDevice, pDevice->pLastContact);
104 
105     spSetComplex(pDevice->matrix);
106     for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
107       pElem = pDevice->elements[eIndex];
108       if (pElem->elemType == SEMICON) {
109 	dxdy = 0.25 * pElem->dx * pElem->dy;
110 	for (index = 0; index <= 3; index++) {
111 	  pNode = pElem->pNodes[index];
112 	  if (pNode->nodeType != CONTACT) {
113 	    if (!OneCarrier) {
114 	      spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -dxdy * omega);
115 	      spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, dxdy * omega);
116 	    } else if (OneCarrier == N_TYPE) {
117 	      spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -dxdy * omega);
118 	    } else if (OneCarrier == P_TYPE) {
119 	      spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, dxdy * omega);
120 	    }
121 	  }
122 	}
123       }
124     }
125     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
126 
127     /* FACTOR */
128     startTime = SPfrontEnd->IFseconds();
129     spFactor(pDevice->matrix);
130     pDevice->pStats->factorTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
131 
132     /* SOLVE */
133     startTime = SPfrontEnd->IFseconds();
134     spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
135     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
136   }
137   /* MISC */
138   startTime = SPfrontEnd->IFseconds();
139   y = contactAdmittance(pDevice, pDevice->pFirstContact, deltaVContact,
140       solnReal, solnImag, &cOmega);
141   CMPLX_ASSIGN_VALUE(yAc, -y->real, -y->imag);
142   CMPLX_ASSIGN(*yd, yAc);
143   CMPLX_MULT_SELF_SCALAR(*yd, GNorm * pDevice->width * LNorm);
144   pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
145 
146   return (AcAnalysisMethod);
147 }
148 
149 
150 int
NBJT2admittance(TWOdevice * pDevice,double omega,SPcomplex * yIeVce,SPcomplex * yIcVce,SPcomplex * yIeVbe,SPcomplex * yIcVbe)151 NBJT2admittance(TWOdevice *pDevice, double omega, SPcomplex *yIeVce,
152                 SPcomplex *yIcVce, SPcomplex *yIeVbe, SPcomplex *yIcVbe)
153 {
154   TWOcontact *pEmitContact = pDevice->pLastContact;
155   TWOcontact *pColContact = pDevice->pFirstContact;
156   TWOcontact *pBaseContact = pDevice->pFirstContact->next;
157   TWOnode *pNode;
158   TWOelem *pElem;
159   int index, eIndex;
160   double width = pDevice->width;
161   double dxdy;
162   double *solnReal, *solnImag;
163   double *rhsReal, *rhsImag;
164   BOOLEAN SORFailed;
165   SPcomplex *y;
166   SPcomplex pIeVce, pIcVce, pIeVbe, pIcVbe;
167   SPcomplex cOmega;
168   double startTime;
169 
170   /* Each time we call this counts as one AC iteration. */
171   pDevice->pStats->numIters[STAT_AC] += 1;
172 
173   pDevice->solverType = SLV_SMSIG;
174   rhsReal = pDevice->rhs;
175   rhsImag = pDevice->rhsImag;
176   solnReal = pDevice->dcDeltaSolution;
177   solnImag = pDevice->copiedSolution;
178 
179   /* use a normalized radian frequency */
180   omega *= TNorm;
181   CMPLX_ASSIGN_VALUE(cOmega, 0.0, omega);
182   CMPLX_ASSIGN_VALUE(pIeVce, NAN, NAN);
183   CMPLX_ASSIGN_VALUE(pIcVce, NAN, NAN);
184 
185   if ((AcAnalysisMethod == SOR) || (AcAnalysisMethod == SOR_ONLY)) {
186     /* LOAD */
187     startTime = SPfrontEnd->IFseconds();
188     /* zero the rhs before loading in the new rhs */
189     for (index = 1; index <= pDevice->numEqns; index++) {
190       rhsImag[index] = 0.0;
191     }
192     storeNewRhs(pDevice, pColContact);
193     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
194 
195     /* SOLVE */
196     startTime = SPfrontEnd->IFseconds();
197     SORFailed = TWOsorSolve(pDevice, solnReal, solnImag, omega);
198     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
199     if (SORFailed && AcAnalysisMethod == SOR) {
200       AcAnalysisMethod = DIRECT;
201       printf("SOR failed at %g Hz, switching to direct-method ac analysis.\n",
202 	  omega / (2 * M_PI * TNorm) );
203     } else if (SORFailed) {	/* Told to only do SOR, so give up. */
204       printf("SOR failed at %g Hz, returning null admittance.\n",
205 	  omega / (2 * M_PI * TNorm) );
206       CMPLX_ASSIGN_VALUE(*yIeVce, 0.0, 0.0);
207       CMPLX_ASSIGN_VALUE(*yIcVce, 0.0, 0.0);
208       CMPLX_ASSIGN_VALUE(*yIeVbe, 0.0, 0.0);
209       CMPLX_ASSIGN_VALUE(*yIcVbe, 0.0, 0.0);
210       return (AcAnalysisMethod);
211     } else {
212       /* MISC */
213       startTime = SPfrontEnd->IFseconds();
214       y = contactAdmittance(pDevice, pEmitContact, FALSE,
215 	  solnReal, solnImag, &cOmega);
216       CMPLX_ASSIGN_VALUE(pIeVce, y->real, y->imag);
217       y = contactAdmittance(pDevice, pColContact, TRUE,
218 	  solnReal, solnImag, &cOmega);
219       CMPLX_ASSIGN_VALUE(pIcVce, y->real, y->imag);
220       pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
221 
222       /* LOAD */
223       startTime = SPfrontEnd->IFseconds();
224       /* load in the base contribution to the rhs */
225       for (index = 1; index <= pDevice->numEqns; index++) {
226 	rhsImag[index] = 0.0;
227       }
228       storeNewRhs(pDevice, pBaseContact);
229       pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
230 
231       /* SOLVE */
232       startTime = SPfrontEnd->IFseconds();
233       SORFailed = TWOsorSolve(pDevice, solnReal, solnImag, omega);
234       pDevice->pStats->solveTime[STAT_AC] +=
235 	  SPfrontEnd->IFseconds() - startTime;
236       if (SORFailed && AcAnalysisMethod == SOR) {
237 	AcAnalysisMethod = DIRECT;
238 	printf("SOR failed at %g Hz, switching to direct-method ac analysis.\n",
239 	    omega / (2 * M_PI * TNorm) );
240       } else if (SORFailed) {	/* Told to only do SOR, so give up. */
241 	printf("SOR failed at %g Hz, returning null admittance.\n",
242 	    omega / (2 * M_PI * TNorm) );
243 	CMPLX_ASSIGN_VALUE(*yIeVce, 0.0, 0.0);
244 	CMPLX_ASSIGN_VALUE(*yIcVce, 0.0, 0.0);
245 	CMPLX_ASSIGN_VALUE(*yIeVbe, 0.0, 0.0);
246 	CMPLX_ASSIGN_VALUE(*yIcVbe, 0.0, 0.0);
247 	return (AcAnalysisMethod);
248       }
249     }
250   }
251   if (AcAnalysisMethod == DIRECT) {
252     /* LOAD */
253     startTime = SPfrontEnd->IFseconds();
254     for (index = 1; index <= pDevice->numEqns; index++) {
255       rhsImag[index] = 0.0;
256     }
257     /* solve the system of equations directly */
258     if (!OneCarrier) {
259       TWO_jacLoad(pDevice);
260     } else if (OneCarrier == N_TYPE) {
261       TWONjacLoad(pDevice);
262     } else if (OneCarrier == P_TYPE) {
263       TWOPjacLoad(pDevice);
264     }
265     storeNewRhs(pDevice, pColContact);
266     spSetComplex(pDevice->matrix);
267     for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
268       pElem = pDevice->elements[eIndex];
269       if (pElem->elemType == SEMICON) {
270 	dxdy = 0.25 * pElem->dx * pElem->dy;
271 	for (index = 0; index <= 3; index++) {
272 	  pNode = pElem->pNodes[index];
273 	  if (pNode->nodeType != CONTACT) {
274 	    if (!OneCarrier) {
275 	      spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -dxdy * omega);
276 	      spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, dxdy * omega);
277 	    } else if (OneCarrier == N_TYPE) {
278 	      spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -dxdy * omega);
279 	    } else if (OneCarrier == P_TYPE) {
280 	      spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, dxdy * omega);
281 	    }
282 	  }
283 	}
284       }
285     }
286     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
287 
288     /* FACTOR */
289     startTime = SPfrontEnd->IFseconds();
290     spFactor(pDevice->matrix);
291     pDevice->pStats->factorTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
292 
293     /* SOLVE */
294     startTime = SPfrontEnd->IFseconds();
295     spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
296     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
297 
298     /* MISC */
299     startTime = SPfrontEnd->IFseconds();
300     y = contactAdmittance(pDevice, pEmitContact, FALSE,
301 	solnReal, solnImag, &cOmega);
302     CMPLX_ASSIGN_VALUE(pIeVce, y->real, y->imag);
303     y = contactAdmittance(pDevice, pColContact, TRUE,
304 	solnReal, solnImag, &cOmega);
305     CMPLX_ASSIGN_VALUE(pIcVce, y->real, y->imag);
306     pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
307 
308     /* LOAD */
309     startTime = SPfrontEnd->IFseconds();
310     for (index = 1; index <= pDevice->numEqns; index++) {
311       rhsImag[index] = 0.0;
312     }
313     storeNewRhs(pDevice, pBaseContact);
314     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
315 
316     /* FACTOR: already done, no need to repeat. */
317 
318     /* SOLVE */
319     startTime = SPfrontEnd->IFseconds();
320     spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
321     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
322   }
323   /* MISC */
324   startTime = SPfrontEnd->IFseconds();
325   y = contactAdmittance(pDevice, pEmitContact, FALSE,
326       solnReal, solnImag, &cOmega);
327   CMPLX_ASSIGN_VALUE(pIeVbe, y->real, y->imag);
328   y = contactAdmittance(pDevice, pColContact, FALSE,
329       solnReal, solnImag, &cOmega);
330   CMPLX_ASSIGN_VALUE(pIcVbe, y->real, y->imag);
331 
332   CMPLX_ASSIGN(*yIeVce, pIeVce);
333   CMPLX_ASSIGN(*yIeVbe, pIeVbe);
334   CMPLX_ASSIGN(*yIcVce, pIcVce);
335   CMPLX_ASSIGN(*yIcVbe, pIcVbe);
336   CMPLX_MULT_SELF_SCALAR(*yIeVce, GNorm * width * LNorm);
337   CMPLX_MULT_SELF_SCALAR(*yIeVbe, GNorm * width * LNorm);
338   CMPLX_MULT_SELF_SCALAR(*yIcVce, GNorm * width * LNorm);
339   CMPLX_MULT_SELF_SCALAR(*yIcVbe, GNorm * width * LNorm);
340   pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
341 
342   return (AcAnalysisMethod);
343 }
344 
345 int
NUMOSadmittance(TWOdevice * pDevice,double omega,struct mosAdmittances * yAc)346 NUMOSadmittance(TWOdevice *pDevice, double omega, struct mosAdmittances *yAc)
347 {
348   TWOcontact *pDContact = pDevice->pFirstContact;
349   TWOcontact *pGContact = pDevice->pFirstContact->next;
350   TWOcontact *pSContact = pDevice->pFirstContact->next->next;
351 /*  TWOcontact *pBContact = pDevice->pLastContact; */
352   TWOnode *pNode;
353   TWOelem *pElem;
354   int index, eIndex;
355   double width = pDevice->width;
356   double dxdy;
357   double *solnReal, *solnImag;
358   double *rhsReal, *rhsImag;
359   BOOLEAN SORFailed;
360   SPcomplex *y, cOmega;
361   double startTime;
362 
363   /* Each time we call this counts as one AC iteration. */
364   pDevice->pStats->numIters[STAT_AC] += 1;
365 
366   pDevice->solverType = SLV_SMSIG;
367   rhsReal = pDevice->rhs;
368   rhsImag = pDevice->rhsImag;
369   solnReal = pDevice->dcDeltaSolution;
370   solnImag = pDevice->copiedSolution;
371 
372   /* use a normalized radian frequency */
373   omega *= TNorm;
374   CMPLX_ASSIGN_VALUE(cOmega, 0.0, omega);
375 
376   if ((AcAnalysisMethod == SOR) || (AcAnalysisMethod == SOR_ONLY)) {
377     /* LOAD */
378     startTime = SPfrontEnd->IFseconds();
379     /* zero the rhs before loading in the new rhs */
380     for (index = 1; index <= pDevice->numEqns; index++) {
381       rhsImag[index] = 0.0;
382     }
383     storeNewRhs(pDevice, pDContact);
384     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
385 
386     /* SOLVE */
387     startTime = SPfrontEnd->IFseconds();
388     SORFailed = TWOsorSolve(pDevice, solnReal, solnImag, omega);
389     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
390     if (SORFailed && AcAnalysisMethod == SOR) {
391       AcAnalysisMethod = DIRECT;
392       printf("SOR failed at %g Hz, switching to direct-method ac analysis.\n",
393 	  omega / (2 * M_PI * TNorm) );
394     } else if (SORFailed) {	/* Told to only do SOR, so give up. */
395       printf("SOR failed at %g Hz, returning null admittance.\n",
396 	  omega / (2 * M_PI * TNorm) );
397       CMPLX_ASSIGN_VALUE(yAc->yIdVdb, 0.0, 0.0);
398       CMPLX_ASSIGN_VALUE(yAc->yIdVsb, 0.0, 0.0);
399       CMPLX_ASSIGN_VALUE(yAc->yIdVgb, 0.0, 0.0);
400       CMPLX_ASSIGN_VALUE(yAc->yIsVdb, 0.0, 0.0);
401       CMPLX_ASSIGN_VALUE(yAc->yIsVsb, 0.0, 0.0);
402       CMPLX_ASSIGN_VALUE(yAc->yIsVgb, 0.0, 0.0);
403       CMPLX_ASSIGN_VALUE(yAc->yIgVdb, 0.0, 0.0);
404       CMPLX_ASSIGN_VALUE(yAc->yIgVsb, 0.0, 0.0);
405       CMPLX_ASSIGN_VALUE(yAc->yIgVgb, 0.0, 0.0);
406       return (AcAnalysisMethod);
407     } else {
408       /* MISC */
409       startTime = SPfrontEnd->IFseconds();
410       y = contactAdmittance(pDevice, pDContact, TRUE,
411 	  solnReal, solnImag, &cOmega);
412       CMPLX_ASSIGN_VALUE(yAc->yIdVdb, y->real, y->imag);
413       y = contactAdmittance(pDevice, pSContact, FALSE,
414 	  solnReal, solnImag, &cOmega);
415       CMPLX_ASSIGN_VALUE(yAc->yIsVdb, y->real, y->imag);
416       y = GateTypeAdmittance(pDevice, pGContact, FALSE,
417 	  solnReal, solnImag, &cOmega);
418       CMPLX_ASSIGN_VALUE(yAc->yIgVdb, y->real, y->imag);
419       pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
420 
421       /* LOAD */
422       startTime = SPfrontEnd->IFseconds();
423       /* load in the source contribution to the rhs */
424       for (index = 1; index <= pDevice->numEqns; index++) {
425 	rhsImag[index] = 0.0;
426       }
427       storeNewRhs(pDevice, pSContact);
428       pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
429 
430       /* SOLVE */
431       startTime = SPfrontEnd->IFseconds();
432       SORFailed = TWOsorSolve(pDevice, solnReal, solnImag, omega);
433       pDevice->pStats->solveTime[STAT_AC] +=
434 	  SPfrontEnd->IFseconds() - startTime;
435       if (SORFailed && AcAnalysisMethod == SOR) {
436 	AcAnalysisMethod = DIRECT;
437 	printf("SOR failed at %g Hz, switching to direct-method ac analysis.\n",
438 	    omega / (2 * M_PI * TNorm) );
439       } else if (SORFailed) {	/* Told to only do SOR, so give up. */
440 	printf("SOR failed at %g Hz, returning null admittance.\n",
441 	    omega / (2 * M_PI * TNorm) );
442 	CMPLX_ASSIGN_VALUE(yAc->yIdVdb, 0.0, 0.0);
443 	CMPLX_ASSIGN_VALUE(yAc->yIdVsb, 0.0, 0.0);
444 	CMPLX_ASSIGN_VALUE(yAc->yIdVgb, 0.0, 0.0);
445 	CMPLX_ASSIGN_VALUE(yAc->yIsVdb, 0.0, 0.0);
446 	CMPLX_ASSIGN_VALUE(yAc->yIsVsb, 0.0, 0.0);
447 	CMPLX_ASSIGN_VALUE(yAc->yIsVgb, 0.0, 0.0);
448 	CMPLX_ASSIGN_VALUE(yAc->yIgVdb, 0.0, 0.0);
449 	CMPLX_ASSIGN_VALUE(yAc->yIgVsb, 0.0, 0.0);
450 	CMPLX_ASSIGN_VALUE(yAc->yIgVgb, 0.0, 0.0);
451 	return (AcAnalysisMethod);
452       } else {
453 	/* MISC */
454 	startTime = SPfrontEnd->IFseconds();
455 	y = contactAdmittance(pDevice, pDContact, FALSE,
456 	    solnReal, solnImag, &cOmega);
457 	CMPLX_ASSIGN_VALUE(yAc->yIdVsb, y->real, y->imag);
458 	y = contactAdmittance(pDevice, pSContact, TRUE,
459 	    solnReal, solnImag, &cOmega);
460 	CMPLX_ASSIGN_VALUE(yAc->yIsVsb, y->real, y->imag);
461 	y = GateTypeAdmittance(pDevice, pGContact, FALSE,
462 	    solnReal, solnImag, &cOmega);
463 	CMPLX_ASSIGN_VALUE(yAc->yIgVsb, y->real, y->imag);
464 	pDevice->pStats->miscTime[STAT_AC] +=
465 	    SPfrontEnd->IFseconds() - startTime;
466 
467 	/* LOAD */
468 	startTime = SPfrontEnd->IFseconds();
469 	/* load in the gate contribution to the rhs */
470 	for (index = 1; index <= pDevice->numEqns; index++) {
471 	  rhsImag[index] = 0.0;
472 	}
473 	storeNewRhs(pDevice, pGContact);
474 	pDevice->pStats->loadTime[STAT_AC] +=
475 	    SPfrontEnd->IFseconds() - startTime;
476 
477 	/* SOLVE */
478 	startTime = SPfrontEnd->IFseconds();
479 	SORFailed = TWOsorSolve(pDevice, solnReal, solnImag, omega);
480 	pDevice->pStats->solveTime[STAT_AC] +=
481 	    SPfrontEnd->IFseconds() - startTime;
482 	if (SORFailed && AcAnalysisMethod == SOR) {
483 	  AcAnalysisMethod = DIRECT;
484 	  printf("SOR failed at %g Hz, switching to direct-method ac analysis.\n",
485 	      omega / (2 * M_PI * TNorm) );
486 	} else if (SORFailed) {	/* Told to only do SOR, so give up. */
487 	  printf("SOR failed at %g Hz, returning null admittance.\n",
488 	      omega / (2 * M_PI * TNorm) );
489 	  CMPLX_ASSIGN_VALUE(yAc->yIdVdb, 0.0, 0.0);
490 	  CMPLX_ASSIGN_VALUE(yAc->yIdVsb, 0.0, 0.0);
491 	  CMPLX_ASSIGN_VALUE(yAc->yIdVgb, 0.0, 0.0);
492 	  CMPLX_ASSIGN_VALUE(yAc->yIsVdb, 0.0, 0.0);
493 	  CMPLX_ASSIGN_VALUE(yAc->yIsVsb, 0.0, 0.0);
494 	  CMPLX_ASSIGN_VALUE(yAc->yIsVgb, 0.0, 0.0);
495 	  CMPLX_ASSIGN_VALUE(yAc->yIgVdb, 0.0, 0.0);
496 	  CMPLX_ASSIGN_VALUE(yAc->yIgVsb, 0.0, 0.0);
497 	  CMPLX_ASSIGN_VALUE(yAc->yIgVgb, 0.0, 0.0);
498 	  return (AcAnalysisMethod);
499 	}
500       }
501     }
502   }
503   if (AcAnalysisMethod == DIRECT) {
504     /* solve the system of equations directly */
505     /* LOAD */
506     startTime = SPfrontEnd->IFseconds();
507     for (index = 1; index <= pDevice->numEqns; index++) {
508       rhsImag[index] = 0.0;
509     }
510     storeNewRhs(pDevice, pDContact);
511 
512     /* Need to load & factor jacobian once. */
513     if (!OneCarrier) {
514       TWO_jacLoad(pDevice);
515     } else if (OneCarrier == N_TYPE) {
516       TWONjacLoad(pDevice);
517     } else if (OneCarrier == P_TYPE) {
518       TWOPjacLoad(pDevice);
519     }
520     spSetComplex(pDevice->matrix);
521     for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
522       pElem = pDevice->elements[eIndex];
523       if (pElem->elemType == SEMICON) {
524 	dxdy = 0.25 * pElem->dx * pElem->dy;
525 	for (index = 0; index <= 3; index++) {
526 	  pNode = pElem->pNodes[index];
527 	  if (pNode->nodeType != CONTACT) {
528 	    if (!OneCarrier) {
529 	      spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -dxdy * omega);
530 	      spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, dxdy * omega);
531 	    } else if (OneCarrier == N_TYPE) {
532 	      spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -dxdy * omega);
533 	    } else if (OneCarrier == P_TYPE) {
534 	      spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, dxdy * omega);
535 	    }
536 	  }
537 	}
538       }
539     }
540     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
541 
542     /* FACTOR */
543     startTime = SPfrontEnd->IFseconds();
544     spFactor(pDevice->matrix);
545     pDevice->pStats->factorTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
546 
547     /* SOLVE */
548     startTime = SPfrontEnd->IFseconds();
549     spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
550     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
551 
552     /* MISC */
553     startTime = SPfrontEnd->IFseconds();
554     y = contactAdmittance(pDevice, pDContact, TRUE,
555 	solnReal, solnImag, &cOmega);
556     CMPLX_ASSIGN_VALUE(yAc->yIdVdb, y->real, y->imag);
557     y = contactAdmittance(pDevice, pSContact, FALSE,
558 	solnReal, solnImag, &cOmega);
559     CMPLX_ASSIGN_VALUE(yAc->yIsVdb, y->real, y->imag);
560     y = GateTypeAdmittance(pDevice, pGContact, FALSE,
561 	solnReal, solnImag, &cOmega);
562     CMPLX_ASSIGN_VALUE(yAc->yIgVdb, y->real, y->imag);
563     pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
564 
565     /* LOAD */
566     startTime = SPfrontEnd->IFseconds();
567     for (index = 1; index <= pDevice->numEqns; index++) {
568       rhsImag[index] = 0.0;
569     }
570     storeNewRhs(pDevice, pSContact);
571     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
572 
573     /* FACTOR: already done, no need to repeat. */
574 
575     /* SOLVE */
576     startTime = SPfrontEnd->IFseconds();
577     spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
578     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
579 
580     /* MISC */
581     startTime = SPfrontEnd->IFseconds();
582     y = contactAdmittance(pDevice, pDContact, FALSE,
583 	solnReal, solnImag, &cOmega);
584     CMPLX_ASSIGN_VALUE(yAc->yIdVsb, y->real, y->imag);
585     y = contactAdmittance(pDevice, pSContact, TRUE,
586 	solnReal, solnImag, &cOmega);
587     CMPLX_ASSIGN_VALUE(yAc->yIsVsb, y->real, y->imag);
588     y = GateTypeAdmittance(pDevice, pGContact, FALSE,
589 	solnReal, solnImag, &cOmega);
590     CMPLX_ASSIGN_VALUE(yAc->yIgVsb, y->real, y->imag);
591     pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
592 
593     /* LOAD */
594     startTime = SPfrontEnd->IFseconds();
595     for (index = 1; index <= pDevice->numEqns; index++) {
596       rhsImag[index] = 0.0;
597     }
598     storeNewRhs(pDevice, pGContact);
599     pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
600 
601     /* FACTOR: already done, no need to repeat. */
602 
603     /* SOLVE */
604     startTime = SPfrontEnd->IFseconds();
605     spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
606     pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
607   }
608   /* MISC */
609   startTime = SPfrontEnd->IFseconds();
610   y = contactAdmittance(pDevice, pDContact, FALSE,
611       solnReal, solnImag, &cOmega);
612   CMPLX_ASSIGN_VALUE(yAc->yIdVgb, y->real, y->imag);
613   y = contactAdmittance(pDevice, pSContact, FALSE,
614       solnReal, solnImag, &cOmega);
615   CMPLX_ASSIGN_VALUE(yAc->yIsVgb, y->real, y->imag);
616   y = GateTypeAdmittance(pDevice, pGContact, TRUE,
617       solnReal, solnImag, &cOmega);
618   CMPLX_ASSIGN_VALUE(yAc->yIgVgb, y->real, y->imag);
619 
620   CMPLX_MULT_SELF_SCALAR(yAc->yIdVdb, GNorm * width * LNorm);
621   CMPLX_MULT_SELF_SCALAR(yAc->yIdVsb, GNorm * width * LNorm);
622   CMPLX_MULT_SELF_SCALAR(yAc->yIdVgb, GNorm * width * LNorm);
623   CMPLX_MULT_SELF_SCALAR(yAc->yIsVdb, GNorm * width * LNorm);
624   CMPLX_MULT_SELF_SCALAR(yAc->yIsVsb, GNorm * width * LNorm);
625   CMPLX_MULT_SELF_SCALAR(yAc->yIsVgb, GNorm * width * LNorm);
626   CMPLX_MULT_SELF_SCALAR(yAc->yIgVdb, GNorm * width * LNorm);
627   CMPLX_MULT_SELF_SCALAR(yAc->yIgVsb, GNorm * width * LNorm);
628   CMPLX_MULT_SELF_SCALAR(yAc->yIgVgb, GNorm * width * LNorm);
629   pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
630 
631   return (AcAnalysisMethod);
632 }
633 
634 BOOLEAN
TWOsorSolve(TWOdevice * pDevice,double * xReal,double * xImag,double omega)635 TWOsorSolve(TWOdevice *pDevice, double *xReal, double *xImag,
636             double omega)
637 {
638   double dxdy;
639   double wRelax = 1.0;		/* SOR relaxation parameter */
640   double *rhsReal = pDevice->rhs;
641   double *rhsSOR = pDevice->rhsImag;
642   BOOLEAN SORConverged = FALSE;
643   BOOLEAN SORFailed = FALSE;
644   int numEqns = pDevice->numEqns;
645   int iterationNum;
646   int indexN, indexP;
647   int index, eIndex;
648   TWOnode *pNode;
649   TWOelem *pElem;
650 
651   /* clear xReal and xImag arrays */
652   for (index = 1; index <= numEqns; index++) {
653     xReal[index] = 0.0;
654     xImag[index] = 0.0;
655   }
656 
657   iterationNum = 1;
658   for (; (!SORConverged) &&(!SORFailed); iterationNum++) {
659     for (index = 1; index <= numEqns; index++) {
660       rhsSOR[index] = 0.0;
661     }
662     for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
663       pElem = pDevice->elements[eIndex];
664       dxdy = 0.25 * pElem->dx * pElem->dy;
665       for (index = 0; index <= 3; index++) {
666 	pNode = pElem->pNodes[index];
667 	if ((pNode->nodeType != CONTACT) && (pElem->elemType == SEMICON)) {
668 	  if (!OneCarrier) {
669 	    indexN = pNode->nEqn;
670 	    indexP = pNode->pEqn;
671 	    rhsSOR[indexN] -= dxdy * omega * xImag[indexN];
672 	    rhsSOR[indexP] += dxdy * omega * xImag[indexP];
673 	  } else if (OneCarrier == N_TYPE) {
674 	    indexN = pNode->nEqn;
675 	    rhsSOR[indexN] -= dxdy * omega * xImag[indexN];
676 	  } else if (OneCarrier == P_TYPE) {
677 	    indexP = pNode->pEqn;
678 	    rhsSOR[indexP] += dxdy * omega * xImag[indexP];
679 	  }
680 	}
681       }
682     }
683 
684     /* now add the terms from rhs to rhsImag */
685     for (index = 1; index <= numEqns; index++) {
686       rhsSOR[index] += rhsReal[index];
687     }
688 
689     /* compute xReal(k+1). solution stored in rhsImag */
690     spSolve(pDevice->matrix, rhsSOR, rhsSOR, NULL, NULL);
691     /* modify solution when wRelax is not 1 */
692     if (wRelax != 1) {
693       for (index = 1; index <= numEqns; index++) {
694 	rhsSOR[index] = (1 - wRelax) * xReal[index] +
695 	    wRelax * rhsSOR[index];
696       }
697     }
698     if (iterationNum > 1) {
699       SORConverged = hasSORConverged(xReal, rhsSOR, numEqns);
700     }
701     /* copy real solution into xReal */
702     for (index = 1; index <= numEqns; index++) {
703       xReal[index] = rhsSOR[index];
704     }
705 
706     /* now compute the imaginary part of the solution xImag */
707     for (index = 1; index <= numEqns; index++) {
708       rhsSOR[index] = 0.0;
709     }
710     for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
711       pElem = pDevice->elements[eIndex];
712       dxdy = 0.25 * pElem->dx * pElem->dy;
713       for (index = 0; index <= 3; index++) {
714 	pNode = pElem->pNodes[index];
715 	if ((pNode->nodeType != CONTACT) && (pElem->elemType == SEMICON)) {
716 	  if (!OneCarrier) {
717 	    indexN = pNode->nEqn;
718 	    indexP = pNode->pEqn;
719 	    rhsSOR[indexN] += dxdy * omega * xReal[indexN];
720 	    rhsSOR[indexP] -= dxdy * omega * xReal[indexP];
721 	  } else if (OneCarrier == N_TYPE) {
722 	    indexN = pNode->nEqn;
723 	    rhsSOR[indexN] += dxdy * omega * xReal[indexN];
724 	  } else if (OneCarrier ==  P_TYPE) {
725 	    indexP = pNode->pEqn;
726 	    rhsSOR[indexP] -= dxdy * omega * xReal[indexP];
727 	  }
728 	}
729       }
730     }
731     /* compute xImag(k+1) */
732     spSolve(pDevice->matrix, rhsSOR, rhsSOR, NULL, NULL);
733     /* modify solution when wRelax is not 1 */
734     if (wRelax != 1) {
735       for (index = 1; index <= numEqns; index++) {
736 	rhsSOR[index] = (1 - wRelax) * xImag[index] +
737 	    wRelax * rhsSOR[index];
738       }
739     }
740     if (iterationNum > 1) {
741       SORConverged = SORConverged && hasSORConverged(xImag, rhsSOR, numEqns);
742     }
743     /* copy imag solution into xImag */
744     for (index = 1; index <= numEqns; index++) {
745       xImag[index] = rhsSOR[index];
746     }
747     if ((iterationNum > 4) && !SORConverged) {
748       SORFailed = TRUE;
749     }
750     if (TWOacDebug)
751       printf("SOR iteration number = %d\n", iterationNum);
752   }
753   return (SORFailed);
754 }
755 
756 
757 SPcomplex *
contactAdmittance(TWOdevice * pDevice,TWOcontact * pContact,BOOLEAN delVContact,double * xReal,double * xImag,SPcomplex * cOmega)758 contactAdmittance(TWOdevice *pDevice, TWOcontact *pContact, BOOLEAN delVContact,
759                   double *xReal, double *xImag, SPcomplex *cOmega)
760 {
761   TWOnode *pNode, *pHNode = NULL, *pVNode = NULL;
762   TWOedge *pHEdge = NULL, *pVEdge = NULL;
763   int index, i, indexPsi, indexN, indexP, numContactNodes;
764   TWOelem *pElem;
765   SPcomplex psiAc, nAc, pAc;
766   SPcomplex prod1, prod2, sum;
767   double temp;
768 
769   NG_IGNORE(pDevice);
770 
771   CMPLX_ASSIGN_VALUE(yTotal, 0.0, 0.0);
772 
773   numContactNodes = pContact->numNodes;
774   for (index = 0; index < numContactNodes; index++) {
775     pNode = pContact->pNodes[index];
776     for (i = 0; i <= 3; i++) {
777       pElem = pNode->pElems[i];
778       if (pElem != NULL) {
779 	switch (i) {
780 	case 0:
781 	  /* the TL element */
782 	  pHNode = pElem->pBLNode;
783 	  pVNode = pElem->pTRNode;
784 	  pHEdge = pElem->pBotEdge;
785 	  pVEdge = pElem->pRightEdge;
786 	  if (pElem->elemType == SEMICON) {
787 	    /* compute the derivatives with n,p */
788 	    if (pHNode->nodeType != CONTACT) {
789 	      indexN = pHNode->nEqn;
790 	      indexP = pHNode->pEqn;
791 	      CMPLX_ASSIGN_VALUE(nAc, xReal[indexN], xImag[indexN]);
792 	      CMPLX_ASSIGN_VALUE(pAc, xReal[indexP], xImag[indexP]);
793 	      CMPLX_MULT_SCALAR(prod1, nAc, pHEdge->dJnDn);
794 	      CMPLX_MULT_SCALAR(prod2, pAc, pHEdge->dJpDp);
795 	      CMPLX_ADD(sum, prod1, prod2);
796 	      CMPLX_MULT_SCALAR(prod1, sum, 0.5 * pElem->dy);
797 	      CMPLX_SUBT_ASSIGN(yTotal, prod1);
798 	    }
799 	    if (pVNode->nodeType != CONTACT) {
800 	      indexN = pVNode->nEqn;
801 	      indexP = pVNode->pEqn;
802 	      CMPLX_ASSIGN_VALUE(nAc, xReal[indexN], xImag[indexN]);
803 	      CMPLX_ASSIGN_VALUE(pAc, xReal[indexP], xImag[indexP]);
804 	      CMPLX_MULT_SCALAR(prod1, nAc, pVEdge->dJnDn);
805 	      CMPLX_MULT_SCALAR(prod2, pAc, pVEdge->dJpDp);
806 	      CMPLX_ADD(sum, prod1, prod2);
807 	      CMPLX_MULT_SCALAR(prod1, sum, 0.5 * pElem->dx);
808 	      CMPLX_SUBT_ASSIGN(yTotal, prod1);
809 	    }
810 	  }
811 	  break;
812 	case 1:
813 	  /* the TR element */
814 	  pHNode = pElem->pBRNode;
815 	  pVNode = pElem->pTLNode;
816 	  pHEdge = pElem->pBotEdge;
817 	  pVEdge = pElem->pLeftEdge;
818 	  if (pElem->elemType == SEMICON) {
819 	    /* compute the derivatives with n,p */
820 	    if (pHNode->nodeType != CONTACT) {
821 	      indexN = pHNode->nEqn;
822 	      indexP = pHNode->pEqn;
823 	      CMPLX_ASSIGN_VALUE(nAc, xReal[indexN], xImag[indexN]);
824 	      CMPLX_ASSIGN_VALUE(pAc, xReal[indexP], xImag[indexP]);
825 	      CMPLX_MULT_SCALAR(prod1, nAc, pHEdge->dJnDnP1);
826 	      CMPLX_MULT_SCALAR(prod2, pAc, pHEdge->dJpDpP1);
827 	      CMPLX_ADD(sum, prod1, prod2);
828 	      CMPLX_MULT_SCALAR(prod1, sum, 0.5 * pElem->dy);
829 	      CMPLX_ADD_ASSIGN(yTotal, prod1);
830 	    }
831 	    if (pVNode->nodeType != CONTACT) {
832 	      indexN = pVNode->nEqn;
833 	      indexP = pVNode->pEqn;
834 	      CMPLX_ASSIGN_VALUE(nAc, xReal[indexN], xImag[indexN]);
835 	      CMPLX_ASSIGN_VALUE(pAc, xReal[indexP], xImag[indexP]);
836 	      CMPLX_MULT_SCALAR(prod1, nAc, pVEdge->dJnDn);
837 	      CMPLX_MULT_SCALAR(prod2, pAc, pVEdge->dJpDp);
838 	      CMPLX_ADD(sum, prod1, prod2);
839 	      CMPLX_MULT_SCALAR(prod1, sum, 0.5 * pElem->dx);
840 	      CMPLX_SUBT_ASSIGN(yTotal, prod1);
841 	    }
842 	  }
843 	  break;
844 	case 2:
845 	  /* the BR element */
846 	  pHNode = pElem->pTRNode;
847 	  pVNode = pElem->pBLNode;
848 	  pHEdge = pElem->pTopEdge;
849 	  pVEdge = pElem->pLeftEdge;
850 	  if (pElem->elemType == SEMICON) {
851 	    /* compute the derivatives with n,p */
852 	    if (pHNode->nodeType != CONTACT) {
853 	      indexN = pHNode->nEqn;
854 	      indexP = pHNode->pEqn;
855 	      CMPLX_ASSIGN_VALUE(nAc, xReal[indexN], xImag[indexN]);
856 	      CMPLX_ASSIGN_VALUE(pAc, xReal[indexP], xImag[indexP]);
857 	      CMPLX_MULT_SCALAR(prod1, nAc, pHEdge->dJnDnP1);
858 	      CMPLX_MULT_SCALAR(prod2, pAc, pHEdge->dJpDpP1);
859 	      CMPLX_ADD(sum, prod1, prod2);
860 	      CMPLX_MULT_SCALAR(prod1, sum, 0.5 * pElem->dy);
861 	      CMPLX_ADD_ASSIGN(yTotal, prod1);
862 	    }
863 	    if (pVNode->nodeType != CONTACT) {
864 	      indexN = pVNode->nEqn;
865 	      indexP = pVNode->pEqn;
866 	      CMPLX_ASSIGN_VALUE(nAc, xReal[indexN], xImag[indexN]);
867 	      CMPLX_ASSIGN_VALUE(pAc, xReal[indexP], xImag[indexP]);
868 	      CMPLX_MULT_SCALAR(prod1, nAc, pVEdge->dJnDnP1);
869 	      CMPLX_MULT_SCALAR(prod2, pAc, pVEdge->dJpDpP1);
870 	      CMPLX_ADD(sum, prod1, prod2);
871 	      CMPLX_MULT_SCALAR(prod1, sum, 0.5 * pElem->dx);
872 	      CMPLX_ADD_ASSIGN(yTotal, prod1);
873 	    }
874 	  }
875 	  break;
876 	case 3:
877 	  /* the BL element */
878 	  pHNode = pElem->pTLNode;
879 	  pVNode = pElem->pBRNode;
880 	  pHEdge = pElem->pTopEdge;
881 	  pVEdge = pElem->pRightEdge;
882 	  if (pElem->elemType == SEMICON) {
883 	    /* compute the derivatives with n,p */
884 	    if (pHNode->nodeType != CONTACT) {
885 	      indexN = pHNode->nEqn;
886 	      indexP = pHNode->pEqn;
887 	      CMPLX_ASSIGN_VALUE(nAc, xReal[indexN], xImag[indexN]);
888 	      CMPLX_ASSIGN_VALUE(pAc, xReal[indexP], xImag[indexP]);
889 	      CMPLX_MULT_SCALAR(prod1, nAc, pHEdge->dJnDn);
890 	      CMPLX_MULT_SCALAR(prod2, pAc, pHEdge->dJpDp);
891 	      CMPLX_ADD(sum, prod1, prod2);
892 	      CMPLX_MULT_SCALAR(prod1, sum, 0.5 * pElem->dy);
893 	      CMPLX_SUBT_ASSIGN(yTotal, prod1);
894 	    }
895 	    if (pVNode->nodeType != CONTACT) {
896 	      indexN = pVNode->nEqn;
897 	      indexP = pVNode->pEqn;
898 	      CMPLX_ASSIGN_VALUE(nAc, xReal[indexN], xImag[indexN]);
899 	      CMPLX_ASSIGN_VALUE(pAc, xReal[indexP], xImag[indexP]);
900 	      CMPLX_MULT_SCALAR(prod1, nAc, pVEdge->dJnDnP1);
901 	      CMPLX_MULT_SCALAR(prod2, pAc, pVEdge->dJpDpP1);
902 	      CMPLX_ADD(sum, prod1, prod2);
903 	      CMPLX_MULT_SCALAR(prod1, sum, 0.5 * pElem->dx);
904 	      CMPLX_ADD_ASSIGN(yTotal, prod1);
905 	    }
906 	  }
907 	  break;
908 	}
909 	if (pElem->elemType == SEMICON) {
910 	  if (pHNode->nodeType != CONTACT) {
911 	    indexPsi = pHNode->psiEqn;
912 	    CMPLX_ASSIGN_VALUE(psiAc, xReal[indexPsi], xImag[indexPsi]);
913 	    temp = 0.5 * pElem->dy * (pHEdge->dJnDpsiP1 + pHEdge->dJpDpsiP1);
914 	    CMPLX_MULT_SCALAR(prod1, psiAc, temp);
915 	    CMPLX_ADD_ASSIGN(yTotal, prod1);
916 	    if (delVContact) {
917 	      CMPLX_ADD_SELF_SCALAR(yTotal, -temp);
918 	    }
919 	  }
920 	  if (pVNode->nodeType != CONTACT) {
921 	    indexPsi = pVNode->psiEqn;
922 	    CMPLX_ASSIGN_VALUE(psiAc, xReal[indexPsi], xImag[indexPsi]);
923 	    temp = 0.5 * pElem->dx * (pVEdge->dJnDpsiP1 + pVEdge->dJpDpsiP1);
924 	    CMPLX_MULT_SCALAR(prod1, psiAc, temp);
925 	    CMPLX_ADD_ASSIGN(yTotal, prod1);
926 	    if (delVContact) {
927 	      CMPLX_ADD_SELF_SCALAR(yTotal, -temp);
928 	    }
929 	  }
930 	}
931 	/* displacement current terms */
932 	if (pHNode->nodeType != CONTACT) {
933 	  indexPsi = pHNode->psiEqn;
934 	  CMPLX_ASSIGN_VALUE(psiAc, xReal[indexPsi], xImag[indexPsi]);
935 	  CMPLX_MULT_SCALAR(prod1, *cOmega, pElem->epsRel * 0.5 * pElem->dyOverDx);
936 	  CMPLX_MULT(prod2, prod1, psiAc);
937 	  CMPLX_SUBT_ASSIGN(yTotal, prod2);
938 	  if (delVContact) {
939 	    CMPLX_ADD_ASSIGN(yTotal, prod1);
940 	  }
941 	}
942 	if (pVNode->nodeType != CONTACT) {
943 	  indexPsi = pVNode->psiEqn;
944 	  CMPLX_ASSIGN_VALUE(psiAc, xReal[indexPsi], xImag[indexPsi]);
945 	  CMPLX_MULT_SCALAR(prod1, *cOmega, pElem->epsRel * 0.5 * pElem->dxOverDy);
946 	  CMPLX_MULT(prod2, prod1, psiAc);
947 	  CMPLX_SUBT_ASSIGN(yTotal, prod2);
948 	  if (delVContact) {
949 	    CMPLX_ADD_ASSIGN(yTotal, prod1);
950 	  }
951 	}
952       }
953     }
954   }
955   return (&yTotal); /* XXX */
956 }
957 
958 
959 SPcomplex *
oxideAdmittance(TWOdevice * pDevice,TWOcontact * pContact,BOOLEAN delVContact,double * xReal,double * xImag,SPcomplex * cOmega)960 oxideAdmittance(TWOdevice *pDevice, TWOcontact *pContact, BOOLEAN delVContact,
961                 double *xReal, double *xImag, SPcomplex *cOmega)
962 {
963   TWOnode *pNode, *pHNode = NULL, *pVNode = NULL;
964   TWOedge *pHEdge, *pVEdge;
965   int index, i, indexPsi, numContactNodes;
966   TWOelem *pElem;
967   SPcomplex psiAc;
968   SPcomplex prod1, prod2;
969 
970   NG_IGNORE(pDevice);
971 
972   CMPLX_ASSIGN_VALUE(yTotal, 0.0, 0.0);
973 
974   numContactNodes = pContact->numNodes;
975   for (index = 0; index < numContactNodes; index++) {
976     pNode = pContact->pNodes[index];
977     for (i = 0; i <= 3; i++) {
978       pElem = pNode->pElems[i];
979       if (pElem != NULL) {
980 	switch (i) {
981 	case 0:
982 	  /* the TL element */
983 	  pHNode = pElem->pBLNode;
984 	  pVNode = pElem->pTRNode;
985 	  pHEdge = pElem->pBotEdge;
986 	  pVEdge = pElem->pRightEdge;
987 	  break;
988 	case 1:
989 	  /* the TR element */
990 	  pHNode = pElem->pBRNode;
991 	  pVNode = pElem->pTLNode;
992 	  pHEdge = pElem->pBotEdge;
993 	  pVEdge = pElem->pLeftEdge;
994 	  break;
995 	case 2:
996 	  /* the BR element */
997 	  pHNode = pElem->pTRNode;
998 	  pVNode = pElem->pBLNode;
999 	  pHEdge = pElem->pTopEdge;
1000 	  pVEdge = pElem->pLeftEdge;
1001 	  break;
1002 	case 3:
1003 	  /* the BL element */
1004 	  pHNode = pElem->pTLNode;
1005 	  pVNode = pElem->pBRNode;
1006 	  pHEdge = pElem->pTopEdge;
1007 	  pVEdge = pElem->pRightEdge;
1008 	  break;
1009 	}
1010 	/* displacement current terms */
1011 	if (pHNode->nodeType != CONTACT) {
1012 	  indexPsi = pHNode->psiEqn;
1013 	  CMPLX_ASSIGN_VALUE(psiAc, xReal[indexPsi], xImag[indexPsi]);
1014 	  CMPLX_MULT_SCALAR(prod1, *cOmega, pElem->epsRel * 0.5 * pElem->dyOverDx);
1015 	  CMPLX_MULT(prod2, prod1, psiAc);
1016 	  CMPLX_SUBT_ASSIGN(yTotal, prod2);
1017 	  if (delVContact) {
1018 	    CMPLX_ADD_ASSIGN(yTotal, prod1);
1019 	  }
1020 	}
1021 	if (pVNode->nodeType != CONTACT) {
1022 	  indexPsi = pVNode->psiEqn;
1023 	  CMPLX_ASSIGN_VALUE(psiAc, xReal[indexPsi], xImag[indexPsi]);
1024 	  CMPLX_MULT_SCALAR(prod1, *cOmega, pElem->epsRel * 0.5 * pElem->dxOverDy);
1025 	  CMPLX_MULT(prod2, prod1, psiAc);
1026 	  CMPLX_SUBT_ASSIGN(yTotal, prod2);
1027 	  if (delVContact) {
1028 	    CMPLX_ADD_ASSIGN(yTotal, prod1);
1029 	  }
1030 	}
1031       }
1032     }
1033   }
1034   return (&yTotal);
1035 }
1036 
1037 void
NUMD2ys(TWOdevice * pDevice,SPcomplex * s,SPcomplex * yIn)1038 NUMD2ys(TWOdevice *pDevice, SPcomplex *s, SPcomplex *yIn)
1039 {
1040   TWOnode *pNode;
1041   TWOelem *pElem;
1042   int index, eIndex;
1043   double dxdy;
1044   double *solnReal, *solnImag;
1045   double *rhsReal, *rhsImag;
1046   SPcomplex yAc, *y;
1047   BOOLEAN deltaVContact = FALSE;
1048   SPcomplex temp, cOmega;
1049 
1050   /*
1051    * change context names of solution vectors for ac analysis dcDeltaSolution
1052    * stores the real part and copiedSolution stores the imaginary part of the
1053    * ac solution vector
1054    */
1055   pDevice->solverType = SLV_SMSIG;
1056   rhsReal = pDevice->rhs;
1057   rhsImag = pDevice->rhsImag;
1058   solnReal = pDevice->dcDeltaSolution;
1059   solnImag = pDevice->copiedSolution;
1060 
1061   /* use a normalized radian frequency */
1062   CMPLX_MULT_SCALAR(cOmega, *s, TNorm);
1063   for (index = 1; index <= pDevice->numEqns; index++) {
1064     rhsImag[index] = 0.0;
1065   }
1066   /* solve the system of equations directly */
1067   if (!OneCarrier) {
1068     TWO_jacLoad(pDevice);
1069   } else if (OneCarrier == N_TYPE) {
1070     TWONjacLoad(pDevice);
1071   } else if (OneCarrier == P_TYPE) {
1072     TWOPjacLoad(pDevice);
1073   }
1074   storeNewRhs(pDevice, pDevice->pLastContact);
1075 
1076   spSetComplex(pDevice->matrix);
1077   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
1078     pElem = pDevice->elements[eIndex];
1079     if (pElem->elemType == SEMICON) {
1080       dxdy = 0.25 * pElem->dx * pElem->dy;
1081       for (index = 0; index <= 3; index++) {
1082 	pNode = pElem->pNodes[index];
1083 	if (pNode->nodeType != CONTACT) {
1084 	  if (!OneCarrier) {
1085 	    CMPLX_MULT_SCALAR(temp, cOmega, dxdy);
1086 	    spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
1087 	    spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
1088 	  } else if (OneCarrier == N_TYPE) {
1089 	    CMPLX_MULT_SCALAR(temp, cOmega, dxdy);
1090 	    spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
1091 	  } else if (OneCarrier == P_TYPE) {
1092 	    CMPLX_MULT_SCALAR(temp, cOmega, dxdy);
1093 	    spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
1094 	  }
1095 	}
1096       }
1097     }
1098   }
1099 
1100   spFactor(pDevice->matrix);
1101   spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
1102   y = contactAdmittance(pDevice, pDevice->pFirstContact, deltaVContact,
1103       solnReal, solnImag, &cOmega);
1104   CMPLX_ASSIGN_VALUE(yAc, y->real, y->imag);
1105   CMPLX_ASSIGN(*yIn, yAc);
1106   CMPLX_NEGATE_SELF(*yIn);
1107   CMPLX_MULT_SELF_SCALAR(*yIn, GNorm * pDevice->width * LNorm);
1108 }
1109 
1110 
1111 void
NBJT2ys(TWOdevice * pDevice,SPcomplex * s,SPcomplex * yIeVce,SPcomplex * yIcVce,SPcomplex * yIeVbe,SPcomplex * yIcVbe)1112 NBJT2ys(TWOdevice *pDevice, SPcomplex *s, SPcomplex *yIeVce, SPcomplex *yIcVce,
1113         SPcomplex *yIeVbe, SPcomplex *yIcVbe)
1114 {
1115   TWOcontact *pEmitContact = pDevice->pLastContact;
1116   TWOcontact *pColContact = pDevice->pFirstContact;
1117   TWOcontact *pBaseContact = pDevice->pFirstContact->next;
1118   TWOnode *pNode;
1119   TWOelem *pElem;
1120   int index, eIndex;
1121   double width = pDevice->width;
1122   double dxdy;
1123   double *solnReal, *solnImag;
1124   double *rhsReal, *rhsImag;
1125   SPcomplex *y;
1126   SPcomplex pIeVce, pIcVce, pIeVbe, pIcVbe;
1127   SPcomplex temp, cOmega;
1128 
1129   pDevice->solverType = SLV_SMSIG;
1130   rhsReal = pDevice->rhs;
1131   rhsImag = pDevice->rhsImag;
1132   solnReal = pDevice->dcDeltaSolution;
1133   solnImag = pDevice->copiedSolution;
1134 
1135   /* use a normalized radian frequency */
1136   CMPLX_MULT_SCALAR(cOmega, *s, TNorm);
1137 
1138   for (index = 1; index <= pDevice->numEqns; index++) {
1139     rhsImag[index] = 0.0;
1140   }
1141   /* solve the system of equations directly */
1142   if (!OneCarrier) {
1143     TWO_jacLoad(pDevice);
1144   } else if (OneCarrier == N_TYPE) {
1145     TWONjacLoad(pDevice);
1146   } else if (OneCarrier == P_TYPE) {
1147     TWOPjacLoad(pDevice);
1148   }
1149   storeNewRhs(pDevice, pColContact);
1150   spSetComplex(pDevice->matrix);
1151   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
1152     pElem = pDevice->elements[eIndex];
1153     if (pElem->elemType == SEMICON) {
1154       dxdy = 0.25 * pElem->dx * pElem->dy;
1155       for (index = 0; index <= 3; index++) {
1156 	pNode = pElem->pNodes[index];
1157 	if (pNode->nodeType != CONTACT) {
1158 	  if (!OneCarrier) {
1159 	    CMPLX_MULT_SCALAR(temp, cOmega, dxdy);
1160 	    spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
1161 	    spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
1162 	  } else if (OneCarrier == N_TYPE) {
1163 	    CMPLX_MULT_SCALAR(temp, cOmega, dxdy);
1164 	    spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
1165 	  } else if (OneCarrier == P_TYPE) {
1166 	    CMPLX_MULT_SCALAR(temp, cOmega, dxdy);
1167 	    spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
1168 	  }
1169 	}
1170       }
1171     }
1172   }
1173   spFactor(pDevice->matrix);
1174   spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
1175 
1176   y = contactAdmittance(pDevice, pEmitContact, FALSE,
1177       solnReal, solnImag, &cOmega);
1178   CMPLX_ASSIGN_VALUE(pIeVce, y->real, y->imag);
1179   y = contactAdmittance(pDevice, pColContact, TRUE,
1180       solnReal, solnImag, &cOmega);
1181   CMPLX_ASSIGN_VALUE(pIcVce, y->real, y->imag);
1182   for (index = 1; index <= pDevice->numEqns; index++) {
1183     rhsImag[index] = 0.0;
1184   }
1185   storeNewRhs(pDevice, pBaseContact);
1186   /* don't need to LU factor the jacobian since it exists */
1187   spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
1188   y = contactAdmittance(pDevice, pEmitContact, FALSE,
1189       solnReal, solnImag, &cOmega);
1190   CMPLX_ASSIGN_VALUE(pIeVbe, y->real, y->imag);
1191   y = contactAdmittance(pDevice, pColContact, FALSE,
1192       solnReal, solnImag, &cOmega);
1193   CMPLX_ASSIGN_VALUE(pIcVbe, y->real, y->imag);
1194 
1195 
1196   CMPLX_ASSIGN(*yIeVce, pIeVce);
1197   CMPLX_ASSIGN(*yIeVbe, pIeVbe);
1198   CMPLX_ASSIGN(*yIcVce, pIcVce);
1199   CMPLX_ASSIGN(*yIcVbe, pIcVbe);
1200   CMPLX_MULT_SELF_SCALAR(*yIeVce, GNorm * width * LNorm);
1201   CMPLX_MULT_SELF_SCALAR(*yIeVbe, GNorm * width * LNorm);
1202   CMPLX_MULT_SELF_SCALAR(*yIcVce, GNorm * width * LNorm);
1203   CMPLX_MULT_SELF_SCALAR(*yIcVbe, GNorm * width * LNorm);
1204 }
1205 
1206 void
NUMOSys(TWOdevice * pDevice,SPcomplex * s,struct mosAdmittances * yAc)1207 NUMOSys(TWOdevice *pDevice, SPcomplex *s, struct mosAdmittances *yAc)
1208 {
1209   TWOcontact *pDContact = pDevice->pFirstContact;
1210   TWOcontact *pGContact = pDevice->pFirstContact->next;
1211   TWOcontact *pSContact = pDevice->pFirstContact->next->next;
1212 /*  TWOcontact *pBContact = pDevice->pLastContact; */
1213   TWOnode *pNode;
1214   TWOelem *pElem;
1215   int index, eIndex;
1216   double width = pDevice->width;
1217   double dxdy;
1218   double *rhsReal, *rhsImag;
1219   double *solnReal, *solnImag;
1220   SPcomplex *y;
1221   SPcomplex temp, cOmega;
1222 
1223   pDevice->solverType = SLV_SMSIG;
1224   rhsReal = pDevice->rhs;
1225   rhsImag = pDevice->rhsImag;
1226   solnReal = pDevice->dcDeltaSolution;
1227   solnImag = pDevice->copiedSolution;
1228 
1229   /* use a normalized radian frequency */
1230   CMPLX_MULT_SCALAR(cOmega, *s, TNorm);
1231   for (index = 1; index <= pDevice->numEqns; index++) {
1232     rhsImag[index] = 0.0;
1233   }
1234   /* solve the system of equations directly */
1235   if (!OneCarrier) {
1236     TWO_jacLoad(pDevice);
1237   } else if (OneCarrier == N_TYPE) {
1238     TWONjacLoad(pDevice);
1239   } else if (OneCarrier == P_TYPE) {
1240     TWOPjacLoad(pDevice);
1241   }
1242   storeNewRhs(pDevice, pDContact);
1243   spSetComplex(pDevice->matrix);
1244 
1245   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
1246     pElem = pDevice->elements[eIndex];
1247     if (pElem->elemType == SEMICON) {
1248       dxdy = 0.25 * pElem->dx * pElem->dy;
1249       for (index = 0; index <= 3; index++) {
1250 	pNode = pElem->pNodes[index];
1251 	if (pNode->nodeType != CONTACT) {
1252 	  if (!OneCarrier) {
1253 	    CMPLX_MULT_SCALAR(temp, cOmega, dxdy);
1254 	    spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
1255 	    spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
1256 	  } else if (OneCarrier == N_TYPE) {
1257 	    CMPLX_MULT_SCALAR(temp, cOmega, dxdy);
1258 	    spADD_COMPLEX_ELEMENT(pNode->fNN, -temp.real, -temp.imag);
1259 	  } else if (OneCarrier == P_TYPE) {
1260 	    CMPLX_MULT_SCALAR(temp, cOmega, dxdy);
1261 	    spADD_COMPLEX_ELEMENT(pNode->fPP, temp.real, temp.imag);
1262 	  }
1263 	}
1264       }
1265     }
1266   }
1267 
1268   spFactor(pDevice->matrix);
1269   spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
1270 
1271   y = contactAdmittance(pDevice, pDContact, TRUE,
1272       solnReal, solnImag, &cOmega);
1273   CMPLX_ASSIGN_VALUE(yAc->yIdVdb, y->real, y->imag);
1274   y = contactAdmittance(pDevice, pSContact, FALSE,
1275       solnReal, solnImag, &cOmega);
1276   CMPLX_ASSIGN_VALUE(yAc->yIsVdb, y->real, y->imag);
1277   y = GateTypeAdmittance(pDevice, pGContact, FALSE,
1278       solnReal, solnImag, &cOmega);
1279   CMPLX_ASSIGN_VALUE(yAc->yIgVdb, y->real, y->imag);
1280 
1281   for (index = 1; index <= pDevice->numEqns; index++) {
1282     rhsImag[index] = 0.0;
1283   }
1284   storeNewRhs(pDevice, pSContact);
1285   /* don't need to LU factor the jacobian since it exists */
1286   spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
1287   y = contactAdmittance(pDevice, pDContact, FALSE,
1288       solnReal, solnImag, &cOmega);
1289   CMPLX_ASSIGN_VALUE(yAc->yIdVsb, y->real, y->imag);
1290   y = contactAdmittance(pDevice, pSContact, TRUE,
1291       solnReal, solnImag, &cOmega);
1292   CMPLX_ASSIGN_VALUE(yAc->yIsVsb, y->real, y->imag);
1293   y = GateTypeAdmittance(pDevice, pGContact, FALSE,
1294       solnReal, solnImag, &cOmega);
1295   CMPLX_ASSIGN_VALUE(yAc->yIgVsb, y->real, y->imag);
1296   for (index = 1; index <= pDevice->numEqns; index++) {
1297     rhsImag[index] = 0.0;
1298   }
1299   storeNewRhs(pDevice, pGContact);
1300   spSolve(pDevice->matrix, rhsReal, solnReal, rhsImag, solnImag);
1301   y = contactAdmittance(pDevice, pDContact, FALSE,
1302       solnReal, solnImag, &cOmega);
1303   CMPLX_ASSIGN_VALUE(yAc->yIdVgb, y->real, y->imag);
1304   y = contactAdmittance(pDevice, pSContact, FALSE,
1305       solnReal, solnImag, &cOmega);
1306   CMPLX_ASSIGN_VALUE(yAc->yIsVgb, y->real, y->imag);
1307   y = GateTypeAdmittance(pDevice, pGContact, TRUE,
1308       solnReal, solnImag, &cOmega);
1309   CMPLX_ASSIGN_VALUE(yAc->yIgVgb, y->real, y->imag);
1310 
1311   CMPLX_MULT_SELF_SCALAR(yAc->yIdVdb, GNorm * width * LNorm);
1312   CMPLX_MULT_SELF_SCALAR(yAc->yIdVsb, GNorm * width * LNorm);
1313   CMPLX_MULT_SELF_SCALAR(yAc->yIdVgb, GNorm * width * LNorm);
1314   CMPLX_MULT_SELF_SCALAR(yAc->yIsVdb, GNorm * width * LNorm);
1315   CMPLX_MULT_SELF_SCALAR(yAc->yIsVsb, GNorm * width * LNorm);
1316   CMPLX_MULT_SELF_SCALAR(yAc->yIsVgb, GNorm * width * LNorm);
1317   CMPLX_MULT_SELF_SCALAR(yAc->yIgVdb, GNorm * width * LNorm);
1318   CMPLX_MULT_SELF_SCALAR(yAc->yIgVsb, GNorm * width * LNorm);
1319   CMPLX_MULT_SELF_SCALAR(yAc->yIgVgb, GNorm * width * LNorm);
1320 }
1321