C语言实现矩阵的LU分解、施密特正交化、Givens分解、Householder分解.doc
- 文档编号:144829
- 上传时间:2023-04-28
- 格式:DOC
- 页数:25
- 大小:45KB
C语言实现矩阵的LU分解、施密特正交化、Givens分解、Householder分解.doc
《C语言实现矩阵的LU分解、施密特正交化、Givens分解、Householder分解.doc》由会员分享,可在线阅读,更多相关《C语言实现矩阵的LU分解、施密特正交化、Givens分解、Householder分解.doc(25页珍藏版)》请在冰点文库上搜索。
C语言实现矩阵的LU分解、施密特正交化、Givens分解、Householder分解
ByKim.Wang,UCAS
#include
#include
#include
#defineHS10
#defineLS10
intn,m;
floata[HS][LS],bc[HS][LS];
voidgivens()
{
floatfm,sc,cos,sin,r[HS][LS],q[HS][LS],swap[HS][LS],p[HS][LS];
intih,jh,i,j,kh,iw;
for(i=0;i for(j=0;j r[i][j]=a[i][j]; for(i=0;i for(j=0;j { if(i==j) q[i][j]=1; else q[i][j]=0; } for(j=0;j for(i=j+1;i { for(ih=0;ih for(jh=0;jh { if(ih==jh) p[ih][jh]=1; else p[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 for(jh=0;jh { sc=0; for(kh=0;kh sc=sc+p[ih][kh]*r[kh][jh]; swap[ih][jh]=sc; } for(ih=0;ih for(jh=0;jh r[ih][jh]=swap[ih][jh]; for(ih=0;ih for(jh=0;jh { sc=0; for(kh=0;kh sc=sc+p[ih][kh]*q[kh][jh]; swap[ih][jh]=sc; } for(ih=0;ih for(jh=0;jh q[ih][jh]=swap[ih][jh]; } printf("\nGivens约减矩阵Q的结果如下: \n"); for(j=0;j { for(i=0;i printf("%0.5f",q[i][j]); printf("\n"); } printf("\nGivens约减矩阵R的结果如下: \n"); for(i=0;i { for(j=0;j printf("%0.5f",r[i][j]); printf("\n"); } } voidhouseholder() { inti,j,k,is,js,ks,iw,ie,je; floatudm,sc,em,nj,q[HS][LS],r[HS][LS],u[LS],swap[HS][LS],dw[HS][LS],p[HS][LS]; for(i=0;i for(j=0;j r[i][j]=a[i][j]; for(i=0;i for(j=0;j { if(i==j) dw[i][j]=q[i][j]=1; else dw[i][j]=q[i][j]=0; } for(k=0;k { for(i=0;i { if(i u[i]=0; elseif(i==k) { em=0; for(iw=k;iw em=em+r[iw][k]*r[iw][k]; em=sqrt(em); u[i]=r[i][k]-em; } else u[i]=r[i][k]; } for(udm=0,ie=k;ie udm=udm+u[ie]*u[ie]; udm=sqrt(udm); if(udm<0.00001) break; for(ie=0;ie u[ie]=u[ie]/udm; for(ie=0;ie for(je=0;je p[ie][je]=dw[ie][je]-2*u[ie]*u[je]; for(is=0;is for(js=0;js { sc=0; for(ks=0;ks sc=sc+p[is][ks]*q[ks][js]; swap[is][js]=sc; } for(is=0;is for(js=0;js q[is][js]=swap[is][js]; for(is=0;is for(js=0;js { sc=0; for(ks=0;ks sc=sc+p[is][ks]*r[ks][js]; swap[is][js]=sc; } for(is=0;is for(js=0;js r[is][js]=swap[is][js]; } printf("\nHouseholder约减矩阵Q的结果如下: \n"); for(j=0;j { for(i=0;i printf("%0.5f",q[i][j]); printf("\n"); } printf("\nHouseholder矩阵R的结果如下: \n"); for(i=0;i { for(j=0;j printf("%0.5f",r[i][j]); printf("\n"); } } voidschmidt() { inti,j,k,l,b,z; floatem,nj,r[HS][LS],q[HS][LS]; for(i=0;i for(j=0;j q[i][j]=a[i][j]; for(k=0;k { for(l=0;l { nj=0; for(i=0;i nj=nj+q[i][l]*a[i][k]; r[l][k]=nj; for(i=0;i q[i][k]=q[i][k]-r[l][k]*q[i][l]; } em=0; for(i=0;i em=em+q[i][l]*q[i][l]; em=sqrt(em); r[l][k]=em; for(i=0;i q[i][k]=q[i][k]/r[l][k]; } for(i=0;i for(j=n;j q[i][j]=bc[i][j]; printf("\nSchmidt正交化矩阵Q的结果如下: \n"); for(i=0;i { for(j=0;j printf("%0.5f",q[i][j]); printf("\n"); } printf("\nSchmidt正交化矩阵R的结果如下: \n"); for(i=0;i { for(j=0;j { if(i>j)r[i][j]=0; printf("%0.5f",r[i][j]); } printf("\n"); } } voidlufj() { inti,j,k,mm,z,o,p[HS][LS],w[HS]; floatx,y,alu[HS][LS],s[HS],l[HS][LS],u[HS][LS];for(i=0;i for(j=0;j { alu[i][j]=a[i][j]; if(i==j) p[i][j]=1; else p[i][j]=0; } for(j=0;j { i=j;x=alu[i][j];k=i; for(;i { mm=i+1; if(abs(x) { k=mm; x=alu[mm][j]; } } if(k! =j) for(o=0,i=j;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 { y=alu[z][j]/alu[j][j]; alu[z][j]=y; for(o=j+1;o alu[z][o]=alu[z][o]-y*alu[j][o]; } } for(i=0;i for(j=0;j if(i>j) { l[i][j]=alu[i][j]; u[i][j]=0; } elseif(i { 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 { for(j=0;j printf("%0.5f",l[i][j]); printf("\n"); } printf("\nLU分解矩阵U的结果如下: \n"); for(i=0;i { for(j=0;j printf("%0.5f",u[i][j]); printf("\n"); } printf("\nLU分解行交换矩阵P的结果如下: \n"); for(i=0;i { for(j=0;j printf("%d",p[i][j]); printf("\n"); } } intmain() { intim,jm,b,z,zl; printf("请输入矩阵行数M(M需为大于1的整数): \n"); scanf("%d",&m); printf("请输入矩阵列数N(N需为大于1的整数): \n"); scanf("%d",&n); for(im=0;im for(jm=0;jm { b=im+1; z=jm+1; printf("请输入第%d行的第%d个数,以enter键结束输入: \n",b,z); scanf("%f",&a[im][jm]); } floatfm,sc,cos,sin,jdz,r[HS][LS],q[HS][LS],swap[HS][LS],p[HS][LS]; intih,jh,i,j,kh,iw; for(i=0;i for(j=0;j r[i][j]=a[i][j]; for(i=0;i for(j=0;j { if(i==j) q[i][j]=1; else q[i][j]=0; } for(j=0;j for(i=j+1;i { for(ih=0;ih for(jh=0;jh { if(ih==jh) p[ih][jh]=1; else p[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 for(jh=0;jh sc=0; for(kh=0;kh sc=sc+p[ih][kh]*r[kh][jh]; swap[ih][jh]=sc; } for(ih=0;ih for(jh=0;jh r[ih][jh]=swap[ih][jh]; for(ih=0;ih for(jh=0;jh { sc=0; for(kh=0;kh sc=sc+p[ih][kh]*q[kh][jh]; swap[ih][jh]=sc; } for(ih=0;ih for(jh=0;jh bc[jh][ih]=q[ih][jh]=swap[ih][jh]; } if(r[n-1][n-1]<0) jdz=-r[n-1][n-1]; else jdz=r[n-1][n-1]; if(m { printf("\n该矩阵列向量线性相关,无法进行分解! \n"); system("pause"); } elseif(jdz<0.0001) { printf("\n该矩阵列向量线性相关,无法进行分解! \n"); system("pause"); } elseif(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) { case0: return0;break; case1: schmidt();break; case2: householder();break; case3: givens();break; case4: lufj();break; } } while(zl! =0); } else do { 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) { case0: return0;break; case1: schmidt();break; case2: householder();break; case3: givens();break; } } while(zl! =0); }
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 语言 实现 矩阵 LU 分解 施密特 正交 Givens Householder