jacobi方法求取_实对称_矩阵的特征值和特征向量.doc

jacobi方法求取_实对称_矩阵的特征值和特征向量
#ifndefksfloat
#defineksfloatdouble
#endif

#ifndefint16
#defineint16int
#endif

#ifndefuint16
#defineuint16unsignedint
#endif

#ifndefksnew
#defineksnew(s,t)((t*)malloc((s)*sizeof(t)))
#endif

#ifndefksfree
#defineksfree(p)if(p!=0){free(p);p=0;}
#endif

/*
ksJacobiGetRealSymMatrixFeatureValue
用jacobi方法求取*实对称*矩阵的特征值和特征向量
aMatrix[in]:n阶实对称阵
side[in]:矩阵阶大小
eps[in]:控制精度
maxTimes[in]:最大迭代次数
返回值：
前side*side个单元存储特征向量，一列为一特征向量；后side个存储单元为特征值，
也就是总大小为size=side*(side+1)存储单元
*/
ksfloat*ksJacobiGetRealSymMatrixFeatureValue(ksfloataMatrix[],uint16side,ksfloateps,uint16maxTimes)
{
ksfloat*t,*s,*s1,*vector,*q1;
floata1,c,s2,t1,maxEps=0,d1,d2,r=1;
uint16i,j,p,q,m,times=0,len=0;

uint16startIndex=0,endIndex;
endIndex=side+startIndex;
side+=startIndex;

len=side*side;
t=ksnew(len,ksfloat);
s=ksnew(len,ksfloat);
s1=ksnew(len,ksfloat);
q1=ksnew(len,ksfloat);
vector=ksnew(len+side,ksfloat);
memset(vector,0,(len+side)*sizeof(ksfloat));
for(i=startIndex;i<endIndex;i++)
{
vector[i*side+i]=1;//将初始Q[i][j]置为单位矩阵
}

while(r>=eps&&times<maxTimes)
{
p=q=startIndex;maxEps=0;
for(i=startIndex;i<endIndex;i++)
{
for(j=startIndex;j<endIndex;j++)
{
if((j!=i)&&fabs(aMatrix[i*side+j])>maxEps)
{
maxEps=(float)fabs(aMatrix[i*side+j]);//找非对角元素绝对值最大的元
p=i;q=j;
}
}
}

r=maxEps;//重置当前误差
a1=(float)((aMatrix[q*side+q]-aMatrix[p*side+p])/(2*aMatrix[p*side+q]));
if(a1>=0)
{
t1=(float)(1/(fabs(a1)+sqrt(1+a1*a1)));
}
else
{
t1=(float)(-1/(fabs(a1)+sqrt(1+a1*a1)));
}
c=(float)(1/sqrt(1+t1*t1));
s2=t1*c;

memset(s,0,len*sizeof(ksfloat));
for(i=startIndex;i<endIndex;i++)
{
s[i*side+i]=1;//将初始Q[i][j]置为单位矩阵
}

s[p*side+p]=c;s[p*side+q]=s2;
s[q*side+p]=-s2;s[q*side+q]=c;
for(i=startIndex;i<endIndex;i++)
{
for(j=startIndex;j<endIndex;j++)
{
s1[i*side+j]=s[j*side+i];//将矩阵s1[i][j]化为s[i][j]的转置
}
}

for(i=startIndex;i<endIndex;i++)
{
for(j=startIndex;j<endIndex;j++)
{
d1=0;
for(m=startIndex;m<endIndex;m++)
{
d1+=(float)(s1[i*side+m]*aMatrix[m*side+j]);//计算s1[i][j]*aMatrix[i][j]
}
t[i*side+j]=d1;
}