Qucs-core  0.0.19
hbsolver.cpp
Go to the documentation of this file.
00001 /*
00002  * hbsolver.cpp - harmonic balance solver class implementation
00003  *
00004  * Copyright (C) 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 #if HAVE_CONFIG_H
00026 # include <config.h>
00027 #endif
00028 
00029 #include<algorithm>
00030 
00031 #include <stdio.h>
00032 
00033 #include "object.h"
00034 #include "logging.h"
00035 #include "complex.h"
00036 #include "vector.h"
00037 #include "node.h"
00038 #include "circuit.h"
00039 #include "component_id.h"
00040 #include "net.h"
00041 #include "netdefs.h"
00042 #include "strlist.h"
00043 #include "ptrlist.h"
00044 #include "tvector.h"
00045 #include "tmatrix.h"
00046 #include "eqnsys.h"
00047 #include "analysis.h"
00048 #include "dataset.h"
00049 #include "fourier.h"
00050 #include "hbsolver.h"
00051 
00052 #define HB_DEBUG 0
00053 
00054 namespace qucs {
00055 
00056 using namespace fourier;
00057 
00058 // Constructor creates an unnamed instance of the hbsolver class.
00059 hbsolver::hbsolver () : analysis () {
00060   type = ANALYSIS_HBALANCE;
00061   frequency = 0;
00062   nlnodes = lnnodes = banodes = nanodes = NULL;
00063   Y = Z = A = NULL;
00064   NA = YV = JQ = JG = JF = NULL;
00065   OM = IR = QR = RH = IG = FQ = VS = VP = FV = IL = IN = IC = IS = NULL;
00066   vs = x = NULL;
00067   runs = 0;
00068   ndfreqs = NULL;
00069 }
00070 
00071 // Constructor creates a named instance of the hbsolver class.
00072 hbsolver::hbsolver (char * n) : analysis (n) {
00073   type = ANALYSIS_HBALANCE;
00074   frequency = 0;
00075   nlnodes = lnnodes = banodes = nanodes = NULL;
00076   Y = Z = A = NULL;
00077   NA = YV = JQ = JG = JF = NULL;
00078   OM = IR = QR = RH = IG = FQ = VS = VP = FV = IL = IN = IC = IS = NULL;
00079   vs = x = NULL;
00080   runs = 0;
00081   ndfreqs = NULL;
00082 }
00083 
00084 // Destructor deletes the hbsolver class object.
00085 hbsolver::~hbsolver () {
00086   // delete nodelists
00087   delete nlnodes;
00088   delete lnnodes;
00089   delete banodes;
00090   delete nanodes;
00091 
00092   // delete temporary matrices
00093   delete A;
00094   delete Z;
00095   delete Y;
00096 
00097   // delete matrices
00098   delete NA;
00099   delete YV;
00100   delete JQ;
00101   delete JG;
00102   delete JF;
00103 
00104   // delete vectors
00105   delete IC;
00106   delete IS;
00107   delete FV;
00108   delete IL;
00109   delete IN;
00110   delete IG;
00111   delete FQ;
00112   delete VS;
00113   delete VP;
00114   delete vs;
00115   delete OM;
00116   delete IR;
00117   delete QR;
00118   delete RH;
00119 
00120   delete x;
00121   delete[] ndfreqs;
00122 }
00123 
00124 /* The copy constructor creates a new instance of the hbsolver class
00125    based on the given hbsolver object. */
00126 hbsolver::hbsolver (hbsolver & o) : analysis (o) {
00127   frequency = o.frequency;
00128   negfreqs = o.negfreqs;
00129   posfreqs = o.posfreqs;
00130   nlnodes = o.nlnodes;
00131   lnnodes = o.lnnodes;
00132   banodes = o.banodes;
00133   nanodes = o.nanodes;
00134   Y = Z = A = NULL;
00135   NA = YV = JQ = JG = JF = NULL;
00136   OM = IR = QR = RH = IG = FQ = VS = VP = FV = IL = IN = IC = IS = NULL;
00137   vs = x = NULL;
00138   runs = o.runs;
00139   ndfreqs = NULL;
00140 }
00141 
00142 #define VS_(r) (*VS) (r)
00143 #define OM_(r) (*OM) (r)
00144 
00145 /* This is the HB netlist solver.  It prepares the circuit list and
00146    solves it then. */
00147 int hbsolver::solve (void) {
00148 
00149   int iterations = 0, done = 0;
00150   int MaxIterations = getPropertyInteger ("MaxIter");
00151 
00152   // collect different parts of the circuit
00153   splitCircuits ();
00154 
00155   // create frequency array
00156   collectFrequencies ();
00157 
00158   // find interconnects between the linear and non-linear subcircuit
00159   getNodeLists ();
00160 
00161   // prepares the linear part --> 0 = IC + [YV] * VS
00162   prepareLinear ();
00163 
00164   runs++;
00165   logprint (LOG_STATUS, "NOTIFY: %s: solving for %d frequencies\n",
00166             getName (), lnfreqs);
00167 
00168   if (nbanodes > 0) {
00169 
00170     // start balancing
00171     logprint (LOG_STATUS, "NOTIFY: %s: balancing at %d nodes\n", getName (),
00172               nbanodes);
00173 
00174     // prepares the non-linear part
00175     prepareNonLinear ();
00176 
00177 #if HB_DEBUG
00178       fprintf (stderr, "YV -- transY in f:\n"); YV->print ();
00179       fprintf (stderr, "IC -- constant current in f:\n"); IC->print ();
00180 #endif
00181 
00182     // start iteration
00183     do {
00184       iterations++;
00185 
00186 #if HB_DEBUG
00187       fprintf (stderr, "\n   -- iteration step: %d\n", iterations);
00188       fprintf (stderr, "vs -- voltage in t:\n"); vs->print ();
00189 #endif
00190 
00191       // evaluate component functionality and fill matrices and vectors
00192       loadMatrices ();
00193 
00194 #if HB_DEBUG
00195       fprintf (stderr, "FQ -- charge in t:\n"); FQ->print ();
00196       fprintf (stderr, "IG -- current in t:\n"); IG->print ();
00197 #endif
00198 
00199       // currents into frequency domain
00200       VectorFFT (IG);
00201 
00202       // charges into frequency domain
00203       VectorFFT (FQ);
00204 
00205       // right hand side currents and charges into the frequency domain
00206       VectorFFT (IR);
00207       VectorFFT (QR);
00208 
00209 #if HB_DEBUG
00210       fprintf (stderr, "VS -- voltage in f:\n"); VS->print ();
00211       fprintf (stderr, "FQ -- charge in f:\n"); FQ->print ();
00212       fprintf (stderr, "IG -- current in f:\n"); IG->print ();
00213       fprintf (stderr, "IR -- corrected Jacobi current in f:\n"); IR->print ();
00214 #endif
00215 
00216       // solve HB equation --> FV = IC + [YV] * VS + j[O] * FQ + IG
00217       solveHB ();
00218 
00219 #if HB_DEBUG
00220       fprintf (stderr, "FV -- error vector F(V) in f:\n"); FV->print ();
00221       fprintf (stderr, "IL -- linear currents in f:\n"); IL->print ();
00222       fprintf (stderr, "IN -- non-linear currents in f:\n"); IN->print ();
00223       fprintf (stderr, "RH -- right-hand side currents in f:\n"); RH->print ();
00224 #endif
00225 
00226       // termination criteria met
00227       if (iterations > 1 && checkBalance ()) {
00228         done = 1;
00229         break;
00230       }
00231 
00232 #if HB_DEBUG
00233       fprintf (stderr, "JG -- G-Jacobian in t:\n"); JG->print ();
00234       fprintf (stderr, "JQ -- C-Jacobian in t:\n"); JQ->print ();
00235 #endif
00236 
00237       // G-Jacobian into frequency domain
00238       MatrixFFT (JG);
00239 
00240       // C-Jacobian into frequency domain
00241       MatrixFFT (JQ);
00242 
00243 #if HB_DEBUG
00244       fprintf (stderr, "JQ -- dQ/dV C-Jacobian in f:\n"); JQ->print ();
00245       fprintf (stderr, "JG -- dI/dV G-Jacobian in f:\n"); JG->print ();
00246 #endif
00247 
00248       // calculate Jacobian --> JF = [YV] + j[O] * JQ + JG
00249       calcJacobian ();
00250 
00251 #if HB_DEBUG
00252       fprintf (stderr, "JF -- full Jacobian in f:\n"); JF->print ();
00253 #endif
00254 
00255       // solve equation system --> JF * VS(n+1) = JF * VS(n) - FV
00256       solveVoltages ();
00257 
00258 #if HB_DEBUG
00259       fprintf (stderr, "VS -- next voltage in f:\n"); VS->print ();
00260 #endif
00261 
00262       // inverse FFT of frequency domain voltage vector VS(n+1)
00263       VectorIFFT (vs);
00264     }
00265     // check termination criteria (balanced frequency domain currents)
00266     while (!done && iterations < MaxIterations);
00267 
00268     if (iterations >= MaxIterations) {
00269       qucs::exception * e = new qucs::exception (EXCEPTION_NO_CONVERGENCE);
00270       e->setText ("no convergence in %s analysis after %d iterations",
00271                   getName (), iterations);
00272       throw_exception (e);
00273       logprint (LOG_ERROR, "%s: no convergence after %d iterations\n",
00274                 getName (), iterations);
00275     }
00276     else {
00277       logprint (LOG_STATUS, "%s: convergence reached after %d iterations\n",
00278                 getName (), iterations);
00279     }
00280   }
00281   else {
00282     // no balancing necessary
00283     logprint (LOG_STATUS, "NOTIFY: %s: no balancing necessary\n", getName ());
00284   }
00285 
00286   // print exception stack
00287   estack.print ();
00288 
00289   // apply AC analysis to the complete network in order to obtain the
00290   // final results
00291   finalSolution ();
00292 
00293   // save results into dataset
00294   saveResults ();
00295   return 0;
00296 }
00297 
00298 /* Goes through the list of circuit objects and runs its calcHB()
00299    function. */
00300 void hbsolver::calc (hbsolver * self) {
00301   circuit * root = self->getNet()->getRoot ();
00302   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00303     c->calcHB (self->frequency);
00304   }
00305 }
00306 
00307 /* Goes through the list of circuit objects and runs its initHB()
00308    function. */
00309 void hbsolver::initHB (void) {
00310   circuit * root = subnet->getRoot ();
00311   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00312     c->initHB ();
00313   }
00314 }
00315 
00316 /* Goes through the list of circuit objects and runs its initDC()
00317    function. */
00318 void hbsolver::initDC (void) {
00319   circuit * root = subnet->getRoot ();
00320   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00321     c->initDC ();
00322   }
00323 }
00324 
00325 // Returns true if circuit is a HB source.
00326 bool hbsolver::isExcitation (circuit * c) {
00327   int type = c->getType ();
00328   if (type == CIR_PAC || type == CIR_VAC || type == CIR_VDC)
00329     return true;
00330   return false;
00331 }
00332 
00333 // Expands the frequency array using the given frequency and the order.
00334 void hbsolver::expandFrequencies (nr_double_t f, int n) {
00335   auto nfreqs = negfreqs;
00336   auto pfreqs = posfreqs;
00337   int i, k, len = nfreqs.size ();
00338   negfreqs.clear ();
00339   posfreqs.clear ();
00340   if (len > 0) {
00341     // frequency expansion for full frequency sets
00342     for (i = 0; i <= n + 1; i++) {
00343       for (k = 0; k < len; k++) {
00344         negfreqs.push_back (i * f + nfreqs[k]);
00345       }
00346     }
00347     for (i = -n; i < 0; i++) {
00348       for (k = 0; k < len; k++) {
00349         negfreqs.push_back (i * f + nfreqs[k]);
00350       }
00351     }
00352     for (i = 0; i <= 2 * n + 1; i++) {
00353       for (k = 0; k < len; k++) {
00354         posfreqs.push_back (i * f + pfreqs[k]);
00355       }
00356     }
00357   }
00358   else {
00359     // first frequency
00360     for (i = 0; i <= n + 1; i++) negfreqs.push_back (i * f);
00361     for (i = -n; i < 0; i++) negfreqs.push_back (i * f);
00362     for (i = 0; i <= 2 * n + 1; i++) posfreqs.push_back (i * f);
00363   }
00364 }
00365 
00366 // Calculates an order fulfilling the "power of two" requirement.
00367 int hbsolver::calcOrder (int n) {
00368   int o, order = n * 2;             // current order + DC + negative freqencies
00369   for (o = 1; o < order; o <<= 1) ; // a power of 2
00370   return o / 2 - 1;
00371 }
00372 
00373 /* The function computes the harmonic frequencies excited in the
00374    circuit list depending on the maximum number of harmonics per
00375    exitation and saves its results into the 'negfreqs' vector. */
00376 void hbsolver::collectFrequencies (void) {
00377 
00378   // initialization
00379   negfreqs.clear ();
00380   posfreqs.clear ();
00381   rfreqs.clear ();
00382   dfreqs.clear ();
00383   delete[] ndfreqs;
00384 
00385   // obtain order
00386   int i, n = calcOrder (getPropertyInteger ("n"));
00387 
00388   // expand frequencies for each exitation
00389   nr_double_t f;
00390   for (auto * c : excitations) {
00391     if (c->getType () != CIR_VDC) { // no extra DC sources
00392       if ((f = c->getPropertyDouble ("f")) != 0.0) {
00393         const auto epsilon = std::numeric_limits<nr_double_t>::epsilon();
00394         auto found =
00395           std::find_if(dfreqs.cbegin(),dfreqs.cend(),
00396                        [f,epsilon](decltype(*dfreqs.cbegin()) x) {
00397                          return (std::abs(x-f) < epsilon);
00398                        })
00399           ;
00400         if (found == dfreqs.cend()) { // no double frequencies
00401           dfreqs.push_back (f);
00402           expandFrequencies (f, n);
00403         }
00404       }
00405     }
00406   }
00407 
00408   // no excitations
00409   if (negfreqs.size () == 0) {
00410     // use specified frequency
00411     f = getPropertyDouble ("f");
00412     dfreqs.push_back (f);
00413     expandFrequencies (f, n);
00414   }
00415 
00416   // build frequency dimension lengths
00417   ndfreqs = new int[dfreqs.size ()];
00418   for (i = 0; i < dfreqs.size (); i++) {
00419     ndfreqs[i] = (n + 1) * 2;
00420   }
00421 
00422 #if HB_DEBUG
00423   fprintf (stderr, "%d frequencies: [ ", negfreqs.getSize ());
00424   for (i = 0; i < negfreqs.getSize (); i++) {
00425     fprintf (stderr, "%g ", (double) real (negfreqs.get (i)));
00426   }
00427   fprintf (stderr, "]\n");
00428 #endif /* HB_DEBUG */
00429 
00430   // build list of positive frequencies including DC
00431   for (n = 0; n < negfreqs.size (); n++) {
00432     if ((f = negfreqs[n]) < 0.0) continue;
00433     rfreqs.push_back (f);
00434   }
00435   lnfreqs = rfreqs.size ();
00436   nlfreqs = negfreqs.size ();
00437 
00438   // pre-calculate the j[O] vector
00439   OM = new tvector<nr_complex_t> (nlfreqs);
00440   for (n = i = 0; n < nlfreqs; n++, i++)
00441     OM_(n) = nr_complex_t (0, 2 * pi * negfreqs[i]);
00442 }
00443 
00444 // Split netlist into excitation, linear and non-linear part.
00445 void hbsolver::splitCircuits (void) {
00446   circuit * root = subnet->getRoot ();
00447   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00448     if (c->isNonLinear ()) {
00449       // non-linear part
00450       nolcircuits.push_front(c);
00451     }
00452     else if (isExcitation (c)) {
00453       // get sinusoidal sources
00454       excitations.push_front(c);
00455     }
00456     else if (c->getType () != CIR_GROUND) {
00457       // linear part
00458       lincircuits.push_front (c);
00459     }
00460   }
00461 }
00462 
00463 // Creates a nodelist for the given circuit list.
00464 strlist * hbsolver::circuitNodes (ptrlist<circuit> circuits) {
00465   strlist * nodes = new strlist ();
00466   for (auto * c : circuits) {
00467     for (int i = 0; i < c->getSize (); i++) {
00468       const char * n = c->getNode(i)->getName ();
00469       if (strcmp (n, "gnd")) {
00470         if (!nodes->contains (n)) nodes->add (n);
00471       }
00472     }
00473   }
00474   return nodes;
00475 }
00476 
00477 // Obtain node lists for linear and non-linear part.
00478 void hbsolver::getNodeLists (void) {
00479   // non-linear nodes
00480   nlnodes = circuitNodes (nolcircuits);
00481   // linear nodes
00482   lnnodes = circuitNodes (lincircuits);
00483   // excitation nodes
00484   exnodes = circuitNodes (excitations);
00485 
00486 #if DEBUG && 0
00487   fprintf (stderr, "nonlinear nodes: [ %s ]\n", nlnodes->toString ());
00488   fprintf (stderr, "   linear nodes: [ %s ]\n", lnnodes->toString ());
00489 #endif
00490 
00491   // organization of the nodes for the MNA:
00492   // --------------------------------------
00493   // 1.)  balanced nodes: all connected to at least one non-linear device
00494   // 2.a) the excitation nodes
00495   // 2.b) all linear nodes not contained in 1. and 2.a
00496   // 2.c) additional gyrators of linear nodes (built-in voltage sources)
00497   // please note: excitation nodes also in 2.b; 1. and 2.a are 'ports'
00498 
00499   nanodes = new strlist (*nlnodes); // list 1.
00500   strlistiterator it;
00501   // add excitation nodes; list 2.a
00502   for (it = strlistiterator (exnodes); *it; ++it)
00503     nanodes->append (*it);
00504   // add linear nodes; list 2.b
00505   for (it = strlistiterator (lnnodes); *it; ++it) {
00506     if (!nanodes->contains (*it))
00507       nanodes->append (*it);
00508   }
00509 
00510   banodes = new strlist (*nlnodes);
00511 #if DEBUG && 0
00512   fprintf (stderr, " balanced nodes: [ %s ]\n", banodes->toString ());
00513   fprintf (stderr, "  exnodes nodes: [ %s ]\n", exnodes->toString ());
00514   fprintf (stderr, "  nanodes nodes: [ %s ]\n", nanodes->toString ());
00515 #endif
00516 }
00517 
00518 /* The function applies unique voltage source identifiers to each
00519    built in internal voltage source in the list of circuits. */
00520 int hbsolver::assignVoltageSources (ptrlist<circuit> circuits) {
00521   int sources = 0;
00522   for (auto *c: circuits) {
00523     if (c->getVoltageSources () > 0) {
00524       c->setVoltageSource (sources);
00525       sources += c->getVoltageSources ();
00526     }
00527   }
00528   return sources;
00529 }
00530 
00531 /* The function assigns unique node identifiers to each circuit node. */
00532 int hbsolver::assignNodes (ptrlist<circuit> circuits, strlist * nodes,
00533                            int offset) {
00534   // through all collected nodes
00535   for (int nr = 0; nr < nodes->length (); nr++) {
00536     char * nn = nodes->get (nr);
00537     // through all circuits
00538     for (auto *c : circuits) {
00539       // assign current identifier if any of the circuit nodes equals
00540       // the current one
00541       for (int i = 0; i < c->getSize (); i++) {
00542         node * n = c->getNode (i);
00543         if (!strcmp (n->getName (), nn)) n->setNode (offset + nr + 1);
00544       }
00545     }
00546   }
00547   return nodes->length ();
00548 }
00549 
00550 // Prepares the linear operations.
00551 void hbsolver::prepareLinear (void) {
00552   for (auto *lc : lincircuits)
00553     lc->initHB ();
00554   nlnvsrcs = assignVoltageSources (lincircuits);
00555   nnlvsrcs = excitations.size ();
00556   nnanodes = nanodes->length ();
00557   nexnodes = exnodes->length ();
00558   nbanodes = banodes->length ();
00559   assignNodes (lincircuits, nanodes);
00560   assignNodes (excitations, nanodes);
00561   createMatrixLinearA ();
00562   createMatrixLinearY ();
00563   calcConstantCurrent ();
00564 }
00565 
00566 /* The function creates the complex linear network MNA matrix.  It
00567    contains the MNA entries for all linear components for each
00568    requested frequency. */
00569 void hbsolver::createMatrixLinearA (void) {
00570   int M = nlnvsrcs;
00571   int N = nnanodes;
00572   int f = 0;
00573   nr_double_t freq;
00574 
00575   // create new MNA matrix
00576   A = new tmatrix<nr_complex_t> ((N + M) * lnfreqs);
00577 
00578   // through each frequency
00579   for (int i = 0; i < rfreqs.size (); i++) {
00580     freq = rfreqs[i];
00581     // calculate components' MNA matrix for the given frequency
00582     for (auto *lc : lincircuits)
00583       lc->calcHB (freq);
00584     // fill in all matrix entries for the given frequency
00585     fillMatrixLinearA (A, f++);
00586   }
00587 
00588   // save a copy of the original MNA matrix
00589   NA = new tmatrix<nr_complex_t> (*A);
00590 }
00591 
00592 // some definitions for the linear matrix filler
00593 #undef  A_
00594 #undef  B_
00595 #define A_(r,c) (*A) ((r)*lnfreqs+f,(c)*lnfreqs+f)
00596 #define G_(r,c) A_(r,c)
00597 #define B_(r,c) A_(r,c+N)
00598 #define C_(r,c) A_(r+N,c)
00599 #define D_(r,c) A_(r+N,c+N)
00600 
00601 /* This function fills in the MNA matrix entries into the A matrix for
00602    a given frequency index. */
00603 void hbsolver::fillMatrixLinearA (tmatrix<nr_complex_t> * A, int f) {
00604   int N = nnanodes;
00605 
00606   // through each linear circuit
00607   for (auto *cir : lincircuits) {
00608     int s = cir->getSize ();
00609     int nr, nc, r, c, v;
00610 
00611     // apply G-matrix entries
00612     for (r = 0; r < s; r++) {
00613       if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue;
00614       for (c = 0; c < s; c++) {
00615         if ((nc = cir->getNode(c)->getNode () - 1) < 0) continue;
00616         G_(nr, nc) += cir->getY (r, c);
00617       }
00618     }
00619 
00620     // augmented part -- built in voltage sources
00621     if ((v = cir->getVoltageSources ()) > 0) {
00622 
00623       // apply B-matrix entries
00624       for (r = 0; r < s; r++) {
00625         if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue;
00626         for (c = 0; c < v; c++) {
00627           nc = cir->getVoltageSource () + c;
00628           B_(nr, nc) += cir->getB (r, nc);
00629         }
00630       }
00631 
00632       // apply C-matrix entries
00633       for (r = 0; r < v; r++) {
00634         nr = cir->getVoltageSource () + r;
00635         for (c = 0; c < s; c++) {
00636           if ((nc = cir->getNode(c)->getNode () - 1) < 0) continue;
00637           C_(nr, nc) += cir->getC (nr, c);
00638         }
00639       }
00640 
00641       // apply D-matrix entries
00642       for (r = 0; r < v; r++) {
00643         nr = cir->getVoltageSource () + r;
00644         for (c = 0; c < v; c++) {
00645           nc = cir->getVoltageSource () + c;
00646           D_(nr, nc) += cir->getD (nr, nc);
00647         }
00648       }
00649     }
00650   }
00651 }
00652 
00653 // The function inverts the given matrix A into the matrix H.
00654 void hbsolver::invertMatrix (tmatrix<nr_complex_t> * A,
00655                              tmatrix<nr_complex_t> * H) {
00656   eqnsys<nr_complex_t> eqns;
00657   int N = A->getCols ();
00658   tvector<nr_complex_t> * x = new tvector<nr_complex_t> (N);
00659   tvector<nr_complex_t> * z = new tvector<nr_complex_t> (N);
00660 
00661   try_running () {
00662     // create LU decomposition of the A matrix
00663     eqns.setAlgo (ALGO_LU_FACTORIZATION_CROUT);
00664     eqns.passEquationSys (A, x, z);
00665     eqns.solve ();
00666   }
00667   // appropriate exception handling
00668   catch_exception () {
00669   case EXCEPTION_PIVOT:
00670   default:
00671     logprint (LOG_ERROR, "WARNING: %s: during TI inversion\n", getName ());
00672     estack.print ();
00673   }
00674 
00675   // use the LU decomposition to obtain the inverse H
00676   eqns.setAlgo (ALGO_LU_SUBSTITUTION_CROUT);
00677   for (int c = 0; c < N; c++) {
00678     z->set (0.0);
00679     z->set (c, 1.0);
00680     eqns.passEquationSys (A, x, z);
00681     eqns.solve ();
00682     for (int r = 0; r < N; r++) H->set (r, c, x->get (r));
00683   }
00684   delete x;
00685   delete z;
00686 }
00687 
00688 // Some defines for matrix element access.
00689 #define V_(r) (*V) (r)
00690 #define I_(r) (*I) (r)
00691 #undef  A_
00692 #define A_(r,c) (*A) (r,c)
00693 
00694 #define Z_(r,c) (*Z) (r,c)
00695 #define Y_(r,c) (*Y) (r,c)
00696 
00697 #define ZVU_(r,c) Z_(r,c)
00698 #define ZVL_(r,c) Z_((r)*lnfreqs+f+sn,c)
00699 #define ZCU_(r,c) Z_(r,(c)*lnfreqs+f+sn)
00700 #define ZCL_(r,c) Z_((r)*lnfreqs+f+sn,(c)*lnfreqs+f+sn)
00701 
00702 #define YV_(r,c) (*YV) (r,c)
00703 #define NA_(r,c) (*NA) (r,c)
00704 #define JF_(r,c) (*JF) (r,c)
00705 
00706 /* The following function performs the following steps:
00707    1. form the MNA matrix A including all nodes (linear, non-linear and
00708       excitations)
00709    2. compute the variable transimpedance matrix entries for the nodes
00710       to be balanced
00711    3. compute the constant transimpedance matrix entries for the constant
00712       current vector caused by the excitations
00713    4. invert this overall transimpedance matrix
00714    5. extract the variable transadmittance matrix entries
00715 */
00716 void hbsolver::createMatrixLinearY (void) {
00717   int M = nlnvsrcs;
00718   int N = nnanodes;
00719   int c, r, f;
00720 
00721   // size of overall MNA matrix
00722   int sa = (N + M) * lnfreqs;
00723   int sv = nbanodes;
00724   int se = nnlvsrcs;
00725   int sy = sv + se;
00726 
00727   // allocate new transimpedance matrix
00728   Z = new tmatrix<nr_complex_t> (sy * lnfreqs);
00729 
00730   // prepare equation system
00731   eqnsys<nr_complex_t> eqns;
00732   tvector<nr_complex_t> * V;
00733   tvector<nr_complex_t> * I;
00734 
00735   // 1. create variable transimpedance matrix entries relating
00736   // voltages at the balanced nodes to the currents through these
00737   // nodes into the non-linear part
00738   int sn = sv * lnfreqs;
00739   V = new tvector<nr_complex_t> (sa);
00740   I = new tvector<nr_complex_t> (sa);
00741 
00742   // connect a 100 Ohm resistor (to ground) to balanced node in the MNA matrix
00743   for (c = 0; c < sv * lnfreqs; c++) A_(c, c) += 0.01;
00744 
00745   // connect a 100 Ohm resistor (in parallel) to each excitation
00746   for (auto *vs : excitations) {
00747     // get positive and negative node
00748     int pnode = vs->getNode(NODE_1)->getNode ();
00749     int nnode = vs->getNode(NODE_2)->getNode ();
00750     for (f = 0; f < lnfreqs; f++) { // for each frequency
00751       int pn = (pnode - 1) * lnfreqs + f;
00752       int nn = (nnode - 1) * lnfreqs + f;
00753       if (pnode) A_(pn, pn) += 0.01;
00754       if (nnode) A_(nn, nn) += 0.01;
00755       if (pnode && nnode) {
00756         A_(pn, nn) -= 0.01;
00757         A_(nn, pn) -= 0.01;
00758       }
00759     }
00760   }
00761 
00762   // LU decompose the MNA matrix
00763   try_running () {
00764     eqns.setAlgo (ALGO_LU_FACTORIZATION_CROUT);
00765     eqns.passEquationSys (A, V, I);
00766     eqns.solve ();
00767   }
00768   // appropriate exception handling
00769   catch_exception () {
00770   case EXCEPTION_PIVOT:
00771   default:
00772     logprint (LOG_ERROR, "WARNING: %s: during A factorization\n", getName ());
00773     estack.print ();
00774   }
00775 
00776   // aquire variable transimpedance matrix entries
00777   eqns.setAlgo (ALGO_LU_SUBSTITUTION_CROUT);
00778   for (c = 0; c < sn; c++) {
00779     I->set (0.0);
00780     I_(c) = 1.0;
00781     eqns.passEquationSys (A, V, I);
00782     eqns.solve ();
00783     // ZV | ..
00784     // ---+---
00785     // .. | ..
00786     for (r = 0; r < sn; r++) ZVU_(r, c) = V_(r);
00787     // .. | ..
00788     // ---+---
00789     // ZV | ..
00790     r = 0;
00791     for (auto ite = excitations.begin(); ite != excitations.end(); ++ite, r++) {
00792       // lower part entries
00793       for (f = 0; f < lnfreqs; f++) {
00794         ZVL_(r, c) = excitationZ (V, *ite, f);
00795       }
00796     }
00797   }
00798 
00799   // create constant transimpedance matrix entries relating the
00800   // source voltages to the interconnection currents
00801   int vsrc = 0;
00802   // aquire constant transadmittance matrix entries
00803   for (auto it = excitations.begin(); it != excitations.end(); ++it, vsrc++) {
00804     circuit * vs = *it;
00805     // get positive and negative node
00806     int pnode = vs->getNode(NODE_1)->getNode ();
00807     int nnode = vs->getNode(NODE_2)->getNode ();
00808     for (f = 0; f < lnfreqs; f++) { // for each frequency
00809       int pn = (pnode - 1) * lnfreqs + f;
00810       int nn = (nnode - 1) * lnfreqs + f;
00811       I->set (0.0);
00812       if (pnode) I_(pn) = +1.0;
00813       if (nnode) I_(nn) = -1.0;
00814       eqns.passEquationSys (A, V, I);
00815       eqns.solve ();
00816       // .. | ZC
00817       // ---+---
00818       // .. | ..
00819       for (r = 0; r < sn; r++) {
00820         // upper part of the entries
00821         ZCU_(r, vsrc) = V_(r);
00822       }
00823       // .. | ..
00824       // ---+---
00825       // .. | ZC
00826       r = 0;
00827       for (auto ite = excitations.begin(); ite != excitations.end(); ++ite, r++) {
00828         // lower part entries
00829         ZCL_(r, vsrc) = excitationZ (V, *ite, f);
00830       }
00831     }
00832   }
00833   delete I;
00834   delete V;
00835 
00836   // allocate new transadmittance matrix
00837   Y = new tmatrix<nr_complex_t> (sy * lnfreqs);
00838 
00839   // invert the Z matrix to a Y matrix
00840   invertMatrix (Z, Y);
00841 
00842   // substract the 100 Ohm resistor
00843   for (c = 0; c < sy * lnfreqs; c++) Y_(c, c) -= 0.01;
00844 
00845   // extract the variable transadmittance matrix
00846   YV = new tmatrix<nr_complex_t> (sv * nlfreqs);
00847 
00848   // variable transadmittance matrix must be continued conjugately
00849   *YV = expandMatrix (*Y, sv);
00850 
00851   // delete overall temporary MNA matrix
00852   delete A; A = NULL;
00853   // delete transimpedance matrix
00854   delete Z; Z = NULL;
00855 }
00856 
00857 /* Little helper function obtaining a transimpedance value for the
00858    given voltage source (excitation) and for a given frequency
00859    index. */
00860 nr_complex_t hbsolver::excitationZ (tvector<nr_complex_t> * V, circuit * vs,
00861                                     int f) {
00862   // get positive and negative node
00863   int pnode = vs->getNode(NODE_1)->getNode ();
00864   int nnode = vs->getNode(NODE_2)->getNode ();
00865   nr_complex_t z = 0.0;
00866   if (pnode) z += V_((pnode - 1) * lnfreqs + f);
00867   if (nnode) z -= V_((nnode - 1) * lnfreqs + f);
00868   return z;
00869 }
00870 
00871 /* This function computes the constant current vectors using the
00872    voltage of the excitations and the transadmittance matrix
00873    entries. */
00874 void hbsolver::calcConstantCurrent (void) {
00875   int se = nnlvsrcs * lnfreqs;
00876   int sn = nbanodes * lnfreqs;
00877   int r, c, vsrc = 0;
00878 
00879   // collect excitation voltages
00880   tvector<nr_complex_t> VC (se);
00881   for (auto it = excitations.begin(); it != excitations.end(); ++it, vsrc++) {
00882     circuit * vs = *it;
00883     vs->initHB ();
00884     vs->setVoltageSource (0);
00885     for (int f = 0; f < rfreqs.size (); f++) { // for each frequency
00886       nr_double_t freq = rfreqs[f];
00887       vs->calcHB (freq);
00888       VC (vsrc * lnfreqs + f) = vs->getE (VSRC_1);
00889     }
00890   }
00891 
00892   // compute constant current vector for balanced nodes
00893   IC = new tvector<nr_complex_t> (sn);
00894   // .. | YC * VC
00895   // ---+---
00896   // .. | ..
00897   for (r = 0; r < sn; r++) {
00898     nr_complex_t i = 0.0;
00899     for (c = 0; c < se; c++) {
00900       i += Y_(r, c + sn) * VC (c);
00901     }
00902     int f = r % lnfreqs;
00903     if (f != 0 && f != lnfreqs - 1) i /= 2;
00904     IC->set (r, i);
00905   }
00906   // expand the constant current conjugate
00907   *IC = expandVector (*IC, nbanodes);
00908 
00909   // compute constant current vector for sources itself
00910   IS = new tvector<nr_complex_t> (se);
00911   // .. | ..
00912   // ---+---
00913   // .. | YC * VC
00914   for (r = 0; r < se; r++) {
00915     nr_complex_t i = 0.0;
00916     for (c = 0; c < se; c++) {
00917       i += Y_(r + sn, c + sn) * VC (c);
00918     }
00919     IS->set (r, i);
00920   }
00921 
00922   // delete overall transadmittance matrix
00923   delete Y; Y = NULL;
00924 }
00925 
00926 /* Checks whether currents through the interconnects of the linear and
00927    non-linear subcircuit (in the frequency domain) are equal. */
00928 int hbsolver::checkBalance (void) {
00929   nr_double_t iabstol = getPropertyDouble ("iabstol");
00930   nr_double_t vabstol = getPropertyDouble ("vabstol");
00931   nr_double_t reltol = getPropertyDouble ("reltol");
00932   int n, len = FV->size ();
00933   for (n = 0; n < len; n++) {
00934     // check iteration voltages
00935     nr_double_t v_abs = abs (VS->get (n) - VP->get (n));
00936     nr_double_t v_rel = abs (VS->get (n));
00937     if (v_abs >= vabstol + reltol * v_rel) return 0;
00938     // check balanced currents
00939     nr_complex_t il = IL->get (n);
00940     nr_complex_t in = IN->get (n);
00941     if (il != in) {
00942       nr_double_t i_abs = abs (il + in);
00943       nr_double_t i_rel = abs ((il + in) / (il - in));
00944       if (i_abs >= iabstol && 2.0 * i_rel >= reltol) return 0;
00945     }
00946   }
00947   return 1;
00948 }
00949 
00950 // some definitions for the non-linear matrix filler
00951 #undef  G_
00952 #undef  C_
00953 #define G_(r,c) (*jg) ((r)*nlfreqs+f,(c)*nlfreqs+f)
00954 #define C_(r,c) (*jq) ((r)*nlfreqs+f,(c)*nlfreqs+f)
00955 #undef  FI_
00956 #undef  FQ_
00957 #define FI_(r) (*ig) ((r)*nlfreqs+f)
00958 #define FQ_(r) (*fq) ((r)*nlfreqs+f)
00959 #define IR_(r) (*ir) ((r)*nlfreqs+f)
00960 #define QR_(r) (*qr) ((r)*nlfreqs+f)
00961 
00962 /* This function fills in the matrix and vector entries for the
00963    non-linear HB equations for a given frequency index. */
00964 void hbsolver::fillMatrixNonLinear (tmatrix<nr_complex_t> * jg,
00965                                     tmatrix<nr_complex_t> * jq,
00966                                     tvector<nr_complex_t> * ig,
00967                                     tvector<nr_complex_t> * fq,
00968                                     tvector<nr_complex_t> * ir,
00969                                     tvector<nr_complex_t> * qr,
00970                                     int f) {
00971   // through each linear circuit
00972   for (auto *cir: nolcircuits) {
00973     int s = cir->getSize ();
00974     int nr, nc, r, c;
00975 
00976     for (r = 0; r < s; r++) {
00977       if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue;
00978       // apply G- and C-matrix entries
00979       for (c = 0; c < s; c++) {
00980         if ((nc = cir->getNode(c)->getNode () - 1) < 0) continue;
00981         G_(nr, nc) += cir->getY (r, c);
00982         C_(nr, nc) += cir->getQV (r, c);
00983       }
00984       // apply I- and Q-vector entries
00985       FI_(nr) -= cir->getI (r);
00986       FQ_(nr) -= cir->getQ (r);
00987       // ThinkME: positive or negative?
00988       IR_(nr) += cir->getGV (r) + cir->getI (r);
00989       QR_(nr) += cir->getCV (r) + cir->getQ (r);
00990     }
00991   }
00992 }
00993 
00994 /* The function initializes the non-linear part of the HB. */
00995 void hbsolver::prepareNonLinear (void) {
00996   int N = nbanodes;
00997 
00998   // allocate matrices and vectors
00999   if (FQ == NULL) {
01000     FQ = new tvector<nr_complex_t> (N * nlfreqs);
01001   }
01002   if (IG == NULL) {
01003     IG = new tvector<nr_complex_t> (N * nlfreqs);
01004   }
01005   if (IR == NULL) {
01006     IR = new tvector<nr_complex_t> (N * nlfreqs);
01007   }
01008   if (QR == NULL) {
01009     QR = new tvector<nr_complex_t> (N * nlfreqs);
01010   }
01011   if (JG == NULL) {
01012     JG = new tmatrix<nr_complex_t> (N * nlfreqs);
01013   }
01014   if (JQ == NULL) {
01015     JQ = new tmatrix<nr_complex_t> (N * nlfreqs);
01016   }
01017   if (JF == NULL) {
01018     JF = new tmatrix<nr_complex_t> (N * nlfreqs);
01019   }
01020 
01021   // voltage vector in frequency and time domain
01022   if (VS == NULL) {
01023     VS = new tvector<nr_complex_t> (N * nlfreqs);
01024   }
01025   if (vs == NULL) {
01026     vs = new tvector<nr_complex_t> (N * nlfreqs);
01027   }
01028   if (VP == NULL) {
01029     VP = new tvector<nr_complex_t> (N * nlfreqs);
01030   }
01031 
01032   // error vector
01033   if (FV == NULL) {
01034     FV = new tvector<nr_complex_t> (N * nlfreqs);
01035   }
01036   // right hand side vector
01037   if (RH == NULL) {
01038     RH = new tvector<nr_complex_t> (N * nlfreqs);
01039   }
01040 
01041   // linear and non-linear current vector
01042   if (IL == NULL) {
01043     IL = new tvector<nr_complex_t> (N * nlfreqs);
01044   }
01045   if (IN == NULL) {
01046     IN = new tvector<nr_complex_t> (N * nlfreqs);
01047   }
01048 
01049   // assign nodes
01050   assignNodes (nolcircuits, nanodes);
01051 
01052   // initialize circuits
01053   for (auto *cir : nolcircuits) {
01054     cir->initHB (nlfreqs);
01055   }
01056 }
01057 
01058 /* Saves the node voltages of the given circuit and for the given
01059    frequency entry into the circuit voltage vector. */
01060 void hbsolver::saveNodeVoltages (circuit * cir, int f) {
01061   int r, nr, s = cir->getSize ();
01062   for (r = 0; r < s; r++) {
01063     if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue;
01064     // apply V-vector entries
01065     cir->setV (r, real (vs->get (nr * nlfreqs + f)));
01066   }
01067 }
01068 
01069 /* The function saves voltages into non-linear circuits, runs each
01070    non-linear components' HB calculator for each frequency and applies
01071    the matrix and vector entries appropriately. */
01072 void hbsolver::loadMatrices (void) {
01073   // clear matrices and vectors before
01074   IG->set (0.0);
01075   FQ->set (0.0);
01076   IR->set (0.0);
01077   QR->set (0.0);
01078   JG->set (0.0);
01079   JQ->set (0.0);
01080   // through each frequency
01081   for (int f = 0; f < nlfreqs; f++) {
01082     // calculate components' HB matrices and vector for the given frequency
01083     for (auto *cir : nolcircuits) {
01084       saveNodeVoltages (cir, f); // node voltages
01085       cir->calcHB (f);         // HB calculator
01086     }
01087     // fill in all matrix entries for the given frequency
01088     fillMatrixNonLinear (JG, JQ, IG, FQ, IR, QR, f);
01089   }
01090 }
01091 
01092 /* The following function transforms a vector using a Fast Fourier
01093    Transformation from the time domain to the frequency domain. 
01094    \todo rewrite ugly sould die
01095 */
01096 void hbsolver::VectorFFT (tvector<nr_complex_t> * V, int isign) {
01097   int i, k, r;
01098   int n = nlfreqs;
01099   int nd = dfreqs.size ();
01100   int nodes = V->size () / n;
01101   nr_double_t * d = (double *)V->getData ();
01102 
01103   if (nd == 1) {
01104     // for each node a single 1d-FFT
01105     for (k = i = 0; i < nodes; i++, k += 2 * n) {
01106       nr_double_t * dst = &d[k];
01107       _fft_1d (dst, n, isign);
01108       if (isign > 0) for (r = 0; r < 2 * n; r++) *dst++ /= n;
01109     }
01110   }
01111   else {
01112     // for each node a single nd-FFT
01113     for (k = i = 0; i < nodes; i++, k += 2 * n) {
01114       nr_double_t * dst = &d[k];
01115       _fft_nd (dst, ndfreqs, nd, isign);
01116       if (isign > 0) for (r = 0; r < 2 * n; r++) *dst++ /= ndfreqs[0];
01117     }
01118   }
01119 }
01120 
01121 /* The following function transforms a vector using an Inverse Fast
01122    Fourier Transformation from the frequency domain to the domain
01123    time. */
01124 void hbsolver::VectorIFFT (tvector<nr_complex_t> * V, int isign) {
01125   VectorFFT (V, -isign);
01126 }
01127 
01128 /* The following function transforms a matrix using a Fast Fourier
01129    Transformation from the time domain to the frequency domain. */
01130 void hbsolver::MatrixFFT (tmatrix<nr_complex_t> * M) {
01131 
01132 #if THE_SLO_ALGO
01133   // each column FFT
01134   for (int c = 0; c < M->getCols (); c++) {
01135     tvector<nr_complex_t> V = M->getCol (c);
01136     VectorFFT (&V);
01137     M->setCol (c, V);
01138   }
01139 
01140   // each row IFFT
01141   for (int r = 0; r < M->getRows (); r++) {
01142     tvector<nr_complex_t> V = M->getRow (r);
01143     VectorIFFT (&V);
01144     M->setRow (r, V);
01145   }
01146 #else
01147   int c, r, nc, nr;
01148   // for each non-linear node block
01149   for (nc = c = 0; c < nbanodes; c++, nc += nlfreqs) {
01150     for (nr = r = 0; r < nbanodes; r++, nr += nlfreqs) {
01151       tvector<nr_complex_t> V (nlfreqs);
01152       int fr, fc, fi;
01153       // transform the sub-diagonal only
01154       for (fc = 0; fc < nlfreqs; fc++) V (fc) = M->get (nr + fc, nc + fc);
01155       VectorFFT (&V);
01156       // fill in resulting sub-matrix for the node
01157       for (fc = 0; fc < nlfreqs; fc++) {
01158         for (fi = nlfreqs - 1 - fc, fr = 0; fr < nlfreqs; fr++) {
01159           if (++fi >= nlfreqs) fi = 0;
01160           M->set (nr + fr, nc + fc, V (fi));
01161         }
01162       }
01163     }
01164   }
01165 #endif
01166 }
01167 
01168 /* This function solves the actual HB equation in the frequency domain.
01169    F(V) = IC + [YV] * VS + j[O] * FQ + IG -> 0
01170    Care must be taken with indexing here: In the frequency domain only
01171    real positive frequencies are used and computed, but in the time
01172    domain we have more values because of the periodic continuation in
01173    the frequency domain.
01174    RHS = j[O] * CV + GV - (IC + j[O] * FQ + IG)
01175    Also the right hand side of the equation system for the new voltage
01176    vector is computed here. */
01177 void hbsolver::solveHB (void) {
01178   // for each non-linear node
01179   for (int r = 0; r < nbanodes * nlfreqs; ) {
01180     // for each frequency
01181     for (int f = 0; f < nlfreqs; f++, r++) {
01182       nr_complex_t il = 0.0, in = 0.0, ir = 0.0;
01183       // constant current vector due to sources
01184       il += IC->get (r);
01185       // part 1 of right hand side vector
01186       ir -= il;
01187       // transadmittance matrix multiplied by voltage vector
01188       for (int c = 0; c < nbanodes * nlfreqs; c++) {
01189         il += YV_(r, c) * VS_(c);
01190       }
01191       // charge vector
01192       in += OM_(f) * FQ->get (r);
01193       // current vector
01194       in += IG->get (r);
01195       // part 2, 3 and 4 of right hand side vector
01196       ir += IR->get (r);
01197       ir += OM_(f) * QR->get (r);
01198       // put values into result vectors
01199       RH->set (r, ir);
01200       FV->set (r, il + in);
01201       IL->set (r, il);
01202       IN->set (r, in);
01203     }
01204   }
01205 }
01206 
01207 /* The function calculates the full Jacobian JF = [YV] + j[O] * JQ + JG */
01208 void hbsolver::calcJacobian (void) {
01209   int c, r, fc, fr, rt, ct;
01210   /* add admittances of capacitance matrix JQ and non-linear
01211      admittances matrix JG into complete Jacobian JF */
01212   for (c = 0; c < nbanodes; c++) {
01213     ct = c * nlfreqs;
01214     for (fc = 0; fc < nlfreqs; fc++, ct++) {
01215       for (r = 0; r < nbanodes; r++) {
01216         rt = r * nlfreqs;
01217         for (fr = 0; fr < nlfreqs; fr++, rt++) {
01218           JF_(rt, ct) = JG->get (rt, ct) + JQ->get (rt, ct) * OM_(fr);
01219         }
01220       }
01221     }
01222   }
01223   *JF += *YV; // add linear admittance matrix
01224 }
01225 
01226 /* The function expands the given vector in the frequency domain to
01227    make it a real valued signal in the time domain. */
01228 tvector<nr_complex_t> hbsolver::expandVector (tvector<nr_complex_t> V,
01229                                               int nodes) {
01230   tvector<nr_complex_t> res (nodes * nlfreqs);
01231   int r, ff, rf, rt;
01232   for (r = 0; r < nodes; r++) {
01233     rt = r * nlfreqs;
01234     rf = r * lnfreqs;
01235     // copy first part of vector
01236     for (ff = 0; ff < lnfreqs; ff++, rf++, rt++) {
01237       res (rt) = V (rf);
01238     }
01239     // continue vector conjugated
01240     for (rf -= 2; ff < nlfreqs; ff++, rf--, rt++) {
01241       res (rt) = conj (V (rf));
01242     }
01243   }
01244   return res;
01245 }
01246 
01247 /* The function expands the given matrix in the frequency domain to
01248    make it a real valued signal in the time domain. */
01249 tmatrix<nr_complex_t> hbsolver::expandMatrix (tmatrix<nr_complex_t> M,
01250                                               int nodes) {
01251   tmatrix<nr_complex_t> res (nodes * nlfreqs);
01252   int r, c, rf, rt, cf, ct, ff;
01253   for (r = 0; r < nodes; r++) {
01254     for (c = 0; c < nodes; c++) {
01255       rf = r * lnfreqs;
01256       rt = r * nlfreqs;
01257       cf = c * lnfreqs;
01258       ct = c * nlfreqs;
01259       // copy first part of diagonal
01260       for (ff = 0; ff < lnfreqs; ff++, cf++, ct++, rf++, rt++) {
01261         res (rt, ct) = M (rf, cf);
01262       }
01263       // continue diagonal conjugated
01264       for (cf -= 2, rf -= 2; ff < nlfreqs; ff++, cf--, ct++, rf--, rt++) {
01265         res (rt, ct) = conj (M (rf, cf));
01266       }
01267     }
01268   }
01269   return res;
01270 }
01271 
01272 /* This function solves the equation system
01273    JF * VS(n+1) = JF * VS(n) - FV
01274    in order to obtains a new voltage vector in the frequency domain. */
01275 void hbsolver::solveVoltages (void) {
01276   // save previous iteration voltage
01277   *VP = *VS;
01278 
01279   // setup equation system
01280   eqnsys<nr_complex_t> eqns;
01281   try_running () {
01282     // use LU decomposition for solving
01283     eqns.setAlgo (ALGO_LU_DECOMPOSITION);
01284     eqns.passEquationSys (JF, VS, RH);
01285     eqns.solve ();
01286   }
01287   // appropriate exception handling
01288   catch_exception () {
01289   case EXCEPTION_PIVOT:
01290   default:
01291     logprint (LOG_ERROR, "WARNING: %s: during NR iteration\n", getName ());
01292     estack.print ();
01293   }
01294 
01295   // save new voltages in time domain vector
01296   *vs = *VS;
01297 }
01298 
01299 /* The following function extends the existing linear MNA matrix to
01300    contain the additional rows and columns for the excitation voltage
01301    sources. */
01302 tmatrix<nr_complex_t> hbsolver::extendMatrixLinear (tmatrix<nr_complex_t> M,
01303                                                     int nodes) {
01304   int no = M.getCols ();
01305   tmatrix<nr_complex_t> res (no + nodes * lnfreqs);
01306   // copy the existing part
01307   for (int r = 0; r < no; r++) {
01308     for (int c = 0; c < no; c++) {
01309       res (r, c) = M (r, c);
01310     }
01311   }
01312   return res;
01313 }
01314 
01315 /* The function fills in the missing MNA entries for the excitation
01316    voltage sources into the extended rows and columns as well as the
01317    actual voltage values into the right hand side vector. */
01318 void hbsolver::fillMatrixLinearExtended (tmatrix<nr_complex_t> * Y,
01319                                          tvector<nr_complex_t> * I) {
01320   // through each excitation source
01321   int of = lnfreqs * (nlnvsrcs + nnanodes);
01322   int sc = of;
01323 
01324   for (auto *vs : excitations) {
01325     // get positive and negative node
01326     int pnode = vs->getNode(NODE_1)->getNode ();
01327     int nnode = vs->getNode(NODE_2)->getNode ();
01328     for (int f = 0; f < lnfreqs; f++, sc++) { // for each frequency
01329       // fill right hand side vector
01330       nr_double_t freq = rfreqs[f];
01331       vs->calcHB (freq);
01332       I_(sc) = vs->getE (VSRC_1);
01333       // fill MNA entries
01334       int pn = (pnode - 1) * lnfreqs + f;
01335       int nn = (nnode - 1) * lnfreqs + f;
01336       if (pnode) {
01337         Y_(pn, sc) = +1.0;
01338         Y_(sc, pn) = +1.0;
01339       }
01340       if (nnode) {
01341         Y_(nn, sc) = -1.0;
01342         Y_(sc, nn) = -1.0;
01343       }
01344     }
01345   }
01346 }
01347 
01348 /* The function calculates and saves the final solution. */
01349 void hbsolver::finalSolution (void) {
01350 
01351   // extend the linear MNA matrix
01352   *NA = extendMatrixLinear (*NA, nnlvsrcs);
01353 
01354   int S = NA->getCols ();
01355   int N = nnanodes * lnfreqs;
01356 
01357   // right hand side vector
01358   tvector<nr_complex_t> * I = new tvector<nr_complex_t> (S);
01359   // temporary solution
01360   tvector<nr_complex_t> * V = new tvector<nr_complex_t> (S);
01361   // final solution
01362   x = new tvector<nr_complex_t> (N);
01363 
01364   // fill in missing MNA entries
01365   fillMatrixLinearExtended (NA, I);
01366 
01367   // put currents through balanced nodes into right hand side
01368   for (int n = 0; n < nbanodes; n++) {
01369     for (int f = 0; f < lnfreqs; f++) {
01370       nr_complex_t i = IL->get (n * nlfreqs + f);
01371       if (f != 0 && f != lnfreqs - 1) i *= 2;
01372       I_(n * lnfreqs + f) = i;
01373     }
01374   }
01375 
01376   // use QR decomposition for the final solution
01377   try_running () {
01378     eqnsys<nr_complex_t> eqns;
01379     eqns.setAlgo (ALGO_LU_DECOMPOSITION);
01380     eqns.passEquationSys (NA, V, I);
01381     eqns.solve ();
01382   }
01383   // appropriate exception handling
01384   catch_exception () {
01385   case EXCEPTION_PIVOT:
01386   default:
01387     logprint (LOG_ERROR, "WARNING: %s: during final AC analysis\n",
01388               getName ());
01389     estack.print ();
01390   }
01391   for (int i = 0; i < N; i++) x->set (i, V_(i));
01392 }
01393 
01394 // Saves simulation results.
01395 void hbsolver::saveResults (void) {
01396   vector * f;
01397   // add current frequency to the dependency of the output dataset
01398   if ((f = data->findDependency ("hbfrequency")) == NULL) {
01399     f = new vector ("hbfrequency");
01400     data->addDependency (f);
01401   }
01402   // save frequency vector
01403   if (runs == 1) {
01404     for (int i = 0; i < lnfreqs; i++) f->add (rfreqs[i]);
01405   }
01406   // save node voltage vectors
01407   int nanode = 0;
01408   for (strlistiterator it (nanodes); *it; ++it, nanode++) {
01409     int l = strlen (*it);
01410     char * n = (char *) malloc (l + 4);
01411     sprintf (n, "%s.Vb", *it);
01412     for (int i = 0; i < lnfreqs; i++) {
01413       saveVariable (n, x->get (i + nanode * lnfreqs), f);
01414     }
01415   }
01416 }
01417 
01418 // properties
01419 PROP_REQ [] = {
01420   { "n", PROP_INT, { 1, PROP_NO_STR }, PROP_MIN_VAL (1) },
01421   PROP_NO_PROP };
01422 PROP_OPT [] = {
01423   { "f", PROP_REAL, { 1e9, PROP_NO_STR }, PROP_POS_RANGEX },
01424   { "iabstol", PROP_REAL, { 1e-12, PROP_NO_STR }, PROP_RNG_X01I },
01425   { "vabstol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
01426   { "reltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
01427   { "MaxIter", PROP_INT, { 150, PROP_NO_STR }, PROP_RNGII (2, 10000) },
01428   PROP_NO_PROP };
01429 struct define_t hbsolver::anadef =
01430   { "HB", 0, PROP_ACTION, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };
01431 
01432 } // namespace qucs
01433 
01434 /* TODO list for HB Solver:
01435    - Take care about nodes with non-linear components only.
01436    - AC Power Sources (extra Z and open loop voltage).
01437    - Current sources.
01438    - Balancing of multi-dimensional non-linear networks.
01439    - Sources directly connected to non-linear components and no other
01440      linear component (insert short).
01441    - Bug: With capacitors at hand there is voltage convergence but no
01442      current balancing.
01443    - Output currents through voltage sources.
01444  */