/*
题目：第2题 对于给定的15阶三对角矩阵，用QR方法求它的全部特征值
时间：2010-7-6
作者：
方法：基于Householder变换的QR分解，分解后Q为正交矩阵，R为上三角矩阵
并且矩阵R的主对角元素即为所求的特征值
*/
#include"iostream.h"
#include"iomanip.h"
#include"stdlib.h"
#include"math.h"
double A[15][15];/*定义全局变量用来记录需要分解的矩阵A*/
double Q[15][15];/*定义全局变量用来记录分解后的正交矩阵*/
/*函数用来实现给矩阵A初始赋值
@param 无——需要赋值的A为全局变量
@return 无——赋值后的A为全局变量
*/
void evaluationMatrixA()
{
for(int i=0;i<15;i++)
{
if(i==0)//第一行没有下次对角元素
{
A[0][0]=4;
A[0][1]=1;
}
else if(i==14)//最后一行没有上次对角元素
{
A[i][i-1]=1;
A[i][i]=4;
}
else //第2行到第14行主对角元素为4，上下次对角元素为1
{
A[i][i-1]=1;
A[i][i]=4;
A[i][i+1]=1;
}
}
}
/*函数用来实现输出矩阵A的全部元素
@param 无——需要处理的A为全局变量
@return 无——处理后的A为全局变量
*/
void printMatrixA()
{
for(int i=0;i<15;i++)
{
cout<<" ";
for(int j=0;j<15;j++)
{
cout<<" "<<A[i][j]<<" ";
}
cout<<endl;
}
}
/*函数用来实现矩阵A的QR分解
@param 无——需要处理的A为全局变量
@return 无——处理后的正交矩阵Q与上三三角矩阵R都是全局变量
步骤：1.做Householder变换，共需要循环进行15步
2.用Q(K)左乘Q,即Q=Q(k)*Q
3.用Q(K)左乘A,即A=Q(k)*A
4.得到上三角矩阵R(即变换后的A)
5.得到Q的转置矩阵
*/
void SplitQR()
{
int i,j,k,l;//定义的循环变量用来记录矩阵的各个元素
double max,sum,w,a,v,p,t,s,temp;//计算公式中的各变量
double u[15];//记录Householder变换后Q(k)中u[i]，i=k,k+1,……，15
/*做Householder变换，共需要进行15步
变换后，Q=|I(k) 0 |
| 0 Q(n-k)|n*n
其中I(k)为单位矩阵，n为矩阵阶数，
Q(n-k)的主对角线为(1-2*u(k)*u(k)),其他元素为（-2*u(i)*u(j)）
u(j)（j=k,k+1,……,n-1）
*/
for(k=0;k<15;k++)
{
for(i=0;i<k;i++)//做单位矩阵I(k)
{
for(j=0;j<k;j++)
{
if(i==j)
{
Q[i][j]=1.0;
}
else
{
Q[i][j]=0.0;
}
}
}//end for
/*计算公式：w=max|a(i)(k)|,(k<=i<=n-1)
k为第k次变换，n为矩阵阶数,
a(i)(k)为矩阵A中的元素
调用math.h中fabs()函数，来实现选最大元素
*/
for(i=k;i<15;i++)
{
w=0.0;
for(j=k;j<15;j++)
{
max=fabs(A[i][j]);//选最大的元素
if(max>w)
w=max;
}//end for
}//end for
/*公式：sum=(a(k)(k)*w)*(a(k)(k)*w)+(a(k+1)(k)*w)*(a(k+1)(k)*w)+…
…+(a(m-1)(k)*w)*(a(m-1)(k)*w)，
a=-sign(a(k)(k))*sqrt(sum)*w,
其中a(i)(k)为矩阵A中的元素，k为第k次变换，n为矩阵阶数,
sign()为符号函数
*/
sum=0.0;
for(i=k;i<15;i++)//计算sum
{
v=A[i][k]/w;
sum=sum+v*v;
}
if(A[k][k]>0)//实现符号函数sign()
{
w=-w;
}
a=w*sqrt(sum);//计算公式a=-sign(a(k)(k))*sqrt(sum)*w
/*
公式：p=sqrt(2*a*(a-a[k][k])),k为第k次变换，sqrt()为开二次方的函数
*/
p=sqrt(2*a*(a-A[k][k]));
if((p+1.0)!=1.0)
{
/*
公式：u[k]=(A[k][k]-a)/p，
u[i]=A[i][k]/p，k+1<=i<=n-1
*/
u[k]=(A[k][k]-a)/p;
for(i=k+1;i<15;i++)
{
u[i]=A[i][k]/p;
}//end for
/*对Q(k)的赋值完成
Q(n-k)的主对角线为(1-2*u(k)*u(k)),其他元素为（-2*u(i)*u(j)）
*/
for(i=k;i<15;i++)
{
for(j=k;j<15;j++)
{
if(i==j)
{
Q[i][j]=1-2*u[i]*u[i];
}
else
{
Q[i][j]=-2*u[i]*u[j];
}//end if=else
}//end for
}//end for
/*用Q(K)左乘Q,即Q=Q(k)*Q,
公式：t=t+A[l][k]*Q[l][j]，k<=l<=n-1,0<=j<=n-1,
Q[i][j]=Q[i][j]-2*A[i][k]*t，k<=i<=n-1,0<=j<=n-1
*/
for(j=0;j<15;j++)
{
t=0.0;
for(l=k;l<15;l++)
{
t=t+A[l][k]*Q[l][j];
}
for(i=k;i<15;i++)
{
Q[i][j]=Q[i][j]-2*A[i][k]*t;
}
}
/*用Q(K)左乘A,即A=Q(k)*A
公式：s=A[l][k]*A[l][j]，k<=l<=n-1,k+1<=j<=n-1
A[i][j]=A[i][j]-2*A[i][k]*s，k<=i<=n-1,k+1<=j<=n-1
*/
for(j=k+1;j<15;j++)
{
s=0.0;
for(l=k;l<15;l++)
{
s=A[l][k]*A[l][j];
}//end for
for(i=k;i<15;i++)
{
A[i][j]=A[i][j]-2*A[i][k]*s;
}
}//end for
/*
公式：A[k][k]=a，并将矩阵A的第k列下三角部分置为0
*/
A[k][k]=a;
for(i=k+1;i<15;i++)
{
A[i][k]=0.0;
}//end for
}//end if
}//end for
/*转置得到矩阵Q*/
for(i=0;i<15;i++)
{
for(j=i+1;j<15;j++)
{
temp=Q[i][j];
Q[i][j]=Q[j][i];
Q[j][i]=temp;
}//end for
}//end for
}//QR分解
/*主函数*/
int main(int argc,char* argv[])
{
evaluationMatrixA();//给矩阵A赋值
cout<<" 需要做QR分解的矩阵A为："<<endl;
printMatrixA();//输出矩阵A
SplitQR();//QR分解
cout<<endl<<" 基于Householder变换的QR分解,根据QR分解 "<<endl<<
" 后得到R矩阵主对角元素为矩阵A的特征值，"<<
endl<<" 即所求的全部的特征值为："<<endl;
for(int i=0;i<15;i++)
{
cout<<" "<<"det("<<i+1<<")="<<A[i][i]<<endl;
}//end for
return 0;
}//end main
