A..0.0.1 chapter1/p13.c

   1 /* 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991)
   2    P. 13 のプログラムをC言語に書き直した。
   3   ** GAUSS ELIMINATION METHOD **
   4   n : NUMBER OF UNKNOWS
   5   a : COEFFICIENT MATRIX
   6   f : KNOWN VECTOR & SOLUTION VECTORS
   7 */
   8 
   9 #include <stdio.h>
  10 #include <stdlib.h> // exit()
  11 #define NDIM (100)
  12 
  13 void gauss(int ndim, double a[ndim][ndim], double f[], int N);
  14 
  15 int main(void)
  16 {
  17   int n, i, j;
  18   double a[NDIM][NDIM], f[NDIM];
  19   printf("INPUT : N=\n"); // 最後の \n は不要かも
  20   scanf("%d", &n);
  21   if (n<=1)
  22     exit(0);
  23   for (i = 0; i < n; i++)
  24     for (j = 0; j < n; j++) {
  25       printf("INPUT : A(%3d,%3d)=\n", i+1, j+1); // 最後の \n は不要かも
  26       scanf("%lf", &a[i][j]);
  27     }
  28   for (i = 0; i < n; i++) {
  29     printf(" INPUT : F(%3d)=\n", i+1); // 最後の \n は不要かも
  30     scanf("%lf", &f[i]);
  31   }
  32   gauss(NDIM,a,f,n);
  33   for (i = 0; i < 5; i++)
  34     printf("  I     X(I)    ");
  35   printf("\n");
  36   for (i = 0; i < n; i++) {
  37     printf("%3d%12.4e ",i+1, f[i]); //I3 は %3d, PE12.4 は %12.4e (1Pは無視)
  38     if ((i+1) % 5 == 0 || i+1 == n)
  39       printf("\n");
  40   }
  41 }
  42 
  43 void gauss(int ndim, double a[ndim][ndim], double f[], int n)
  44 {
  45   int i, j, k;
  46   double aa;
  47   // FORWARD ELIMITNATION
  48   for (i=0; i<n-1; i++)
  49     for (j=i+1; j<n; j++) {
  50       aa = a[j][i] / a[i][i];
  51       for (k=i+1; k<n; k++)
  52     a[j][k] -= aa * a[i][k];
  53       f[j] -= aa * f[i];
  54     }
  55   // BACKWARD SUBSTITUTION
  56   f[n-1] /= a[n-1][n-1];
  57   for (i=n-2;i>=0;i--) {
  58     for (j=i+1; j<n; j++)
  59       f[i] -= a[i][j] * f[j];
  60     f[i] /= a[i][i];
  61   }
  62 }