//===========================================
//矩陣上hessenberg化及QR方法（帶原點位移）
//===========================================
#include <stdio.h>
#include <math.h>
#include <iostream.h>
//--------------------------------------------
void main(){
//變量聲明-----------------------------------
    int n;                                 //矩陣階數
    double a1[20][20],a2[20][20];              
    double u[20][20]={0};
    double h[20][20];
    int i,j,k,l,m;
    double alpha;
    double b;
//輸入---------------------------------------
	cout<<"請輸入矩陣階數:"<<endl;
	cin>>n;
	cout<<"請輸入矩陣元素:"<<endl;
	for(i=1;i<=n;i++){
		for(j=1;j<=n;j++){
			cin>>a1[i][j];
		}
	}
//--------------------------------------------
	for(k=1;k<=n-2;k++){
	//計算alpha
		alpha=0;
		for(i=k+1;i<=n;i++){
			alpha=alpha+a1[i][k]*a1[i][k];
		}
		alpha=sqrt(alpha);
		if(a1[k+1][k]>=0)alpha=-alpha;
	//計算第k次的u向量
		u[k][k+1]=a1[k+1][k]-alpha;
		for(j=k+2;j<=n;j++){
			u[k][j]=a1[j][k];
		}
	//計算b
		b=alpha*(alpha-a1[k+1][k]);
	//計算內部householder矩陣
		for(i=k+1;i<=n;i++){
			for(j=k+1;j<=n;j++){
				if(i==j)h[i][j]=1-u[k][i]*u[k][j]/b;
				else h[i][j]=-u[k][i]*u[k][j]/b;
			}
		}
	//計算經過一次變換后的矩陣
		//i<=k&&j>k時
        for(i=1;i<=k;i++){
			for(j=k+1;j<=n;j++){
				a2[i][j]=0;
				for(l=k+1;l<=n;l++){
					a2[i][j]=a2[i][j]+a1[i][l]*h[l][j];
				}
			}
		}
		//i>k&&j<=k時
		for(i=k+1;i<=n;i++){
			for(j=1;j<=k;j++){
				a2[i][j]=0;
				for(m=k+1;m<=n;m++){
					a2[i][j]=a2[i][j]+h[i][m]*a1[m][j];
				}
			}
		}
        //i>k&&j>k時
		for(i=k+1;i<=n;i++){
			for(j=k+1;j<=n;j++){
				a2[i][j]=0;
				for(l=k+1;l<=n;l++){
					for(m=k+1;m<=n;m++){
						a2[i][j]=a2[i][j]+h[i][m]*a1[m][l]*h[l][j];
					}
				}
			}
		}
		for(i=1;i<=n;i++){
			for(j=1;j<=n;j++){
				if(i<=k&&j<=k)continue;
				a1[i][j]=a2[i][j];
			}
		}
	}
//輸出最終變換后的hessenberg矩陣
	cout<<"Hessenberg矩陣為:"<<endl;
	for(i=1;i<=n;i++){
		for(j=1;j<=n;j++){
			if(fabs(a1[i][j])<0.000000001)a1[i][j]=0;
			cout<<a1[i][j]<<"\t";
		}
		cout<<endl;
	}
//QR方法（帶原點平移）
	{
		double s[20],c[20],d;
		double delta=0.000000001;
		double t;
		for(m=1;fabs(a1[n][n-1])>delta*(fabs(a1[n][n])+fabs(a1[n-1][n-1]));m++){
			t=a1[n][n];
			a1[1][1]=a1[1][1]-t;
			for(k=1;k<=n;k++){
				if(k!=n){
				//確定R(k,k+1)
					if(a1[k+1][k]==0){
						c[k]=1;s[k]=0;d=a1[k][k];
					}
					else{
						d=sqrt(a1[k][k]*a1[k][k]+a1[k+1][k]*a1[k+1][k]);
						c[k]=a1[k][k]/d;
						s[k]=a1[k+1][k]/d;
					}
				//用R(k,k+1)左乘A-tI
					a1[k][k]=d;
					a1[k+1][k]=0;
					a1[k+1][k+1]=a1[k+1][k+1]-t;
					for(j=k+1;j<=n;j++){
						b=a1[k][j]*c[k]+a1[k+1][j]*s[k];
						a1[k+1][j]=-a1[k][j]*s[k]+a1[k+1][j]*c[k];
						a1[k][j]=b;
					}
				}
				if(k!=1){
					for(i=1;i<=k;i++){
						b=a1[i][k-1]*c[k-1]+a1[i][k]*s[k-1];
						a1[i][k]=-a1[i][k-1]*s[k-1]+a1[i][k]*c[k-1];
						a1[i][k-1]=b;
					}
					a1[k-1][k-1]=a1[k-1][k-1]+t;
				}
			}
			a1[n][n]=a1[n][n]+t;
		}
	//輸出矩陣特征值
		cout<<"矩陣特征值為："<<endl;
		for(i=1;i<=n;i++){
			cout<<a1[i][i]<<"\t";
		}
		cout<<endl;
	}
}