関数 trid() の使用例

   1 /*
   2  * test-lu.c --- 三重対角行列係数の連立1次方程式を解く trilu(), trisol()
   3  *               のテスト
   4  *   コンパイル＆リンク: gcc -o test-lu test-lu.c trid-lu.c
   5  *   trid-lu.c, trid-lu.h の入手は以下の URL で
   6  *   http://www.math.meiji.ac.jp/~mk/program/linear/trid-lu.c
   7  *   http://www.math.meiji.ac.jp/~mk/program/linear/trid-lu.h
   8  */
   9 
  10 #define	NDIM 100
  11 
  12 #include <stdio.h>
  13 #include <math.h>
  14 #include "trid-lu.h"
  15 
  16 int main()
  17 {
  18   int i, n, nm1;
  19   double al[NDIM],ad[NDIM],au[NDIM],b[NDIM],x[NDIM];
  20   n = 10; nm1 = n - 1;
  21   /* 係数行列 A */
  22   ad[0] = 2.0; au[0] = -1.0;
  23   for (i = 1; i < nm1; i++) {
  24     al[i] = -1.0; ad[i] = 2.0; au[i] = - 1.0;
  25   }	
  26   al[nm1] = - 1.0; ad[nm1] = 2.0;
  27   /* 解 x */
  28   for (i = 0; i < n; i++)
  29     x[i] = i;
  30   /* 右辺 */
  31   b[0] = ad[0] * x[0] + au[0] * x[1];
  32   for (i = 1; i < nm1; i++)
  33     b[i] = al[i] * x[i-1] + ad[i] * x[i] + au[i] * x[i+1];
  34   b[nm1] = al[nm1] * x[nm1-1] + ad[nm1] * x[nm1];
  35   /* A を LU分解する */
  36   trilu(n, al, ad, au);
  37   /* 後退代入して連立1次方程式 A x=b を解く */
  38   trisol(n, al, ad, au, b);
  39   /* 結果を表示する */
  40   for (i = 0; i < n; i++)
  41     printf("%f\n", b[i]);
  42   return 0;
  43 }

oyabun% make
gcc -O -W -Wall   -c -o test-trid-lu.o test-trid-lu.c
gcc -O -W -Wall -o test-trid-lu test-trid-lu.o trid-lu.o
oyabun% ./test-trid-lu 
-0.000000
1.000000
2.000000
3.000000
4.000000
5.000000
6.000000
7.000000
8.000000
9.000000
oyabun%