C语言实现矩阵的LU分解、施密特正交化、Givens分解、Householder分解By Kim.Wang，UCAS#include<stdio.h>#include<math.h>#include<windows.h>#define HS 10#define LS 10int n, m;float a[HS][LS],bc[HS][LS];void givens(){float fm,sc,cos,sin,r[HS][LS],q[HS][LS],swap[HS][LS],p[HS][LS];int ih,jh,i, j,kh,iw;for (i = 0; i < m; i++)for (j = 0; j < n; j++)r[i][j]=a[i][j];for (i = 0; i < m; i++)for (j = 0; j < m; j++){if(i==j)q[i][j]=1;elseq[i][j]=0;}for(j=0;j<n;j++)for (i = j+1; i < m; i++){for (ih = 0; ih < m; ih++)for (jh = 0; jh < m; jh++){if(ih==jh)p[ih][jh]=1;elsep[ih][jh]=0;}fm=sqrt(r[i][j]*r[i][j]+r[j][j]*r[j][j]);cos=r[j][j]/fm;sin=r[i][j]/fm;p[i][i]=p[j][j]=cos;p[i][j]=-sin;p[j][i]=sin;for (ih = 0; ih < m; ih++)for (jh = 0; jh < n; jh++){sc=0;for (kh = 0; kh < m; kh++)sc=sc+p[ih][kh]*r[kh][jh];swap[ih][jh]=sc;}for (ih = 0; ih < m; ih++)for (jh = 0; jh < n; jh++)r[ih][jh]=swap[ih][jh];for (ih = 0; ih < m; ih++)for (jh = 0; jh < m; jh++){sc=0;for (kh = 0; kh < m; kh++)sc=sc+p[ih][kh]*q[kh][jh];swap[ih][jh]=sc;}for (ih = 0; ih < m; ih++)for (jh = 0; jh < m; jh++)q[ih][jh]=swap[ih][jh];}printf("\n Givens约减矩阵Q的结果如下：\n");for (j = 0; j < m; j++){for (i = 0; i < m; i++)printf("%0.5f ", q[i][j]);printf("\n");}printf("\n Givens约减矩阵R的结果如下：\n");for (i = 0; i < m; i++){for (j = 0; j < n; j++)printf("%0.5f ", r[i][j]);printf("\n");}}void householder(){int i, j, k, is,js,ks,iw, ie,je;float udm,sc,em,nj,q[HS][LS],r[HS][LS],u[LS],swap[HS][LS],dw[HS][LS],p[HS][LS];for (i = 0; i < m; i++)for (j = 0; j < n; j++)r[i][j]=a[i][j];for (i = 0; i < m; i++)for (j = 0; j < m; j++){if(i==j)dw[i][j]=q[i][j]=1;elsedw[i][j]=q[i][j]=0;}for(k=0;k<n;k++){for (i = 0; i < m; i++){if(i<k)u[i]=0;else if(i==k){em=0;for(iw=k;iw<m;iw++)em=em+r[iw][k]*r[iw][k];em=sqrt(em);u[i]=r[i][k]-em;}elseu[i]=r[i][k];}for (udm=0,ie = k; ie < m; ie++)udm=udm+u[ie]*u[ie];udm=sqrt(udm);if(udm<0.00001)break;for (ie = 0; ie < m; ie++)u[ie]=u[ie]/udm;for (ie = 0; ie < m; ie++)for (je = 0; je < m; je++)p[ie][je]=dw[ie][je]-2*u[ie]*u[je];for (is = 0; is < m; is++)for (js = 0; js < m; js++){sc=0;for (ks = 0; ks < m; ks++)sc=sc+p[is][ks]*q[ks][js];swap[is][js]=sc;}for (is = 0; is < m; is++)for (js = 0; js < m; js++)q[is][js]=swap[is][js];for (is = 0; is < m; is++)for (js = 0; js < m; js++){sc=0;for (ks = 0; ks < m; ks++)sc=sc+p[is][ks]*r[ks][js];swap[is][js]=sc;}for (is = 0; is < m; is++)for (js = 0; js < m; js++)r[is][js]=swap[is][js];}printf("\n Householder约减矩阵Q的结果如下：\n");for (j = 0; j < m; j++){for (i = 0; i < m; i++)printf("%0.5f ", q[i][j]);printf("\n");}printf("\n Householder矩阵R的结果如下：\n");for (i = 0; i < m; i++){for (j = 0; j < n; j++)printf("%0.5f ", r[i][j]);printf("\n");}}void schmidt(){int i, j, k, l, b, z;float em,nj, r[HS][LS],q[HS][LS];for (i = 0; i < m; i++)for (j = 0; j < n; j++)q[i][j]=a[i][j];for(k=0;k<n;k++){for(l=0;l<k;l++){nj=0;for(i=0;i<m;i++)nj=nj+q[i][l]*a[i][k];r[l][k]=nj;for(i=0;i<m;i++)q[i][k]=q[i][k]-r[l][k]*q[i][l];}em=0;for(i=0;i<m;i++)em=em+q[i][l]*q[i][l];em=sqrt(em);r[l][k]=em;for(i=0;i<m;i++)q[i][k]=q[i][k]/r[l][k];}for(i=0;i<m;i++)for(j=n;j<m;j++)q[i][j]=bc[i][j];printf("\n Schmidt正交化矩阵Q的结果如下：\n");for (i = 0; i < m; i++){for (j = 0; j < m; j++)printf("%0.5f ", q[i][j]);printf("\n");}printf("\n Schmidt正交化矩阵R的结果如下：\n");for (i = 0; i < m; i++){for (j = 0; j < n; j++){if(i>j) r[i][j]=0;printf("%0.5f ", r[i][j]);}printf("\n");}}void lufj(){int i, j, k,mm, z, o, p[HS][LS], w[HS];float x, y,alu[HS][LS], s[HS],l[HS][LS],u[HS][LS]; for (i = 0; i < n; i++)for (j = 0; j < n; j++){alu[i][j]=a[i][j];if (i == j)p[i][j] = 1;elsep[i][j] = 0;}for (j = 0; j < n; j++){i = j; x = alu[i][j]; k = i;for (; i < n - 1; i++){mm = i + 1;if (abs(x) < abs(alu[mm][j])){k = mm;x = alu[mm][j];}}if (k != j)for (o = 0,i = j; o < n; o++){s[o] = alu[i][o];alu[i][o] = alu[k][o];alu[k][o] = s[o];w[o] = p[i][o];p[i][o] = p[k][o];p[k][o] = w[o];}for (z = j + 1; z < n; z++){y = alu[z][j] / alu[j][j];alu[z][j] = y;for (o = j+1; o < n; o++)alu[z][o] = alu[z][o] - y*alu[j][o];}}for (i = 0; i < n; i++)for (j = 0; j < n; j++)if(i>j){l[i][j]=alu[i][j];u[i][j]=0;}else if(i<j){l[i][j]=0;u[i][j]=alu[i][j];}else{l[i][j]=1;u[i][j]=alu[i][j];}printf(" LU分解矩阵L的结果如下：\n");for (i = 0; i < n; i++){for (j = 0; j < n; j++)printf("%0.5f ", l[i][j]);printf("\n");}printf("\n LU分解矩阵U的结果如下：\n");for (i = 0; i < n; i++){for (j = 0; j < n; j++)printf("%0.5f ", u[i][j]);printf("\n");}printf("\n LU分解行交换矩阵P的结果如下：\n");for (i = 0; i < n; i++){for (j = 0; j < n; j++)printf("%d ", p[i][j]);printf("\n");}}int main(){int im,jm,b,z,zl;printf("请输入矩阵行数M(M需为大于1的整数)：\n");scanf("%d",&m);printf("请输入矩阵列数N(N需为大于1的整数)：\n");scanf("%d",&n);for (im = 0; im < m; im++)for (jm = 0; jm < n; jm++){b = im + 1;z = jm + 1;printf("请输入第%d行的第%d个数，以enter键结束输入：\n", b, z);scanf("%f", &a[im][jm]);}float fm,sc,cos,sin,jdz,r[HS][LS],q[HS][LS],swap[HS][LS],p[HS][LS];int ih,jh,i, j,kh,iw;for (i = 0; i < m; i++)for (j = 0; j < n; j++)r[i][j]=a[i][j];for (i = 0; i < m; i++)for (j = 0; j < m; j++){if(i==j)q[i][j]=1;elseq[i][j]=0;}for(j=0;j<n;j++)for (i = j+1; i < m; i++){for (ih = 0; ih < m; ih++)for (jh = 0; jh < m; jh++){if(ih==jh)p[ih][jh]=1;elsep[ih][jh]=0;}fm=sqrt(r[i][j]*r[i][j]+r[j][j]*r[j][j]);cos=r[j][j]/fm;sin=r[i][j]/fm;p[i][i]=p[j][j]=cos;p[i][j]=-sin;p[j][i]=sin;for (ih = 0; ih < m; ih++)for (jh = 0; jh < n; jh++)sc=0;for (kh = 0; kh < m; kh++)sc=sc+p[ih][kh]*r[kh][jh];swap[ih][jh]=sc;}for (ih = 0; ih < m; ih++)for (jh = 0; jh < n; jh++)r[ih][jh]=swap[ih][jh];for (ih = 0; ih < m; ih++)for (jh = 0; jh < m; jh++){sc=0;for (kh = 0; kh < m; kh++)sc=sc+p[ih][kh]*q[kh][jh];swap[ih][jh]=sc;}for (ih = 0; ih < m; ih++)for (jh = 0; jh < m; jh++)bc[jh][ih]=q[ih][jh]=swap[ih][jh];}if(r[n-1][n-1]<0)jdz=-r[n-1][n-1];elsejdz=r[n-1][n-1];if(m<n){printf("\n 该矩阵列向量线性相关，无法进行分解! \n");system("pause");}else if(jdz<0.0001){printf("\n 该矩阵列向量线性相关，无法进行分解! \n");system("pause");}else if(m==n){do{printf("\n");printf("该矩阵可进行四种分解，请输入选项前的数字选择！\n");printf("0：退出程序\n");printf("1：Gram-Schmidt正交化\n");printf("2：Householder分解\n");printf("3：Givens分解\n");printf("4：LU分解\n\n");scanf("%d", &zl);switch(zl){case 0:return 0; break;case 1:schmidt(); break;case 2:householder(); break;case 3:givens(); break;case 4:lufj(); break;}}while(zl!=0);}elsedo{printf("\n");printf("该矩阵可进行除LU分解外的三种分解，请输入选项前的数字选择！\n");printf("0：退出程序\n");printf("1：Gram-Schmidt正交化\n");printf("2：Householder分解\n");printf("3：Givens分解\n\n");scanf("%d", &zl);switch(zl){case 0:return 0; break;case 1:schmidt(); break;case 2:householder(); break;case 3:givens(); break;}}while(zl!=0);}。