関数 trid() の使用例

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