1 /***************************************************************************
2 JSPICE3 adaptation of Spice3f2 - Copyright (c) Stephen R. Whiteley 1992
3 Copyright 1990 Regents of the University of California. All rights reserved.
4 Authors: UCB CAD Group
5 1993 Stephen R. Whiteley
6 ****************************************************************************/
7
8 #include "spice.h"
9 #include "spmatrix.h"
10 #include "util.h"
11 #include "srcdefs.h"
12 #include "outdata.h"
13 #include "sensdefs.h"
14 #include "sensgen.h"
15 #include "uflags.h"
16
17
18 #define NEWN(TYPE,COUNT) ((TYPE *) MALLOC(sizeof(TYPE) * (COUNT)))
19 #ifndef M_LOG2E
20 #define M_LOG2E 1.4426950408889634074
21 #endif
22
23 static double Sens_Delta = 0.000001;
24 static double Sens_Abs_Delta = 0.000001;
25
26 #ifdef __STDC__
27 static int sens_acoperation(CKTcircuit*,struct sOUTdata*,int);
28 static int sens_dcoperation(CKTcircuit*,struct sOUTdata*,int);
29 static int sens_loadnew(sgen*,CKTcircuit*,int);
30 static int sens_setp(sgen*,CKTcircuit*,IFvalue*);
31 #else
32 static int sens_acoperation();
33 static int sens_dcoperation();
34 static int sens_loadnew();
35 static int sens_setp();
36 #endif
37
38 extern SPICEdev *DEVices[];
39
40 /*
41 * Procedure:
42 *
43 * Determine operating point (call CKTop)
44 *
45 * For each frequency point:
46 * (for AC) call NIacIter to get base node voltages
47 * For each element/parameter in the test list:
48 * construct the perturbation matrix
49 * Solve for the sensitivities:
50 * delta_E = Y^-1 (delta_I - delta_Y E)
51 * save results
52 */
53
54
55
56 int
SENSan(cktp,restart)57 SENSan(cktp, restart)
58
59 GENERIC *cktp;
60 int restart;
61 {
62 CKTcircuit *ckt = (CKTcircuit *)cktp;
63 SENSAN *job = (SENSAN *) ckt->CKTcurJob;
64 struct sSENSint *st = &job->ST;
65 sgen *sg;
66 static struct sOUTdata outd;
67 int i;
68 int is_dc;
69 int code;
70 int error;
71 char namebuf[513];
72
73 if (ckt->CKTjjPresent) {
74 (*(SPfrontEnd->IFerror))(ERR_FATAL,
75 "Sensitivity analysis not possible with Josephson junctions",
76 NULL);
77 return (OK);
78 }
79 if (restart) {
80
81 /* make sure sources are present */
82 code = CKTtypelook("Source");
83
84 if (job->SENSoutSrc) {
85 error = CKTfndDev((GENERIC*)ckt,&code,&job->SENSoutSrcDev,
86 job->SENSoutSrc, (GENERIC *)NULL,(IFuid)NULL);
87 if (error != 0) {
88 (*(SPfrontEnd->IFerror))(ERR_FATAL,
89 "Sensitivity Analysis source %s not in circuit",
90 &job->SENSoutSrc);
91 return (E_NOTFOUND);
92 }
93 }
94 if (job->DC.eltName[0] != NULL) {
95 /* DC source was given */
96
97 for (i = 0; i <= job->DC.nestLevel; i++) {
98 SRCinstance *here;
99
100 here = (SRCinstance*)NULL;
101 error = CKTfndDev((GENERIC*)ckt,&code,(GENERIC**)&here,
102 job->DC.eltName[i], (GENERIC *)NULL,(IFuid)NULL);
103 if (error) {
104 (*(SPfrontEnd->IFerror))(ERR_FATAL,
105 "DCtrCurv: source %s not in circuit",
106 &(job->DC.eltName[i]));
107 return (E_NOTFOUND);
108 }
109 job->DC.elt[i] = (GENinstance*)here;
110 job->DC.vsave[i] = here->SRCdcValue;
111 }
112 }
113
114 /* in case these never got freed */
115 if (st->dY) {
116 spDestroy(st->dY);
117 st->dY = NULL;
118 }
119 FREE(st->dIr);
120 FREE(st->dIi);
121 FREE(st->dIdYr);
122 FREE(st->dIdYi);
123 FREE(st->o_values);
124 FREE(st->o_cvalues);
125
126 is_dc = (job->AC.stepType == DCSTEP);
127
128 st->size = spGetSize(ckt->CKTmatrix, 1);
129
130 /* Create the perturbation matrix */
131 st->dY = spCreate(st->size, 1, &error);
132 if (error)
133 return (error);
134
135 /* Create extra rhs */
136 st->dIr = NEWN(double, st->size + 1);
137 st->dIi = NEWN(double, st->size + 1);
138 st->dIdYr = NEWN(double, st->size + 1);
139 st->dIdYi = NEWN(double, st->size + 1);
140
141 outd.numNames = 0;
142 for (sg = sgen_init(ckt, is_dc); sg; sgen_next(&sg)) {
143 outd.numNames++;
144 }
145 if (!outd.numNames)
146 return (E_NOTFOUND);
147
148 outd.dataNames = NEWN(IFuid,outd.numNames);
149 for (i = 0,sg = sgen_init(ckt, is_dc); sg; i++,sgen_next(&sg)) {
150 if (!sg->is_instparam) {
151 sprintf(namebuf, "%s:%s",
152 sg->instance->GENname,
153 sg->ptable[sg->param].keyword);
154 }
155 else if ((sg->ptable[sg->param].dataType
156 & IF_PRINCIPAL) && sg->is_principle == 1) {
157 sprintf(namebuf, "%s", sg->instance->GENname);
158 }
159 else {
160 sprintf(namebuf, "%s.%s",
161 sg->instance->GENname,
162 sg->ptable[sg->param].keyword);
163 }
164
165 (*SPfrontEnd->IFnewUid)((GENERIC *) ckt,
166 outd.dataNames + i, NULL,
167 namebuf, UID_OTHER, NULL);
168 }
169
170 if (is_dc) {
171 outd.dataType = IF_REAL;
172 (*(SPfrontEnd->IFnewUid))((GENERIC *)ckt,
173 &outd.refName,(IFuid )NULL,
174 "sweep", UID_OTHER, NULL);
175 }
176 else {
177 outd.dataType = IF_COMPLEX;
178 (*SPfrontEnd->IFnewUid)((GENERIC *) ckt,
179 &outd.refName, NULL,
180 "frequency", UID_OTHER, NULL);
181 }
182
183 outd.circuitPtr = (GENERIC *)ckt;
184 outd.analysisPtr = (GENERIC*)ckt->CKTcurJob;
185 outd.analName = ckt->CKTcurJob->JOBname;
186 outd.refType = IF_REAL;
187 outd.plotPtr = &job->SENSplot;
188 outd.initValue = 0;
189 outd.finalValue = 0;
190 outd.step = 0;
191 outd.numPts = 1;
192 outd.count = 0;
193 (*(SPfrontEnd->OUTbeginPlot))((GENERIC*)&outd);
194
195 FREE(outd.dataNames);
196
197 if (is_dc) {
198 st->o_values = NEWN(double, outd.numNames);
199 st->o_cvalues = NULL;
200 }
201 else {
202 st->o_values = NULL;
203 st->o_cvalues = NEWN(IFcomplex, outd.numNames);
204 }
205 }
206 error = DCTloop(sens_dcoperation,ckt,restart,&job->DC,&outd);
207
208 if (error == E_PAUSE)
209 return (error);
210
211 (*(SPfrontEnd->OUTendPlot))(*outd.plotPtr);
212
213 spDestroy(st->dY);
214 st->dY = NULL;
215 FREE(st->dIr);
216 FREE(st->dIi);
217 FREE(st->dIdYr);
218 FREE(st->dIdYi);
219 FREE(st->o_values);
220 FREE(st->o_cvalues);
221 return (error);
222 }
223
224
225 /* ARGSUSED */
226 static int
sens_dcoperation(ckt,outd,restart)227 sens_dcoperation(ckt,outd,restart)
228
229 CKTcircuit *ckt;
230 struct sOUTdata *outd;
231 int restart;
232 {
233 SENSAN *job = (SENSAN *) ckt->CKTcurJob;
234 struct sSENSint *st = &job->ST;
235 sgen *sg;
236 IFvalue value, nvalue;
237 double *E;
238 SMPmatrix *Y;
239 double delta_var;
240 int i, j;
241 int bypass;
242 int error;
243
244 /* find the operating point */
245
246 error = CKTic(ckt);
247 if (error)
248 return (error);
249
250 error = CKTop(ckt,
251 MODEDCOP | MODEINITJCT,
252 MODEDCOP | MODEINITFLOAT,
253 ckt->CKTdcMaxIter);
254 if (error)
255 return (error);
256
257 if (job->AC.stepType != DCSTEP) {
258 error = ACloop(sens_acoperation,ckt,restart,&job->AC,outd);
259 return (error);
260 }
261
262 bypass = ckt->CKTbypass;
263 ckt->CKTbypass = 0;
264
265 /* The unknown vector of node voltages overwrites rhs */
266 E = ckt->CKTrhs;
267 ckt->CKTrhsOld = E;
268 Y = ckt->CKTmatrix;
269
270 /* Use a different vector & matrix */
271 ckt->CKTrhs = st->dIr;
272 ckt->CKTmatrix = st->dY;
273
274 /* calc. effect of each param */
275 for (i = 0,sg = sgen_init(ckt, TRUE); sg; i++,sgen_next(&sg)) {
276
277 /* clear CKTmatrix, CKTrhs */
278 spSetReal(st->dY);
279 spClear(st->dY);
280 for (j = 0; j <= st->size; j++) {
281 st->dIr[j] = 0.0;
282 }
283
284 error = sens_loadnew(sg,ckt,TRUE);
285 if (error)
286 return (error);
287
288 /* Alter the parameter */
289
290 if (sg->value != 0.0)
291 delta_var = sg->value * Sens_Delta;
292 else
293 delta_var = Sens_Abs_Delta;
294
295 nvalue.rValue = sg->value + delta_var;
296
297 sens_setp(sg, ckt, &nvalue);
298 if (error)
299 return (error);
300
301 /* change sign of CKTmatrix, CKTrhs */
302 spConstMult(st->dY, -1.0);
303 for (j = 0; j <= st->size; j++) {
304 st->dIr[j] = -st->dIr[j];
305 }
306
307 error = sens_loadnew(sg,ckt,TRUE);
308 if (error)
309 return (error);
310
311 /* now have delta_Y = CKTmatrix, delta_I = CKTrhs */
312
313 /* reset parameter */
314 value.rValue = sg->value;
315 sens_setp(sg, ckt, &value);
316
317 /* delta_Y E */
318 spMultiply(st->dY, st->dIdYr, E, NULL, NULL);
319
320 /* delta_I - delta_Y E */
321 for (j = 0; j <= st->size; j++) {
322 st->dIr[j] -= st->dIdYr[j];
323 }
324
325 /* Solve; Y already factored */
326 spSolve(Y, st->dIr, st->dIr, NULL, NULL);
327
328 /* delta_I is now equal to delta_E */
329
330 st->dIr[0] = 0.0;
331 if (job->SENSoutName)
332 st->o_values[i] =
333 st->dIr[job->SENSoutPos->number] -
334 st->dIr[job->SENSoutNeg->number];
335 else {
336 st->o_values[i] =
337 st->dIr
338 [((SRCinstance*)job->SENSoutSrcDev)->SRCbranch];
339 }
340 st->o_values[i] /= delta_var;
341 }
342 nvalue.v.vec.rVec = st->o_values;
343
344 if (job->DC.elt[0])
345 value.rValue = ((SRCinstance*)job->DC.elt[0])->SRCdcValue;
346 else
347 value.rValue = 0;
348 OUTdata(job->SENSplot, &value, &nvalue);
349 outd->count++;
350
351 ckt->CKTrhs = E;
352 ckt->CKTmatrix = Y;
353 ckt->CKTbypass = bypass;
354
355 return (OK);
356 }
357
358
359 /* ARGSUSED */
360 static int
sens_acoperation(ckt,outd,restart)361 sens_acoperation(ckt,outd,restart)
362
363 CKTcircuit *ckt;
364 struct sOUTdata *outd;
365 int restart;
366 {
367 SENSAN *job = (SENSAN *) ckt->CKTcurJob;
368 struct sSENSint *st = &job->ST;
369 sgen *sg;
370 IFvalue value, nvalue;
371 double *E, *iE;
372 SMPmatrix *Y;
373 double delta_var;
374 int i, j;
375 int bypass;
376 int error;
377 bypass = ckt->CKTbypass;
378
379 ckt->CKTbypass = 0;
380
381 /* The unknown vector of node voltages overwrites rhs */
382 E = ckt->CKTrhs;
383 iE = ckt->CKTirhs;
384 ckt->CKTrhsOld = E;
385 ckt->CKTirhsOld = iE;
386 Y = ckt->CKTmatrix;
387
388 for (i = 0; i <= st->size; i++) {
389 st->dIr[i] = 0.0;
390 st->dIi[i] = 0.0;
391 }
392
393 ckt->CKTrhs = E;
394 ckt->CKTirhs = iE;
395 ckt->CKTmatrix = Y;
396
397 /* This generates Y in LU form */
398
399 /* Yes, all this has to be re-done */
400 ckt->CKTpreload = 0;
401 error = CKTsetup(ckt);
402 if (error)
403 return (error);
404 error = CKTtemp(ckt);
405 if (error)
406 return (error);
407 ckt->CKTmode = MODEDCOP | MODEINITSMSIG;
408 error = CKTload(ckt);
409 if (error)
410 return (error);
411 ckt->CKTmode = MODEAC;
412 error = NIacIter(ckt);
413 if (error)
414 return (error);
415
416 /* Use a different vector & matrix */
417 ckt->CKTrhs = st->dIr;
418 ckt->CKTirhs = st->dIi;
419 ckt->CKTmatrix = st->dY;
420
421 /* calc. effect of each param */
422 for (i = 0,sg = sgen_init(ckt, FALSE); sg; i++,sgen_next(&sg)) {
423
424 /* clear CKTmatrix, CKTrhs */
425 spSetComplex(st->dY);
426 spClear(st->dY);
427 for (j = 0; j <= st->size; j++) {
428 st->dIr[j] = 0.0;
429 st->dIi[j] = 0.0;
430 }
431
432 error = sens_loadnew(sg,ckt,FALSE);
433 if (error)
434 return (error);
435
436 /* Alter the parameter */
437
438 if (sg->value != 0.0)
439 delta_var = sg->value * Sens_Delta;
440 else
441 delta_var = Sens_Abs_Delta;
442
443 nvalue.rValue = sg->value + delta_var;
444
445 sens_setp(sg, ckt, &nvalue);
446 if (error)
447 return (error);
448
449 /* change sign of CKTmatrix, CKTrhs */
450 spConstMult(st->dY, -1.0);
451 for (j = 0; j <= st->size; j++) {
452 st->dIr[j] = -st->dIr[j];
453 st->dIi[j] = -st->dIi[j];
454 }
455
456 error = sens_loadnew(sg,ckt,FALSE);
457 if (error)
458 return (error);
459
460 /* now have delta_Y = CKTmatrix, delta_I = CKTrhs */
461
462 /* reset parameter */
463 value.rValue = sg->value;
464 sens_setp(sg, ckt, &value);
465
466 /* delta_Y E */
467 spMultiply(st->dY, st->dIdYr, E, st->dIdYi, iE);
468
469 /* delta_I - delta_Y E */
470 for (j = 0; j <= st->size; j++) {
471 st->dIr[j] -= st->dIdYr[j];
472 st->dIi[j] -= st->dIdYi[j];
473 }
474
475 /* Solve; Y already factored */
476 spSolve(Y, st->dIr, st->dIr, st->dIi, st->dIi);
477
478 /* delta_I is now equal to delta_E */
479
480 st->dIr[0] = 0.0;
481 st->dIi[0] = 0.0;
482 if (job->SENSoutName) {
483 st->o_cvalues[i].real =
484 st->dIr[job->SENSoutPos->number] -
485 st->dIr[job->SENSoutNeg->number];
486 st->o_cvalues[i].imag =
487 st->dIi[job->SENSoutPos->number] -
488 st->dIi[job->SENSoutNeg->number];
489 }
490 else {
491 st->o_cvalues[i].real =
492 st->dIr
493 [((SRCinstance*)job->SENSoutSrcDev)->SRCbranch];
494 st->o_cvalues[i].imag =
495 st->dIi
496 [((SRCinstance*)job->SENSoutSrcDev)->SRCbranch];
497 }
498 st->o_cvalues[i].real /= delta_var;
499 st->o_cvalues[i].imag /= delta_var;
500 }
501 nvalue.v.vec.cVec = st->o_cvalues;
502
503 value.rValue = ckt->CKTomega/(2*M_PI);
504 OUTdata(job->SENSplot, &value, &nvalue);
505 outd->count++;
506
507 ckt->CKTrhs = E;
508 ckt->CKTirhs = iE;
509 ckt->CKTmatrix = Y;
510 ckt->CKTbypass = bypass;
511
512 return (OK);
513 }
514
515
516 /* Get parameter value */
517 int
sens_getp(sg,ckt,val)518 sens_getp(sg, ckt, val)
519
520 sgen *sg;
521 CKTcircuit *ckt;
522 IFvalue *val;
523 {
524 int (*fn)( );
525 int pid;
526 int error = 0;
527
528 if (sg->is_instparam) {
529 fn = DEVices[sg->dev]->DEVask;
530 pid = DEVices[sg->dev]->DEVpublic.instanceParms[sg->param].id;
531 if (fn)
532 error = (*fn)(ckt, sg->instance, pid, val, NULL);
533 else
534 return (1);
535 }
536 else {
537 fn = DEVices[sg->dev]->DEVmodAsk;
538 pid = DEVices[sg->dev]->DEVpublic.modelParms[sg->param].id;
539 if (fn)
540 error = (*fn)(ckt, sg->model, pid, val);
541 else
542 return (1);
543 }
544
545 if (error) {
546 if (sg->is_instparam)
547 printf("GET ERROR: %s:%s:%s -> param %s (%d)\n",
548 DEVices[sg->dev]->DEVpublic.name,
549 sg->model->GENmodName,
550 sg->instance->GENname,
551 sg->ptable[sg->param].keyword, pid);
552 else
553 printf("GET ERROR: %s:%s:%s -> mparam %s (%d)\n",
554 DEVices[sg->dev]->DEVpublic.name,
555 sg->model->GENmodName,
556 sg->instance->GENname,
557 sg->ptable[sg->param].keyword, pid);
558 }
559
560 return (error);
561 }
562
563
564 /* Set parameter value */
565 static int
sens_setp(sg,ckt,val)566 sens_setp(sg, ckt, val)
567
568 sgen *sg;
569 CKTcircuit *ckt;
570 IFvalue *val;
571 {
572 int (*fn)( );
573 int pid;
574 int error = 0;
575
576 if (sg->is_instparam) {
577 fn = DEVices[sg->dev]->DEVparam;
578 pid = DEVices[sg->dev]->DEVpublic.instanceParms[sg->param].id;
579 if (fn)
580 error = (*fn)(ckt, pid, val, sg->instance, NULL);
581 else
582 return (1);
583 }
584 else {
585 fn = DEVices[sg->dev]->DEVmodParam;
586 pid = DEVices[sg->dev]->DEVpublic.modelParms[sg->param].id;
587 if (fn)
588 error = (*fn)(pid, val, sg->model);
589 else
590 return (1);
591 }
592
593 if (error) {
594 if (sg->is_instparam)
595 printf("SET ERROR: %s:%s:%s -> param %s (%d)\n",
596 DEVices[sg->dev]->DEVpublic.name,
597 sg->model->GENmodName,
598 sg->instance->GENname,
599 sg->ptable[sg->param].keyword, pid);
600 else
601 printf("SET ERROR: %s:%s:%s -> mparam %s (%d)\n",
602 DEVices[sg->dev]->DEVpublic.name,
603 sg->model->GENmodName,
604 sg->instance->GENname,
605 sg->ptable[sg->param].keyword, pid);
606 }
607 return (error);
608 }
609
610
611 static int
sens_loadnew(sg,ckt,is_dc)612 sens_loadnew(sg, ckt, is_dc)
613
614 sgen *sg;
615 CKTcircuit *ckt;
616 int is_dc;
617 {
618 int error = 0;
619 int (*fn)();
620
621 if (!is_dc)
622 ckt->CKTpreload = 0;
623 else
624 ckt->CKTpreload = 1;
625
626 /* call setup */
627 ckt->CKTnumStates = sg->istate;
628 fn = DEVices[sg->dev]->DEVsetup;
629 if (fn)
630 (*fn)(ckt->CKTmatrix, sg->model, ckt, &ckt->CKTnumStates);
631
632 /* call temp */
633 fn = DEVices[sg->dev]->DEVtemperature;
634 if (fn)
635 (*fn)(sg->model, ckt);
636
637 /* call load */
638 if (!is_dc) {
639 fn = DEVices[sg->dev]->DEVacLoad;
640 }
641 else {
642 fn = DEVices[sg->dev]->DEVload;
643 }
644 if (fn)
645 error = (*fn)(sg->model, ckt);
646
647 return (error);
648 }
649