Jacobi法プログラム
#include <stdio.h>
#include <math.h>
#define M 10
void PrintMatrix(double a[M][M], int n);
int main(void){
double a[M][M],v[M][M];
double eps=0.0001,div,r,t,s,c,apj,aqj,aip,aiq,vip,viq;
int i,j,n,iter,count,iterMax=100,p,q;
//------------------------------------------------------------------
// Input
//------------------------------------------------------------------
scanf("%d",&n);
for(i=1;i<=n;i++){
for(j=1;j<=n;j++){
scanf("%lf",&a[i][j]);
}
}
PrintMatrix(a,n);
//------------------------------------------------------------------
// Jacobi
//------------------------------------------------------------------
for(i=1;i<=n;i++){
for(j=1;j<=n;j++){
v[i][j]=0.;
}
v[i][i]=1.;
}
for(iter=1;iter<=iterMax;iter++){
count=0;
for(p=1;p<=n-1;p++){
for(q=p+1;q<=n;q++){
if(fabs(a[p][q])<eps) continue;
count++;
div=a[p][p]-a[q][q];
if (div != 0.0){
r=2.0*a[p][q]/div;
t=0.5*atan(r);
} else {
t=0.78539818;
}
s=sin(t);
c=cos(t);
for(j=1;j<=n;j++){
apj=a[p][j];
aqj=a[q][j];
a[p][j]=apj*c+aqj*s;
a[q][j]=-apj*s+aqj*c;
}
for(i=1;i<=n;i++){
aip=a[i][p];
aiq=a[i][q];
a[i][p]=aip*c+aiq*s;
a[i][q]=-aip*s+aiq*c;
vip=v[i][p];
viq=v[i][q];
v[i][p]=vip*c+viq*s;
v[i][q]=-vip*s+viq*c;
}
printf("p,q=%3d,%3d\n",p,q);
PrintMatrix(a,n);
}
}
if (count==0) break;
}
//------------------------------------------------------------------
//------------------------------------------------------------------
// Output
//------------------------------------------------------------------
printf("Eigen values:\n");
for(i=1;i<=n;i++) printf("%6.2f",a[i][i]);
printf("\nEigen vectors:\n");
PrintMatrix(v,n);
//------------------------------------------------------------------
return 0;
}
void PrintMatrix(double a[M][M], int n){
int i,j;
for(i=1;i<=n;i++){
for(j=1;j<=n;j++) printf("%6.2f",a[i][j]);
printf("\n");
}
printf("\n");
}