A..0.0.9 chapter4/p88.c

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