//QR方法求矩阵特征值
/******
数据
-1 2 1
2 4 -1
1 1 -6
******/
#include<iostream>
#include<fstream>
#include<iomanip>
using namespace std;
#include<math.h>
#define N 3 //矩阵的维数
#define NUM 15 //QR分解次数
#define SNUM 13 //用来控制输出格式
void Print(long double A[N][N],long double Q[N][N],long double R[N][N]);                  //矩阵输出
void QR(long double A[N][N],long double Q[N][N],long double R[N][N]);                    //QR分解
void Multiplicate(long double A[N][N],long double R[N][N],long double Q[N][N]);        //迭代，获得下一次矩阵A=QR
int main()
{
    int i,j;
    long double A[N][N];
    long double Q[N][N]={0};
    long double R[N][N]={0};
    cout<<"*************************************************************"<<endl;
    cout<<"输入初始矩阵:\n";
    cout<<"-------------------------------------------------------------"<<endl;
    for(i=0;i<N;i++)
        for(j=0;j<N;j++)
            cin>>A[i][j];
    cout<<"*************************************************************"<<endl;
    cout<<"输出每一步迭代过程: \n";
    cout<<"*************************************************************"<<endl;
    for(i=1;i<=NUM;i++)
   {
    QR(A,Q,R);
    cout<<"第"<<i<<"步QR分解:\n";
    cout<<"-------------------------------------------------------------"<<endl;
    Print(A,Q,R);
    Multiplicate(A,R,Q);
   }
    cout<<"矩阵特征值为:\n";
    cout<<"-------------------------------------------------------------"<<endl;
    for(i=0;i<N;i++) //输出特征值
    cout<<"r["<<(i+1)<<"]="<<R[i][i]<<endl;
    cout<<"*************************************************************"<<endl;
    return 0;
}
void QR(long double A[N][N],long double Q[N][N],long double R[N][N])
{
    int i,j,k,m;
    long double temp;
    long double a[N],b[N];
    for(j=0;j<N;j++)
   {
        for(i=0;i<N;i++)
       {
          a[i]=A[i][j];
          b[i]=A[i][j];
        }
       for(k=0;k<j;k++)
      {
         R[k][j]=0;
         for(m=0;m<N;m++)
         R[k][j]+=a[m]*Q[m][k];
         for(m=0;m<N;m++)
         b[m]-=R[k][j]*Q[m][k];
      }
      temp=0;
      for(i=0;i<N;i++)
      temp+=b[i]*b[i];
      R[j][j]=sqrt(temp);
      for(i=0;i<N;i++)
      Q[i][j]=b[i]/sqrt(temp);
   }
}
void Multiplicate(long double A[N][N],long double R[N][N],long double Q[N][N])
{
    int i,j,k;
    long double temp;
    for(i=0;i<N;i++)
    for(j=0;j<N;j++)
    {
        temp=0;
        for(k=0;k<N;k++)
        temp+=R[i][k]*Q[k][j];
        A[i][j]=temp;
    }
}
void Print(long double A[N][N],long double Q[N][N],long double R[N][N])
{
    int i,j;
    cout<<left;
    cout<<"矩阵A:\n";
    for(i=0;i<N;i++){
        for(j=0;j<N;j++)
           cout<<setw(SNUM)<<A[i][j];
        cout<<endl;
    }
    cout<<"-------------------------------------------------------------"<<endl;
    cout<<"矩阵Q:\n";
    for(i=0;i<N;i++){
         for(j=0;j<N;j++)
            cout<<setw(SNUM)<<Q[i][j];
    cout<<endl;
    }
    cout<<"-------------------------------------------------------------"<<endl;
    cout<<"矩阵R:\n";
    for(i=0;i<N;i++){
        for(j=0;j<N;j++)
            cout<<setw(SNUM)<<R[i][j];
    cout<<endl;
    }
    cout<<"*************************************************************"<<endl;
}