ベキ乗法を用いて行列の絶対値最大の固有値を求める。

以下は全２分木の隣接行列をベキ乗法で解くためのプログラム。Mは行列のサイズにより適宜書き直す必要あり。

全2分木ではｎを3,7,15,31,63,127,255,511,1023,2047,4095,8191,16383とした時のデータを取る。

---

#include <stdio.h>
#include <math.h>
#define M 383                    //入力するnよりも大きくしておくこと
double a[M][M],x[M]={0},y[M];
int n;

void Beki(void);
void Throw(void);

int main(){

  int i,j;
  printf("行列の大きさ：");    scanf("%d", &n);  //解きたい行列のサイズ

  Throw();           //行列の要素を代入。
  Beki();             //ベキ乗法のcode
  return 0;
}

void Beki(void){
   double s,xmax,t,b;
   int i,j,imax,l,m,triger,number=1000;
   x[0] = 1.0;

   imax = 50000000;   l = 0;   xmax = 0.0;   s = 0.0;

   for(i=0;i<n;i++)  s += x[i]*x[i];
   for(i=0;i<n;i++)  x[i] /= sqrt(s);

   for(m=1;m<=imax;m++){
      for(i=0;i<n;i++){
         y[i] = 0.0;
         for(j=0;j<n;j++)  y[i] += a[i][j]*x[j];
      }

      s = 0.0;
      for(i=0;i<n;i++)  s += y[i]*y[i];

      t = 0.0;
      for(i=0;i<n;i++)  t += y[i]*x[i];

      b = s/t;  l = m;

      if(m >= number){
         printf("xmax = %15.12f  now = %d  \n", xmax,m);
         printf("continue ? yes--追加反復回数 no--1:");  //収束が悪ければ追加で反復
         scanf("%d", &triger);
         number += triger;
         if(triger == 1)  break;
      }

      xmax = b;
      for(i=0;i<n;i++)  x[i] = y[i]/sqrt(s);
   }
   imax = l;

   printf("最大固有値　%15.12f \t 反復回数　%d\n",xmax,imax);
}

void Throw(void){
  int i,j;
  for(i=0;i<n;i++)   a[i][i] = 1;


   for(i=0;i<=2;i++){
     a[0][i] = 1;
     a[i][0] = 1;
   }
   j = 3;
   for(i=1; i<n;i++){
      a[i][j] = 1;
      a[j][i] = 1;
      j++;
      a[i][j] = 1;
      a[j][i] = 1;
      if(j>=n)  break;
      j++;
   }
}