Qucs-core  0.0.19
nasolver.cpp
Go to the documentation of this file.
00001 /*
00002  * nasolver.cpp - nodal analysis solver class implementation
00003  *
00004  * Copyright (C) 2004, 2005, 2006, 2007, 2008 Stefan Jahn <stefan@lkcc.org>
00005  *
00006  * This is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2, or (at your option)
00009  * any later version.
00010  *
00011  * This software is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this package; see the file COPYING.  If not, write to
00018  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00019  * Boston, MA 02110-1301, USA.
00020  *
00021  * $Id$
00022  *
00023  */
00024 
00025 // the types required for qucs library files are defined
00026 // in qucs_typedefs.h, created by configure from
00027 // qucs_typedefs.h.in
00028 #include "qucs_typedefs.h"
00029 
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include <string.h>
00033 #include <cmath>
00034 #include <float.h>
00035 #include <assert.h>
00036 #include <limits>
00037 
00038 #include "logging.h"
00039 #include "complex.h"
00040 #include "object.h"
00041 #include "node.h"
00042 #include "circuit.h"
00043 #include "vector.h"
00044 #include "dataset.h"
00045 #include "net.h"
00046 #include "analysis.h"
00047 #include "nodelist.h"
00048 #include "nodeset.h"
00049 #include "strlist.h"
00050 #include "tvector.h"
00051 #include "tmatrix.h"
00052 #include "eqnsys.h"
00053 #include "precision.h"
00054 #include "operatingpoint.h"
00055 #include "exception.h"
00056 #include "exceptionstack.h"
00057 #include "nasolver.h"
00058 #include "constants.h"
00059 
00060 namespace qucs {
00061 
00062 // Constructor creates an unnamed instance of the nasolver class.
00063 template <class nr_type_t>
00064 nasolver<nr_type_t>::nasolver () : analysis ()
00065 {
00066     nlist = NULL;
00067     A = C = NULL;
00068     z = x = xprev = zprev = NULL;
00069     reltol = abstol = vntol = 0;
00070     calculate_func = NULL;
00071     convHelper = fixpoint = 0;
00072     eqnAlgo = ALGO_LU_DECOMPOSITION;
00073     updateMatrix = 1;
00074     gMin = srcFactor = 0;
00075     eqns = new eqnsys<nr_type_t> ();
00076 }
00077 
00078 // Constructor creates a named instance of the nasolver class.
00079 template <class nr_type_t>
00080 nasolver<nr_type_t>::nasolver (const std::string &n) : analysis (n)
00081 {
00082     nlist = NULL;
00083     A = C = NULL;
00084     z = x = xprev = zprev = NULL;
00085     reltol = abstol = vntol = 0;
00086     calculate_func = NULL;
00087     convHelper = fixpoint = 0;
00088     eqnAlgo = ALGO_LU_DECOMPOSITION;
00089     updateMatrix = 1;
00090     gMin = srcFactor = 0;
00091     eqns = new eqnsys<nr_type_t> ();
00092 }
00093 
00094 // Destructor deletes the nasolver class object.
00095 template <class nr_type_t>
00096 nasolver<nr_type_t>::~nasolver ()
00097 {
00098     delete nlist;
00099     delete C;
00100     delete A;
00101     delete z;
00102     delete x;
00103     delete xprev;
00104     delete zprev;
00105     delete eqns;
00106 }
00107 
00108 /* The copy constructor creates a new instance of the nasolver class
00109    based on the given nasolver object. */
00110 template <class nr_type_t>
00111 nasolver<nr_type_t>::nasolver (nasolver & o) : analysis (o)
00112 {
00113     nlist = o.nlist ? new nodelist (*(o.nlist)) : NULL;
00114     A = o.A ? new tmatrix<nr_type_t> (*(o.A)) : NULL;
00115     C = o.C ? new tmatrix<nr_type_t> (*(o.C)) : NULL;
00116     z = o.z ? new tvector<nr_type_t> (*(o.z)) : NULL;
00117     x = o.x ? new tvector<nr_type_t> (*(o.x)) : NULL;
00118     xprev = zprev = NULL;
00119     reltol = o.reltol;
00120     abstol = o.abstol;
00121     vntol = o.vntol;
00122     desc = o.desc;
00123     calculate_func = o.calculate_func;
00124     convHelper = o.convHelper;
00125     eqnAlgo = o.eqnAlgo;
00126     updateMatrix = o.updateMatrix;
00127     fixpoint = o.fixpoint;
00128     gMin = o.gMin;
00129     srcFactor = o.srcFactor;
00130     eqns = new eqnsys<nr_type_t> (*(o.eqns));
00131     solution = nasolution<nr_type_t> (o.solution);
00132 }
00133 
00134 /* The function runs the nodal analysis solver once, reports errors if
00135    any and save the results into each circuit. */
00136 template <class nr_type_t>
00137 int nasolver<nr_type_t>::solve_once (void)
00138 {
00139     qucs::exception * e;
00140     int error = 0, d;
00141 
00142     // run the calculation function for each circuit
00143     calculate ();
00144 
00145     // generate A matrix and z vector
00146     createMatrix ();
00147 
00148     // solve equation system
00149     try_running ()
00150     {
00151         runMNA ();
00152     }
00153     // appropriate exception handling
00154     catch_exception ()
00155     {
00156     case EXCEPTION_PIVOT:
00157     case EXCEPTION_WRONG_VOLTAGE:
00158         e = new qucs::exception (EXCEPTION_NA_FAILED);
00159         d = top_exception()->getData ();
00160         pop_exception ();
00161         if (d >= countNodes ())
00162         {
00163             d -= countNodes ();
00164             e->setText ("voltage source `%s' conflicts with some other voltage "
00165                         "source", findVoltageSource(d)->getName ());
00166         }
00167         else
00168         {
00169             e->setText ("circuit admittance matrix in %s solver is singular at "
00170                         "node `%s' connected to [%s]", desc.c_str(), nlist->get (d).c_str(),
00171                         nlist->getNodeString (d).c_str());
00172         }
00173         throw_exception (e);
00174         error++;
00175         break;
00176     case EXCEPTION_SINGULAR:
00177         do
00178         {
00179             d = top_exception()->getData ();
00180             pop_exception ();
00181             if (d < countNodes ())
00182             {
00183                 logprint (LOG_ERROR, "WARNING: %s: inserted virtual resistance at "
00184                           "node `%s' connected to [%s]\n", getName (), nlist->get (d).c_str(),
00185                           nlist->getNodeString (d).c_str());
00186             }
00187         }
00188         while (top_exception() != NULL &&
00189                 top_exception()->getCode () == EXCEPTION_SINGULAR);
00190         break;
00191     default:
00192         estack.print ();
00193         break;
00194     }
00195 
00196     // save results into circuits
00197     if (!error) saveSolution ();
00198     return error;
00199 }
00200 
00201 /* Run this function after the actual solver run and before evaluating
00202    the results. */
00203 template <class nr_type_t>
00204 void nasolver<nr_type_t>::solve_post (void)
00205 {
00206     delete nlist;
00207     nlist = NULL;
00208 }
00209 
00210 /* Run this function before the actual solver. */
00211 template <class nr_type_t>
00212 void nasolver<nr_type_t>::solve_pre (void)
00213 {
00214     // create node list, enumerate nodes and voltage sources
00215 #if DEBUG
00216     logprint (LOG_STATUS, "NOTIFY: %s: creating node list for %s analysis\n",
00217               getName (), desc.c_str());
00218 #endif
00219     nlist = new nodelist (subnet);
00220     nlist->assignNodes ();
00221     assignVoltageSources ();
00222 #if DEBUG && 0
00223     nlist->print ();
00224 #endif
00225 
00226     // create matrix, solution vector and right hand side vector
00227     int M = countVoltageSources ();
00228     int N = countNodes ();
00229     delete A;
00230     A = new tmatrix<nr_type_t> (M + N);
00231     delete z;
00232     z = new tvector<nr_type_t> (N + M);
00233     delete x;
00234     x = new tvector<nr_type_t> (N + M);
00235 
00236 #if DEBUG
00237     logprint (LOG_STATUS, "NOTIFY: %s: solving %s netlist\n", getName (), desc.c_str());
00238 #endif
00239 }
00240 
00241 /* This function goes through the nodeset list of the current netlist
00242    and applies the stored values to the current solution vector.  Then
00243    the function saves the solution vector back into the actual
00244    component nodes. */
00245 template <class nr_type_t>
00246 void nasolver<nr_type_t>::applyNodeset (bool nokeep)
00247 {
00248     if (x == NULL || nlist == NULL) return;
00249 
00250     // set each solution to zero
00251     if (nokeep) for (int i = 0; i < x->size (); i++) x->set (i, 0);
00252 
00253     // then apply the nodeset itself
00254     for (nodeset * n = subnet->getNodeset (); n; n = n->getNext ())
00255     {
00256         struct nodelist_t * nl = nlist->getNode (n->getName ());
00257         if (nl != NULL)
00258         {
00259             x->set (nl->n, n->getValue ());
00260         }
00261         else
00262         {
00263             logprint (LOG_ERROR, "WARNING: %s: no such node `%s' found, cannot "
00264                       "initialize node\n", getName (), n->getName ());
00265         }
00266     }
00267     if (xprev != NULL) *xprev = *x;
00268     saveSolution ();
00269 }
00270 
00271 /* The following function uses the gMin-stepping algorithm in order to
00272    solve the given non-linear netlist by continuous iterations. */
00273 template <class nr_type_t>
00274 int nasolver<nr_type_t>::solve_nonlinear_continuation_gMin (void)
00275 {
00276     qucs::exception * e;
00277     int convergence, run = 0, MaxIterations, error = 0;
00278     nr_double_t gStep, gPrev;
00279 
00280     // fetch simulation properties
00281     MaxIterations = getPropertyInteger ("MaxIter") / 4 + 1;
00282     updateMatrix = 1;
00283     fixpoint = 0;
00284 
00285     // initialize the stepper
00286     gPrev = gMin = 0.01;
00287     gStep = gMin / 100;
00288     gMin -= gStep;
00289 
00290     do
00291     {
00292         // run solving loop until convergence is reached
00293         run = 0;
00294         do
00295         {
00296             error = solve_once ();
00297             if (!error)
00298             {
00299                 // convergence check
00300                 convergence = (run > 0) ? checkConvergence () : 0;
00301                 savePreviousIteration ();
00302                 run++;
00303             }
00304             else break;
00305         }
00306         while (!convergence && run < MaxIterations);
00307         iterations += run;
00308 
00309         // not yet converged, so decreased the gMin-step
00310         if (run >= MaxIterations || error)
00311         {
00312             gStep /= 2;
00313             // here the absolute minimum step checker
00314             if (gStep < std::numeric_limits<nr_double_t>::epsilon())
00315             {
00316                 error = 1;
00317                 e = new qucs::exception (EXCEPTION_NO_CONVERGENCE);
00318                 e->setText ("no convergence in %s analysis after %d gMinStepping "
00319                             "iterations", desc.c_str(), iterations);
00320                 throw_exception (e);
00321                 break;
00322             }
00323             gMin = MAX (gPrev - gStep, 0);
00324         }
00325         // converged, increased the gMin-step
00326         else
00327         {
00328             gPrev = gMin;
00329             gMin = MAX (gMin - gStep, 0);
00330             gStep *= 2;
00331         }
00332     }
00333     // continue until no additional resistances is necessary
00334     while (gPrev > 0);
00335 
00336     return error;
00337 }
00338 
00339 /* The following function uses the source-stepping algorithm in order
00340    to solve the given non-linear netlist by continuous iterations. */
00341 template <class nr_type_t>
00342 int nasolver<nr_type_t>::solve_nonlinear_continuation_Source (void)
00343 {
00344     qucs::exception * e;
00345     int convergence, run = 0, MaxIterations, error = 0;
00346     nr_double_t sStep, sPrev;
00347 
00348     // fetch simulation properties
00349     MaxIterations = getPropertyInteger ("MaxIter") / 4 + 1;
00350     updateMatrix = 1;
00351     fixpoint = 0;
00352 
00353     // initialize the stepper
00354     sPrev = srcFactor = 0;
00355     sStep = 0.01;
00356     srcFactor += sStep;
00357 
00358     do
00359     {
00360         // run solving loop until convergence is reached
00361         run = 0;
00362         do
00363         {
00364             subnet->setSrcFactor (srcFactor);
00365             error = solve_once ();
00366             if (!error)
00367             {
00368                 // convergence check
00369                 convergence = (run > 0) ? checkConvergence () : 0;
00370                 savePreviousIteration ();
00371                 run++;
00372             }
00373             else break;
00374         }
00375         while (!convergence && run < MaxIterations);
00376         iterations += run;
00377 
00378         // not yet converged, so decreased the source-step
00379         if (run >= MaxIterations || error)
00380         {
00381             if (error)
00382                 sStep *= 0.1;
00383             else
00384                 sStep *= 0.5;
00385             restorePreviousIteration ();
00386             saveSolution ();
00387             // here the absolute minimum step checker
00388             if (sStep < std::numeric_limits<nr_double_t>::epsilon())
00389             {
00390                 error = 1;
00391                 e = new qucs::exception (EXCEPTION_NO_CONVERGENCE);
00392                 e->setText ("no convergence in %s analysis after %d sourceStepping "
00393                             "iterations", desc.c_str(), iterations);
00394                 throw_exception (e);
00395                 break;
00396             }
00397             srcFactor = std::min (sPrev + sStep, 1.0);
00398         }
00399         // converged, increased the source-step
00400         else if (run < MaxIterations / 4)
00401         {
00402             sPrev = srcFactor;
00403             srcFactor = std::min (srcFactor + sStep, 1.0);
00404             sStep *= 1.5;
00405         }
00406         else
00407         {
00408             srcFactor = std::min (srcFactor + sStep, 1.0);
00409         }
00410     }
00411     // continue until no source factor is necessary
00412     while (sPrev < 1);
00413 
00414     subnet->setSrcFactor (1);
00415     return error;
00416 }
00417 
00418 /* The function returns an appropriate text representation for the
00419    currently used convergence helper algorithm. */
00420 template <class nr_type_t>
00421 const char * nasolver<nr_type_t>::getHelperDescription (void)
00422 {
00423     if (convHelper == CONV_Attenuation)
00424     {
00425         return "RHS attenuation";
00426     }
00427     else if  (convHelper == CONV_LineSearch)
00428     {
00429         return "line search";
00430     }
00431     else if  (convHelper == CONV_SteepestDescent)
00432     {
00433         return "steepest descent";
00434     }
00435     else if  (convHelper == CONV_GMinStepping)
00436     {
00437         return "gMin stepping";
00438     }
00439     else if  (convHelper == CONV_SourceStepping)
00440     {
00441         return "source stepping";
00442     }
00443     return "none";
00444 }
00445 
00446 /* This is the non-linear iterative nodal analysis netlist solver. */
00447 template <class nr_type_t>
00448 int nasolver<nr_type_t>::solve_nonlinear (void)
00449 {
00450     qucs::exception * e;
00451     int convergence, run = 0, MaxIterations, error = 0;
00452 
00453     // fetch simulation properties
00454     MaxIterations = getPropertyInteger ("MaxIter");
00455     reltol = getPropertyDouble ("reltol");
00456     abstol = getPropertyDouble ("abstol");
00457     vntol = getPropertyDouble ("vntol");
00458     updateMatrix = 1;
00459 
00460     if (convHelper == CONV_GMinStepping)
00461     {
00462         // use the alternative non-linear solver solve_nonlinear_continuation_gMin
00463         // instead of the basic solver provided by this function
00464         iterations = 0;
00465         error = solve_nonlinear_continuation_gMin ();
00466         return error;
00467     }
00468     else if (convHelper == CONV_SourceStepping)
00469     {
00470         // use the alternative non-linear solver solve_nonlinear_continuation_Source
00471         // instead of the basic solver provided by this function
00472         iterations = 0;
00473         error = solve_nonlinear_continuation_Source ();
00474         return error;
00475     }
00476 
00477     // run solving loop until convergence is reached
00478     do
00479     {
00480         error = solve_once ();
00481         if (!error)
00482         {
00483             // convergence check
00484             convergence = (run > 0) ? checkConvergence () : 0;
00485             savePreviousIteration ();
00486             run++;
00487             // control fixpoint iterations
00488             if (fixpoint)
00489             {
00490                 if (convergence && !updateMatrix)
00491                 {
00492                     updateMatrix = 1;
00493                     convergence = 0;
00494                 }
00495                 else
00496                 {
00497                     updateMatrix = 0;
00498                 }
00499             }
00500         }
00501         else
00502         {
00503             break;
00504         }
00505     }
00506     while (!convergence &&
00507             run < MaxIterations * (1 + convHelper ? 1 : 0));
00508 
00509     if (run >= MaxIterations || error)
00510     {
00511         e = new qucs::exception (EXCEPTION_NO_CONVERGENCE);
00512         e->setText ("no convergence in %s analysis after %d iterations",
00513                     desc.c_str(), run);
00514         throw_exception (e);
00515         error++;
00516     }
00517 
00518     iterations = run;
00519     return error;
00520 }
00521 
00522 /* This is the linear nodal analysis netlist solver. */
00523 template <class nr_type_t>
00524 int nasolver<nr_type_t>::solve_linear (void)
00525 {
00526     updateMatrix = 1;
00527     return solve_once ();
00528 }
00529 
00530 /* Applying the MNA (Modified Nodal Analysis) to a circuit with
00531    passive elements and independent current and voltage sources
00532    results in a matrix equation of the form Ax = z.  This function
00533    generates the A and z matrix. */
00534 template <class nr_type_t>
00535 void nasolver<nr_type_t>::createMatrix (void)
00536 {
00537 
00538     /* Generate the A matrix.  The A matrix consists of four (4) minor
00539        matrices in the form     +-   -+
00540                             A = | G B |
00541                                 | C D |
00542                       +-   -+.
00543        Each of these minor matrices is going to be generated here. */
00544     if (updateMatrix)
00545     {
00546         createGMatrix ();
00547         createBMatrix ();
00548         createCMatrix ();
00549         createDMatrix ();
00550     }
00551 
00552     /* Adjust G matrix if requested. */
00553     if (convHelper == CONV_GMinStepping)
00554     {
00555         int N = countNodes ();
00556         int M = countVoltageSources ();
00557         for (int n = 0; n < N + M; n++)
00558         {
00559             A->set (n, n, A->get (n, n) + gMin);
00560         }
00561     }
00562 
00563     /* Generate the z Matrix.  The z Matrix consists of two (2) minor
00564        matrices in the form     +- -+
00565                             z = | i |
00566                                 | e |
00567                       +- -+.
00568        Each of these minor matrices is going to be generated here. */
00569     createZVector ();
00570 }
00571 
00572 /* This MatVal() functionality is just helper to get the correct
00573    values from the circuit's matrices.  The additional (unused)
00574    argument is used to differentiate between the two possible
00575    types. */
00576 #define MatVal(x) MatValX (x, (nr_type_t *) 0)
00577 
00578 template <class nr_type_t>
00579 nr_type_t nasolver<nr_type_t>::MatValX (nr_complex_t z, nr_complex_t *)
00580 {
00581     return z;
00582 }
00583 
00584 template <class nr_type_t>
00585 nr_type_t nasolver<nr_type_t>::MatValX (nr_complex_t z, nr_double_t *)
00586 {
00587     return real (z);
00588 }
00589 
00590 /* The B matrix is an MxN matrix with only 0, 1 and -1 elements.  Each
00591    location in the matrix corresponds to a particular voltage source
00592    (first dimension) or a node (second dimension).  If the positive
00593    terminal of the ith voltage source is connected to node k, then the
00594    element (i,k) in the B matrix is a 1.  If the negative terminal of
00595    the ith voltage source is connected to node k, then the element
00596    (i,k) in the B matrix is a -1.  Otherwise, elements of the B matrix
00597    are zero. */
00598 template <class nr_type_t>
00599 void nasolver<nr_type_t>::createBMatrix (void)
00600 {
00601     int N = countNodes ();
00602     int M = countVoltageSources ();
00603     circuit * vs;
00604     struct nodelist_t * n;
00605     nr_type_t val;
00606 
00607     // go through each voltage sources (first dimension)
00608     for (int c = 0; c < M; c++)
00609     {
00610         vs = findVoltageSource (c);
00611         // go through each node (second dimension)
00612         for (int r = 0; r < N; r++)
00613         {
00614             val = 0.0;
00615             n = nlist->getNode (r);
00616             for (auto &current : *n)
00617             {
00618                 // is voltage source connected to node ?
00619               if (current->getCircuit () == vs)
00620                 {
00621                   val += MatVal (vs->getB (current->getPort (), c));
00622                 }
00623             }
00624             // put value into B matrix
00625             A->set (r, c + N, val);
00626         }
00627     }
00628 }
00629 
00630 /* The C matrix is an NxM matrix with only 0, 1 and -1 elements.  Each
00631    location in the matrix corresponds to a particular node (first
00632    dimension) or a voltage source (first dimension).  If the positive
00633    terminal of the ith voltage source is connected to node k, then the
00634    element (k,i) in the C matrix is a 1.  If the negative terminal of
00635    the ith voltage source is connected to node k, then the element
00636    (k,i) in the C matrix is a -1.  Otherwise, elements of the C matrix
00637    are zero. */
00638 template <class nr_type_t>
00639 void nasolver<nr_type_t>::createCMatrix (void)
00640 {
00641     int N = countNodes ();
00642     int M = countVoltageSources ();
00643     circuit * vs;
00644     struct nodelist_t * n;
00645     nr_type_t val;
00646 
00647     // go through each voltage sources (second dimension)
00648     for (int r = 0; r < M; r++)
00649     {
00650         vs = findVoltageSource (r);
00651         // go through each node (first dimension)
00652         for (int c = 0; c < N; c++)
00653         {
00654             val = 0.0;
00655             n = nlist->getNode (c);
00656             for (auto &current: *n)
00657             {
00658                 // is voltage source connected to node ?
00659               if (current->getCircuit () == vs)
00660                 {
00661                   val += MatVal (vs->getC (r, current->getPort ()));
00662                 }
00663             }
00664             // put value into C matrix
00665             A->set (r + N, c, val);
00666         }
00667     }
00668 }
00669 
00670 /* The D matrix is an MxM matrix that is composed entirely of zeros.
00671    It can be non-zero if dependent sources are considered. */
00672 template <class nr_type_t>
00673 void nasolver<nr_type_t>::createDMatrix (void)
00674 {
00675     int M = countVoltageSources ();
00676     int N = countNodes ();
00677     circuit * vsr, * vsc;
00678     nr_type_t val;
00679     for (int r = 0; r < M; r++)
00680     {
00681         vsr = findVoltageSource (r);
00682         for (int c = 0; c < M; c++)
00683         {
00684             vsc = findVoltageSource (c);
00685             val = 0.0;
00686             if (vsr == vsc)
00687             {
00688                 val = MatVal (vsr->getD (r, c));
00689             }
00690             A->set (r + N, c + N, val);
00691         }
00692     }
00693 }
00694 
00695 /* The G matrix is an NxN matrix formed in two steps.
00696    1. Each element in the diagonal matrix is equal to the sum of the
00697    conductance of each element connected to the corresponding node.
00698    2. The off diagonal elements are the negative conductance of the
00699    element connected to the pair of corresponding nodes.  Therefore a
00700    resistor between nodes 1 and 2 goes into the G matrix at location
00701    (1,2) and location (2,1).  If an element is grounded, it will only
00702    have contribute to one entry in the G matrix -- at the appropriate
00703    location on the diagonal. */
00704 template <class nr_type_t>
00705 void nasolver<nr_type_t>::createGMatrix (void)
00706 {
00707     int pr, pc, N = countNodes ();
00708     nr_type_t g;
00709     struct nodelist_t * nr, * nc;
00710     circuit * ct;
00711 
00712     // go through each column of the G matrix
00713     for (int c = 0; c < N; c++)
00714     {
00715         nc = nlist->getNode (c);
00716         // go through each row of the G matrix
00717         for (int r = 0; r < N; r++)
00718         {
00719             nr = nlist->getNode (r);
00720             g = 0.0;
00721             // sum up the conductance of each connected circuit
00722             for (auto & currentnc  : *nc)
00723               for (auto & currentnr: *nr)
00724                 if (currentnc->getCircuit () == currentnr->getCircuit ())
00725                   {
00726                     ct = currentnc->getCircuit ();
00727                     pc = currentnc->getPort ();
00728                     pr = currentnr->getPort ();
00729                     g += MatVal (ct->getY (pr, pc));
00730                   }
00731             // put value into G matrix
00732             A->set (r, c, g);
00733         }
00734     }
00735 }
00736 
00737 /* The following function creates the (N+M)x(N+M) noise current
00738    correlation matrix used during the AC noise computations.  */
00739 template <class nr_type_t>
00740 void nasolver<nr_type_t>::createNoiseMatrix (void)
00741 {
00742     int pr, pc, N = countNodes ();
00743     int M = countVoltageSources ();
00744     struct nodelist_t * n;
00745     nr_type_t val;
00746     int r, c, ri, ci;
00747     struct nodelist_t * nr, * nc;
00748     circuit * ct;
00749 
00750     // create new Cy matrix if necessary
00751     delete C;
00752     C = new tmatrix<nr_type_t> (N + M);
00753 
00754     // go through each column of the Cy matrix
00755     for (c = 0; c < N; c++)
00756     {
00757         nc = nlist->getNode (c);
00758         // go through each row of the Cy matrix
00759         for (r = 0; r < N; r++)
00760         {
00761             nr = nlist->getNode (r);
00762             val = 0.0;
00763             // sum up the noise-correlation of each connected circuit
00764             for (auto & currentnc: *nc)
00765                 /* a = 0; a < nc->size(); a++ */
00766               for (auto &currentnr : *nr)
00767                 /* b = 0; b < nr->size(); b++) */
00768                     if (currentnc->getCircuit () == currentnr->getCircuit ())
00769                     {
00770                         ct = currentnc->getCircuit ();
00771                         pc = currentnc->getPort ();
00772                         pr = currentnr->getPort ();
00773                         val += MatVal (ct->getN (pr, pc));
00774                     }
00775             // put value into Cy matrix
00776             C->set (r, c, val);
00777         }
00778     }
00779 
00780     // go through each additional voltage source and put coefficients into
00781     // the noise current correlation matrix
00782     circuit * vsr, * vsc;
00783     for (r = 0; r < M; r++)
00784     {
00785         vsr = findVoltageSource (r);
00786         for (c = 0; c < M; c++)
00787         {
00788             vsc = findVoltageSource (c);
00789             val = 0.0;
00790             if (vsr == vsc)
00791             {
00792                 ri = vsr->getSize () + r - vsr->getVoltageSource ();
00793                 ci = vsc->getSize () + c - vsc->getVoltageSource ();
00794                 val = MatVal (vsr->getN (ri, ci));
00795             }
00796             C->set (r + N, c + N, val);
00797         }
00798     }
00799 
00800     // go through each additional voltage source
00801     for (r = 0; r < M; r++)
00802     {
00803         vsr = findVoltageSource (r);
00804         // go through each node
00805         for (c = 0; c < N; c++)
00806         {
00807             val = 0.0;
00808             n = nlist->getNode (c);
00809             for (auto &currentn: *n)
00810               /*i = 0; i < n->size(); i++ )*/
00811             {
00812                 // is voltage source connected to node ?
00813                 if (currentn->getCircuit () == vsr)
00814                 {
00815                     ri = vsr->getSize () + r - vsr->getVoltageSource ();
00816                     ci = currentn->getPort ();
00817                     val += MatVal (vsr->getN (ri, ci));
00818                 }
00819             }
00820             // put value into Cy matrix
00821             C->set (r + N, c, val);
00822         }
00823     }
00824 
00825     // go through each voltage source
00826     for (c = 0; c < M; c++)
00827     {
00828         vsc = findVoltageSource (c);
00829         // go through each node
00830         for (r = 0; r < N; r++)
00831         {
00832             val = 0.0;
00833             n = nlist->getNode (r);
00834             for (auto & currentn: *n)/*i = 0; i < n->size(); i++)*/
00835             {
00836                 // is voltage source connected to node ?
00837                 if (currentn->getCircuit () == vsc)
00838                 {
00839                     ci = vsc->getSize () + c - vsc->getVoltageSource ();
00840                     ri = currentn->getPort ();
00841                     val += MatVal (vsc->getN (ri, ci));
00842                 }
00843             }
00844             // put value into Cy matrix
00845             C->set (r, c + N, val);
00846         }
00847     }
00848 
00849 }
00850 
00851 /* The i matrix is an 1xN matrix with each element of the matrix
00852    corresponding to a particular node.  The value of each element of i
00853    is determined by the sum of current sources into the corresponding
00854    node.  If there are no current sources connected to the node, the
00855    value is zero. */
00856 template <class nr_type_t>
00857 void nasolver<nr_type_t>::createIVector (void)
00858 {
00859     int N = countNodes ();
00860     nr_type_t val;
00861     struct nodelist_t * n;
00862     circuit * is;
00863 
00864     // go through each node
00865     for (int r = 0; r < N; r++)
00866     {
00867         val = 0.0;
00868         n = nlist->getNode (r);
00869         // go through each circuit connected to the node
00870         for (auto &currentn: *n)/* int i = 0; i < n->size(); i++)*/
00871         {
00872           is = currentn->getCircuit ();
00873           // is this a current source ?
00874           if (is->isISource () || is->isNonLinear ())
00875             {
00876               val += MatVal (is->getI (currentn->getPort ()));
00877             }
00878         }
00879         // put value into i vector
00880         z->set (r, val);
00881     }
00882 }
00883 
00884 /* The e matrix is a 1xM matrix with each element of the matrix equal
00885    in value to the corresponding independent voltage source. */
00886 template <class nr_type_t>
00887 void nasolver<nr_type_t>::createEVector (void)
00888 {
00889     int N = countNodes ();
00890     int M = countVoltageSources ();
00891     nr_type_t val;
00892     circuit * vs;
00893 
00894     // go through each voltage source
00895     for (int r = 0; r < M; r++)
00896     {
00897         vs = findVoltageSource (r);
00898         val = MatVal (vs->getE (r));
00899         // put value into e vector
00900         z->set (r + N, val);
00901     }
00902 }
00903 
00904 // The function loads the right hand side vector.
00905 template <class nr_type_t>
00906 void nasolver<nr_type_t>::createZVector (void)
00907 {
00908     createIVector ();
00909     createEVector ();
00910 }
00911 
00912 // Returns the number of nodes in the nodelist, excluding the ground node.
00913 template <class nr_type_t>
00914 int nasolver<nr_type_t>::countNodes (void)
00915 {
00916     return nlist->length () - 1;
00917 }
00918 
00919 // Returns the node number of the give node name.
00920 template <class nr_type_t>
00921 int nasolver<nr_type_t>::getNodeNr (const std::string &str)
00922 {
00923     return nlist->getNodeNr (str);
00924 }
00925 
00926 /* The function returns the assigned node number for the port of the
00927    given circuits.  It returns -1 if there is no such node. */
00928 template <class nr_type_t>
00929 int nasolver<nr_type_t>::findAssignedNode (circuit * c, int port)
00930 {
00931     int N = countNodes ();
00932     struct nodelist_t * n;
00933     for (int r = 0; r < N; r++)
00934     {
00935         n = nlist->getNode (r);
00936         for (auto &currentn : *n) /*int i = 0; i < n->size(); i++)*/
00937             if (c == currentn->getCircuit ())
00938                 if (port == currentn->getPort ())
00939                     return r;
00940     }
00941     return -1;
00942 }
00943 
00944 // Returns the number of voltage sources in the nodelist.
00945 template <class nr_type_t>
00946 int nasolver<nr_type_t>::countVoltageSources (void)
00947 {
00948     return subnet->getVoltageSources ();
00949 }
00950 
00951 /* The function returns the voltage source circuit object
00952    corresponding to the given number.  If there is no such voltage
00953    source it returns NULL. */
00954 template <class nr_type_t>
00955 circuit * nasolver<nr_type_t>::findVoltageSource (int n)
00956 {
00957     circuit * root = subnet->getRoot ();
00958     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00959     {
00960         if (n >= c->getVoltageSource () &&
00961                 n <= c->getVoltageSource () + c->getVoltageSources () - 1)
00962             return c;
00963     }
00964     return NULL;
00965 }
00966 
00967 /* The function applies unique voltage source identifiers to each
00968    voltage source (explicit and built in internal ones) in the list of
00969    registered circuits. */
00970 template <class nr_type_t>
00971 void nasolver<nr_type_t>::assignVoltageSources (void)
00972 {
00973     circuit * root = subnet->getRoot ();
00974     int nSources = 0;
00975     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00976     {
00977         if (c->getVoltageSources () > 0)
00978         {
00979             c->setVoltageSource (nSources);
00980             nSources += c->getVoltageSources ();
00981         }
00982     }
00983     subnet->setVoltageSources (nSources);
00984 }
00985 
00986 /* The matrix equation Ax = z is solved by x = A^-1*z.  The function
00987    applies the operation to the previously generated matrices. */
00988 template <class nr_type_t>
00989 void nasolver<nr_type_t>::runMNA (void)
00990 {
00991 
00992     // just solve the equation system here
00993     eqns->setAlgo (eqnAlgo);
00994     eqns->passEquationSys (updateMatrix ? A : NULL, x, z);
00995     eqns->solve ();
00996 
00997     // if damped Newton-Raphson is requested
00998     if (xprev != NULL && top_exception () == NULL)
00999     {
01000         if (convHelper == CONV_Attenuation)
01001         {
01002             applyAttenuation ();
01003         }
01004         else if (convHelper == CONV_LineSearch)
01005         {
01006             lineSearch ();
01007         }
01008         else if (convHelper == CONV_SteepestDescent)
01009         {
01010             steepestDescent ();
01011         }
01012     }
01013 }
01014 
01015 /* This function applies a damped Newton-Raphson (limiting scheme) to
01016    the current solution vector in the form x1 = x0 + a * (x1 - x0).  This
01017    convergence helper is heuristic and does not ensure global convergence. */
01018 template <class nr_type_t>
01019 void nasolver<nr_type_t>::applyAttenuation (void)
01020 {
01021     nr_double_t alpha = 1.0, nMax;
01022 
01023     // create solution difference vector and find maximum deviation
01024     tvector<nr_type_t> dx = *x - *xprev;
01025     nMax = maxnorm (dx);
01026 
01027     // compute appropriate damping factor
01028     if (nMax > 0.0)
01029     {
01030         nr_double_t g = 1.0;
01031         alpha = std::min (0.9, g / nMax);
01032         if (alpha < 0.1) alpha = 0.1;
01033     }
01034 
01035     // apply damped solution vector
01036     *x = *xprev + alpha * dx;
01037 }
01038 
01039 /* This is damped Newton-Raphson using nested iterations in order to
01040    find a better damping factor.  It identifies a damping factor in
01041    the interval [0,1] which minimizes the right hand side vector.  The
01042    algorithm actually ensures global convergence but pushes the
01043    solution to local minimums, i.e. where the Jacobian matrix A may be
01044    singular. */
01045 template <class nr_type_t>
01046 void nasolver<nr_type_t>::lineSearch (void)
01047 {
01048     nr_double_t alpha = 0.5, n, nMin, aprev = 1.0, astep = 0.5, adiff;
01049     int dir = -1;
01050 
01051     // compute solution deviation vector
01052     tvector<nr_type_t> dx = *x - *xprev;
01053     nMin = std::numeric_limits<nr_double_t>::max();
01054 
01055     do
01056     {
01057         // apply current damping factor and see what happens
01058         *x = *xprev + alpha * dx;
01059 
01060         // recalculate Jacobian and right hand side
01061         saveSolution ();
01062         calculate ();
01063         createZVector ();
01064 
01065         // calculate norm of right hand side vector
01066         n = norm (*z);
01067 
01068         // TODO: this is not perfect, but usable
01069         astep /= 2;
01070         adiff = fabs (alpha - aprev);
01071         if (adiff > 0.005)
01072         {
01073             aprev = alpha;
01074             if (n < nMin)
01075             {
01076                 nMin = n;
01077                 if (alpha == 1) dir = -dir;
01078                 alpha += astep * dir;
01079             }
01080             else
01081             {
01082                 dir = -dir;
01083                 alpha += 1.5 * astep * dir;
01084             }
01085         }
01086     }
01087     while (adiff > 0.005);
01088 
01089     // apply final damping factor
01090     assert (alpha > 0 && alpha <= 1);
01091     *x = *xprev + alpha * dx;
01092 }
01093 
01094 /* The function looks for the optimal gradient for the right hand side
01095    vector using the so-called 'steepest descent' method.  Though
01096    better than the one-dimensional linesearch (it doesn't push
01097    iterations into local minimums) it converges painfully slow. */
01098 template <class nr_type_t>
01099 void nasolver<nr_type_t>::steepestDescent (void)
01100 {
01101     nr_double_t alpha = 1.0, sl, n;
01102 
01103     // compute solution deviation vector
01104     tvector<nr_type_t> dx = *x - *xprev;
01105     tvector<nr_type_t> dz = *z - *zprev;
01106     n = norm (*zprev);
01107 
01108     do
01109     {
01110         // apply current damping factor and see what happens
01111         *x = *xprev + alpha * dx;
01112 
01113         // recalculate Jacobian and right hand side
01114         saveSolution ();
01115         calculate ();
01116         createZVector ();
01117 
01118         // check gradient criteria, ThinkME: Is this correct?
01119         dz = *z - *zprev;
01120         sl = real (sum (dz * -dz));
01121         if (norm (*z) < n + alpha * sl) break;
01122         alpha *= 0.7;
01123     }
01124     while (alpha > 0.001);
01125 
01126     // apply final damping factor
01127     *x = *xprev + alpha * dx;
01128 }
01129 
01130 /* The function checks whether the iterative algorithm for linearizing
01131    the non-linear components in the network shows convergence.  It
01132    returns non-zero if it converges and zero otherwise. */
01133 template <class nr_type_t>
01134 int nasolver<nr_type_t>::checkConvergence (void)
01135 {
01136 
01137     int N = countNodes ();
01138     int M = countVoltageSources ();
01139     nr_double_t v_abs, v_rel, i_abs, i_rel;
01140     int r;
01141 
01142     // check the nodal voltage changes against the allowed absolute
01143     // and relative tolerance values
01144     for (r = 0; r < N; r++)
01145     {
01146         v_abs = abs (x->get (r) - xprev->get (r));
01147         v_rel = abs (x->get (r));
01148         if (v_abs >= vntol + reltol * v_rel) return 0;
01149         if (!convHelper)
01150         {
01151             i_abs = abs (z->get (r) - zprev->get (r));
01152             i_rel = abs (z->get (r));
01153             if (i_abs >= abstol + reltol * i_rel) return 0;
01154         }
01155     }
01156 
01157     for (r = 0; r < M; r++)
01158     {
01159         i_abs = abs (x->get (r + N) - xprev->get (r + N));
01160         i_rel = abs (x->get (r + N));
01161         if (i_abs >= abstol + reltol * i_rel) return 0;
01162         if (!convHelper)
01163         {
01164             v_abs = abs (z->get (r + N) - zprev->get (r + N));
01165             v_rel = abs (z->get (r + N));
01166             if (v_abs >= vntol + reltol * v_rel) return 0;
01167         }
01168     }
01169     return 1;
01170 }
01171 
01172 /* The function saves the solution and right hand vector of the previous
01173    iteration. */
01174 template <class nr_type_t>
01175 void nasolver<nr_type_t>::savePreviousIteration (void)
01176 {
01177     if (xprev != NULL)
01178         *xprev = *x;
01179     else
01180         xprev = new tvector<nr_type_t> (*x);
01181     if (zprev != NULL)
01182         *zprev = *z;
01183     else
01184         zprev = new tvector<nr_type_t> (*z);
01185 }
01186 
01187 /* The function restores the solution and right hand vector of the
01188    previous (successful) iteration. */
01189 template <class nr_type_t>
01190 void nasolver<nr_type_t>::restorePreviousIteration (void)
01191 {
01192     if (xprev != NULL) *x = *xprev;
01193     if (zprev != NULL) *z = *zprev;
01194 }
01195 
01196 /* The function restarts the NR iteration for each non-linear
01197    circuit. */
01198 template <class nr_type_t>
01199 void nasolver<nr_type_t>::restartNR (void)
01200 {
01201     circuit * root = subnet->getRoot ();
01202     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
01203     {
01204         if (c->isNonLinear ()) c->restartDC ();
01205     }
01206 }
01207 
01208 /* This function goes through solution (the x vector) and saves the
01209    node voltages of the last iteration into each non-linear
01210    circuit. */
01211 template <class nr_type_t>
01212 void nasolver<nr_type_t>::saveNodeVoltages (void)
01213 {
01214     int N = countNodes ();
01215     struct nodelist_t * n;
01216     // save all nodes except reference node
01217     for (int r = 0; r < N; r++)
01218     {
01219         n = nlist->getNode (r);
01220         /* for (int i = 0; i < n->size(); i++)*/
01221         for(auto &currentn: *n)
01222         {
01223           currentn->getCircuit()->setV (currentn->getPort (), x->get (r));
01224         }
01225     }
01226     // save reference node
01227     n = nlist->getNode (-1);
01228     for(auto &currentn: *n)
01229       currentn->getCircuit()->setV (currentn->getPort (), 0.0);
01230 }
01231 
01232 /* This function goes through solution (the x vector) and saves the
01233    branch currents through the voltage sources of the last iteration
01234    into each circuit. */
01235 template <class nr_type_t>
01236 void nasolver<nr_type_t>::saveBranchCurrents (void)
01237 {
01238     int N = countNodes ();
01239     int M = countVoltageSources ();
01240     circuit * vs;
01241     // save all branch currents of voltage sources
01242     for (int r = 0; r < M; r++)
01243     {
01244         vs = findVoltageSource (r);
01245         vs->setJ (r, x->get (r + N));
01246     }
01247 }
01248 
01249 // The function saves the solution vector into each circuit.
01250 template <class nr_type_t>
01251 void nasolver<nr_type_t>::saveSolution (void)
01252 {
01253     saveNodeVoltages ();
01254     saveBranchCurrents ();
01255 }
01256 
01257 // This function stores the solution (node voltages and branch currents).
01258 template <class nr_type_t>
01259 void nasolver<nr_type_t>::storeSolution (void)
01260 {
01261     // cleanup solution previously
01262     solution.clear ();
01263     int r;
01264     int N = countNodes ();
01265     int M = countVoltageSources ();
01266     // store all nodes except reference node
01267     for (r = 0; r < N; r++)
01268     {
01269         struct nodelist_t * n = nlist->getNode (r);
01270         nr_type_t gr = x->get (r);
01271         qucs::naentry<nr_type_t> entry(gr, 0);
01272         solution.insert({{n->name, entry }});
01273     }
01274     // store all branch currents of voltage sources
01275     for (r = 0; r < M; r++)
01276     {
01277         circuit * vs = findVoltageSource (r);
01278         int vn = r - vs->getVoltageSource () + 1;
01279         nr_type_t xg = x->get (r + N);
01280         qucs::naentry<nr_type_t> entry(xg, vn);
01281         solution.insert({{vs->getName (), entry}});
01282     }
01283 }
01284 
01285 // This function recalls the solution (node voltages and branch currents).
01286 template <class nr_type_t>
01287 void nasolver<nr_type_t>::recallSolution (void)
01288 {
01289     int r;
01290     int N = countNodes ();
01291     int M = countVoltageSources ();
01292     // store all nodes except reference node
01293     for (r = 0; r < N; r++)
01294     {
01295         struct nodelist_t * n = nlist->getNode (r);
01296         auto na = solution.find(n->name);
01297         if (na != solution.end())
01298           if ((*na).second.current == 0)
01299             x->set (r, (*na).second.value);
01300     }
01301     // store all branch currents of voltage sources
01302     for (r = 0; r < M; r++)
01303     {
01304         circuit * vs = findVoltageSource (r);
01305         int vn = r - vs->getVoltageSource () + 1;
01306         auto na = solution.find(vs->getName ());
01307         if (na != solution.end())
01308           if ((*na).second.current == vn)
01309             x->set (r + N, (*na).second.value);
01310     }
01311 }
01312 
01313 /* This function saves the results of a single solve() functionality
01314    into the output dataset. */
01315 template <class nr_type_t>
01316 void nasolver<nr_type_t>::saveResults (const std::string &volts, const std::string &amps,
01317                                        int saveOPs, qucs::vector * f)
01318 {
01319     int N = countNodes ();
01320     int M = countVoltageSources ();
01321 
01322     // add node voltage variables
01323     if (!volts.empty())
01324     {
01325         for (int r = 0; r < N; r++)
01326         {
01327           std::string n = createV (r, volts, saveOPs);
01328           if(!n.empty())
01329             {
01330                 saveVariable (n, x->get (r), f);
01331             }
01332         }
01333     }
01334 
01335     // add branch current variables
01336     if (!amps.empty())
01337     {
01338         for (int r = 0; r < M; r++)
01339         {
01340           std::string n = createI (r, amps, saveOPs);
01341           if (!n.empty())
01342             {
01343                 saveVariable (n, x->get (r + N), f);
01344             }
01345         }
01346     }
01347 
01348     // add voltage probe data
01349     if (!volts.empty())
01350     {
01351         circuit * root = subnet->getRoot ();
01352         for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
01353         {
01354             if (!c->isProbe ()) continue;
01355             if (!c->getSubcircuit().empty() && !(saveOPs & SAVE_ALL)) continue;
01356             if (volts != "vn")
01357                 c->saveOperatingPoints ();          
01358             std::string n = createOP (c->getName (), volts);
01359             saveVariable (n, nr_complex_t (c->getOperatingPoint ("Vr"),
01360                                    c->getOperatingPoint ("Vi")), f);
01361         }
01362     }
01363 
01364     // save operating points of non-linear circuits if requested
01365     if (saveOPs & SAVE_OPS)
01366     {
01367         circuit * root = subnet->getRoot ();
01368         for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
01369         {
01370             if (!c->isNonLinear ()) continue;
01371             if (!c->getSubcircuit ().empty() && !(saveOPs & SAVE_ALL)) continue;
01372             c->calcOperatingPoints ();
01373             for (auto ops: c->getOperatingPoints ())
01374             {
01375                 operatingpoint &p = ops.second;
01376                 std::string n = createOP (c->getName (), p.getName ());
01377                 saveVariable (n, p.getValue (), f);
01378             }
01379         }
01380     }
01381 }
01382 
01383 /* Create an appropriate variable name for operating points.  The
01384    caller is responsible to free() the returned string. */
01385 template <class nr_type_t>
01386 std::string nasolver<nr_type_t>::createOP (const std::string &c, const std::string &n)
01387 {
01388     return c+"."+n;
01389 }
01390 
01391 /* Creates an appropriate variable name for voltages.  The caller is
01392    responsible to free() the returned string. */
01393 template <class nr_type_t>
01394 std::string nasolver<nr_type_t>::createV (int n, const std::string &volts, int saveOPs)
01395 {
01396     if (nlist->isInternal (n))
01397       return std::string();
01398     std::string node = nlist->get (n);
01399     if(node.find('.')!=std::string::npos && !(saveOPs & SAVE_ALL))
01400       return std::string();
01401     std::string ret = node+"."+volts;
01402     return ret;
01403 }
01404 
01405 /* Create an appropriate variable name for currents.  The caller is
01406    responsible to free() the returned string. */
01407 template <class nr_type_t>
01408 std::string nasolver<nr_type_t>::createI (int n, const std::string &amps, int saveOPs)
01409 {
01410     circuit * vs = findVoltageSource (n);
01411 
01412     // don't output internal (helper) voltage sources
01413     if (vs->isInternalVoltageSource ())
01414       return std::string();
01415 
01416     /* save only current through real voltage sources and explicit
01417        current probes */
01418     if (!vs->isVSource () && !(saveOPs & SAVE_OPS))
01419       return std::string();
01420 
01421     // don't output subcircuit components if not requested
01422     if (!vs->getSubcircuit ().empty() && !(saveOPs & SAVE_ALL))
01423       return std::string();
01424 
01425     // create appropriate current name for single/multiple voltage sources
01426     std::string name = vs->getName ();
01427     if (vs->getVoltageSources () > 1)
01428       return name+"."+amps+std::to_string(n - vs->getVoltageSource () + 1);
01429     else
01430       return name+"."+amps;
01431 }
01432 
01433 
01434 /* Alternaive to countNodes () */
01435 template <class nr_type_t>
01436 int nasolver<nr_type_t>::getN()
01437 {
01438     return countNodes ();
01439 }
01440 
01441 /* Alternative to countVoltageSources () */
01442 template <class nr_type_t>
01443 int nasolver<nr_type_t>::getM()
01444 {
01445     return countVoltageSources ();
01446 }
01447 
01448 } // namespace qucs