A..0.0.11 chapter5/p98.c

   1 // 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 98
   2 /* ** INVERSE ITERATION METHOD **
   3    N : ORDER OF MATRIX
   4    A : COEFFICIENT MATRIX
   5    EPS : THRESHOLD VALUE FOR ERROR
   6    M : MAXIMUM NO. OF ITERATIONS
   7    EVAL : EIGENVALUE
   8    X : EIGENVECTOR
   9    NITER : NO. OF ITERATIONS
  10 */
  11 
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include <math.h>
  15 
  16 void gauss(int ndim, double a[][ndim], double f[], int N, int ID);
  17 void inpr(double x[], double y[], int n, double *s);
  18 
  19 #define NDIM (100)
  20 
  21 int main(void)
  22 {
  23   int n, m, i, j, niter;
  24   double eps, eval, s, t, enew, error;
  25   double a[NDIM][NDIM], x[NDIM], y[NDIM];
  26   // INPUT
  27   printf("INPUT : N= \n");
  28   scanf("%d", &n);
  29   if (n <= 1) exit(1);
  30   printf("INPUT: EPS=\n");
  31   scanf("%lf", &eps);
  32   printf("INPUT : M= \n");
  33   scanf("%d", &m);
  34   for (i = 0; i < n; i++) {
  35     for (j=0; j < n; j++) {
  36       printf("INPUT : A(%3d,%3d)=\n", i+1, j+1);
  37       scanf("%lf", &a[i][j]);
  38     }
  39   }
  40   for (i=0; i<n; i++) {
  41     printf("INPUT:X(%3d)=\n",i+1);
  42     scanf("%lf", &x[i]);
  43   }
  44   // ITERATION
  45   niter = 0;
  46   eval = 0.0;
  47   do {
  48     for (i=0; i<n; i++)
  49       y[i] = x[i];
  50     gauss(NDIM,a,y,n,niter);
  51     inpr(x,y,n,&s);
  52     inpr(y,y,n,&t);
  53     enew = s / t;
  54     t = sqrt(t);
  55     for (i=0; i<n; i++)
  56       x[i] = y[i] / t;
  57     error = fabs((enew-eval)/enew);
  58     eval = enew;
  59     niter++;
  60     printf("(NITER,EVAL,ERROR)=%d %g %g\n", niter, eval, error);
  61     if (niter == m) break;
  62   }
  63   while (error > eps);
  64   // OUTPUT
  65   printf("NO. OF ITERATIONS = %3d   EPS=%12.4e\n", niter, eps);
  66   printf("EIGENVALUE=%12.4e\n", eval);
  67   printf("**EIGENVECTOR**\n");
  68   for (i = 0; i < 5; i++)
  69     printf("  I     X(I)    ");
  70   printf("\n");
  71   for (i=0; i<n; i++) {
  72     printf("%3d%12.4e ", i+1, x[i]);
  73     if ((i+1) % 5 == 0 || (i+1) == n)
  74       printf("\n");
  75   }
  76 }
  77 
  78 void gauss(int ndim, double a[][ndim], double f[], int n, int id)
  79 {
  80   int i, j, k;
  81   double aa;
  82   // FORWARD ELIMINATION
  83   for (i=0; i<n-1; i++) {
  84     for (j=i+1; j<n; j++) {
  85       aa = a[j][i] / a[i][i];
  86       if (id == 0) {
  87     for (k=i+1; k<n; k++)
  88       a[j][k] -= aa * a[i][k];
  89       }
  90       f[j] -= aa * f[i];
  91     }
  92   }
  93   // BACKWARD SUBSTITUTION
  94   f[n-1] /= a[n-1][n-1];
  95   for (i=n-2; i>= 0; i--) {
  96     for (j=i+1; j<n; j++)
  97       f[i] -= a[i][j] * f[j];
  98     f[i] /= a[i][i];
  99   }
 100 }
 101 
 102 void inpr(double x[], double y[], int n, double *s)
 103 {
 104   int i;
 105   // INNER PRODUCT
 106   *s = 0.0;
 107   for (i=0; i<n; i++)
 108     *s += x[i] * y[i];
 109 }