Qucs-core  0.0.19
tridiag.cpp
Go to the documentation of this file.
00001 /*
00002  * tridiag.cpp - tridiagonal matrix template class implementation
00003  *
00004  * Copyright (C) 2005, 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 <assert.h>
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include <string.h>
00033 #include <cmath>
00034 #include <vector>
00035 
00036 #include "compat.h"
00037 #include "complex.h"
00038 //#include "tvector.h"
00039 
00040 namespace qucs {
00041 
00042 // Constructor creates an instance of the tridiag class.
00043 template <class nr_type_t>
00044 tridiag<nr_type_t>::tridiag () {
00045   abov = belo = diag = offdiag = rhs = NULL;
00046   type = TRIDIAG_UNKNOWN;
00047 }
00048 
00049 /* The copy constructor creates a new instance based on the given
00050    tridiag object. */
00051 template <class nr_type_t>
00052 tridiag<nr_type_t>::tridiag (const tridiag & t) {
00053   abov = t.abov;
00054   belo = t.belo;
00055   diag = t.diag;
00056   offdiag = t.offdiag;
00057   rhs = t.rhs;
00058   type = t.type;
00059 }
00060 
00061 /* The assignment copy constructor creates a new instance based on the
00062    given tridiag object. */
00063 template <class nr_type_t>
00064 const tridiag<nr_type_t>&
00065 tridiag<nr_type_t>::operator=(const tridiag<nr_type_t> & t) {
00066   if (&t != this) {
00067     abov = t.abov;
00068     belo = t.belo;
00069     diag = t.diag;
00070     offdiag = t.offdiag;
00071     rhs = t.rhs;
00072     type = t.type;
00073   }
00074   return *this;
00075 }
00076 
00077 // Destructor deletes a tridiag object.
00078 template <class nr_type_t>
00079 tridiag<nr_type_t>::~tridiag () {
00080 }
00081 
00082 // Set the diagonal vector of the tridiagonal matrix.
00083 template <class nr_type_t>
00084 void tridiag<nr_type_t>::setDiagonal (::std::vector<nr_type_t> * v) {
00085   diag = v;
00086 }
00087 
00088 // Set the off-diagonal vector of the tridiagonal matrix.
00089 template <class nr_type_t>
00090 void tridiag<nr_type_t>::setOffDiagonal (::std::vector<nr_type_t> * v) {
00091   abov = belo = offdiag = v;
00092 }
00093 
00094 // Set the above off-diagonal vector of the tridiagonal matrix.
00095 template <class nr_type_t>
00096 void tridiag<nr_type_t>::setA (::std::vector<nr_type_t> * v) {
00097   abov = v;
00098 }
00099 
00100 // Set the below off-diagonal vector of the tridiagonal matrix.
00101 template <class nr_type_t>
00102 void tridiag<nr_type_t>::setB (::std::vector<nr_type_t> * v) {
00103   belo = v;
00104 }
00105 
00106 // Set the right hand side vector of the equation system.
00107 template <class nr_type_t>
00108 void tridiag<nr_type_t>::setRHS (::std::vector<nr_type_t> * v) {
00109   rhs = v;
00110 }
00111 
00112 /* Depending on the type of tridiagonal matrix the function runs one
00113    of the solvers specialized on this type.  The solvers are in-place
00114    algorithms meaning the right hand side is replaced by the solution
00115    and the original matrix entries are destroyed during solving the
00116    system. */
00117 template <class nr_type_t>
00118 void tridiag<nr_type_t>::solve (void) {
00119   switch (type) {
00120   case TRIDIAG_NONSYM:
00121     solve_ns (); break;
00122   case TRIDIAG_SYM:
00123     solve_s (); break;
00124   case TRIDIAG_NONSYM_CYCLIC:
00125     solve_ns_cyc (); break;
00126   case TRIDIAG_SYM_CYCLIC:
00127     solve_s_cyc (); break;
00128   }
00129 }
00130 
00131 /* This function solves a tridiagonal equation system given
00132    by
00133        diag[0]  abov[0]        0   .....            0
00134        belo[1]  diag[1]  abov[1]   .....            0
00135              0  belo[2]  diag[2]
00136              0        0  belo[3]   .....    abov[n-2]
00137            ...         ...
00138              0         ...       belo[n-1]  diag[n-1]
00139 */
00140 template <class nr_type_t>
00141 void tridiag<nr_type_t>::solve_ns (void) {
00142   d = al = &diag->front(); //diag->getData ();
00143   f = ga = &abov->front (); //abov->getData ();
00144   e = &belo->front(); // belo->getData ();
00145   b = c = x = &rhs->front(); //rhs->getData ();
00146   int i, n = (int)diag->size(); //diag->getSize ();
00147 
00148   // factorize A = LU
00149   al[0] = d[0];
00150   ga[0] = f[0] / al[0];
00151   for (i = 1; i < n - 1; i++) {
00152     al[i] = d[i] - e[i] * ga[i-1];
00153     ga[i] = f[i] / al[i];
00154   }
00155   al[n-1] = d[n-1] - e[n-1] * ga[n-2];
00156 
00157   // update b = Lc
00158   c[0] = b[0] / d[0];
00159   for (i = 1; i < n; i++) {
00160     c[i] = (b[i] - e[i] * c[i-1]) / al[i];
00161   }
00162 
00163   // back substitution Rx = c
00164   x[n-1] = c[n-1];
00165   for (i = n - 2; i >= 0; i--) {
00166     x[i] = c[i] - ga[i] * x[i+1];
00167   }
00168 }
00169 
00170 /* This function solves a cyclically tridiagonal equation system given
00171    by
00172        diag[0]  abov[0]        0   .....      belo[0]
00173        belo[1]  diag[1]  abov[1]   .....           0
00174              0  belo[2]  diag[2]
00175              0        0  belo[3]   .....    abov[n-2]
00176            ...         ...
00177       abov[n-1]        ...       belo[n-1]  diag[n-1]
00178 */
00179 template <class nr_type_t>
00180 void tridiag<nr_type_t>::solve_ns_cyc (void) {
00181   d = al = &diag->front(); //diag->getData ();
00182   f = ga = &abov->front (); //abov->getData ();
00183   e = be = &belo->front(); // belo->getData ();
00184   b = x = c = &rhs->front(); //rhs->getData ();
00185   int i, n = (int)diag->size(); //diag->getSize ();
00186   de = new nr_type_t[n];
00187   ep = new nr_type_t[n];
00188 
00189   // factorize A = LU
00190   al[0] = d[0];
00191   ga[0] = f[0] / al[0];
00192   de[0] = e[0] / al[0];
00193   for (i = 1; i < n - 2; i++) {
00194     al[i] = d[i] - e[i] * ga[i-1];
00195     ga[i] = f[i] / al[i];
00196     be[i] = e[i];
00197     de[i] = -be[i] * de[i-1] / al[i];
00198   }
00199   al[n-2] = d[n-2] - e[n-2] * ga[n-3];
00200   be[n-2] = e[n-2];
00201   ep[2] = f[n-1];
00202   for (i = 3; i < n; i++) {
00203     ep[i] = -ep[i-1] * ga[i-3];
00204   }
00205   ga[n-2] = (f[n-2] - be[n-2] * de[n-3]) / al[n-2];
00206   be[n-1] = e[n-1] - ep[n-1] * ga[n-3];
00207   al[n-1] = d[n-1] - be[n-1] * ga[n-2];
00208   for (i = 2; i < n; i++) {
00209     al[n-1] -= ep[i] * de[i-2];
00210   }
00211 
00212   // update Lc = b
00213   c[0] = b[0] / al[0];
00214   for (i = 1; i < n - 1; i++) {
00215     c[i] = (b[i] - c[i-1] * be[i]) / al[i];
00216   }
00217   c[n-1] = b[n-1] - be[n-1] * c[n-2];
00218   for (i = 2; i < n; i++) {
00219     c[n-1] -= ep[i] * c[i-2];
00220   }
00221   c[n-1] /= al[n-1];
00222 
00223   // back substitution Rx = c
00224   x[n-1] = c[n-1];
00225   x[n-2] = c[n-2] - ga[n-2] * x[n-1];
00226   for (i = n - 3; i >= 0; i--) {
00227     x[i] = c[i] - ga[i] * x[i+1] - de[i] * x[n-1];
00228   }
00229 
00230   delete[] de;
00231   delete[] ep;
00232 }
00233 
00234 /* This function solves a symmetric tridiagonal strongly nonsingular
00235    equation system given by
00236        diag[0]  offdiag[0]             0   .....                   0
00237     offdiag[0]     diag[1]    offdiag[1]   .....                   0
00238              0  offdiag[1]       diag[2]
00239              0           0    offdiag[2]   .....        offdiag[n-2]
00240            ...         ...
00241              0         ...                offdiag[n-2]     diag[n-1]
00242 */
00243 template <class nr_type_t>
00244 void tridiag<nr_type_t>::solve_s (void) {
00245   d = al = &diag->front(); //diag->getData ();
00246   f = ga = &offdiag->front (); //offdiag->getData ();
00247   b = z = x = c = &rhs->front(); //rhs->getData ();
00248   nr_type_t t;
00249   int i, n = (int)diag->size(); //diag->getSize ();
00250   de = new nr_type_t[n];
00251 
00252   // factorize A = LDL'
00253   al[0] = d[0];
00254   t = f[0];
00255   ga[0] = t / al[0];
00256   for (i = 1; i < n - 1; i++) {
00257     al[i] = d[i] - t * ga[i-1];
00258     t = f[i];
00259     ga[i] = t / al[i];
00260   }
00261   al[n-1] = d[n-1] - t * ga[n-2];
00262 
00263   // update L'z = b and Dc = z
00264   z[0] = b[0];
00265   for (i = 1; i < n; i++) {
00266     z[i] = b[i] - ga[i-1] * z[i-1];
00267   }
00268   for (i = 0; i < n; i++) {
00269     c[i] = z[i] / al[i];
00270   }
00271 
00272   // back substitution L'x = c
00273   x[n-1] = c[n-1];
00274   for (i = n-2; i >= 0; i--) {
00275     x[i] = c[i] - ga[i] * x[i+1];
00276   }
00277 
00278   delete[] de;
00279 }
00280 
00281 /* This function solves a symmetric cyclically tridiagonal strongly
00282    nonsingular equation system given by
00283        diag[0]  offdiag[0]             0   .....        offdiag[n-1]
00284     offdiag[0]     diag[1]    offdiag[1]   .....                   0
00285              0  offdiag[1]       diag[2]
00286              0           0    offdiag[2]   .....        offdiag[n-2]
00287            ...         ...
00288   offdiag[n-1]         ...                offdiag[n-2]     diag[n-1]
00289 */
00290 template <class nr_type_t>
00291 void tridiag<nr_type_t>::solve_s_cyc (void) {
00292   d = al = &diag->front(); //diag->getData ();
00293   f = ga = &offdiag->front (); //offdiag->getData ();
00294   b = c = z  = x = &rhs->front(); //rhs->getData ();
00295   nr_type_t t;
00296   int i, n = (int)diag->size(); //diag->getSize ();
00297   de = new nr_type_t[n];
00298 
00299   // factorize A = LDL'
00300   al[0] = d[0];
00301   t = f[0];
00302   ga[0] = t / al[0];
00303   de[0] = f[n-1] / al[0];
00304   for (i = 1; i < n - 2; i++) {
00305     al[i] = d[i] - t * ga[i-1];
00306     de[i] = -de[i-1] * t / al[i];
00307     t = f[i];
00308     ga[i] = t / al[i];
00309   }
00310   al[n-2] = d[n-2] - t * ga[n-3];
00311   ga[n-2] = (f[n-2] - t * de[n-3]) / al[n-2];
00312   al[n-1] = d[n-1] - al[n-2] * ga[n-2] * ga[n-2];
00313   for (i = 0; i < n - 2; i++) {
00314     al[n-1] -= al[i] * de[i] * de[i];
00315   }
00316 
00317   // update Lz = b and Dc = z
00318   z[0] = b[0];
00319   for (i = 1; i < n - 1; i++) {
00320     z[i] = b[i] - ga[i-1] * z[i-1];
00321   }
00322   z[n-1] = b[n-1] - ga[n-2] * z[n-2];
00323   for (i = 0; i < n - 2; i++) {
00324     z[n-1] -= de[i] * z[i];
00325   }
00326   for (i = 0; i < n; i++) {
00327     c[i] = z[i] / al[i];
00328   }
00329 
00330   // back substitution L'x = c
00331   x[n-1] = c[n-1];
00332   x[n-2] = c[n-2] - ga[n-2] * x[n-1];
00333   for (i = n - 3; i >= 0; i--) {
00334     x[i] = c[i] - ga[i] * x[i+1] - de[i] * x[n-1];
00335   }
00336 
00337   delete[] de;
00338 }
00339 
00340 } // namespace qucs