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

7707次浏览 |255点赞 | 34评论 | 2020-08-01 15:38 更新

来源：互联网

最佳文档

本文由作者推荐

吐成语-清角

#ifndef ksfloat
#define ksfloat double
#endif

#ifndef int16
#define int16 int
#endif

#ifndef uint16
#define uint16 unsigned int
#endif

#ifndef ksnew
#define ksnew(s,t) ((t*)malloc((s)*sizeof(t)))
#endif

#ifndef ksfree
#define ksfree(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( ksfloat aMatrix[] , uint16 side , ksfloat eps , uint16 maxTimes )
{
ksfloat *t, *s, *s1, *vector, *q1
float a1 , c , s2 , t1 , maxEps = 0 , d1 , d2, r = 1;
uint16 i , j , p , q , m , times = 0 , len = 0;

uint16 startIndex = 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 = startIn
dex; 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;
}
}

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

memcpy ( vector , q1 , len*sizeof(ksfloat));
times += 1;
}

// 对特征值与特征向量进行整合
for(i = startIndex; i < endIndex; i++)
{
vector[len+i] = aMatrix[i*side+i];
}

ksfree ( t );
ksfree ( s );
ksfree ( s1 );
ksfree ( q1 );

if(times > maxTimes)
{
ksfree ( vector );
return NULL;
}

return vector;
}