实对称矩阵特征值求解算法：Jacobi行循环法

做过程故障诊断，需要计算含小于100个变量的数据的主成分，下面这个代码还是合适的，特此备份。经测试，当矩阵A阶数n=100时，耗费时间在0.2s左右。

int jacobi_loop(double** A, double** V, double* eigsv, double epsl, int maxt, int n){
/************************************************************************
* 作者: taoshaohui from Qingdao University of Science & Technology, 2010/08/17
* 函数名称: jacobi_loop
* 函数功能: Jacobi行循环法求取n阶实对称矩阵A的所有特征值和特征向量V(:,i)
* 输入参数:

A: n阶实对称矩阵, 算法结束后其第j个对角线元素为A的第j个特征值

n: A的阶数
V: n阶单位矩阵, 算法结束后其第j列存储对应第j个特征值的特征向量
eigsv:按从大到小顺序排列的特征值向量
epsl: 迭代收敛标准, 当A所有非对角线元素绝对值皆小于epsl时算法成功
maxt: 算法最大迭代次数

* 输出参数: V, eigsv
* 返回值 : 整数success, 0迭代未收敛/1迭代成功
* 注意 : 本程序中，2×2Givens旋转矩阵为: [cn, sn; -sn, cn]
* 文献 : 《矩阵计算》, G. H. Golub & C. F. Van Loan, 2002, 袁亚湘 译, P 494.
************************************************************************/
int success=0; // 函数返回值
double tao, t, cn, sn; // 临时变量
double maxa;
for (int it=0; it<maxt; it++)
{
maxa=0;
for (int p=0; p<n-1; p++)
{
for (int q=p+1; q<n; q++)
{
if (fabs(A[p][q])>maxa) // 记录非对角线元素最大值
{
maxa=fabs(A[p][q]);
}

if (fabs(A[p][q])>epsl) // 非对角线元素非0时才执行Jacobi变换
{
// 计算Givens旋转矩阵的重要元素:cos(theta), 即cn, sin(theta), 即sn
tao=0.5*(A[q][q]-A[p][p])/A[p][q];

if (tao>=0) // t为方程 t^2 + 2*t*tao - 1 = 0的根, 取绝对值较小的根为t
{
t=-tao+sqrt(1+tao*tao);
}
else
{
t=-tao-sqrt(1+tao*tao);
}
cn=1/sqrt(1+t*t);
sn=t*cn;

// Givens旋转矩阵之转置左乘A, 即更新A的p行和q行
for (int j=0; j<n; j++)
{
double apj=A[p][j];
double aqj=A[q][j];
A[p][j]=cn*apj-sn*aqj;
A[q][j]=sn*apj+cn*aqj;
}

// Givens旋转矩阵右乘A, 即更新A的p列和q列
for (int i=0; i<n; i++)
{
double aip=A[i][p];
double aiq=A[i][q];
A[i][p]=cn*aip-sn*aiq;
A[i][q]=sn*aip+cn*aiq;
}

// 更新特征向量存储矩阵V, V=J0×J1×J2...×Jit, 所以每次只更新V的p, q两列
for (int i=0; i<n; i++)
{
double vip=V[i][p];
double viq=V[i][q];
V[i][p]=cn*vip-sn*viq;
V[i][q]=sn*vip+cn*viq;
}
}

}

}

if (maxa<epsl) // 非对角线元素已小于收敛标准，迭代结束
{
// 特征值向量
for (int j=0; j<n; j++)
{
eigsv[j]=A[j][j];
}
// 对特征值向量从大到小进行排序, 并调整特征向量顺序 (直接插入法)

double* tmp=new double
;
for (int j=1; j<n; j++)
{
int i=j;
double a=eigsv[j];
for (int k=0; k<n; k++)
{
tmp[k]=V[k][j];
}
while(i>0 && eigsv[i-1]<a){
eigsv[i]=eigsv[i-1];
for (int k=0; k<n; k++)
{
V[k][i]=V[k][i-1];
}
i--;
}
eigsv[i]=a;
for (int k=0; k<n; k++)
{
V[k][i]=tmp[k];
}
}

delete[] tmp;

return success=1;
}
}
return success;
}