//QR方法求矩阵特征值

#include<iostream>
#include<fstream>
#include<iomanip>
#include<cmath>
using namespace std;
#define N 3   //矩阵的维数
#define NUM 50 //QR分解次数
#define SNUM 16

ofstream fout("QR_out.txt");

void QR(double (*A)[N],double (*Q)[N],double (*R)[N]);               //QR分解
void Multiplicate(double (*A)[N],double (*R)[N],double (*Q)[N]);     //迭代，获得下一次矩阵A=QR
void Print(double (*A)[N],double (*Q)[N],double (*R)[N]);            //矩阵输出
int main()
{
 int i,j,k;
 double a[N][N]={-1,2,1,2,4,-1,1,1,-6};
 double (*A)[N]=a;
 double qq[N][N]={0};
 double (*Q)[N]=qq;
 double rr[N][N]={0};
 double (*R)[N]=rr;
//cout<<"*************************************************************"<<endl;
 //cout<<"输入初始矩阵:\n";
 //cout<<"-------------------------------------------------------------"<<endl;
 //for(i=0;i<N;i++)
  //for(j=0;j<N;j++)
   //fin>>A[i][j];               
 //cout<<"*************************************************************"<<endl;
 fout<<"输出每一步迭代过程: \n";
 fout<<"*************************************************************"<<endl;
 for(k=1;k<=NUM;k++)
 {
  
  QR(A,Q,R);
  fout<<"第"<<k<<"步QR分解:\n";
  fout<<"-------------------------------------------------------------"<<endl;
  Print(A,Q,R);
  Multiplicate(A,R,Q);
 }

 fout<<"矩阵特征值为:\n";
 fout<<"-------------------------------------------------------------"<<endl;
 for(i=0;i<N;i++)   //输出特征值
  fout<<"r["<<(i+1)<<"]="<<A[i][i]<<endl;
 fout<<"-------------------------------------------------------------"<<endl;
 fout<<"矩阵特征向量为:\n";
 fout<<left;
 for(i=0;i<N;i++){
  for(j=0;j<N;j++)
   fout<<setw(SNUM)<<Q[i][j];
  fout<<endl;
 }
 fout<<"*************************************************************"<<endl;
 return 0;
}

void QR(double (*A)[N],double (*Q)[N],double (*R)[N])
{
 int i,j,k,m;
 double temp;
 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);
 }
 return;
}

void Multiplicate(double (*A)[N],double (*R)[N],double (*Q)[N])
{
 int i,j,k;
 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;
  }
 return;
}

void Print(double (*A)[N],double (*Q)[N],double (*R)[N])
{
 int i,j;
  fout<<left;
 fout<<"矩阵A:\n";
 for(i=0;i<N;i++){
  //if(i!=0) cout<<setw(SNUM)<<" ";
  for(j=0;j<N;j++)
   fout<<setw(SNUM)<<A[i][j];
  fout<<endl;
 }
 fout<<"-------------------------------------------------------------"<<endl;
 fout<<"矩阵Q:\n";
 for(i=0;i<N;i++){
  //if(i!=0) cout<<setw(SNUM)<<" ";
  for(j=0;j<N;j++)
   fout<<setw(SNUM)<<Q[i][j];
  fout<<endl;
 }
 fout<<"-------------------------------------------------------------"<<endl;
 fout<<"矩阵R:\n";
 for(i=0;i<N;i++){
  //if(i!=0) cout<<setw(SNUM)<<" ";
  for(j=0;j<N;j++)
   fout<<setw(SNUM)<<R[i][j];
  fout<<endl;
 }
 fout<<"*************************************************************"<<endl;
 return;
}