関数 trid(), trilu(), trisol()

   1 /* trid-lu.c -- 3項方程式を Gauss の消去法で解く */
   2 
   3 #include "trid-lu.h"
   4 
   5 /* 3項方程式 (係数行列が三重対角である連立1次方程式のこと) Ax=b を解く
   6  *
   7  *   入力
   8  *     n: 未知数の個数
   9  *     al,ad,au: 連立1次方程式の係数行列
  10  *       (al: 対角線の下側 i.e. 下三角部分 (lower part)
  11  *        ad: 対角線       i.e. 対角部分 (diagonal part)
  12  *        au: 対角線の上側 i.e. 上三角部分 (upper part)
  13  *        つまり
  14  *
  15  *         ad[0] au[0]   0   .................. 0
  16  *         al[1] ad[1] au[1]   0  ............. 0
  17  *           0   al[2] ad[2] au[2]  0 ......... 0
  18  *                         ....................
  19  *                         al[n-2] ad[n-2] au[n-2]
  20  *                           0     al[n-1] ad[n-1]
  21 
  22  *        al[i] = A_{i,i-1}, ad[i] = A_{i,i}, au[i] = A_{i,i+1},
  23  *        al[0], au[n-1] は意味がない)
  24  *
  25  *     b: 連立1次方程式の右辺の既知ベクトル
  26  *        (添字は 0 から。i.e. b[0],b[1],...,b[n-1] にデータが入っている。)
  27  *   出力
  28  *     al,ad,au: 入力した係数行列を LU 分解したもの
  29  *     b: 連立1次方程式の解
  30  *   能書き
  31  *     一度 call すると係数行列を LU 分解したものが返されるので、
  32  *     以後は同じ係数行列に関する連立1次方程式を解くために、
  33  *     関数 trisol() が使える。
  34  *   注意
  35  *     Gauss の消去法を用いているが、ピボットの選択等はしていな
  36  *     いので、ピボットの選択をしていないので、係数行列が正定値である
  37  *     などの適切な条件がない場合は結果が保証できない。
  38  */
  39 
  40 void trid(int n, double *al, double *ad, double *au, double *b)
  41 {
  42   trilu(n,al,ad,au);
  43   trisol(n,al,ad,au,b);
  44 }
  45 
  46 /* 三重対角行列の LU 分解 (pivoting なし) */
  47 void trilu(int n, double *al, double *ad, double *au)
  48 {
  49   int i, nm1 = n - 1;
  50   /* 前進消去 (forward elimination) */
  51   for (i = 0; i < nm1; i++) {
  52     al[i + 1] /= ad[i];
  53     ad[i + 1] -= au[i] * al[i + 1];
  54   }
  55 }
  56 
  57 /* LU 分解済みの三重対角行列を係数に持つ3項方程式を解く */
  58 void trisol(int n, double *al, double *ad, double *au, double *b)
  59 {
  60   int i, nm1 = n - 1;
  61   /* 前進消去 (forward elimination) */
  62   for (i = 0; i < nm1; i++) b[i + 1] -= b[i] * al[i + 1];
  63   /* 後退代入 (backward substitution) */
  64   b[nm1] /= ad[nm1];
  65   for (i = n - 2; i >= 0; i--) b[i] = (b[i] - au[i] * b[i + 1]) / ad[i];
  66 }