Branch data Line data Source code
1 : : /*
2 : : * nasolver.cpp - nodal analysis solver class implementation
3 : : *
4 : : * Copyright (C) 2004, 2005, 2006, 2007, 2008 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : // the types required for qucs library files are defined
26 : : // in qucs_typedefs.h, created by configure from
27 : : // qucs_typedefs.h.in
28 : : #include "qucs_typedefs.h"
29 : :
30 : : #include <stdio.h>
31 : : #include <stdlib.h>
32 : : #include <string.h>
33 : : #include <cmath>
34 : : #include <float.h>
35 : : #include <assert.h>
36 : : #include <limits>
37 : :
38 : : #include "logging.h"
39 : : #include "complex.h"
40 : : #include "object.h"
41 : : #include "node.h"
42 : : #include "circuit.h"
43 : : #include "vector.h"
44 : : #include "dataset.h"
45 : : #include "net.h"
46 : : #include "analysis.h"
47 : : #include "nodelist.h"
48 : : #include "nodeset.h"
49 : : #include "strlist.h"
50 : : #include "tvector.h"
51 : : #include "tmatrix.h"
52 : : #include "eqnsys.h"
53 : : #include "precision.h"
54 : : #include "operatingpoint.h"
55 : : #include "exception.h"
56 : : #include "exceptionstack.h"
57 : : #include "nasolver.h"
58 : : #include "constants.h"
59 : :
60 : : namespace qucs {
61 : :
62 : : // Constructor creates an unnamed instance of the nasolver class.
63 : : template <class nr_type_t>
64 [ + - ]: 103 : nasolver<nr_type_t>::nasolver () : analysis ()
65 : : {
66 : 103 : nlist = NULL;
67 : 103 : A = C = NULL;
68 : 103 : z = x = xprev = zprev = NULL;
69 : 103 : reltol = abstol = vntol = 0;
70 : 103 : desc = NULL;
71 : 103 : calculate_func = NULL;
72 : 103 : convHelper = fixpoint = 0;
73 : 103 : eqnAlgo = ALGO_LU_DECOMPOSITION;
74 : 103 : updateMatrix = 1;
75 : 103 : gMin = srcFactor = 0;
76 [ + - ]: 103 : eqns = new eqnsys<nr_type_t> ();
77 : 103 : }
78 : :
79 : : // Constructor creates a named instance of the nasolver class.
80 : : template <class nr_type_t>
81 [ # # ]: 0 : nasolver<nr_type_t>::nasolver (char * n) : analysis (n)
82 : : {
83 : 0 : nlist = NULL;
84 : 0 : A = C = NULL;
85 : 0 : z = x = xprev = zprev = NULL;
86 : 0 : reltol = abstol = vntol = 0;
87 : 0 : desc = NULL;
88 : 0 : calculate_func = NULL;
89 : 0 : convHelper = fixpoint = 0;
90 : 0 : eqnAlgo = ALGO_LU_DECOMPOSITION;
91 : 0 : updateMatrix = 1;
92 : 0 : gMin = srcFactor = 0;
93 [ # # ]: 0 : eqns = new eqnsys<nr_type_t> ();
94 : 0 : }
95 : :
96 : : // Destructor deletes the nasolver class object.
97 : : template <class nr_type_t>
98 : 103 : nasolver<nr_type_t>::~nasolver ()
99 : : {
100 [ # # ][ # # ]: 103 : if (nlist) delete nlist;
[ # # ][ - + ]
[ # # ][ # # ]
101 [ # # ][ # # ]: 103 : if (C) delete C;
[ + + ][ + - ]
102 [ # # ][ + - ]: 103 : delete A;
103 [ # # ][ + - ]: 103 : delete z;
104 [ # # ][ + - ]: 103 : delete x;
105 [ # # ][ + + ]: 103 : delete xprev;
106 [ # # ][ + + ]: 103 : delete zprev;
107 [ # # ][ # # ]: 206 : delete eqns;
[ + - ][ - + ]
108 [ # # ][ + - ]: 206 : }
109 : :
110 : : /* The copy constructor creates a new instance of the nasolver class
111 : : based on the given nasolver object. */
112 : : template <class nr_type_t>
113 [ # # ]: 0 : nasolver<nr_type_t>::nasolver (nasolver & o) : analysis (o)
114 : : {
115 [ # # ][ # # ]: 0 : nlist = o.nlist ? new nodelist (*(o.nlist)) : NULL;
[ # # ]
116 [ # # ][ # # ]: 0 : A = o.A ? new tmatrix<nr_type_t> (*(o.A)) : NULL;
[ # # ]
117 [ # # ][ # # ]: 0 : C = o.C ? new tmatrix<nr_type_t> (*(o.C)) : NULL;
[ # # ]
118 [ # # ][ # # ]: 0 : z = o.z ? new tvector<nr_type_t> (*(o.z)) : NULL;
[ # # ]
119 [ # # ][ # # ]: 0 : x = o.x ? new tvector<nr_type_t> (*(o.x)) : NULL;
[ # # ]
120 : 0 : xprev = zprev = NULL;
121 : 0 : reltol = o.reltol;
122 : 0 : abstol = o.abstol;
123 : 0 : vntol = o.vntol;
124 : 0 : desc = o.desc;
125 : 0 : calculate_func = o.calculate_func;
126 : 0 : convHelper = o.convHelper;
127 : 0 : eqnAlgo = o.eqnAlgo;
128 : 0 : updateMatrix = o.updateMatrix;
129 : 0 : fixpoint = o.fixpoint;
130 : 0 : gMin = o.gMin;
131 : 0 : srcFactor = o.srcFactor;
132 [ # # ][ # # ]: 0 : eqns = new eqnsys<nr_type_t> (*(o.eqns));
133 [ # # ][ # # ]: 0 : solution = nasolution<nr_type_t> (o.solution);
134 : 0 : }
135 : :
136 : : /* The function runs the nodal analysis solver once, reports errors if
137 : : any and save the results into each circuit. */
138 : : template <class nr_type_t>
139 : 658511 : int nasolver<nr_type_t>::solve_once (void)
140 : : {
141 : : qucs::exception * e;
142 : 658511 : int error = 0, d;
143 : :
144 : : // run the calculation function for each circuit
145 : 658511 : calculate ();
146 : :
147 : : // generate A matrix and z vector
148 : 658511 : createMatrix ();
149 : :
150 : : // solve equation system
151 : : try_running ()
152 : : {
153 : 658511 : runMNA ();
154 : : }
155 : : // appropriate exception handling
156 [ + + ]: 658511 : catch_exception ()
[ - + - ]
157 : : {
158 : : case EXCEPTION_PIVOT:
159 : : case EXCEPTION_WRONG_VOLTAGE:
160 [ # # ]: 0 : e = new qucs::exception (EXCEPTION_NA_FAILED);
161 : 0 : d = top_exception()->getData ();
162 : 0 : pop_exception ();
163 [ # # ]: 0 : if (d >= countNodes ())
164 : : {
165 : 0 : d -= countNodes ();
166 : 0 : e->setText ("voltage source `%s' conflicts with some other voltage "
167 : : "source", findVoltageSource(d)->getName ());
168 : : }
169 : : else
170 : : {
171 : 0 : e->setText ("circuit admittance matrix in %s solver is singular at "
172 : : "node `%s' connected to [%s]", desc, nlist->get (d),
173 : : nlist->getNodeString (d));
174 : : }
175 : 0 : throw_exception (e);
176 : 0 : error++;
177 : 0 : break;
178 : : case EXCEPTION_SINGULAR:
179 [ - + # # ]: 17 : do
[ - + ]
180 : : {
181 : 17 : d = top_exception()->getData ();
182 : 17 : pop_exception ();
183 [ - + ]: 17 : if (d < countNodes ())
184 : : {
185 : 17 : logprint (LOG_ERROR, "WARNING: %s: inserted virtual resistance at "
186 : : "node `%s' connected to [%s]\n", getName (), nlist->get (d),
187 : : nlist->getNodeString (d));
188 : : }
189 : : }
190 : : while (top_exception() != NULL &&
191 : : top_exception()->getCode () == EXCEPTION_SINGULAR);
192 : 17 : break;
193 : : default:
194 : 0 : estack.print ();
195 : 17 : break;
196 : : }
197 : :
198 : : // save results into circuits
199 [ + - ]: 658511 : if (!error) saveSolution ();
200 : 658511 : return error;
201 : : }
202 : :
203 : : /* Run this function after the actual solver run and before evaluating
204 : : the results. */
205 : : template <class nr_type_t>
206 : 22664 : void nasolver<nr_type_t>::solve_post (void)
207 : : {
208 [ + - ]: 22664 : delete nlist;
209 : 22664 : nlist = NULL;
210 : 22664 : }
211 : :
212 : : /* Run this function before the actual solver. */
213 : : template <class nr_type_t>
214 : 22664 : void nasolver<nr_type_t>::solve_pre (void)
215 : : {
216 : : // create node list, enumerate nodes and voltage sources
217 : : #if DEBUG
218 : 22664 : logprint (LOG_STATUS, "NOTIFY: %s: creating node list for %s analysis\n",
219 : : getName (), desc);
220 : : #endif
221 [ + - ]: 22664 : nlist = new nodelist (subnet);
222 : 22664 : nlist->assignNodes ();
223 : 22664 : assignVoltageSources ();
224 : : #if DEBUG && 0
225 : : nlist->print ();
226 : : #endif
227 : :
228 : : // create matrix, solution vector and right hand side vector
229 : 22664 : int M = countVoltageSources ();
230 : 22664 : int N = countNodes ();
231 [ + + ][ + - ]: 22664 : if (A != NULL) delete A;
232 [ + - ]: 22664 : A = new tmatrix<nr_type_t> (M + N);
233 [ + + ][ + - ]: 22664 : if (z != NULL) delete z;
234 [ + - ]: 22664 : z = new tvector<nr_type_t> (N + M);
235 [ + + ][ + - ]: 22664 : if (x != NULL) delete x;
236 [ + - ]: 22664 : x = new tvector<nr_type_t> (N + M);
237 : :
238 : : #if DEBUG
239 : 22664 : logprint (LOG_STATUS, "NOTIFY: %s: solving %s netlist\n", getName (), desc);
240 : : #endif
241 : 22664 : }
242 : :
243 : : /* This function goes through the nodeset list of the current netlist
244 : : and applies the stored values to the current solution vector. Then
245 : : the function saves the solution vector back into the actual
246 : : component nodes. */
247 : : template <class nr_type_t>
248 : 20880 : void nasolver<nr_type_t>::applyNodeset (bool nokeep)
249 : : {
250 [ + - ][ - + ]: 41760 : if (x == NULL || nlist == NULL) return;
251 : :
252 : : // set each solution to zero
253 [ + + ][ + + ]: 133584 : if (nokeep) for (int i = 0; i < x->getSize (); i++) x->set (i, 0);
254 : :
255 : : // then apply the nodeset itself
256 [ + + ]: 20882 : for (nodeset * n = subnet->getNodeset (); n; n = n->getNext ())
257 : : {
258 : 2 : struct nodelist_t * nl = nlist->getNode (n->getName ());
259 [ + - ]: 2 : if (nl != NULL)
260 : : {
261 : 2 : x->set (nl->n, n->getValue ());
262 : : }
263 : : else
264 : : {
265 : 0 : logprint (LOG_ERROR, "WARNING: %s: no such node `%s' found, cannot "
266 : : "initialize node\n", getName (), n->getName ());
267 : : }
268 : : }
269 [ + + ]: 20880 : if (xprev != NULL) *xprev = *x;
270 : 20880 : saveSolution ();
271 : : }
272 : :
273 : : /* The following function uses the gMin-stepping algorithm in order to
274 : : solve the given non-linear netlist by continuous iterations. */
275 : : template <class nr_type_t>
276 : 0 : int nasolver<nr_type_t>::solve_nonlinear_continuation_gMin (void)
277 : : {
278 : : qucs::exception * e;
279 : 0 : int convergence, run = 0, MaxIterations, error = 0;
280 : : nr_double_t gStep, gPrev;
281 : :
282 : : // fetch simulation properties
283 : 0 : MaxIterations = getPropertyInteger ("MaxIter") / 4 + 1;
284 : 0 : updateMatrix = 1;
285 : 0 : fixpoint = 0;
286 : :
287 : : // initialize the stepper
288 : 0 : gPrev = gMin = 0.01;
289 : 0 : gStep = gMin / 100;
290 : 0 : gMin -= gStep;
291 : :
292 [ # # ]: 0 : do
293 : : {
294 : : // run solving loop until convergence is reached
295 : 0 : run = 0;
296 [ # # ][ # # ]: 0 : do
[ # # ]
297 : : {
298 : 0 : error = solve_once ();
299 [ # # ]: 0 : if (!error)
300 : : {
301 : : // convergence check
302 [ # # ]: 0 : convergence = (run > 0) ? checkConvergence () : 0;
303 : 0 : savePreviousIteration ();
304 : 0 : run++;
305 : : }
306 : 0 : else break;
307 : : }
308 : : while (!convergence && run < MaxIterations);
309 : 0 : iterations += run;
310 : :
311 : : // not yet converged, so decreased the gMin-step
312 [ # # ][ # # ]: 0 : if (run >= MaxIterations || error)
313 : : {
314 : 0 : gStep /= 2;
315 : : // here the absolute minimum step checker
316 [ # # ]: 0 : if (gStep < std::numeric_limits<nr_double_t>::epsilon())
317 : : {
318 : 0 : error = 1;
319 [ # # ]: 0 : e = new qucs::exception (EXCEPTION_NO_CONVERGENCE);
320 : 0 : e->setText ("no convergence in %s analysis after %d gMinStepping "
321 : : "iterations", desc, iterations);
322 : 0 : throw_exception (e);
323 : 0 : break;
324 : : }
325 [ # # ]: 0 : gMin = MAX (gPrev - gStep, 0);
326 : : }
327 : : // converged, increased the gMin-step
328 : : else
329 : : {
330 : 0 : gPrev = gMin;
331 [ # # ]: 0 : gMin = MAX (gMin - gStep, 0);
332 : 0 : gStep *= 2;
333 : : }
334 : : }
335 : : // continue until no additional resistances is necessary
336 : : while (gPrev > 0);
337 : :
338 : 0 : return error;
339 : : }
340 : :
341 : : /* The following function uses the source-stepping algorithm in order
342 : : to solve the given non-linear netlist by continuous iterations. */
343 : : template <class nr_type_t>
344 : 0 : int nasolver<nr_type_t>::solve_nonlinear_continuation_Source (void)
345 : : {
346 : : qucs::exception * e;
347 : 0 : int convergence, run = 0, MaxIterations, error = 0;
348 : : nr_double_t sStep, sPrev;
349 : :
350 : : // fetch simulation properties
351 : 0 : MaxIterations = getPropertyInteger ("MaxIter") / 4 + 1;
352 : 0 : updateMatrix = 1;
353 : 0 : fixpoint = 0;
354 : :
355 : : // initialize the stepper
356 : 0 : sPrev = srcFactor = 0;
357 : 0 : sStep = 0.01;
358 : 0 : srcFactor += sStep;
359 : :
360 [ # # ]: 0 : do
361 : : {
362 : : // run solving loop until convergence is reached
363 : 0 : run = 0;
364 [ # # ][ # # ]: 0 : do
[ # # ]
365 : : {
366 : 0 : subnet->setSrcFactor (srcFactor);
367 : 0 : error = solve_once ();
368 [ # # ]: 0 : if (!error)
369 : : {
370 : : // convergence check
371 [ # # ]: 0 : convergence = (run > 0) ? checkConvergence () : 0;
372 : 0 : savePreviousIteration ();
373 : 0 : run++;
374 : : }
375 : 0 : else break;
376 : : }
377 : : while (!convergence && run < MaxIterations);
378 : 0 : iterations += run;
379 : :
380 : : // not yet converged, so decreased the source-step
381 [ # # ][ # # ]: 0 : if (run >= MaxIterations || error)
382 : : {
383 [ # # ]: 0 : if (error)
384 : 0 : sStep *= 0.1;
385 : : else
386 : 0 : sStep *= 0.5;
387 : 0 : restorePreviousIteration ();
388 : 0 : saveSolution ();
389 : : // here the absolute minimum step checker
390 [ # # ]: 0 : if (sStep < std::numeric_limits<nr_double_t>::epsilon())
391 : : {
392 : 0 : error = 1;
393 [ # # ]: 0 : e = new qucs::exception (EXCEPTION_NO_CONVERGENCE);
394 : 0 : e->setText ("no convergence in %s analysis after %d sourceStepping "
395 : : "iterations", desc, iterations);
396 : 0 : throw_exception (e);
397 : 0 : break;
398 : : }
399 [ # # ]: 0 : srcFactor = MIN (sPrev + sStep, 1);
400 : : }
401 : : // converged, increased the source-step
402 [ # # ]: 0 : else if (run < MaxIterations / 4)
403 : : {
404 : 0 : sPrev = srcFactor;
405 [ # # ]: 0 : srcFactor = MIN (srcFactor + sStep, 1);
406 : 0 : sStep *= 1.5;
407 : : }
408 : : else
409 : : {
410 [ # # ]: 0 : srcFactor = MIN (srcFactor + sStep, 1);
411 : : }
412 : : }
413 : : // continue until no source factor is necessary
414 : : while (sPrev < 1);
415 : :
416 : 0 : subnet->setSrcFactor (1);
417 : 0 : return error;
418 : : }
419 : :
420 : : /* The function returns an appropriate text representation for the
421 : : currently used convergence helper algorithm. */
422 : : template <class nr_type_t>
423 : 0 : const char * nasolver<nr_type_t>::getHelperDescription (void)
424 : : {
425 [ # # ]: 0 : if (convHelper == CONV_Attenuation)
426 : : {
427 : 0 : return "RHS attenuation";
428 : : }
429 [ # # ]: 0 : else if (convHelper == CONV_LineSearch)
430 : : {
431 : 0 : return "line search";
432 : : }
433 [ # # ]: 0 : else if (convHelper == CONV_SteepestDescent)
434 : : {
435 : 0 : return "steepest descent";
436 : : }
437 [ # # ]: 0 : else if (convHelper == CONV_GMinStepping)
438 : : {
439 : 0 : return "gMin stepping";
440 : : }
441 [ # # ]: 0 : else if (convHelper == CONV_SourceStepping)
442 : : {
443 : 0 : return "source stepping";
444 : : }
445 : 0 : return "none";
446 : : }
447 : :
448 : : /* This is the non-linear iterative nodal analysis netlist solver. */
449 : : template <class nr_type_t>
450 : 235300 : int nasolver<nr_type_t>::solve_nonlinear (void)
451 : : {
452 : : qucs::exception * e;
453 : 235300 : int convergence, run = 0, MaxIterations, error = 0;
454 : :
455 : : // fetch simulation properties
456 : 235300 : MaxIterations = getPropertyInteger ("MaxIter");
457 : 235300 : reltol = getPropertyDouble ("reltol");
458 : 235300 : abstol = getPropertyDouble ("abstol");
459 : 235300 : vntol = getPropertyDouble ("vntol");
460 : 235300 : updateMatrix = 1;
461 : :
462 [ - + ]: 235300 : if (convHelper == CONV_GMinStepping)
463 : : {
464 : : // use the alternative non-linear solver solve_nonlinear_continuation_gMin
465 : : // instead of the basic solver provided by this function
466 : 0 : iterations = 0;
467 : 0 : error = solve_nonlinear_continuation_gMin ();
468 : 0 : return error;
469 : : }
470 [ - + ]: 235300 : else if (convHelper == CONV_SourceStepping)
471 : : {
472 : : // use the alternative non-linear solver solve_nonlinear_continuation_Source
473 : : // instead of the basic solver provided by this function
474 : 0 : iterations = 0;
475 : 0 : error = solve_nonlinear_continuation_Source ();
476 : 0 : return error;
477 : : }
478 : :
479 : : // run solving loop until convergence is reached
480 [ + + ][ + - ]: 650214 : do
[ + + ][ + + ]
481 : : {
482 : 650214 : error = solve_once ();
483 [ + - ]: 650214 : if (!error)
484 : : {
485 : : // convergence check
486 [ + + ]: 650214 : convergence = (run > 0) ? checkConvergence () : 0;
487 : 650214 : savePreviousIteration ();
488 : 650214 : run++;
489 : : // control fixpoint iterations
490 [ - + ]: 650214 : if (fixpoint)
491 : : {
492 [ # # ][ # # ]: 0 : if (convergence && !updateMatrix)
493 : : {
494 : 0 : updateMatrix = 1;
495 : 0 : convergence = 0;
496 : : }
497 : : else
498 : : {
499 : 0 : updateMatrix = 0;
500 : : }
501 : : }
502 : : }
503 : : else
504 : : {
505 : 0 : break;
506 : : }
507 : : }
508 : : while (!convergence &&
509 : : run < MaxIterations * (1 + convHelper ? 1 : 0));
510 : :
511 [ + + ][ - + ]: 235300 : if (run >= MaxIterations || error)
512 : : {
513 [ + - ]: 65 : e = new qucs::exception (EXCEPTION_NO_CONVERGENCE);
514 : 65 : e->setText ("no convergence in %s analysis after %d iterations",
515 : : desc, run);
516 : 65 : throw_exception (e);
517 : 65 : error++;
518 : : }
519 : :
520 : 235300 : iterations = run;
521 : 235300 : return error;
522 : : }
523 : :
524 : : /* This is the linear nodal analysis netlist solver. */
525 : : template <class nr_type_t>
526 : 8297 : int nasolver<nr_type_t>::solve_linear (void)
527 : : {
528 : 8297 : updateMatrix = 1;
529 : 8297 : return solve_once ();
530 : : }
531 : :
532 : : /* Applying the MNA (Modified Nodal Analysis) to a circuit with
533 : : passive elements and independent current and voltage sources
534 : : results in a matrix equation of the form Ax = z. This function
535 : : generates the A and z matrix. */
536 : : template <class nr_type_t>
537 : 661514 : void nasolver<nr_type_t>::createMatrix (void)
538 : : {
539 : :
540 : : /* Generate the A matrix. The A matrix consists of four (4) minor
541 : : matrices in the form +- -+
542 : : A = | G B |
543 : : | C D |
544 : : +- -+.
545 : : Each of these minor matrices is going to be generated here. */
546 [ + - ]: 661514 : if (updateMatrix)
547 : : {
548 : 661514 : createGMatrix ();
549 : 661514 : createBMatrix ();
550 : 661514 : createCMatrix ();
551 : 661514 : createDMatrix ();
552 : : }
553 : :
554 : : /* Adjust G matrix if requested. */
555 [ - + ]: 661514 : if (convHelper == CONV_GMinStepping)
556 : : {
557 : 0 : int N = countNodes ();
558 : 0 : int M = countVoltageSources ();
559 [ # # ][ # # ]: 0 : for (int n = 0; n < N + M; n++)
560 : : {
561 [ # # ]: 0 : A->set (n, n, A->get (n, n) + gMin);
562 : : }
563 : : }
564 : :
565 : : /* Generate the z Matrix. The z Matrix consists of two (2) minor
566 : : matrices in the form +- -+
567 : : z = | i |
568 : : | e |
569 : : +- -+.
570 : : Each of these minor matrices is going to be generated here. */
571 : 661514 : createZVector ();
572 : 661514 : }
573 : :
574 : : /* This MatVal() functionality is just helper to get the correct
575 : : values from the circuit's matrices. The additional (unused)
576 : : argument is used to differentiate between the two possible
577 : : types. */
578 : : #define MatVal(x) MatValX (x, (nr_type_t *) 0)
579 : :
580 : : template <class nr_type_t>
581 : 708477 : nr_type_t nasolver<nr_type_t>::MatValX (nr_complex_t z, nr_complex_t *)
582 : : {
583 : 708477 : return z;
584 : : }
585 : :
586 : : template <class nr_type_t>
587 : 32994149 : nr_type_t nasolver<nr_type_t>::MatValX (nr_complex_t z, nr_double_t *)
588 : : {
589 : 32994149 : return real (z);
590 : : }
591 : :
592 : : /* The B matrix is an MxN matrix with only 0, 1 and -1 elements. Each
593 : : location in the matrix corresponds to a particular voltage source
594 : : (first dimension) or a node (second dimension). If the positive
595 : : terminal of the ith voltage source is connected to node k, then the
596 : : element (i,k) in the B matrix is a 1. If the negative terminal of
597 : : the ith voltage source is connected to node k, then the element
598 : : (i,k) in the B matrix is a -1. Otherwise, elements of the B matrix
599 : : are zero. */
600 : : template <class nr_type_t>
601 : 661514 : void nasolver<nr_type_t>::createBMatrix (void)
602 : : {
603 [ + - ]: 661514 : int N = countNodes ();
604 : 661514 : int M = countVoltageSources ();
605 : : circuit * vs;
606 : : struct nodelist_t * n;
607 : 9665 : nr_type_t val;
608 : :
609 : : // go through each voltage sources (first dimension)
610 [ + + ]: 2665803 : for (int c = 0; c < M; c++)
611 : : {
612 [ + - ]: 2004289 : vs = findVoltageSource (c);
613 : : // go through each node (second dimension)
614 [ + + ][ + + ]: 14807769 : for (int r = 0; r < N; r++)
615 : : {
616 : 12803480 : val = 0.0;
617 [ + - ]: 12803480 : n = nlist->getNode (r);
618 [ + + ][ + + ]: 47375063 : for (int i = 0; i < n->nNodes; i++)
619 : : {
620 : : // is voltage source connected to node ?
621 [ + + ]: 34571583 : if (n->nodes[i]->getCircuit () == vs)
622 : : {
623 [ + - ]: 3726536 : val += MatVal (vs->getB (n->nodes[i]->getPort (), c));
624 : : }
625 : : }
626 : : // put value into B matrix
627 [ + - ]: 12803480 : A->set (r, c + N, val);
628 : : }
629 : : }
630 : 661514 : }
631 : :
632 : : /* The C matrix is an NxM matrix with only 0, 1 and -1 elements. Each
633 : : location in the matrix corresponds to a particular node (first
634 : : dimension) or a voltage source (first dimension). If the positive
635 : : terminal of the ith voltage source is connected to node k, then the
636 : : element (k,i) in the C matrix is a 1. If the negative terminal of
637 : : the ith voltage source is connected to node k, then the element
638 : : (k,i) in the C matrix is a -1. Otherwise, elements of the C matrix
639 : : are zero. */
640 : : template <class nr_type_t>
641 : 661514 : void nasolver<nr_type_t>::createCMatrix (void)
642 : : {
643 [ + - ]: 661514 : int N = countNodes ();
644 : 661514 : int M = countVoltageSources ();
645 : : circuit * vs;
646 : : struct nodelist_t * n;
647 : 9665 : nr_type_t val;
648 : :
649 : : // go through each voltage sources (second dimension)
650 [ + + ]: 2665803 : for (int r = 0; r < M; r++)
651 : : {
652 [ + - ]: 2004289 : vs = findVoltageSource (r);
653 : : // go through each node (first dimension)
654 [ + + ][ + + ]: 14807769 : for (int c = 0; c < N; c++)
655 : : {
656 : 12803480 : val = 0.0;
657 [ + - ]: 12803480 : n = nlist->getNode (c);
658 [ + + ][ + + ]: 47375063 : for (int i = 0; i < n->nNodes; i++)
659 : : {
660 : : // is voltage source connected to node ?
661 [ + + ]: 34571583 : if (n->nodes[i]->getCircuit () == vs)
662 : : {
663 [ + - ]: 3726536 : val += MatVal (vs->getC (r, n->nodes[i]->getPort ()));
664 : : }
665 : : }
666 : : // put value into C matrix
667 [ + - ]: 12803480 : A->set (r + N, c, val);
668 : : }
669 : : }
670 : 661514 : }
671 : :
672 : : /* The D matrix is an MxM matrix that is composed entirely of zeros.
673 : : It can be non-zero if dependent sources are considered. */
674 : : template <class nr_type_t>
675 : 661514 : void nasolver<nr_type_t>::createDMatrix (void)
676 : : {
677 : 661514 : int M = countVoltageSources ();
678 [ + - ]: 661514 : int N = countNodes ();
679 : : circuit * vsr, * vsc;
680 : 9665 : nr_type_t val;
681 [ + + ][ + + ]: 2665803 : for (int r = 0; r < M; r++)
682 : : {
683 [ + - ]: 2004289 : vsr = findVoltageSource (r);
684 [ + + ][ + + ]: 10507840 : for (int c = 0; c < M; c++)
685 : : {
686 [ + - ]: 8503551 : vsc = findVoltageSource (c);
687 : 8503551 : val = 0.0;
688 [ + + ]: 8503551 : if (vsr == vsc)
689 : : {
690 [ + - ]: 2433325 : val = MatVal (vsr->getD (r, c));
691 : : }
692 [ + - ]: 8503551 : A->set (r + N, c + N, val);
693 : : }
694 : : }
695 : 661514 : }
696 : :
697 : : /* The G matrix is an NxN matrix formed in two steps.
698 : : 1. Each element in the diagonal matrix is equal to the sum of the
699 : : conductance of each element connected to the corresponding node.
700 : : 2. The off diagonal elements are the negative conductance of the
701 : : element connected to the pair of corresponding nodes. Therefore a
702 : : resistor between nodes 1 and 2 goes into the G matrix at location
703 : : (1,2) and location (2,1). If an element is grounded, it will only
704 : : have contribute to one entry in the G matrix -- at the appropriate
705 : : location on the diagonal. */
706 : : template <class nr_type_t>
707 : 661514 : void nasolver<nr_type_t>::createGMatrix (void)
708 : : {
709 [ + - ]: 661514 : int pr, pc, N = countNodes ();
710 : 9665 : nr_type_t g;
711 : : struct nodelist_t * nr, * nc;
712 : : circuit * ct;
713 : :
714 : : // go through each column of the G matrix
715 [ + + ][ + + ]: 3864187 : for (int c = 0; c < N; c++)
716 : : {
717 [ + - ]: 3202673 : nc = nlist->getNode (c);
718 : : // go through each row of the G matrix
719 [ + + ][ + + ]: 26557656 : for (int r = 0; r < N; r++)
720 : : {
721 [ + - ]: 23354983 : nr = nlist->getNode (r);
722 : 23354983 : g = 0.0;
723 : : // sum up the conductance of each connected circuit
724 [ + + ]: 86160963 : for (int a = 0; a < nc->nNodes; a++)
725 [ + + ]: 237021650 : for (int b = 0; b < nr->nNodes; b++)
726 [ + + ]: 174215670 : if (nc->nodes[a]->getCircuit () == nr->nodes[b]->getCircuit ())
727 : : {
728 : 18419608 : ct = nc->nodes[a]->getCircuit ();
729 : 18419608 : pc = nc->nodes[a]->getPort ();
730 : 18419608 : pr = nr->nodes[b]->getPort ();
731 [ + - ]: 18419608 : g += MatVal (ct->getY (pr, pc));
732 : : }
733 : : // put value into G matrix
734 [ + - ]: 23354983 : A->set (r, c, g);
735 : : }
736 : : }
737 : 661514 : }
738 : :
739 : : /* The following function creates the (N+M)x(N+M) noise current
740 : : correlation matrix used during the AC noise computations. */
741 : : template <class nr_type_t>
742 : 3003 : void nasolver<nr_type_t>::createNoiseMatrix (void)
743 : : {
744 [ + - ]: 3003 : int pr, pc, N = countNodes ();
745 : 3003 : int M = countVoltageSources ();
746 : : struct nodelist_t * n;
747 : 3003 : nr_type_t val;
748 : : int r, c, a, b, ri, ci, i;
749 : : struct nodelist_t * nr, * nc;
750 : : circuit * ct;
751 : :
752 : : // create new Cy matrix if necessary
753 [ + + ][ + - ]: 3003 : if (C != NULL) delete C;
754 [ + - ][ + - ]: 3003 : C = new tmatrix<nr_type_t> (N + M);
755 : :
756 : : // go through each column of the Cy matrix
757 [ + + ]: 17199 : for (c = 0; c < N; c++)
758 : : {
759 [ + - ]: 14196 : nc = nlist->getNode (c);
760 : : // go through each row of the Cy matrix
761 [ + + ]: 83538 : for (r = 0; r < N; r++)
762 : : {
763 [ + - ]: 69342 : nr = nlist->getNode (r);
764 : 69342 : val = 0.0;
765 : : // sum up the noise-correlation of each connected circuit
766 [ + + ]: 262626 : for (a = 0; a < nc->nNodes; a++)
767 [ + + ]: 732732 : for (b = 0; b < nr->nNodes; b++)
768 [ + + ]: 539448 : if (nc->nodes[a]->getCircuit () == nr->nodes[b]->getCircuit ())
769 : : {
770 : 78078 : ct = nc->nodes[a]->getCircuit ();
771 : 78078 : pc = nc->nodes[a]->getPort ();
772 : 78078 : pr = nr->nodes[b]->getPort ();
773 [ + - ]: 78078 : val += MatVal (ct->getN (pr, pc));
774 : : }
775 : : // put value into Cy matrix
776 [ + - ]: 69342 : C->set (r, c, val);
777 : : }
778 : : }
779 : :
780 : : // go through each additional voltage source and put coefficients into
781 : : // the noise current correlation matrix
782 : : circuit * vsr, * vsc;
783 [ + + ]: 14469 : for (r = 0; r < M; r++)
784 : : {
785 [ + - ]: 11466 : vsr = findVoltageSource (r);
786 [ + + ]: 56238 : for (c = 0; c < M; c++)
787 : : {
788 [ + - ]: 44772 : vsc = findVoltageSource (c);
789 : 44772 : val = 0.0;
790 [ + + ]: 44772 : if (vsr == vsc)
791 : : {
792 : 16926 : ri = vsr->getSize () + r - vsr->getVoltageSource ();
793 : 16926 : ci = vsc->getSize () + c - vsc->getVoltageSource ();
794 [ + - ]: 16926 : val = MatVal (vsr->getN (ri, ci));
795 : : }
796 [ + - ]: 44772 : C->set (r + N, c + N, val);
797 : : }
798 : : }
799 : :
800 : : // go through each additional voltage source
801 [ + + ]: 14469 : for (r = 0; r < M; r++)
802 : : {
803 [ + - ]: 11466 : vsr = findVoltageSource (r);
804 : : // go through each node
805 [ + + ]: 67158 : for (c = 0; c < N; c++)
806 : : {
807 : 55692 : val = 0.0;
808 [ + - ]: 55692 : n = nlist->getNode (c);
809 [ + + ]: 210756 : for (i = 0; i < n->nNodes; i++)
810 : : {
811 : : // is voltage source connected to node ?
812 [ + + ]: 155064 : if (n->nodes[i]->getCircuit () == vsr)
813 : : {
814 : 25116 : ri = vsr->getSize () + r - vsr->getVoltageSource ();
815 : 25116 : ci = n->nodes[i]->getPort ();
816 [ + - ]: 25116 : val += MatVal (vsr->getN (ri, ci));
817 : : }
818 : : }
819 : : // put value into Cy matrix
820 [ + - ]: 55692 : C->set (r + N, c, val);
821 : : }
822 : : }
823 : :
824 : : // go through each voltage source
825 [ + + ]: 14469 : for (c = 0; c < M; c++)
826 : : {
827 [ + - ]: 11466 : vsc = findVoltageSource (c);
828 : : // go through each node
829 [ + + ]: 67158 : for (r = 0; r < N; r++)
830 : : {
831 : 55692 : val = 0.0;
832 [ + - ]: 55692 : n = nlist->getNode (r);
833 [ + + ]: 210756 : for (i = 0; i < n->nNodes; i++)
834 : : {
835 : : // is voltage source connected to node ?
836 [ + + ]: 155064 : if (n->nodes[i]->getCircuit () == vsc)
837 : : {
838 : 25116 : ci = vsc->getSize () + c - vsc->getVoltageSource ();
839 : 25116 : ri = n->nodes[i]->getPort ();
840 [ + - ]: 25116 : val += MatVal (vsc->getN (ri, ci));
841 : : }
842 : : }
843 : : // put value into Cy matrix
844 [ + - ]: 55692 : C->set (r, c + N, val);
845 : : }
846 : : }
847 : :
848 : 3003 : }
849 : :
850 : : /* The i matrix is an 1xN matrix with each element of the matrix
851 : : corresponding to a particular node. The value of each element of i
852 : : is determined by the sum of current sources into the corresponding
853 : : node. If there are no current sources connected to the node, the
854 : : value is zero. */
855 : : template <class nr_type_t>
856 : 670763 : void nasolver<nr_type_t>::createIVector (void)
857 : : {
858 [ + - ]: 670763 : int N = countNodes ();
859 : 9665 : nr_type_t val;
860 : : struct nodelist_t * n;
861 : : circuit * is;
862 : :
863 : : // go through each node
864 [ + + ][ + + ]: 3955709 : for (int r = 0; r < N; r++)
865 : : {
866 : 3284946 : val = 0.0;
867 [ + - ]: 3284946 : n = nlist->getNode (r);
868 : : // go through each circuit connected to the node
869 [ + + ][ + + ]: 11999427 : for (int i = 0; i < n->nNodes; i++)
870 : : {
871 : 8714481 : is = n->nodes[i]->getCircuit ();
872 : : // is this a current source ?
873 [ + + ][ + + ]: 8714481 : if (is->isISource () || is->isNonLinear ())
[ + + ]
874 : : {
875 [ + - ]: 3201749 : val += MatVal (is->getI (n->nodes[i]->getPort ()));
876 : : }
877 : : }
878 : : // put value into i vector
879 [ + - ]: 3284946 : z->set (r, val);
880 : : }
881 : 670763 : }
882 : :
883 : : /* The e matrix is a 1xM matrix with each element of the matrix equal
884 : : in value to the corresponding independent voltage source. */
885 : : template <class nr_type_t>
886 : 670763 : void nasolver<nr_type_t>::createEVector (void)
887 : : {
888 [ + - ]: 670763 : int N = countNodes ();
889 : 670763 : int M = countVoltageSources ();
890 : 9665 : nr_type_t val;
891 : : circuit * vs;
892 : :
893 : : // go through each voltage source
894 [ + + ][ + + ]: 2720399 : for (int r = 0; r < M; r++)
895 : : {
896 [ + - ]: 2049636 : vs = findVoltageSource (r);
897 [ + - ]: 2049636 : val = MatVal (vs->getE (r));
898 : : // put value into e vector
899 [ + - ]: 2049636 : z->set (r + N, val);
900 : : }
901 : 670763 : }
902 : :
903 : : // The function loads the right hand side vector.
904 : : template <class nr_type_t>
905 : 670763 : void nasolver<nr_type_t>::createZVector (void)
906 : : {
907 : 670763 : createIVector ();
908 : 670763 : createEVector ();
909 : 670763 : }
910 : :
911 : : // Returns the number of nodes in the nodelist, excluding the ground node.
912 : : template <class nr_type_t>
913 : 6744568 : int nasolver<nr_type_t>::countNodes (void)
914 : : {
915 : 6744568 : return nlist->length () - 1;
916 : : }
917 : :
918 : : // Returns the node number of the give node name.
919 : : template <class nr_type_t>
920 : 5460 : int nasolver<nr_type_t>::getNodeNr (char * str)
921 : : {
922 : 5460 : return nlist->getNodeNr (str);
923 : : }
924 : :
925 : : /* The function returns the assigned node number for the port of the
926 : : given circuits. It returns -1 if there is no such node. */
927 : : template <class nr_type_t>
928 : 0 : int nasolver<nr_type_t>::findAssignedNode (circuit * c, int port)
929 : : {
930 : 0 : int N = countNodes ();
931 : : struct nodelist_t * n;
932 [ # # ]: 0 : for (int r = 0; r < N; r++)
933 : : {
934 : 0 : n = nlist->getNode (r);
935 [ # # ]: 0 : for (int i = 0; i < n->nNodes; i++)
936 [ # # ]: 0 : if (c == n->nodes[i]->getCircuit ())
937 [ # # ]: 0 : if (port == n->nodes[i]->getPort ())
938 : 0 : return r;
939 : : }
940 : 0 : return -1;
941 : : }
942 : :
943 : : // Returns the number of voltage sources in the nodelist.
944 : : template <class nr_type_t>
945 : 4509149 : int nasolver<nr_type_t>::countVoltageSources (void)
946 : : {
947 : 4509149 : return subnet->getVoltageSources ();
948 : : }
949 : :
950 : : /* The function returns the voltage source circuit object
951 : : corresponding to the given number. If there is no such voltage
952 : : source it returns NULL. */
953 : : template <class nr_type_t>
954 : 20542191 : circuit * nasolver<nr_type_t>::findVoltageSource (int n)
955 : : {
956 : 20542191 : circuit * root = subnet->getRoot ();
957 [ + - ]: 155556807 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
958 : : {
959 [ + + ][ + + ]: 155556807 : if (n >= c->getVoltageSource () &&
[ + + ]
960 : : n <= c->getVoltageSource () + c->getVoltageSources () - 1)
961 : 20542191 : return c;
962 : : }
963 : 20542191 : return NULL;
964 : : }
965 : :
966 : : /* The function applies unique voltage source identifiers to each
967 : : voltage source (explicit and built in internal ones) in the list of
968 : : registered circuits. */
969 : : template <class nr_type_t>
970 : 22664 : void nasolver<nr_type_t>::assignVoltageSources (void)
971 : : {
972 : 22664 : circuit * root = subnet->getRoot ();
973 : 22664 : int nSources = 0;
974 [ + + ]: 141291 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
975 : : {
976 [ + + ]: 118627 : if (c->getVoltageSources () > 0)
977 : : {
978 : 30974 : c->setVoltageSource (nSources);
979 : 30974 : nSources += c->getVoltageSources ();
980 : : }
981 : : }
982 : 22664 : subnet->setVoltageSources (nSources);
983 : 22664 : }
984 : :
985 : : /* The matrix equation Ax = z is solved by x = A^-1*z. The function
986 : : applies the operation to the previously generated matrices. */
987 : : template <class nr_type_t>
988 : 687176 : void nasolver<nr_type_t>::runMNA (void)
989 : : {
990 : :
991 : : // just solve the equation system here
992 : 687176 : eqns->setAlgo (eqnAlgo);
993 [ + + ]: 687176 : eqns->passEquationSys (updateMatrix ? A : NULL, x, z);
994 : 687176 : eqns->solve ();
995 : :
996 : : // if damped Newton-Raphson is requested
997 [ + + ][ + + ]: 687176 : if (xprev != NULL && top_exception () == NULL)
[ + + ]
998 : : {
999 [ - + ]: 650130 : if (convHelper == CONV_Attenuation)
1000 : : {
1001 : 0 : applyAttenuation ();
1002 : : }
1003 [ - + ]: 650130 : else if (convHelper == CONV_LineSearch)
1004 : : {
1005 : 0 : lineSearch ();
1006 : : }
1007 [ + + ]: 650130 : else if (convHelper == CONV_SteepestDescent)
1008 : : {
1009 : 654 : steepestDescent ();
1010 : : }
1011 : : }
1012 : 687176 : }
1013 : :
1014 : : /* This function applies a damped Newton-Raphson (limiting scheme) to
1015 : : the current solution vector in the form x1 = x0 + a * (x1 - x0). This
1016 : : convergence helper is heuristic and does not ensure global convergence. */
1017 : : template <class nr_type_t>
1018 : 0 : void nasolver<nr_type_t>::applyAttenuation (void)
1019 : : {
1020 : 0 : nr_double_t alpha = 1.0, nMax;
1021 : :
1022 : : // create solution difference vector and find maximum deviation
1023 [ # # ][ # # ]: 0 : tvector<nr_type_t> dx = *x - *xprev;
[ # # ]
1024 [ # # ][ # # ]: 0 : nMax = maxnorm (dx);
1025 : :
1026 : : // compute appropriate damping factor
1027 [ # # ]: 0 : if (nMax > 0.0)
1028 : : {
1029 : 0 : nr_double_t g = 1.0;
1030 [ # # ]: 0 : alpha = MIN (0.9, g / nMax);
1031 [ # # ]: 0 : if (alpha < 0.1) alpha = 0.1;
1032 : : }
1033 : :
1034 : : // apply damped solution vector
1035 [ # # ][ # # ]: 0 : *x = *xprev + alpha * dx;
[ # # ][ # # ]
[ # # ]
1036 : 0 : }
1037 : :
1038 : : /* This is damped Newton-Raphson using nested iterations in order to
1039 : : find a better damping factor. It identifies a damping factor in
1040 : : the interval [0,1] which minimizes the right hand side vector. The
1041 : : algorithm actually ensures global convergence but pushes the
1042 : : solution to local minimums, i.e. where the Jacobian matrix A may be
1043 : : singular. */
1044 : : template <class nr_type_t>
1045 : 0 : void nasolver<nr_type_t>::lineSearch (void)
1046 : : {
1047 : 0 : nr_double_t alpha = 0.5, n, nMin, aprev = 1.0, astep = 0.5, adiff;
1048 : 0 : int dir = -1;
1049 : :
1050 : : // compute solution deviation vector
1051 [ # # ][ # # ]: 0 : tvector<nr_type_t> dx = *x - *xprev;
[ # # ]
1052 : 0 : nMin = std::numeric_limits<nr_double_t>::max();
1053 : :
1054 [ # # ]: 0 : do
1055 : : {
1056 : : // apply current damping factor and see what happens
1057 [ # # ][ # # ]: 0 : *x = *xprev + alpha * dx;
[ # # ][ # # ]
[ # # ]
1058 : :
1059 : : // recalculate Jacobian and right hand side
1060 [ # # ]: 0 : saveSolution ();
1061 [ # # ]: 0 : calculate ();
1062 [ # # ]: 0 : createZVector ();
1063 : :
1064 : : // calculate norm of right hand side vector
1065 [ # # ][ # # ]: 0 : n = norm (*z);
1066 : :
1067 : : // TODO: this is not perfect, but usable
1068 : 0 : astep /= 2;
1069 : 0 : adiff = fabs (alpha - aprev);
1070 [ # # ]: 0 : if (adiff > 0.005)
1071 : : {
1072 : 0 : aprev = alpha;
1073 [ # # ]: 0 : if (n < nMin)
1074 : : {
1075 : 0 : nMin = n;
1076 [ # # ]: 0 : if (alpha == 1) dir = -dir;
1077 : 0 : alpha += astep * dir;
1078 : : }
1079 : : else
1080 : : {
1081 : 0 : dir = -dir;
1082 : 0 : alpha += 1.5 * astep * dir;
1083 : : }
1084 : : }
1085 : : }
1086 : : while (adiff > 0.005);
1087 : :
1088 : : // apply final damping factor
1089 [ # # ][ # # ]: 0 : assert (alpha > 0 && alpha <= 1);
[ # # ]
1090 [ # # ][ # # ]: 0 : *x = *xprev + alpha * dx;
[ # # ][ # # ]
[ # # ]
1091 : 0 : }
1092 : :
1093 : : /* The function looks for the optimal gradient for the right hand side
1094 : : vector using the so-called 'steepest descent' method. Though
1095 : : better than the one-dimensional linesearch (it doesn't push
1096 : : iterations into local minimums) it converges painfully slow. */
1097 : : template <class nr_type_t>
1098 : 654 : void nasolver<nr_type_t>::steepestDescent (void)
1099 : : {
1100 : 654 : nr_double_t alpha = 1.0, sl, n;
1101 : :
1102 : : // compute solution deviation vector
1103 [ + - ][ + - ]: 654 : tvector<nr_type_t> dx = *x - *xprev;
[ + - ]
1104 [ + - ][ + - ]: 654 : tvector<nr_type_t> dz = *z - *zprev;
[ + - ]
1105 [ + - ][ + - ]: 654 : n = norm (*zprev);
1106 : :
1107 [ + + ]: 9027 : do
1108 : : {
1109 : : // apply current damping factor and see what happens
1110 [ + - ][ + - ]: 9249 : *x = *xprev + alpha * dx;
[ + - ][ + - ]
[ + - ]
1111 : :
1112 : : // recalculate Jacobian and right hand side
1113 [ + - ]: 9249 : saveSolution ();
1114 [ + - ]: 9249 : calculate ();
1115 [ + - ]: 9249 : createZVector ();
1116 : :
1117 : : // check gradient criteria, ThinkME: Is this correct?
1118 [ + - ][ + - ]: 9249 : dz = *z - *zprev;
[ + - ][ + - ]
1119 [ + - ][ + - ]: 9249 : sl = real (sum (dz * -dz));
[ + - ][ + - ]
[ + - ][ + - ]
1120 [ + - ][ + - ]: 9249 : if (norm (*z) < n + alpha * sl) break;
[ + + ]
1121 : 9027 : alpha *= 0.7;
1122 : : }
1123 : : while (alpha > 0.001);
1124 : :
1125 : : // apply final damping factor
1126 [ + - ][ + - ]: 654 : *x = *xprev + alpha * dx;
[ + - ][ + - ]
[ + - ]
1127 : 654 : }
1128 : :
1129 : : /* The function checks whether the iterative algorithm for linearizing
1130 : : the non-linear components in the network shows convergence. It
1131 : : returns non-zero if it converges and zero otherwise. */
1132 : : template <class nr_type_t>
1133 : 414914 : int nasolver<nr_type_t>::checkConvergence (void)
1134 : : {
1135 : :
1136 : 414914 : int N = countNodes ();
1137 : 414914 : int M = countVoltageSources ();
1138 : : nr_double_t v_abs, v_rel, i_abs, i_rel;
1139 : : int r;
1140 : :
1141 : : // check the nodal voltage changes against the allowed absolute
1142 : : // and relative tolerance values
1143 [ + + ]: 1773355 : for (r = 0; r < N; r++)
1144 : : {
1145 : 1536976 : v_abs = abs (x->get (r) - xprev->get (r));
1146 : 1536976 : v_rel = abs (x->get (r));
1147 [ + + ]: 1536976 : if (v_abs >= vntol + reltol * v_rel) return 0;
1148 [ + + ]: 1455200 : if (!convHelper)
1149 : : {
1150 : 1452940 : i_abs = abs (z->get (r) - zprev->get (r));
1151 : 1452940 : i_rel = abs (z->get (r));
1152 [ + + ]: 1452940 : if (i_abs >= abstol + reltol * i_rel) return 0;
1153 : : }
1154 : : }
1155 : :
1156 [ + + ]: 1062739 : for (r = 0; r < M; r++)
1157 : : {
1158 : 827504 : i_abs = abs (x->get (r + N) - xprev->get (r + N));
1159 : 827504 : i_rel = abs (x->get (r + N));
1160 [ + + ]: 827504 : if (i_abs >= abstol + reltol * i_rel) return 0;
1161 [ + + ]: 826462 : if (!convHelper)
1162 : : {
1163 : 825523 : v_abs = abs (z->get (r + N) - zprev->get (r + N));
1164 : 825523 : v_rel = abs (z->get (r + N));
1165 [ + + ]: 825523 : if (v_abs >= vntol + reltol * v_rel) return 0;
1166 : : }
1167 : : }
1168 : 414914 : return 1;
1169 : : }
1170 : :
1171 : : /* The function saves the solution and right hand vector of the previous
1172 : : iteration. */
1173 : : template <class nr_type_t>
1174 : 650214 : void nasolver<nr_type_t>::savePreviousIteration (void)
1175 : : {
1176 [ + + ]: 650214 : if (xprev != NULL)
1177 : 650141 : *xprev = *x;
1178 : : else
1179 [ + - ]: 73 : xprev = new tvector<nr_type_t> (*x);
1180 [ + + ]: 650214 : if (zprev != NULL)
1181 : 650141 : *zprev = *z;
1182 : : else
1183 [ + - ]: 73 : zprev = new tvector<nr_type_t> (*z);
1184 : 650214 : }
1185 : :
1186 : : /* The function restores the solution and right hand vector of the
1187 : : previous (successful) iteration. */
1188 : : template <class nr_type_t>
1189 : 0 : void nasolver<nr_type_t>::restorePreviousIteration (void)
1190 : : {
1191 [ # # ]: 0 : if (xprev != NULL) *x = *xprev;
1192 [ # # ]: 0 : if (zprev != NULL) *z = *zprev;
1193 : 0 : }
1194 : :
1195 : : /* The function restarts the NR iteration for each non-linear
1196 : : circuit. */
1197 : : template <class nr_type_t>
1198 : : void nasolver<nr_type_t>::restartNR (void)
1199 : : {
1200 : : circuit * root = subnet->getRoot ();
1201 : : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
1202 : : {
1203 : : if (c->isNonLinear ()) c->restartDC ();
1204 : : }
1205 : : }
1206 : :
1207 : : /* This function goes through solution (the x vector) and saves the
1208 : : node voltages of the last iteration into each non-linear
1209 : : circuit. */
1210 : : template <class nr_type_t>
1211 : 903125 : void nasolver<nr_type_t>::saveNodeVoltages (void)
1212 : : {
1213 : 903125 : int N = countNodes ();
1214 : : struct nodelist_t * n;
1215 : : // save all nodes except reference node
1216 [ + + ]: 5370952 : for (int r = 0; r < N; r++)
1217 : : {
1218 : 4467827 : n = nlist->getNode (r);
1219 [ + + ][ + + ]: 16366309 : for (int i = 0; i < n->nNodes; i++)
1220 : : {
1221 [ + - ]: 11898482 : n->nodes[i]->getCircuit()->setV (n->nodes[i]->getPort (), x->get (r));
1222 : : }
1223 : : }
1224 : : // save reference node
1225 : 903125 : n = nlist->getNode (-1);
1226 [ + + ]: 5452855 : for (int i = 0; i < n->nNodes; i++)
1227 : : {
1228 [ + - ]: 4549730 : n->nodes[i]->getCircuit()->setV (n->nodes[i]->getPort (), 0.0);
1229 : : }
1230 : 903125 : }
1231 : :
1232 : : /* This function goes through solution (the x vector) and saves the
1233 : : branch currents through the voltage sources of the last iteration
1234 : : into each circuit. */
1235 : : template <class nr_type_t>
1236 : 903125 : void nasolver<nr_type_t>::saveBranchCurrents (void)
1237 : : {
1238 : 903125 : int N = countNodes ();
1239 : 903125 : int M = countVoltageSources ();
1240 : : circuit * vs;
1241 : : // save all branch currents of voltage sources
1242 [ + + ][ + + ]: 3767883 : for (int r = 0; r < M; r++)
1243 : : {
1244 : 2864758 : vs = findVoltageSource (r);
1245 [ + - ]: 2864758 : vs->setJ (r, x->get (r + N));
1246 : : }
1247 : 903125 : }
1248 : :
1249 : : // The function saves the solution vector into each circuit.
1250 : : template <class nr_type_t>
1251 : 903125 : void nasolver<nr_type_t>::saveSolution (void)
1252 : : {
1253 : 903125 : saveNodeVoltages ();
1254 : 903125 : saveBranchCurrents ();
1255 : 903125 : }
1256 : :
1257 : : // This function stores the solution (node voltages and branch currents).
1258 : : template <class nr_type_t>
1259 : 64 : void nasolver<nr_type_t>::storeSolution (void)
1260 : : {
1261 : : // cleanup solution previously
1262 : 64 : solution.clear ();
1263 : : int r;
1264 : 64 : int N = countNodes ();
1265 : 64 : int M = countVoltageSources ();
1266 : : // store all nodes except reference node
1267 [ + + ]: 496 : for (r = 0; r < N; r++)
1268 : : {
1269 : 432 : struct nodelist_t * n = nlist->getNode (r);
1270 : 432 : solution.add (n->name, x->get (r), 0);
1271 : : }
1272 : : // store all branch currents of voltage sources
1273 [ + + ]: 338 : for (r = 0; r < M; r++)
1274 : : {
1275 : 274 : circuit * vs = findVoltageSource (r);
1276 : 274 : int vn = r - vs->getVoltageSource () + 1;
1277 : 274 : solution.add (vs->getName (), x->get (r + N), vn);
1278 : : }
1279 : 64 : }
1280 : :
1281 : : // This function recalls the solution (node voltages and branch currents).
1282 : : template <class nr_type_t>
1283 : 65 : void nasolver<nr_type_t>::recallSolution (void)
1284 : : {
1285 : : int r;
1286 : 65 : int N = countNodes ();
1287 : 65 : int M = countVoltageSources ();
1288 : : naentry<nr_type_t> * na;
1289 : : // store all nodes except reference node
1290 [ + + ]: 506 : for (r = 0; r < N; r++)
1291 : : {
1292 : 441 : struct nodelist_t * n = nlist->getNode (r);
1293 [ + + ]: 441 : if ((na = solution.find (n->name, 0)) != NULL)
1294 : 432 : x->set (r, na->value);
1295 : : }
1296 : : // store all branch currents of voltage sources
1297 [ + + ]: 335 : for (r = 0; r < M; r++)
1298 : : {
1299 : 270 : circuit * vs = findVoltageSource (r);
1300 : 270 : int vn = r - vs->getVoltageSource () + 1;
1301 [ + + ]: 270 : if ((na = solution.find (vs->getName (), vn)) != NULL)
1302 : 265 : x->set (r + N, na->value);
1303 : : }
1304 : 65 : }
1305 : :
1306 : : /* This function saves the results of a single solve() functionality
1307 : : into the output dataset. */
1308 : : template <class nr_type_t>
1309 : 75228 : void nasolver<nr_type_t>::saveResults (const char * volts, const char * amps,
1310 : : int saveOPs, qucs::vector * f)
1311 : : {
1312 : 75228 : int N = countNodes ();
1313 : 75228 : int M = countVoltageSources ();
1314 : : char * n;
1315 : :
1316 : : // add node voltage variables
1317 [ + - ]: 75228 : if (volts)
1318 : : {
1319 [ + + ]: 459042 : for (int r = 0; r < N; r++)
1320 : : {
1321 [ + + ]: 383814 : if ((n = createV (r, volts, saveOPs)) != NULL)
1322 : : {
1323 [ + - ]: 319473 : saveVariable (n, x->get (r), f);
1324 : 319473 : free (n);
1325 : : }
1326 : : }
1327 : : }
1328 : :
1329 : : // add branch current variables
1330 [ + - ]: 75228 : if (amps)
1331 : : {
1332 [ + + ]: 308460 : for (int r = 0; r < M; r++)
1333 : : {
1334 [ + + ]: 233232 : if ((n = createI (r, amps, saveOPs)) != NULL)
1335 : : {
1336 [ + - ]: 133643 : saveVariable (n, x->get (r + N), f);
1337 : 133643 : free (n);
1338 : : }
1339 : : }
1340 : : }
1341 : :
1342 : : // add voltage probe data
1343 [ + - ]: 75228 : if (volts)
1344 : : {
1345 : 75228 : circuit * root = subnet->getRoot ();
1346 [ + + ]: 701826 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
1347 : : {
1348 [ + + ]: 626598 : if (!c->isProbe ()) continue;
1349 [ - + ][ # # ]: 36629 : if (c->getSubcircuit () && !(saveOPs & SAVE_ALL)) continue;
[ - + ]
1350 [ + + ]: 36629 : if (strcmp (volts, "vn"))
1351 : 33899 : c->saveOperatingPoints ();
1352 : 36629 : n = createOP (c->getName (), volts);
1353 [ + - ]: 36629 : saveVariable (n, nr_complex_t (c->getOperatingPoint ("Vr"),
1354 : : c->getOperatingPoint ("Vi")), f);
1355 : 36629 : free (n);
1356 : : }
1357 : : }
1358 : :
1359 : : // save operating points of non-linear circuits if requested
1360 [ - + ]: 75228 : if (saveOPs & SAVE_OPS)
1361 : : {
1362 : 0 : circuit * root = subnet->getRoot ();
1363 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
1364 : : {
1365 [ # # ]: 0 : if (!c->isNonLinear ()) continue;
1366 [ # # ][ # # ]: 0 : if (c->getSubcircuit () && !(saveOPs & SAVE_ALL)) continue;
[ # # ]
1367 [ # # ]: 0 : c->calcOperatingPoints ();
1368 [ # # ]: 0 : valuelistiterator<operatingpoint> it (c->getOperatingPoints ());
1369 [ # # ][ # # ]: 0 : for (; *it; ++it)
1370 : : {
1371 : 0 : operatingpoint * p = it.currentVal ();
1372 : 0 : n = createOP (c->getName (), p->getName ());
1373 [ # # ]: 0 : saveVariable (n, p->getValue (), f);
1374 : 0 : free (n);
1375 : : }
1376 : : }
1377 : : }
1378 : 75228 : }
1379 : :
1380 : : /* Create an appropriate variable name for operating points. The
1381 : : caller is responsible to free() the returned string. */
1382 : : template <class nr_type_t>
1383 : 36629 : char * nasolver<nr_type_t>::createOP (const char * c, const char * n)
1384 : : {
1385 : 36629 : char * text = (char *) malloc (strlen (c) + strlen (n) + 2);
1386 : 36629 : sprintf (text, "%s.%s", c, n);
1387 : 36629 : return text;
1388 : : }
1389 : :
1390 : : /* Creates an appropriate variable name for voltages. The caller is
1391 : : responsible to free() the returned string. */
1392 : : template <class nr_type_t>
1393 : 383814 : char * nasolver<nr_type_t>::createV (int n, const char * volts, int saveOPs)
1394 : : {
1395 [ + + ]: 383814 : if (nlist->isInternal (n)) return NULL;
1396 : 335347 : char * node = nlist->get (n);
1397 [ + + ][ + - ]: 335347 : if (strchr (node, '.') && !(saveOPs & SAVE_ALL)) return NULL;
1398 : 319473 : char * text = (char *) malloc (strlen (node) + 2 + strlen (volts));
1399 : 319473 : sprintf (text, "%s.%s", node, volts);
1400 : 383814 : return text;
1401 : : }
1402 : :
1403 : : /* Create an appropriate variable name for currents. The caller is
1404 : : responsible to free() the returned string. */
1405 : : template <class nr_type_t>
1406 : 233232 : char * nasolver<nr_type_t>::createI (int n, const char * amps, int saveOPs)
1407 : : {
1408 : 233232 : circuit * vs = findVoltageSource (n);
1409 : :
1410 : : // don't output internal (helper) voltage sources
1411 [ + + ]: 233232 : if (vs->isInternalVoltageSource ())
1412 : 808 : return NULL;
1413 : :
1414 : : /* save only current through real voltage sources and explicit
1415 : : current probes */
1416 [ + + ][ + - ]: 232424 : if (!vs->isVSource () && !(saveOPs & SAVE_OPS))
[ + + ]
1417 : 98280 : return NULL;
1418 : :
1419 : : // don't output subcircuit components if not requested
1420 [ + + ][ + - ]: 134144 : if (vs->getSubcircuit () && !(saveOPs & SAVE_ALL))
[ + + ]
1421 : 501 : return NULL;
1422 : :
1423 : : // create appropriate current name for single/multiple voltage sources
1424 : 133643 : char * name = vs->getName ();
1425 : 133643 : char * text = (char *) malloc (strlen (name) + 4 + strlen (amps));
1426 [ - + ]: 133643 : if (vs->getVoltageSources () > 1)
1427 : 0 : sprintf (text, "%s.%s%d", name, amps, n - vs->getVoltageSource () + 1);
1428 : : else
1429 : 133643 : sprintf (text, "%s.%s", name, amps);
1430 : 233232 : return text;
1431 : : }
1432 : :
1433 : :
1434 : : /* Alternaive to countNodes () */
1435 : : template <class nr_type_t>
1436 : 0 : int nasolver<nr_type_t>::getN()
1437 : : {
1438 : 0 : return countNodes ();
1439 : : }
1440 : :
1441 : : /* Alternative to countVoltageSources () */
1442 : : template <class nr_type_t>
1443 : 0 : int nasolver<nr_type_t>::getM()
1444 : : {
1445 : 0 : return countVoltageSources ();
1446 : : }
1447 : :
1448 : : } // namespace qucs
|
