熟悉雅可比法的帮忙了!!!!
程序如下:(运行出现死循环,不知哪儿出了毛病?请指点!)
#define EPSL 1.0e-12 //迭代的精度水平
#define PI 3.1415926
#include <stdio.h>
#include <math.h>
void main() //主程序开始;
{
int size,//协方差矩阵的阶数
row,col,k,
p,q;//记录maxa的行与列
float a[99][99],//存放协方差矩阵
R[99][99]={0},
RT[99][99],//矩阵R的转置
S[99][99],//存放特征向量,赋初值为单位矩阵
tan,//旋转角度的正切
alfa,//旋转角度
maxa=0;//非主对角线上绝对值最大元素
FILE *fp;
fp=fopen("sj.txt","r");
fscanf(fp,"%d",&size);//读入协方差矩阵的阶数
for(row=0;row<size;row++)
for(col=0;col<size;col++)
fscanf(fp,"%f",&a[row][col]);
fclose(fp);
for(row=0;row<size;row++)//给矩阵S赋初值为单位矩阵
for(col=0;col<size;col++)
if(row==col)
S[row][col]=1;
else
S[row][col]=0;
while(maxa==0||maxa>EPSL)
{
for(row=0;row<size;row++)
for(col=row+1;col<size;col++)
if(fabs(a[row][col])>maxa)
{
maxa=a[row][col];//找出非主对角线上绝对值最大的元素
p=row;q=col;//记录maxa的行与列数
}
if(a[p][p]=a[q][q])
tan=PI/2.0;
else
tan=2*a[p][q]/(a[p][p]-a[q][q]);
alfa=atan(tan);//求出旋转角度
R[p][p]=cos(alfa);
R[q][q]=cos(alfa);
R[p][q]=-sin(alfa);
R[q][p]=sin(alfa);//形成矩阵R
for(row=0;row<size;row++)
for(col=0;col<size;col++)
RT[col][row]=R[row][col];//形成R矩阵的转置矩阵RT
for(row=0;row<size;row++)
for(col=0;col<size;col++)
for(k=0;k<size;k++)
a[row][col]+=RT[row][k]*a[k][col]; //计算RT与a的乘积,并赋给a
for(row=0;row<size;row++)
for(col=0;col<size;col++)
for(k=0;k<size;k++)
a[row][col]+=a[row][k]*R[k][col]; //计算a与R的乘积,并赋给a
for(row=0;row<size;row++)
for(col=0;col<size;col++)
for(k=0;k<size;k++)
S[row][col]+=S[row][k]*R[k][col]; //计算S与R的乘积,并赋给S
}
fp=fopen("sc1.txt","w");
for(row=0;row<size;row++)
for(col=0;col<size;col++)
fscanf(fp,"%f",a[row][col]);
fclose(fp);
fp=fopen("sc2.txt","w");
for(row=0;row<size;row++)
for(col=0;col<size;col++)
fscanf(fp,"%f",S[row][col]);
fclose(fp);
}