Qucs-core
0.0.19
|
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