QR.rar_qr_qr decomposition_qr分解_qr分解C

#include<iostream.h> #include<math.h> #include<iomanip.h> #define Max 10 void Gongshi1(double x[Max][Max],double y[Max][Max],double z[Max][Max],int m,int n) { int i,j,k; for(i=1;i<=m;i++) for(j=1;j<=m;j++) { z[i][j]=0; for(k=1;k<=n;k++) z[i][j]+=x[i][k]*y[k][j]; } } void main() { int i,j,k,n,p; double s,m,b[Max]={0},X[Max]={0},Y[Max]={0},w[Max]={0},I[Max][Max]={0},H[Max][Max]={0},a[Max][Max]={0},Q1[Max][Max][Max]={0},Q[Max][Max]={0},R[Max][Max]={0},w1[Max][Max]={0}; cout<<"请输入矩阵的阶次:"<<endl; cin>>n; cout<<"请输入系数矩阵a:"<<endl; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { cin>>a[i][j]; } } for(i=1;i<=n-1;i++) { for(j=1;j<=n;j++) { for(k=1;k<=n;k++) { if(j==k) { Q[j][j]=1; I[j][j]=1; Q1[i][j][j]=1; } else { I[j][k]=0; Q[j][k]=0; Q1[i][j][k]=0; } } } s=0; for(j=i;j<=n;j++) { s+=pow(a[j][i],2); } cout<<"输出w:"<<endl; for(j=i;j<=n;j++) { if(j==i) { w[j]=a[j][i]-pow(s,0.5); } else w[j]=a[j][i]; cout<<setiosflags(ios::fixed)<<setprecision(6)<<w[j]<<'\t'; } cout<<endl; m=0; for(j=i;j<=n;j++) { m+=pow(w[j],2); } for(j=i;j<=n;j++) { for(k=i;k<=n;k++) { w1[j][k]=w[j]*w[k]; } } cout<<"输出变换矩阵H(w"<<i<<"):"<<endl; for(j=i;j<=n;j++) { for(k=i;k<=n;k++) { H[j][k]=I[j][k]-(2/m)*w1[j][k]; cout<<setiosflags(ios::fixed)<<setprecision(6)<<H[j][k]<<'\t'; Q1[i][j][k]=H[j][k]; Q[j][k]=H[j][k]; } cout<<endl; } cout<<"输出每矩阵Q"<<i<<":"<<endl; for(j=1;j<=n;j++) { for(k=1;k<=n;k++) { cout<<setiosflags(ios::fixed)<<setprecision(6)<<Q[j][k]<<'\t'; } cout<<endl; } Gongshi1(Q,a,R,n,n); for(k=i+1;k<=n;k++) { R[k][i]=0; } cout<<"输出矩阵R("<<i<<"):"<<endl; for(j=1;j<=n;j++) { for(k=1;k<=n;k++) { a[j][k]=R[j][k]; cout<<setiosflags(ios::fixed)<<setprecision(6)<<R[j][k]<<'\t'; } cout<<endl; } } for(i=1;i<n-1;i++) { for(j=1;j<=n;j++) { for(k=1;k<=n;k++) { Q[j][k]=0; for(p=1;p<=n;p++) { Q[j][k]+=Q1[i+1][j][p]*Q1[i][p][k]; } } } for(j=1;j<=n;j++) { for(k=1;k<=n;k++) { Q1[i+1][j][k]=Q[j][k]; } } } cout<<"输出矩阵Q的转置矩阵QT:"<<endl; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { cout<<setiosflags(ios::fixed)<<setprecision(6)<<Q[i][j]<<'\t'; } cout<<endl; } cout<<"请输入常数矩阵b:"<<endl; for(i=1;i<=n;i++) { cin>>b[i]; } for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { Y[i]+=Q[i][j]*b[j]; } } for(i=n;i>=1;i--) { s=0; for(j=n;j>=i+1;j--) { s+=R[i][j]*X[j]; } X[i]=(Y[i]-s)/R[i][i]; } cout<<"输出矩阵的解x:"<<endl; for(i=1;i<=n;i++) cout<<setiosflags(ios::fixed)<<setprecision(6)<<X[i]<<endl; }