00001 /* 00002 ------------------------------------------------------------------- 00003 00004 Copyright (C) 2006, 2007, 2008, Andrew W. Steiner 00005 00006 This file is part of O2scl. It has been modified from the 00007 original GSL version written by Gerard Jungman, 00008 Brian Gough, and David Necas. 00009 00010 O2scl is free software; you can redistribute it and/or modify 00011 it under the terms of the GNU General Public License as published by 00012 the Free Software Foundation; either version 3 of the License, or 00013 (at your option) any later version. 00014 00015 O2scl is distributed in the hope that it will be useful, 00016 but WITHOUT ANY WARRANTY; without even the implied warranty of 00017 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00018 GNU General Public License for more details. 00019 00020 You should have received a copy of the GNU General Public License 00021 along with O2scl. If not, see <http://www.gnu.org/licenses/>. 00022 00023 ------------------------------------------------------------------- 00024 */ 00025 00026 /** \file tridiag_base.h 00027 \brief File for solving tridiagonal systems 00028 */ 00029 00030 #ifdef DOXYGENP 00031 namespace o2scl { 00032 #endif 00033 00034 /** 00035 \brief Solve a symmetric tridiagonal linear system 00036 00037 \future Convert into class for memory managment and combine 00038 with other functions below 00039 00040 For description of method see [Engeln-Mullges + Uhlig, p. 92] 00041 * 00042 * diag[0] offdiag[0] 0 ..... 00043 * offdiag[0] diag[1] offdiag[1] ..... 00044 * 0 offdiag[1] diag[2] 00045 * 0 0 offdiag[2] ..... 00046 00047 */ 00048 template<class vec_t, class vec2_t> 00049 int solve_tridiag_sym(const vec_t &diag, const vec2_t &offdiag, 00050 const vec_t &b, vec_t &x, size_t N) { 00051 int status; 00052 double *gamma = (double *) malloc (N * sizeof (double)); 00053 double *alpha = (double *) malloc (N * sizeof (double)); 00054 double *c = (double *) malloc (N * sizeof (double)); 00055 double *z = (double *) malloc (N * sizeof (double)); 00056 00057 if (gamma == 0 || alpha == 0 || c == 0 || z == 0) { 00058 status = GSL_ENOMEM; 00059 } else { 00060 size_t i, j; 00061 00062 /* Cholesky decomposition 00063 A = L.D.L^t 00064 lower_diag(L) = gamma 00065 diag(D) = alpha 00066 */ 00067 alpha[0] = O2SCL_IX(diag,0); 00068 gamma[0] = O2SCL_IX(offdiag,0) / alpha[0]; 00069 00070 for (i = 1; i < N - 1; i++) { 00071 alpha[i] = O2SCL_IX(diag,i) - O2SCL_IX(offdiag,i-1) * gamma[i - 1]; 00072 gamma[i] = O2SCL_IX(offdiag,i) / alpha[i]; 00073 } 00074 00075 if (N > 1) { 00076 alpha[N-1]=O2SCL_IX(diag,N-1)-O2SCL_IX(offdiag,N-2)*gamma[N-2]; 00077 } 00078 00079 /* update RHS */ 00080 z[0] = b[0]; 00081 for (i = 1; i < N; i++) { 00082 z[i] = O2SCL_IX(b,i) - gamma[i - 1] * z[i - 1]; 00083 } for (i = 0; i < N; i++) { 00084 c[i] = z[i] / alpha[i]; 00085 } 00086 00087 /* backsubstitution */ 00088 O2SCL_IX(x,N-1)=c[N-1]; 00089 if (N >= 2) { 00090 for (i = N - 2, j = 0; j <= N - 2; j++, i--) { 00091 O2SCL_IX(x,i) = c[i] - gamma[i] * O2SCL_IX(x,i+1); 00092 } 00093 } 00094 00095 status = GSL_SUCCESS; 00096 } 00097 00098 if (z != 0) 00099 free (z); 00100 if (c != 0) 00101 free (c); 00102 if (alpha != 0) 00103 free (alpha); 00104 if (gamma != 0) 00105 free (gamma); 00106 00107 return status; 00108 } 00109 00110 /** 00111 \brief Solve an asymmetric tridiagonal linear system 00112 00113 plain gauss elimination, only not bothering with the zeroes 00114 * 00115 * diag[0] abovediag[0] 0 ..... 00116 * belowdiag[0] diag[1] abovediag[1] ..... 00117 * 0 belowdiag[1] diag[2] 00118 * 0 0 belowdiag[2] ..... 00119 00120 */ 00121 template<class vec_t, class vec2_t> 00122 int solve_tridiag_nonsym(const vec_t &diag, const vec2_t &abovediag, 00123 const vec2_t &belowdiag, const vec_t &rhs, 00124 vec_t &x, size_t N) { 00125 int status; 00126 double *alpha = (double *) malloc (N * sizeof (double)); 00127 double *z = (double *) malloc (N * sizeof (double)); 00128 00129 if (alpha == 0 || z == 0) { 00130 status = GSL_ENOMEM; 00131 } else { 00132 size_t i, j; 00133 00134 /* Bidiagonalization (eliminating belowdiag) 00135 & rhs update 00136 diag' = alpha 00137 rhs' = z 00138 */ 00139 alpha[0] = O2SCL_IX(diag,0); 00140 z[0] = O2SCL_IX(rhs,0); 00141 00142 for (i = 1; i < N; i++) { 00143 const double t = O2SCL_IX(belowdiag,i-1)/alpha[i-1]; 00144 alpha[i] = O2SCL_IX(diag,i) - t*O2SCL_IX(abovediag,i-1); 00145 z[i] = O2SCL_IX(rhs,i) - t*z[i-1]; 00146 /* FIXME!!! */ 00147 if (alpha[i] == 0) { 00148 status = GSL_EZERODIV; 00149 goto solve_tridiag_nonsym_END; 00150 } 00151 } 00152 00153 /* backsubstitution */ 00154 O2SCL_IX(x,N-1) = z[N - 1]/alpha[N - 1]; 00155 if (N >= 2) { 00156 for (i = N - 2, j = 0; j <= N - 2; j++, i--) { 00157 O2SCL_IX(x,i) = (z[i] - O2SCL_IX(abovediag,i) * 00158 O2SCL_IX(x,i+1))/alpha[i]; 00159 } 00160 } 00161 00162 status = GSL_SUCCESS; 00163 } 00164 00165 solve_tridiag_nonsym_END: 00166 00167 if (z != 0) 00168 free (z); 00169 if (alpha != 0) 00170 free (alpha); 00171 00172 return status; 00173 } 00174 00175 /** 00176 \brief Solve a symmetric cyclic tridiagonal linear system 00177 00178 for description of method see [Engeln-Mullges + Uhlig, p. 96] 00179 * 00180 * diag[0] offdiag[0] 0 ..... offdiag[N-1] 00181 * offdiag[0] diag[1] offdiag[1] ..... 00182 * 0 offdiag[1] diag[2] 00183 * 0 0 offdiag[2] ..... 00184 * ... ... 00185 * offdiag[N-1] ... 00186 00187 */ 00188 template<class vec_t> 00189 int solve_cyc_tridiag_sym(const vec_t &diag, const vec_t &offdiag, 00190 const vec_t &b, vec_t &x, size_t N) { 00191 int status; 00192 double * delta = (double *) malloc (N * sizeof (double)); 00193 double * gamma = (double *) malloc (N * sizeof (double)); 00194 double * alpha = (double *) malloc (N * sizeof (double)); 00195 double * c = (double *) malloc (N * sizeof (double)); 00196 double * z = (double *) malloc (N * sizeof (double)); 00197 00198 if (delta == 0 || gamma == 0 || alpha == 0 || c == 0 || z == 0) { 00199 status = GSL_ENOMEM; 00200 } else { 00201 size_t i, j; 00202 double sum = 0.0; 00203 00204 /* factor */ 00205 00206 if (N == 1) { 00207 x[0] = b[0] / O2SCL_IX(diag,0); 00208 return GSL_SUCCESS; 00209 } 00210 00211 alpha[0] = O2SCL_IX(diag,0); 00212 gamma[0] = O2SCL_IX(offdiag,0) / alpha[0]; 00213 delta[0] = O2SCL_IX(offdiag,N-1) / alpha[0]; 00214 00215 for (i = 1; i < N - 2; i++) { 00216 alpha[i] = O2SCL_IX(diag,i) - O2SCL_IX(offdiag,i-1) * gamma[i - 1]; 00217 gamma[i] = O2SCL_IX(offdiag,i) / alpha[i]; 00218 delta[i] = -delta[i - 1] * O2SCL_IX(offdiag,i-1) / alpha[i]; 00219 } 00220 00221 for (i = 0; i < N - 2; i++) { 00222 sum += alpha[i] * delta[i] * delta[i]; 00223 } 00224 00225 alpha[N - 2] = diag[ (N - 2)] - O2SCL_IX(offdiag,N-3) * gamma[N - 3]; 00226 00227 gamma[N - 2] = (offdiag[(N - 2)] - offdiag[(N - 3)] * 00228 delta[N - 3]) / alpha[N - 2]; 00229 00230 alpha[N - 1] = diag[(N - 1)] - sum - alpha[(N - 2)] * 00231 gamma[N - 2] * gamma[N - 2]; 00232 00233 /* update */ 00234 z[0] = b[0]; 00235 for (i = 1; i < N - 1; i++) { 00236 z[i] = O2SCL_IX(b,i) - z[i - 1] * gamma[i - 1]; 00237 } 00238 sum = 0.0; 00239 for (i = 0; i < N - 2; i++) { 00240 sum += delta[i] * z[i]; 00241 } 00242 z[N - 1] = b[(N - 1)] - sum - gamma[N - 2] * z[N - 2]; 00243 for (i = 0; i < N; i++) { 00244 c[i] = z[i] / alpha[i]; 00245 } 00246 00247 /* backsubstitution */ 00248 O2SCL_IX(x,N-1) = c[N - 1]; 00249 x[(N - 2)] = c[N - 2] - gamma[N - 2] * O2SCL_IX(x,N-1); 00250 if (N >= 3) { 00251 for (i = N - 3, j = 0; j <= N - 3; j++, i--) { 00252 O2SCL_IX(x,i) = c[i] - gamma[i] * x[(i + 1)] - 00253 delta[i] * O2SCL_IX(x,N-1); 00254 } 00255 } 00256 00257 status = GSL_SUCCESS; 00258 } 00259 00260 if (z != 0) 00261 free (z); 00262 if (c != 0) 00263 free (c); 00264 if (alpha != 0) 00265 free (alpha); 00266 if (gamma != 0) 00267 free (gamma); 00268 if (delta != 0) 00269 free (delta); 00270 00271 return status; 00272 } 00273 00274 /** 00275 \brief Solve an asymmetric cyclic tridiagonal linear system 00276 00277 solve following system w/o the corner elements and then use 00278 Sherman-Morrison formula to compensate for them 00279 * 00280 * diag[0] abovediag[0] 0 ..... belowdiag[N-1] 00281 * belowdiag[0] diag[1] abovediag[1] ..... 00282 * 0 belowdiag[1] diag[2] 00283 * 0 0 belowdiag[2] ..... 00284 * ... ... 00285 * abovediag[N-1] ... 00286 00287 */ 00288 template<class vec_t> 00289 int solve_cyc_tridiag_nonsym(const vec_t &diag, const vec_t &abovediag, 00290 const vec_t &belowdiag, const vec_t &rhs, 00291 vec_t &x, size_t N) { 00292 int status; 00293 double *alpha = (double *) malloc (N * sizeof (double)); 00294 double *zb = (double *) malloc (N * sizeof (double)); 00295 double *zu = (double *) malloc (N * sizeof (double)); 00296 double *w = (double *) malloc (N * sizeof (double)); 00297 double beta; 00298 00299 if (alpha == 0 || zb == 0 || zu == 0 || w == 0) { 00300 status = GSL_ENOMEM; 00301 } else { 00302 /* Bidiagonalization (eliminating belowdiag) 00303 & rhs update 00304 diag' = alpha 00305 rhs' = zb 00306 rhs' for Aq=u is zu 00307 */ 00308 zb[0] = O2SCL_IX(rhs,0); 00309 if (O2SCL_IX(diag,0) != 0) { 00310 beta = -O2SCL_IX(diag,0); 00311 } else { 00312 beta = 1; 00313 } 00314 const double q = 1 - O2SCL_IX(abovediag,0)*O2SCL_IX(belowdiag,0)/ 00315 (O2SCL_IX(diag,0)*diag[1]); 00316 if (fabs(q/beta) > 0.5 && fabs(q/beta) < 2) { 00317 beta *= (fabs(q/beta) < 1) ? 0.5 : 2; 00318 } 00319 zu[0] = beta; 00320 alpha[0] = O2SCL_IX(diag,0) - beta; 00321 00322 { 00323 size_t i; 00324 for (i = 1; i+1 < N; i++) { 00325 const double t = O2SCL_IX(belowdiag,i-1)/alpha[i-1]; 00326 alpha[i] = O2SCL_IX(diag,i) - t*O2SCL_IX(abovediag,i-1); 00327 zb[i] = O2SCL_IX(rhs,i) - t*zb[i-1]; 00328 zu[i] = -t*zu[i-1]; 00329 /* FIXME!!! */ 00330 if (alpha[i] == 0) { 00331 status = GSL_EZERODIV; 00332 goto solve_cyc_tridiag_nonsym_END; 00333 } 00334 } 00335 } 00336 00337 { 00338 const size_t i = N-1; 00339 const double t = O2SCL_IX(belowdiag,i-1)/alpha[i-1]; 00340 alpha[i]=O2SCL_IX(diag,i)-O2SCL_IX(abovediag,i)* 00341 O2SCL_IX(belowdiag,i)/beta- 00342 t*O2SCL_IX(abovediag,i-1); 00343 zb[i] = O2SCL_IX(rhs,i) - t*zb[i-1]; 00344 zu[i] = O2SCL_IX(abovediag,i) - t*zu[i-1]; 00345 00346 /* FIXME!!! */ 00347 if (alpha[i] == 0) { 00348 status = GSL_EZERODIV; 00349 goto solve_cyc_tridiag_nonsym_END; 00350 } 00351 } 00352 00353 { 00354 /* backsubstitution */ 00355 size_t i, j; 00356 w[N-1] = zu[N-1]/alpha[N-1]; 00357 O2SCL_IX(x,N-1) = zb[N-1]/alpha[N-1]; 00358 for (i = N - 2, j = 0; j <= N - 2; j++, i--) { 00359 w[i] = (zu[i] - O2SCL_IX(abovediag,i) * w[i+1])/alpha[i]; 00360 O2SCL_IX(x,i) = (zb[i] - O2SCL_IX(abovediag,i) * 00361 O2SCL_IX(x,i + 1))/alpha[i]; 00362 } 00363 } 00364 00365 /* Sherman-Morrison */ 00366 const double vw = w[0] + O2SCL_IX(belowdiag,N-1)/beta * w[N-1]; 00367 const double vx = O2SCL_IX(x,0) + 00368 O2SCL_IX(belowdiag,N-1)/beta * O2SCL_IX(x,N-1); 00369 00370 /* FIXME!!! */ 00371 if (vw + 1 == 0) { 00372 status = GSL_EZERODIV; 00373 goto solve_cyc_tridiag_nonsym_END; 00374 } 00375 00376 { 00377 size_t i; 00378 for (i = 0; i < N; i++) 00379 O2SCL_IX(x,i) -= vx/(1 + vw)*w[i]; 00380 } 00381 00382 status = GSL_SUCCESS; 00383 } 00384 00385 solve_cyc_tridiag_nonsym_END: 00386 00387 if (zb != 0) 00388 free (zb); 00389 if (zu != 0) 00390 free (zu); 00391 if (w != 0) 00392 free (w); 00393 if (alpha != 0) 00394 free (alpha); 00395 00396 return status; 00397 } 00398 00399 #ifdef DOXYGENP 00400 } 00401 #endif
Documentation generated with Doxygen and provided under the GNU Free Documentation License. See License Information for details.
Project hosting provided by
,
O2scl Sourceforge Project Page