include/m_matrix.h

   1 /*$Id: m_matrix.h,v 26.131 2009/11/20 08:22:10 al Exp $ -*- C++ -*-
   2  * Copyright (C) 2001 Albert Davis
   3  * Author: Albert Davis <aldavis@gnu.org>
   4  *
   5  * This file is part of "Gnucap", the Gnu Circuit Analysis Package
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  20  * 02110-1301, USA.
  21  *------------------------------------------------------------------
  22  * Sparse matrix package
  23  * Bump and spike - bordered block diagonal pattern
  24  * -------------------------------------------------
  25  * To use it (simple method) .....
  26  * 1. Declare it ...
  27  *      BSMATRIX<double> a;
  28  * 2. Then tell it what slots to allocate, in a loop ...
  29  *      for (all pairs i,j that might exist)
  30  *         a.iwant(i,j);
  31  * 3. Then do the actual allocation ...
  32  *      a.allocate();
  33  * 4. If you want, you can add an offset to the diagonal ...
  34  *      a.dezero(gmin);
  35  * 5. You probably should set a minimum pivot value ...
  36  *      a.setminpivot(gmin);
  37  * 6. Fill in the matrix ...
  38  *      for (all pairs i,j you want to fill in)
  39  *         a.m(i,j) += data;
  40  * 7. Then do the LU decomposition ...
  41  *      a.lu_decomp();
  42  * 8. Then get the solution by applying the right side vector (v)
  43  * and doing the fwd and back substitution.
  44  *      a.fbsub(v);
  45  * The solution is now in v.
  46  * -------------------------------------------------
  47  * To reuse it .....
  48  * Get rid of the old data ...
  49  *      a.zero().
  50  * Restore the offset, if you want ...
  51  *      a.dezero(gmin);
  52  * Then fill and solve as before. (steps 6-8)
  53  * -------------------------------------------------
  54  * In the above case, LU replaced the matrix and the solution replaced
  55  * the right side, losing the data.
  56  * To keep the matrix, and keep the right side ...
  57  * 1. Declare two matrices ...
  58  *      BSMATRIX<double> a;
  59  *      BSMATRIX<double> lu;
  60  * 2. as before
  61  * 2a. Say the lu has the same structure as a.
  62  *      lu.clone(a);
  63  * 3-5. Allocate both.
  64  * 6. Fill in "a" only.
  65  * 7. Do the LU decomposition, keeping a, with result in lu ...
  66  *      lu.lu_decomp(a, false);
  67  * 8. Do the f&B sub, keeping b (the right side), with result in x ...
  68  *      lu.fbsub(x, b);
  69  * -------------------------------------------------
  70  * To make a change to the matrix and re-solve ...
  71  * Apply the change to a ...
  72  *      for (all changes you want to make)
  73  *         a.m(i,j) += delta;  a.set_changed(i).set_changed(j)
  74  * Then solve again .. step 7 above for a full solution,
  75  * or for an update ...
  76  *      lu.lu_decomp(a, true);
  77  * Then do the same step 8 as above.
  78  * -------------------------------------------------
  79  * some other functions that might be useful ....
  80  * reinit(newsize) -- change the size (and lose the contents and structure)
  81  * size()const -- return the size (# of rows & cols)
  82  * density()const -- return the matrix density, as a number between 0 and 1
  83  * sparsity()const -- 1-density
  84  * -------------------------------------------------
  85  * individual element access ...
  86  * 5 access functions are provided.
  87  * All return lvalues (references to the actual entry).
  88  * All take 2 args, row and column.  They differ in speed and generality.
  89  * Since they are used in equations, they all have 1 letter names.
  90  * s(r,c) -- "safe" -- most general case, with bounds checking
  91  *      if completely out of bounds, returns trash
  92  *      if in the matrix but out of band, returns a reference to zero
  93  *      Changing it will corrupt the matrix.
  94  *      All others have no bounds checking.
  95  * m(r,c) -- "matrix" -- known to be within the allocated band
  96  * u(r,c) -- "upper" -- known to be in the upper triangle. (c>=r)
  97  * l(r,c) -- "lower" -- known to be in the lower triangle. (r>=c)
  98  * d(r,c) -- "diagonal" -- known to be on the diagonal. (r==c)
  99  * Using s() will always work, but will be slower than the other methods.
 100  * You should use the most restricted function that you know is correct.
 101  * -------------------------------------------------
 102  * Other notes ...
 103  * The numbering starts at 1 (Fortran style).
 104  * When using "s" access, it is ok to load row and column 0, then ignore it.
 105  *   This may simplify load functions, at the expense of run speed.
 106  * "s" will let you change the value of zero,
 107  *   but you will find out about it later.
 108  */
 109 //testing=script 2008.09.19
 110 #ifndef M_MATRIX_H
 111 #define M_MATRIX_H
 112 /*--------------------------------------------------------------------------*/
 113 #include "l_stlextra.h"
 114 /*--------------------------------------------------------------------------*/
 115 template <class T>
 116 class BSMATRIX {
 117 private:
 118   mutable bool* _changed;// flag: this node changed value
 119   int*  _lownode;       // lowest node connecting to this one
 120   T*    _space;         // ptr to actual memory space used
 121   T**   _rowptr;        // ptrs to col 0 of every row
 122   T**   _colptr;        // ptrs to row 0 of every col
 123   T**   _diaptr;        // ptrs to diagonal
 124   int   _nzcount;       // count of non-zero elements
 125   int   _size;          // # of rows and columns
 126   T     _zero;          // always 0 but not const
 127   T     _trash;         // depository for row and col 0, write only
 128   T     _min_pivot;     // minimum pivot value
 129 private:
 130   explicit      BSMATRIX(const BSMATRIX<T>&) {incomplete();unreachable();}
 131   void          uninit();
 132   void          init(int s=0);
 133   T&            subtract_dot_product(int r, int c, int d);
 134   T&            subtract_dot_product(int r, int c, int d, const T& in);
 135   int           lownode(int i)const     {return _lownode[i];}
 136   bool          is_changed(int n)const  {return _changed[n];}
 137   void          set_changed(int n, bool x = true)const {_changed[n] = x;}
 138 public:
 139   explicit      BSMATRIX(int ss=0);
 140                 ~BSMATRIX()             {uninit();}
 141   void          reinit(int ss=0)        {uninit(); init(ss);}
 142   //void        clone(const BSMATRIX<T>&);
 143   void          iwant(int, int);
 144   void          unallocate();
 145   void          allocate();
 146   void          reallocate()            {unallocate(); allocate();}
 147   void          set_min_pivot(double x) {_min_pivot = x;}
 148   void          zero();
 149   void          dezero(T& o);
 150   int           size()const             {return _size;}
 151   double        density();
 152   T     d(int r, int  )const    {return *(_diaptr[r]);}
 153 private:
 154   T     u(int r, int c)const    {return _colptr[c][r];}
 155   T     l(int r, int c)const    {return _rowptr[r][-c];}
 156   T&    d(int r, int c);
 157   T&    u(int r, int c);
 158   T&    l(int r, int c);
 159   T&    m(int r, int c);
 160   //T&  s(int r, int c);
 161 public:
 162   void          load_diagonal_point(int i, T value);
 163   void          load_point(int i, int j, T value);
 164   void          load_couple(int i, int j, T value);
 165   void          load_symmetric(int i, int j, T value);
 166   void          load_asymmetric(int r1, int r2, int c1, int c2, T value);
 167
 168   void          lu_decomp(const BSMATRIX<T>&, bool do_partial);
 169   void          lu_decomp();
 170   void          fbsub(T* v) const;
 171   void          fbsub(T* x, const T* b, T* c = NULL) const;
 172 };
 173 /*--------------------------------------------------------------------------*/
 174 // private implementations
 175 /*--------------------------------------------------------------------------*/
 176 template <class T>
 177 void BSMATRIX<T>::uninit()
 178 {
 179   unallocate();
 180   delete [] _lownode;
 181   _lownode = NULL;
 182   delete [] _changed;
 183   _changed = NULL;
 184 }
 185 /*--------------------------------------------------------------------------*/
 186 template <class T>
 187 void BSMATRIX<T>::init(int ss)
 188 {
 189   assert(!_lownode);
 190   assert(!_colptr);
 191   assert(!_rowptr);
 192   assert(!_diaptr);
 193   assert(!_space);
 194
 195   assert(_zero == 0.);
 196   _min_pivot = _trash = 0.;
 197   _nzcount = 0;
 198   _size = ss;
 199   _lownode = new int[size()+1];
 200   assert(_lownode);
 201   for (int ii = 0;  ii <= size();  ++ii) {
 202     _lownode[ii] = ii;
 203   }
 204   _changed = new bool[size()+1];
 205   assert(_changed);
 206   for (int ii = 0;  ii <= size();  ++ii) {
 207     set_changed(ii, false);
 208   }
 209 }
 210 /*--------------------------------------------------------------------------*/
 211 template <class T>
 212 T& BSMATRIX<T>::subtract_dot_product(int rr, int cc, int dd)
 213 {
 214   assert(_lownode);
 215   int kk = std::max(_lownode[rr], _lownode[cc]);
 216   int len = dd - kk;
 217   T& dot = m(rr, cc);
 218   if (len > 0) {
 219     T* row = &(l(rr,kk));
 220     T* col = &(u(kk,cc));
 221     /* for (ii = kk;   ii < dd;   ++ii) */
 222     for (int ii = 0;   ii < len;   ++ii) {
 223       dot -= row[-ii] * col[ii];
 224     }
 225   }else{
 226   }
 227   return dot;
 228 }
 229 /*--------------------------------------------------------------------------*/
 230 template <class T>
 231 T& BSMATRIX<T>::subtract_dot_product(int rr, int cc, int dd, const T& in)
 232 {
 233   assert(_lownode);
 234   int kk = std::max(_lownode[rr], _lownode[cc]);
 235   int len = dd - kk;
 236   T& dot = m(rr, cc);
 237   dot = in;
 238   if (len > 0) {
 239     T* row = &(l(rr,kk));
 240     T* col = &(u(kk,cc));
 241     /* for (ii = kk;   ii < dd;   ++ii) */
 242     for (int ii = 0;   ii < len;   ++ii) {
 243       dot -= row[-ii] * col[ii];
 244     }
 245   }else{
 246   }
 247   return dot;
 248 }
 249 /*--------------------------------------------------------------------------*/
 250 // public implementations
 251 /*--------------------------------------------------------------------------*/
 252 template <class T>
 253 BSMATRIX<T>::BSMATRIX(int ss)
 254   :_changed(NULL),
 255    _lownode(NULL),
 256    _space(NULL),
 257    _rowptr(NULL),
 258    _colptr(NULL),
 259    _diaptr(NULL),
 260    _nzcount(0),
 261    _size(ss),
 262    _zero(0.),
 263    _trash(0.),
 264    _min_pivot(0.)
 265 {
 266   init(ss);
 267 }
 268 /*--------------------------------------------------------------------------*/
 269 #if 0
 270 /* clone: copy to self the structure of another BSMATRIX
 271  * this does not copy the values stored in the matrix
 272  */
 273 template <class T>
 274 void BSMATRIX<T>::clone(const BSMATRIX<T> & aa)
 275 {untested();
 276   reinit(aa.size());
 277   for (int ii = 0;  ii <= size();  ++ii) {untested();
 278     _lownode[ii] = aa.lownode(ii);
 279   }
 280 }
 281 #endif
 282 /*--------------------------------------------------------------------------*/
 283 /* iwant: indicate that "iwant" to allocate this spot in the matrix
 284  */
 285 template <class T>
 286 void BSMATRIX<T>::iwant(int node1, int node2)
 287 {
 288   assert(_lownode);
 289   assert(node1 <= size());
 290   assert(node2 <= size());
 291
 292   if (node1 <= 0  ||  node2 <= 0) {
 293     // node 0 is ground, and doesn't count as a connection
 294     // negative is invalid, not used but still may be in a node list
 295   }else if (node1 < _lownode[node2]) {
 296     _lownode[node2]=node1;
 297   }else if (node2 < _lownode[node1]) {
 298     _lownode[node1]=node2;
 299   }else{
 300   }
 301 }
 302 /*--------------------------------------------------------------------------*/
 303 template <class T>
 304 void BSMATRIX<T>::unallocate()
 305 {
 306   assert (_zero == 0.);
 307   delete [] _rowptr;
 308   delete [] _colptr;
 309   delete [] _diaptr;
 310   delete [] _space;
 311
 312   _rowptr = _colptr = _diaptr = NULL;
 313   _space = NULL;
 314 }
 315 /*--------------------------------------------------------------------------*/
 316 /* allocate: really get the space to work
 317  */
 318 template <class T>
 319 void BSMATRIX<T>::allocate()
 320 {
 321   assert(_lownode);
 322   assert(!_colptr);
 323   assert(!_rowptr);
 324   assert(!_diaptr);
 325   assert(!_space);
 326
 327   _nzcount = 0;
 328   for (int ii = 0;   ii <= size();   ++ii) {
 329     _nzcount += 2 * (ii - _lownode[ii]) + 1;
 330   }
 331
 332   _colptr = new T*[size()+1];
 333   _rowptr = new T*[size()+1];
 334   _diaptr = new T*[size()+1];
 335   _space  = new T[_nzcount];
 336
 337   assert(_colptr);
 338   assert(_rowptr);
 339   assert(_diaptr);
 340
 341   zero();
 342
 343   {
 344     T* point = _space;
 345     for (int ii = 0;   ii <= size();   ++ii) {
 346       _colptr[ii] = point - _lownode[ii];
 347       _rowptr[ii] = _colptr[ii] + 2*ii;
 348       _diaptr[ii] = _colptr[ii] + ii;
 349       point += 2 * (ii - _lownode[ii]) + 1;
 350     }
 351   }
 352 }
 353 /*--------------------------------------------------------------------------*/
 354 /* zero: wipe the whole array
 355  */
 356 template <class T>
 357 void BSMATRIX<T>::zero()
 358 {
 359   assert(_space);
 360   assert(_zero == 0.);
 361   _trash = 0.;
 362   std::fill_n(_space, _nzcount, 0.);
 363 }
 364 /*--------------------------------------------------------------------------*/
 365 /* dezero: make sure(?) the diagonal is non-zero
 366  */
 367 template <class T>
 368 void BSMATRIX<T>::dezero(T& offset)
 369 {
 370   for (int ii = 1;  ii <= size();  ++ii) {
 371     d(ii,ii) += offset;
 372   }
 373 }
 374 /*--------------------------------------------------------------------------*/
 375 template <class T>
 376 double BSMATRIX<T>::density()
 377 {
 378   if (size() > 0) {
 379     assert(_lownode);
 380     _nzcount = 0;
 381     for (int ii = 0;   ii <= size();   ++ii) {
 382       _nzcount += 2 * (ii - _lownode[ii]) + 1;
 383     }
 384     return static_cast<double>(_nzcount-1)/(static_cast<double>(size())*size());
 385   }else{
 386     return 0;
 387   }
 388 }
 389 /*--------------------------------------------------------------------------*/
 390 /* d: fast matrix entry access, lvalue
 391  * It is known that the entry is valid and on the diagonal
 392  */
 393 template <class T>
 394 T& BSMATRIX<T>::d(int r, int c)
 395 {
 396   assert(_diaptr);
 397   assert(r == c);
 398   assert(0 <= r);
 399   assert(r <= _size);
 400
 401   return *(_diaptr[r]);
 402 }
 403 /*--------------------------------------------------------------------------*/
 404 /* u: fast matrix entry access
 405  * It is known that the entry is valid and in the upper triangle
 406  */
 407 template <class T>
 408 T& BSMATRIX<T>::u(int r, int c)
 409 {
 410   assert(_colptr);
 411   assert(_lownode);
 412   assert(0 < r);
 413   assert(r <= c);
 414   assert(c <= _size);
 415   assert(1 <= _lownode[c]);
 416   assert(_lownode[c] <= r);
 417
 418   return _colptr[c][r];
 419 }
 420 /*--------------------------------------------------------------------------*/
 421 /* l: fast matrix entry access
 422  * It is known that the entry is valid and in the lower triangle
 423  */
 424 template <class T>
 425 T& BSMATRIX<T>::l(int r, int c)
 426 {
 427   assert(_rowptr);
 428   assert(_lownode);
 429   assert(0 < c);
 430   assert(c <= r);
 431   assert(r <= _size);
 432   assert(1 <= _lownode[r]);
 433   assert(_lownode[r] <= c);
 434
 435   return _rowptr[r][-c];
 436 }
 437 /*--------------------------------------------------------------------------*/
 438 /* m: semi-fast matrix entry access
 439  * It is known that the entry refers to a valid location
 440  * but it is not known whether lower, upper, or diagonal
 441  */
 442 template <class T>
 443 T& BSMATRIX<T>::m(int r, int c)
 444 {
 445   return (c>=r) ? u(r,c) : l(r,c);
 446 }
 447 /*--------------------------------------------------------------------------*/
 448 /* s: general matrix entry access.
 449  * It is known that the location is strictly in bounds,
 450  *   but it is not known whether the location actually exists.
 451  * If access is attempted to a non-allocated location,
 452  *   it returns a reference to a shared zero variable.
 453  *   Writing to this zero is not prohibited,
 454  *   but will corrupt the matrix in a known and testable way.
 455  * If access is attempted to row 0 or column 0,
 456  *   it returns a reference to a shared trash variable.
 457  *   Writing to trash is allowed and encouraged,
 458  *   but reading it gives a number not useful for anything.
 459  */
 460 #if 0
 461 template <class T>
 462 T& BSMATRIX<T>::s(int row, int col)
 463 {untested();
 464   assert(_lownode);
 465   assert(0 <= col);
 466   assert(col <= size());
 467   assert(0 <= row);
 468   assert(row <= size());
 469   assert(_zero == 0.);
 470
 471   if (col == row) {untested();
 472     return d(row,col);
 473   }else if (col > row) {untested();     /* above the diagonal */
 474     if (row == 0) {untested();
 475       return _trash;
 476     }else if (row < _lownode[col]) {untested();
 477       return _zero;
 478     }else{untested();
 479       return u(row,col);
 480     }
 481   }else{untested();                     /* below the diagonal */
 482     assert(col < row);
 483     if (col == 0) {untested();
 484       return _trash;
 485     }else if (col < _lownode[row]) {untested();
 486       return _zero;
 487     }else{untested();
 488       return l(row,col);
 489     }
 490   }
 491   unreachable();
 492 }
 493 #endif
 494 /*--------------------------------------------------------------------------*/
 495 template <class T>
 496 void BSMATRIX<T>::load_point(int i, int j, T value)
 497 {
 498   if (i > 0 && j > 0) {
 499     set_changed(j);
 500     set_changed(i);
 501     m(i,j) += value;
 502   }else{
 503   }
 504 }
 505 /*--------------------------------------------------------------------------*/
 506 // load_point(i, i, value);
 507 template <class T>
 508 void BSMATRIX<T>::load_diagonal_point(int i, T value)
 509 {
 510   if (i > 0) {
 511     set_changed(i);
 512     d(i,i) += value;
 513   }else{untested();
 514   }
 515 }
 516 /*--------------------------------------------------------------------------*/
 517 // load_point(i, j, -value);
 518 // load_point(j, i, -value);
 519 template <class T>
 520 void BSMATRIX<T>::load_couple(int i, int j, T value)
 521 {
 522   if (j > 0) {
 523     set_changed(j);
 524     if (i > 0) {
 525       set_changed(i);
 526       m(i,j) -= value;
 527       m(j,i) -= value;
 528     }else{
 529     }
 530   }else{untested();
 531   }
 532 }
 533 /*--------------------------------------------------------------------------*/
 534 // load_point(i, i, value); or load_diagonal_point(i, value);
 535 // load_point(j, j, value); or load_diagonal_point(j, value);
 536 // load_point(i, j, -value);
 537 // load_point(j, i, -value);
 538 template <class T>
 539 void BSMATRIX<T>::load_symmetric(int i, int j, T value)
 540 {
 541   if (j > 0) {
 542     set_changed(j);
 543     d(j,j) += value;
 544     if (i > 0) {
 545       set_changed(i);
 546       d(i,i) += value;
 547       m(i,j) -= value;
 548       m(j,i) -= value;
 549     }else{
 550     }
 551   }else if (i > 0) {
 552     set_changed(i);
 553     d(i,i) += value;
 554   }else{
 555   }
 556 }
 557 /*--------------------------------------------------------------------------*/
 558 // load_point(r1, c1, value);
 559 // load_point(r2, c2, value);
 560 // load_point(r1, c2, -value);
 561 // load_point(r2, c1, -value);
 562 template <class T>
 563 void BSMATRIX<T>::load_asymmetric(int r1,int r2,int c1,int c2,T value)
 564 {
 565   set_changed(c1);
 566   set_changed(c2);
 567   if (r1 > 0) {
 568     set_changed(r1);
 569     if (c1 > 0) {
 570       m(r1,c1) += value;
 571     }else{
 572     }
 573     if (c2 > 0) {
 574       m(r1,c2) -= value;
 575     }else{
 576     }
 577   }else{
 578   }
 579   if (r2 > 0) {
 580     set_changed(r2);
 581     if (c1 > 0) {
 582       m(r2,c1) -= value;
 583     }else{
 584     }
 585     if (c2 > 0) {
 586       m(r2,c2) += value;
 587     }else{
 588     }
 589   }else{
 590   }
 591 }
 592 /*--------------------------------------------------------------------------*/
 593 template <class T>
 594 void BSMATRIX<T>::lu_decomp(const BSMATRIX<T>& aa, bool do_partial)
 595 {
 596   int prop = 0;   /* change propagation indicator */
 597   assert(_lownode);
 598   assert(aa._lownode);
 599   assert(aa.size() == size());
 600   for (int mm = 1;   mm <= size();   ++mm) {
 601     assert(aa.lownode(mm) == _lownode[mm]);
 602     int bn = _lownode[mm];
 603     if (!do_partial  ||  aa.is_changed(mm)  ||  bn <= prop) {
 604       aa.set_changed(mm, false);
 605       prop = mm;
 606       if (bn < mm) {
 607         prop = mm;
 608         u(bn,mm) = aa.u(bn,mm) / d(bn,bn);
 609         for (int ii = bn+1;  ii<mm;  ii++) {
 610           /* u(ii,mm) = (aa.u(ii,mm) - dot(ii,mm,ii)) / d(ii,ii); */
 611           subtract_dot_product(ii,mm,ii,aa.u(ii,mm)) /= d(ii,ii);
 612         }
 613         l(mm,bn) = aa.l(mm,bn);
 614         for (int jj = bn+1;  jj<mm;  jj++) {
 615           /* l(mm,jj) = aa.l(mm,jj) - dot(mm,jj,jj); */
 616           subtract_dot_product(mm,jj,jj,aa.l(mm,jj));
 617         }
 618         { /* jj == mm */
 619           /* d(mm,mm) = aa.d(mm,mm) - dot(mm,mm,mm); then test */
 620           if (subtract_dot_product(mm,mm,mm,aa.d(mm,mm)) == 0.) {untested();
 621             error(bWARNING, "open circuit: internal node %u\n", mm);
 622             d(mm,mm) = _min_pivot;
 623           }else{
 624           }
 625         }
 626       }else{    /* bn == mm */
 627         d(mm,mm) = aa.d(mm,mm);
 628         if (d(mm,mm)==0.) {untested();
 629           d(mm,mm) = _min_pivot;
 630         }else{
 631         }
 632       }
 633     }else{
 634     }
 635   }
 636 }
 637 /*--------------------------------------------------------------------------*/
 638 template <class T>
 639 void BSMATRIX<T>::lu_decomp()
 640 {
 641   assert(_lownode);
 642   for (int mm = 1;   mm <= size();   ++mm) {
 643     int bn = _lownode[mm];
 644     if (bn < mm) {
 645       u(bn,mm) /= d(bn,bn);
 646       for (int ii =bn+1;  ii<mm;  ii++) {
 647         /* (m(ii,mm) -= dot(ii,mm,ii)) /= d(ii,ii); */
 648         subtract_dot_product(ii,mm,ii) /= d(ii,ii);
 649       }
 650       for (int jj = bn+1;  jj<mm;  jj++) {
 651         /* m(mm,jj) -= dot(mm,jj,jj); */
 652         subtract_dot_product(mm,jj,jj);
 653       }
 654       { /* jj == mm */
 655         /* m(mm,mm) -= dot(mm,mm,mm); then test */
 656         if (subtract_dot_product(mm,mm,mm) == 0.) {untested();
 657           error(bWARNING, "open circuit: internal node %u\n", mm);
 658           d(mm,mm) = _min_pivot;
 659         }else{
 660         }
 661       }
 662     }else{    /* bn == mm */
 663       if (d(mm,mm)==0.) {
 664         d(mm,mm) = _min_pivot;
 665       }else{
 666       }
 667     }
 668   }
 669 }
 670 /*--------------------------------------------------------------------------*/
 671 /* fbsub: forward and back sub, shared storage
 672  * v = right side vector, changed in place to solution vector
 673  */
 674 template <class T>
 675 void BSMATRIX<T>::fbsub(T* v) const
 676 {
 677   assert(_lownode);
 678   assert(v);
 679
 680   for (int ii = 1; ii <= size(); ++ii) {        /* forward substitution */
 681     for (int jj = _lownode[ii]; jj < ii; ++jj) {
 682       v[ii] -= l(ii,jj) * v[jj];
 683     }
 684     v[ii] /= d(ii,ii);
 685   }
 686
 687   for (int jj = size(); jj > 1; --jj) {         /* back substitution    */
 688     for (int ii = _lownode[jj]; ii < jj; ++ii) {
 689       v[ii] -= u(ii,jj) * v[jj];
 690     }
 691   }
 692 }
 693 /*--------------------------------------------------------------------------*/
 694 /* fbsub: forward and back sub, separate storage
 695  * b = right side vector
 696  * c = intermediate vector after fwd sub
 697  * x = solution vector
 698  */
 699 template <class T>
 700 void BSMATRIX<T>::fbsub(T* x, const T* b, T* c) const
 701 {
 702   assert(_lownode);
 703   assert(x);
 704   assert(b);
 705   assert(c);
 706
 707   {
 708     int ii = 1;
 709     for (   ; ii <= size(); ++ii) {
 710       if (b[ii] != 0.) {
 711         break;
 712       }else{
 713       }
 714       c[ii] = 0.;
 715     }
 716
 717     int first_nz = ii;
 718     for (   ; ii <= size(); ++ii) {             /* forward substitution */
 719       int low_node = std::max(_lownode[ii], first_nz);
 720       c[ii] = b[ii];
 721       for (int jj = low_node; jj < ii; ++jj) {
 722         c[ii] -= l(ii,jj) * c[jj];
 723       }
 724       c[ii] /= d(ii,ii);
 725     }
 726   }
 727
 728   notstd::copy_n(c, size()+1, x);
 729
 730   for (int jj = size(); jj > 1; --jj) {         /* back substitution    */
 731     for (int ii = _lownode[jj]; ii < jj; ++ii) {
 732       x[ii] -= u(ii,jj) * x[jj];
 733     }
 734   }
 735   x[0] = 0.;
 736   //BUG// some things don't work unless there is a zero here.
 737 }
 738 /*--------------------------------------------------------------------------*/
 739 /*--------------------------------------------------------------------------*/
 740 #endif
 741
 742 // vim:ts=8:sw=2:noet: