tridiag_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /* linalg/tridiag.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2004,
26  * 2007 Gerard Jungman, Brian Gough, David Necas
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 3 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41  * 02110-1301, USA.
42  */
43 /** \file tridiag_base.h
44  \brief File for solving tridiagonal systems
45 */
46 
47 #ifdef DOXYGEN
48 namespace o2scl_linalg {
49 #endif
50 
51  /** \brief Allocation object for 2 arrays of equal size
52 
53  This class is used in solve_tridiag_nonsym().
54  */
56  public:
58  ubvector v1, v2;
59  void resize(size_t n) {
60  v1.resize(n);
61  v2.resize(n);
62  }
63  };
64 
65  /** \brief Allocation object for 4 arrays of equal size
66 
67  This class is used in \ref solve_tridiag_sym() and
68  \ref solve_cyc_tridiag_nonsym().
69  */
71  public:
73  ubvector v1, v2, v3, v4;
74  void resize(size_t n) {
75  v1.resize(n);
76  v2.resize(n);
77  v3.resize(n);
78  v4.resize(n);
79  }
80  };
81 
82  /** \brief Allocation object for 5 arrays of equal size
83 
84  This class is used in solve_cyc_tridiag_sym().
85  */
87  public:
89  ubvector v1, v2, v3, v4, v5;
90  void resize(size_t n) {
91  v1.resize(n);
92  v2.resize(n);
93  v3.resize(n);
94  v4.resize(n);
95  v5.resize(n);
96  }
97  };
98 
99  /** \brief Solve a symmetric tridiagonal linear system with
100  user-specified memory
101 
102  This function solves the system \f$ A x = b \f$ where
103  \f$ A \f$ is a matrix of the form
104  \verbatim
105  *
106  * diag[0] offdiag[0] 0 .....
107  * offdiag[0] diag[1] offdiag[1] .....
108  * 0 offdiag[1] diag[2]
109  * 0 0 offdiag[2] .....
110  \endverbatim
111  given the \c N diagonal elements in \c diag, \c N-1 diagonal
112  elements in \c offdiag, and the \c N elements \c b from the RHS.
113 
114  See \ref EngelnMullges96 .
115  */
116  template<class vec_t, class vec2_t, class vec3_t,
117  class vec4_t, class mem_t, class mem_vec_t>
118  void solve_tridiag_sym(const vec_t &diag, const vec2_t &offdiag,
119  const vec3_t &b, vec4_t &x, size_t N, mem_t &m) {
120 
121  mem_vec_t &gamma=m.v1;
122  mem_vec_t &alpha=m.v2;
123  mem_vec_t &c=m.v3;
124  mem_vec_t &z=m.v4;
125 
126  size_t i, j;
127 
128  /* [GSL] Cholesky decomposition
129  A = L.D.L^t
130  lower_diag(L) = gamma
131  diag(D) = alpha
132  */
133  alpha[0] = O2SCL_IX(diag,0);
134  gamma[0] = O2SCL_IX(offdiag,0) / alpha[0];
135 
136  for (i = 1; i < N - 1; i++) {
137  alpha[i] = O2SCL_IX(diag,i) - O2SCL_IX(offdiag,i-1) * gamma[i - 1];
138  gamma[i] = O2SCL_IX(offdiag,i) / alpha[i];
139  }
140 
141  if (N > 1) {
142  alpha[N-1]=O2SCL_IX(diag,N-1)-O2SCL_IX(offdiag,N-2)*gamma[N-2];
143  }
144 
145  // [GSL] update RHS
146  z[0] = b[0];
147  for (i = 1; i < N; i++) {
148  z[i] = O2SCL_IX(b,i) - gamma[i - 1] * z[i - 1];
149  } for (i = 0; i < N; i++) {
150  c[i] = z[i] / alpha[i];
151  }
152 
153  // [GSL] backsubstitution
154  O2SCL_IX(x,N-1)=c[N-1];
155  if (N >= 2) {
156  for (i = N - 2, j = 0; j <= N - 2; j++, i--) {
157  O2SCL_IX(x,i) = c[i] - gamma[i] * O2SCL_IX(x,i+1);
158  }
159  }
160 
161  return;
162  }
163 
164  /** \brief Solve an asymmetric tridiagonal linear system with
165  user-specified memory
166 
167  This function solves the system \f$ A x = b \f$ where
168  \f$ A \f$ is a matrix of the form
169  \verbatim
170  *
171  * diag[0] abovediag[0] 0 .....
172  * belowdiag[0] diag[1] abovediag[1] .....
173  * 0 belowdiag[1] diag[2]
174  * 0 0 belowdiag[2] .....
175  \endverbatim
176  This function uses plain Gauss elimination, not bothering
177  with the zeroes.
178 
179  \future Offer an option to avoid throwing on divide by zero?
180  */
181  template<class vec_t, class vec2_t, class vec3_t, class vec4_t,
182  class vec5_t, class mem_t, class mem_vec_t>
183  void solve_tridiag_nonsym(const vec_t &diag, const vec2_t &abovediag,
184  const vec3_t &belowdiag, const vec4_t &rhs,
185  vec5_t &x, size_t N, mem_t &m) {
186  mem_vec_t &alpha=m.v1;
187  mem_vec_t &z=m.v2;
188 
189  size_t i, j;
190 
191  /* [GSL] Bidiagonalization (eliminating belowdiag)
192  & rhs update
193  diag' = alpha
194  rhs' = z
195  */
196  alpha[0] = O2SCL_IX(diag,0);
197  z[0] = O2SCL_IX(rhs,0);
198 
199  for (i = 1; i < N; i++) {
200  const double t = O2SCL_IX(belowdiag,i-1)/alpha[i-1];
201  alpha[i] = O2SCL_IX(diag,i) - t*O2SCL_IX(abovediag,i-1);
202  z[i] = O2SCL_IX(rhs,i) - t*z[i-1];
203  if (alpha[i] == 0) {
204  O2SCL_ERR2("Divide by zero in ",
205  "solve_tridiag_nonsym().",o2scl::exc_ezerodiv);
206  return;
207  }
208  }
209 
210  // [GSL] backsubstitution
211  O2SCL_IX(x,N-1) = z[N - 1]/alpha[N - 1];
212  if (N >= 2) {
213  for (i = N - 2, j = 0; j <= N - 2; j++, i--) {
214  O2SCL_IX(x,i) = (z[i] - O2SCL_IX(abovediag,i) *
215  O2SCL_IX(x,i+1))/alpha[i];
216  }
217  }
218 
219  return;
220  }
221 
222  /** \brief Solve a symmetric cyclic tridiagonal linear system with
223  user specified memory
224 
225  This function solves the system \f$ A x = b \f$ where
226  \f$ A \f$ is a matrix of the form
227  \verbatim
228  *
229  * diag[0] offdiag[0] 0 ..... offdiag[N-1]
230  * offdiag[0] diag[1] offdiag[1] .....
231  * 0 offdiag[1] diag[2]
232  * 0 0 offdiag[2] .....
233  * ... ...
234  * offdiag[N-1] ...
235  \endverbatim
236 
237  See \ref EngelnMullges96 .
238  */
239  template<class vec_t, class vec2_t, class vec3_t, class vec4_t,
240  class mem_t, class mem_vec_t>
241  void solve_cyc_tridiag_sym(const vec_t &diag, const vec2_t &offdiag,
242  const vec3_t &b, vec4_t &x, size_t N,
243  mem_t &m) {
244 
245  mem_vec_t &delta=m.v1;
246  mem_vec_t &gamma=m.v2;
247  mem_vec_t &alpha=m.v3;
248  mem_vec_t &c=m.v4;
249  mem_vec_t &z=m.v5;
250 
251  size_t i, j;
252  double sum = 0.0;
253 
254  // [GSL] factor
255 
256  if (N == 1) {
257  x[0] = b[0] / O2SCL_IX(diag,0);
258  return;
259  }
260 
261  alpha[0] = O2SCL_IX(diag,0);
262  gamma[0] = O2SCL_IX(offdiag,0) / alpha[0];
263  delta[0] = O2SCL_IX(offdiag,N-1) / alpha[0];
264 
265  for (i = 1; i < N - 2; i++) {
266  alpha[i] = O2SCL_IX(diag,i) - O2SCL_IX(offdiag,i-1) * gamma[i - 1];
267  gamma[i] = O2SCL_IX(offdiag,i) / alpha[i];
268  delta[i] = -delta[i - 1] * O2SCL_IX(offdiag,i-1) / alpha[i];
269  }
270 
271  for (i = 0; i < N - 2; i++) {
272  sum += alpha[i] * delta[i] * delta[i];
273  }
274 
275  alpha[N - 2] = diag[ (N - 2)] - O2SCL_IX(offdiag,N-3) * gamma[N - 3];
276 
277  gamma[N - 2] = (offdiag[(N - 2)] - offdiag[(N - 3)] *
278  delta[N - 3]) / alpha[N - 2];
279 
280  alpha[N - 1] = diag[(N - 1)] - sum - alpha[(N - 2)] *
281  gamma[N - 2] * gamma[N - 2];
282 
283  // [GSL] update
284 
285  z[0] = b[0];
286  for (i = 1; i < N - 1; i++) {
287  z[i] = O2SCL_IX(b,i) - z[i - 1] * gamma[i - 1];
288  }
289  sum = 0.0;
290  for (i = 0; i < N - 2; i++) {
291  sum += delta[i] * z[i];
292  }
293  z[N - 1] = b[(N - 1)] - sum - gamma[N - 2] * z[N - 2];
294  for (i = 0; i < N; i++) {
295  c[i] = z[i] / alpha[i];
296  }
297 
298  // [GSL] backsubstitution
299  O2SCL_IX(x,N-1) = c[N - 1];
300  x[(N - 2)] = c[N - 2] - gamma[N - 2] * O2SCL_IX(x,N-1);
301  if (N >= 3) {
302  for (i = N - 3, j = 0; j <= N - 3; j++, i--) {
303  O2SCL_IX(x,i) = c[i] - gamma[i] * x[(i + 1)] -
304  delta[i] * O2SCL_IX(x,N-1);
305  }
306  }
307 
308  return;
309  }
310 
311  /** \brief Solve an asymmetric cyclic tridiagonal linear system
312  with user-specified memory
313 
314  This function solves the system \f$ A x = b \f$ where
315  \f$ A \f$ is a matrix of the form
316  \verbatim
317  *
318  * diag[0] abovediag[0] 0 ..... belowdiag[N-1]
319  * belowdiag[0] diag[1] abovediag[1] .....
320  * 0 belowdiag[1] diag[2]
321  * 0 0 belowdiag[2] .....
322  * ... ...
323  * abovediag[N-1] ...
324  \endverbatim
325 
326  This function solves the following system without the corner
327  elements and then use Sherman-Morrison formula to compensate for
328  them.
329 
330  \comment
331  Note that the three FIXME!!! entries are from the GSL original.
332  \endcomment
333 
334  \future Offer an option to avoid throwing on divide by zero?
335  */
336  template<class vec_t, class vec2_t, class vec3_t, class vec4_t,
337  class vec5_t, class mem_t, class mem_vec_t>
338  void solve_cyc_tridiag_nonsym(const vec_t &diag, const vec2_t &abovediag,
339  const vec3_t &belowdiag, const vec4_t &rhs,
340  vec5_t &x, size_t N, mem_t &m) {
341 
342  double beta;
343  mem_vec_t &alpha=m.v1;
344  mem_vec_t &zb=m.v2;
345  mem_vec_t &zu=m.v3;
346  mem_vec_t &w=m.v4;
347 
348  /* [GSL] Bidiagonalization (eliminating belowdiag)
349  & rhs update
350  diag' = alpha
351  rhs' = zb
352  rhs' for Aq=u is zu
353  */
354  zb[0] = O2SCL_IX(rhs,0);
355  if (O2SCL_IX(diag,0) != 0) {
356  beta = -O2SCL_IX(diag,0);
357  } else {
358  beta = 1;
359  }
360  const double q = 1 - O2SCL_IX(abovediag,0)*O2SCL_IX(belowdiag,0)/
361  (O2SCL_IX(diag,0)*diag[1]);
362  if (fabs(q/beta) > 0.5 && fabs(q/beta) < 2) {
363  beta *= (fabs(q/beta) < 1) ? 0.5 : 2;
364  }
365  zu[0] = beta;
366  alpha[0] = O2SCL_IX(diag,0) - beta;
367 
368  {
369  size_t i;
370  for (i = 1; i+1 < N; i++) {
371  const double t = O2SCL_IX(belowdiag,i-1)/alpha[i-1];
372  alpha[i] = O2SCL_IX(diag,i) - t*O2SCL_IX(abovediag,i-1);
373  zb[i] = O2SCL_IX(rhs,i) - t*zb[i-1];
374  zu[i] = -t*zu[i-1];
375  // [GSL] FIXME!!!
376  if (alpha[i] == 0) {
377  O2SCL_ERR2("Divide by zero (1) in ",
378  "solve_cyc_tridiag_nonsym().",o2scl::exc_ezerodiv);
379  }
380  }
381  }
382 
383  {
384  const size_t i = N-1;
385  const double t = O2SCL_IX(belowdiag,i-1)/alpha[i-1];
386  alpha[i]=O2SCL_IX(diag,i)-O2SCL_IX(abovediag,i)*
387  O2SCL_IX(belowdiag,i)/beta-
388  t*O2SCL_IX(abovediag,i-1);
389  zb[i] = O2SCL_IX(rhs,i) - t*zb[i-1];
390  zu[i] = O2SCL_IX(abovediag,i) - t*zu[i-1];
391 
392  // [GSL] FIXME!!!
393  if (alpha[i] == 0) {
394  O2SCL_ERR2("Divide by zero (2) in ",
395  "solve_cyc_tridiag_nonsym().",o2scl::exc_ezerodiv);
396  }
397  }
398 
399  {
400  // [GSL] backsubstitution
401  size_t i, j;
402  w[N-1] = zu[N-1]/alpha[N-1];
403  O2SCL_IX(x,N-1) = zb[N-1]/alpha[N-1];
404  for (i = N - 2, j = 0; j <= N - 2; j++, i--) {
405  w[i] = (zu[i] - O2SCL_IX(abovediag,i) * w[i+1])/alpha[i];
406  O2SCL_IX(x,i) = (zb[i] - O2SCL_IX(abovediag,i) *
407  O2SCL_IX(x,i + 1))/alpha[i];
408  }
409  }
410 
411  // [GSL] Sherman-Morrison
412  const double vw = w[0] + O2SCL_IX(belowdiag,N-1)/beta * w[N-1];
413  const double vx = O2SCL_IX(x,0) +
414  O2SCL_IX(belowdiag,N-1)/beta * O2SCL_IX(x,N-1);
415 
416  // [GSL] FIXME!!!
417  if (vw + 1 == 0) {
418  O2SCL_ERR2("Divide by zero (3) in ",
419  "solve_cyc_tridiag_nonsym().",o2scl::exc_ezerodiv);
420  }
421 
422  {
423  size_t i;
424  for (i = 0; i < N; i++)
425  O2SCL_IX(x,i) -= vx/(1 + vw)*w[i];
426  }
427 
428  return;
429  }
430 
431  /** \brief Solve a symmetric tridiagonal linear system
432  */
433  template<class vec_t, class vec2_t, class vec3_t, class vec4_t>
434  void solve_tridiag_sym(const vec_t &diag, const vec2_t &offdiag,
435  const vec3_t &b, vec4_t &x, size_t N) {
436  typedef boost::numeric::ublas::vector<double> ubvector;
437  ubvector_4_mem v4m;
438  v4m.resize(N);
439  solve_tridiag_sym<vec_t,vec2_t,vec3_t,vec4_t,ubvector_4_mem,
440  ubvector>(diag,offdiag,b,x,N,v4m);
441  return;
442  }
443 
444  /** \brief Solve an asymmetric tridiagonal linear system
445  */
446  template<class vec_t, class vec2_t, class vec3_t, class vec4_t,
447  class vec5_t>
448  void solve_tridiag_nonsym(const vec_t &diag, const vec2_t &abovediag,
449  const vec3_t &belowdiag, const vec4_t &rhs,
450  vec5_t &x, size_t N) {
451  typedef boost::numeric::ublas::vector<double> ubvector;
452  ubvector_2_mem v2m;
453  v2m.resize(N);
454  solve_tridiag_nonsym<vec_t,vec2_t,vec3_t,vec4_t,vec5_t,ubvector_2_mem,
455  ubvector>(diag,abovediag,belowdiag,rhs,x,N,v2m);
456  return;
457  }
458 
459  /** \brief Solve a symmetric cyclic tridiagonal linear system
460  */
461  template<class vec_t, class vec2_t, class vec3_t, class vec4_t>
462  void solve_cyc_tridiag_sym(const vec_t &diag, const vec2_t &offdiag,
463  const vec3_t &b, vec4_t &x, size_t N) {
464  typedef boost::numeric::ublas::vector<double> ubvector;
465  ubvector_5_mem v5m;
466  v5m.resize(N);
467  solve_cyc_tridiag_sym<vec_t,vec2_t,vec3_t,vec4_t,ubvector_5_mem,
468  ubvector>(diag,offdiag,b,x,N,v5m);
469  return;
470  }
471 
472  /** \brief Solve an asymmetric cyclic tridiagonal linear system
473  */
474  template<class vec_t, class vec2_t, class vec3_t, class vec4_t,
475  class vec5_t>
476  void solve_cyc_tridiag_nonsym(const vec_t &diag, const vec2_t &abovediag,
477  const vec3_t &belowdiag, const vec4_t &rhs,
478  vec5_t &x, size_t N) {
479  typedef boost::numeric::ublas::vector<double> ubvector;
480  ubvector_4_mem v4m;
481  v4m.resize(N);
482  solve_cyc_tridiag_nonsym<vec_t,vec2_t,vec3_t,vec4_t,vec5_t,
483  ubvector_4_mem,ubvector>(diag,abovediag,belowdiag,rhs,x,N,v4m);
484  return;
485  }
486 
487 #ifdef DOXYGEN
488 }
489 #endif
tried to divide by zero
Definition: err_hnd.h:75
void solve_tridiag_sym(const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
Solve a symmetric tridiagonal linear system with user-specified memory.
Definition: tridiag_base.h:118
Allocation object for 5 arrays of equal size.
Definition: tridiag_base.h:86
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
Allocation object for 4 arrays of equal size.
Definition: tridiag_base.h:70
void solve_cyc_tridiag_nonsym(const vec_t &diag, const vec2_t &abovediag, const vec3_t &belowdiag, const vec4_t &rhs, vec5_t &x, size_t N, mem_t &m)
Solve an asymmetric cyclic tridiagonal linear system with user-specified memory.
Definition: tridiag_base.h:338
void solve_tridiag_nonsym(const vec_t &diag, const vec2_t &abovediag, const vec3_t &belowdiag, const vec4_t &rhs, vec5_t &x, size_t N, mem_t &m)
Solve an asymmetric tridiagonal linear system with user-specified memory.
Definition: tridiag_base.h:183
Allocation object for 2 arrays of equal size.
Definition: tridiag_base.h:55
void solve_cyc_tridiag_sym(const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
Solve a symmetric cyclic tridiagonal linear system with user specified memory.
Definition: tridiag_base.h:241

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).