1 /**********
2 Copyright 1990 Regents of the University of California. All rights reserved.
3 Author: 1985 Thomas L. Quarles
4 Modified: 2001 AlansFixes
5 **********/
6
7 /*
8 * NIiter(ckt,maxIter)
9 *
10 * This subroutine performs the actual numerical iteration.
11 * It uses the sparse matrix stored in the circuit struct
12 * along with the matrix loading program, the load data, the
13 * convergence test function, and the convergence parameters
14 */
15
16 #include "ngspice/ngspice.h"
17 #include "ngspice/trandefs.h"
18 #include "ngspice/cktdefs.h"
19 #include "ngspice/smpdefs.h"
20 #include "ngspice/sperror.h"
21
22
23 /* NIiter() - return value is non-zero for convergence failure */
24
25 int
NIiter(CKTcircuit * ckt,int maxIter)26 NIiter(CKTcircuit *ckt, int maxIter)
27 {
28 double startTime, *OldCKTstate0 = NULL;
29 int error, i, j;
30
31 int iterno = 0;
32 int ipass = 0;
33
34 /* some convergence issues that get resolved by increasing max iter */
35 if (maxIter < 100)
36 maxIter = 100;
37
38 if ((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC)) {
39 SWAP(double *, ckt->CKTrhs, ckt->CKTrhsOld);
40 error = CKTload(ckt);
41 if (error)
42 return(error);
43 return(OK);
44 }
45
46 #ifdef WANT_SENSE2
47 if (ckt->CKTsenInfo) {
48 error = NIsenReinit(ckt);
49 if (error)
50 return(error);
51 }
52 #endif
53
54 if (ckt->CKTniState & NIUNINITIALIZED) {
55 error = NIreinit(ckt);
56 if (error) {
57 #ifdef STEPDEBUG
58 printf("re-init returned error \n");
59 #endif
60 return(error);
61 }
62 }
63
64 /* OldCKTstate0 = TMALLOC(double, ckt->CKTnumStates + 1); */
65
66 for (;;) {
67
68 ckt->CKTnoncon = 0;
69
70 #ifdef NEWPRED
71 if (!(ckt->CKTmode & MODEINITPRED))
72 #endif
73 {
74
75 error = CKTload(ckt);
76 /* printf("loaded, noncon is %d\n", ckt->CKTnoncon); */
77 /* fflush(stdout); */
78 iterno++;
79 if (error) {
80 ckt->CKTstat->STATnumIter += iterno;
81 #ifdef STEPDEBUG
82 printf("load returned error \n");
83 #endif
84 FREE(OldCKTstate0);
85 return (error);
86 }
87
88 /* printf("after loading, before solving\n"); */
89 /* CKTdump(ckt); */
90
91 if (!(ckt->CKTniState & NIDIDPREORDER)) {
92 error = SMPpreOrder(ckt->CKTmatrix);
93 if (error) {
94 ckt->CKTstat->STATnumIter += iterno;
95 #ifdef STEPDEBUG
96 printf("pre-order returned error \n");
97 #endif
98 FREE(OldCKTstate0);
99 return(error); /* badly formed matrix */
100 }
101 ckt->CKTniState |= NIDIDPREORDER;
102 }
103
104 if ((ckt->CKTmode & MODEINITJCT) ||
105 ((ckt->CKTmode & MODEINITTRAN) && (iterno == 1)))
106 {
107 ckt->CKTniState |= NISHOULDREORDER;
108 }
109
110 if (ckt->CKTniState & NISHOULDREORDER) {
111 startTime = SPfrontEnd->IFseconds();
112 error = SMPreorder(ckt->CKTmatrix, ckt->CKTpivotAbsTol,
113 ckt->CKTpivotRelTol, ckt->CKTdiagGmin);
114 ckt->CKTstat->STATreorderTime +=
115 SPfrontEnd->IFseconds() - startTime;
116 if (error) {
117 /* new feature - we can now find out something about what is
118 * wrong - so we ask for the troublesome entry
119 */
120 SMPgetError(ckt->CKTmatrix, &i, &j);
121 SPfrontEnd->IFerrorf (ERR_WARNING, "singular matrix: check nodes %s and %s\n", NODENAME(ckt, i), NODENAME(ckt, j));
122 ckt->CKTstat->STATnumIter += iterno;
123 #ifdef STEPDEBUG
124 printf("reorder returned error \n");
125 #endif
126 FREE(OldCKTstate0);
127 return(error); /* can't handle these errors - pass up! */
128 }
129 ckt->CKTniState &= ~NISHOULDREORDER;
130 } else {
131 startTime = SPfrontEnd->IFseconds();
132 error = SMPluFac(ckt->CKTmatrix, ckt->CKTpivotAbsTol,
133 ckt->CKTdiagGmin);
134 ckt->CKTstat->STATdecompTime +=
135 SPfrontEnd->IFseconds() - startTime;
136 if (error) {
137 if (error == E_SINGULAR) {
138 ckt->CKTniState |= NISHOULDREORDER;
139 DEBUGMSG(" forced reordering....\n");
140 continue;
141 }
142 /* CKTload(ckt); */
143 /* SMPprint(ckt->CKTmatrix, stdout); */
144 /* seems to be singular - pass the bad news up */
145 ckt->CKTstat->STATnumIter += iterno;
146 #ifdef STEPDEBUG
147 printf("lufac returned error \n");
148 #endif
149 FREE(OldCKTstate0);
150 return(error);
151 }
152 }
153
154 /* moved it to here as if xspice is included then CKTload changes
155 CKTnumStates the first time it is run */
156 if (!OldCKTstate0)
157 OldCKTstate0 = TMALLOC(double, ckt->CKTnumStates + 1);
158 memcpy(OldCKTstate0, ckt->CKTstate0,
159 (size_t) ckt->CKTnumStates * sizeof(double));
160
161 startTime = SPfrontEnd->IFseconds();
162 SMPsolve(ckt->CKTmatrix, ckt->CKTrhs, ckt->CKTrhsSpare);
163 ckt->CKTstat->STATsolveTime +=
164 SPfrontEnd->IFseconds() - startTime;
165 #ifdef STEPDEBUG
166 /*XXXX*/
167 if (ckt->CKTrhs[0] != 0.0)
168 printf("NIiter: CKTrhs[0] = %g\n", ckt->CKTrhs[0]);
169 if (ckt->CKTrhsSpare[0] != 0.0)
170 printf("NIiter: CKTrhsSpare[0] = %g\n", ckt->CKTrhsSpare[0]);
171 if (ckt->CKTrhsOld[0] != 0.0)
172 printf("NIiter: CKTrhsOld[0] = %g\n", ckt->CKTrhsOld[0]);
173 /*XXXX*/
174 #endif
175 ckt->CKTrhs[0] = 0;
176 ckt->CKTrhsSpare[0] = 0;
177 ckt->CKTrhsOld[0] = 0;
178
179 if (iterno > maxIter) {
180 ckt->CKTstat->STATnumIter += iterno;
181 /* we don't use this info during transient analysis */
182 if (ckt->CKTcurrentAnalysis != DOING_TRAN) {
183 FREE(errMsg);
184 errMsg = copy("Too many iterations without convergence");
185 #ifdef STEPDEBUG
186 fprintf(stderr, "too many iterations without convergence: %d iter's (max iter == %d)\n",
187 iterno, maxIter);
188 #endif
189 }
190 FREE(OldCKTstate0);
191 return(E_ITERLIM);
192 }
193
194 if ((ckt->CKTnoncon == 0) && (iterno != 1))
195 ckt->CKTnoncon = NIconvTest(ckt);
196 else
197 ckt->CKTnoncon = 1;
198
199 #ifdef STEPDEBUG
200 printf("noncon is %d\n", ckt->CKTnoncon);
201 #endif
202 }
203
204 if ((ckt->CKTnodeDamping != 0) && (ckt->CKTnoncon != 0) &&
205 ((ckt->CKTmode & MODETRANOP) || (ckt->CKTmode & MODEDCOP)) &&
206 (iterno > 1))
207 {
208 CKTnode *node;
209 double diff, maxdiff = 0;
210 for (node = ckt->CKTnodes->next; node; node = node->next)
211 if (node->type == SP_VOLTAGE) {
212 diff = fabs(ckt->CKTrhs[node->number] - ckt->CKTrhsOld[node->number]);
213 if (maxdiff < diff)
214 maxdiff = diff;
215 }
216
217 if (maxdiff > 10) {
218 double damp_factor = 10 / maxdiff;
219 if (damp_factor < 0.1)
220 damp_factor = 0.1;
221 for (node = ckt->CKTnodes->next; node; node = node->next) {
222 diff = ckt->CKTrhs[node->number] - ckt->CKTrhsOld[node->number];
223 ckt->CKTrhs[node->number] =
224 ckt->CKTrhsOld[node->number] + (damp_factor * diff);
225 }
226 for (i = 0; i < ckt->CKTnumStates; i++) {
227 diff = ckt->CKTstate0[i] - OldCKTstate0[i];
228 ckt->CKTstate0[i] = OldCKTstate0[i] + (damp_factor * diff);
229 }
230 }
231 }
232
233 if (ckt->CKTmode & MODEINITFLOAT) {
234 if ((ckt->CKTmode & MODEDC) && ckt->CKThadNodeset) {
235 if (ipass)
236 ckt->CKTnoncon = ipass;
237 ipass = 0;
238 }
239 if (ckt->CKTnoncon == 0) {
240 ckt->CKTstat->STATnumIter += iterno;
241 FREE(OldCKTstate0);
242 return(OK);
243 }
244 } else if (ckt->CKTmode & MODEINITJCT) {
245 ckt->CKTmode = (ckt->CKTmode & ~INITF) | MODEINITFIX;
246 ckt->CKTniState |= NISHOULDREORDER;
247 } else if (ckt->CKTmode & MODEINITFIX) {
248 if (ckt->CKTnoncon == 0)
249 ckt->CKTmode = (ckt->CKTmode & ~INITF) | MODEINITFLOAT;
250 ipass = 1;
251 } else if (ckt->CKTmode & MODEINITSMSIG) {
252 ckt->CKTmode = (ckt->CKTmode & ~INITF) | MODEINITFLOAT;
253 } else if (ckt->CKTmode & MODEINITTRAN) {
254 if (iterno <= 1)
255 ckt->CKTniState |= NISHOULDREORDER;
256 ckt->CKTmode = (ckt->CKTmode & ~INITF) | MODEINITFLOAT;
257 } else if (ckt->CKTmode & MODEINITPRED) {
258 ckt->CKTmode = (ckt->CKTmode & ~INITF) | MODEINITFLOAT;
259 } else {
260 ckt->CKTstat->STATnumIter += iterno;
261 #ifdef STEPDEBUG
262 printf("bad initf state \n");
263 #endif
264 FREE(OldCKTstate0);
265 return(E_INTERN);
266 /* impossible - no such INITF flag! */
267 }
268
269 /* build up the lvnim1 array from the lvn array */
270 SWAP(double *, ckt->CKTrhs, ckt->CKTrhsOld);
271 /* printf("after loading, after solving\n"); */
272 /* CKTdump(ckt); */
273 }
274 /*NOTREACHED*/
275 }
276