A..0.0.10 chapter5/p104.c

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