矩阵四种分解:PALU,QR,Household,Givens的c程序.docx
- 文档编号:260293
- 上传时间:2022-10-08
- 格式:DOCX
- 页数:10
- 大小:11.07KB
矩阵四种分解:PALU,QR,Household,Givens的c程序.docx
《矩阵四种分解:PALU,QR,Household,Givens的c程序.docx》由会员分享,可在线阅读,更多相关《矩阵四种分解:PALU,QR,Household,Givens的c程序.docx(10页珍藏版)》请在冰豆网上搜索。
#include"stdio.h"#include"malloc.h"#include"math.h"
voidScanMatrix(float*arr,intm,intn); //输入矩阵voidmatrixmul(float*arr1,float*arr2,float*arr3,intm,intn,intk); //矩阵数乘函数
voidprintm(charc,float*arr,intm,intn); //矩阵打印函数voidPALU(float*arr,intm,intn); //PA=LU分解
函数
voidGram_Schmint(float*arr,intm,intn); //QR分解函数
voidHousehold(float*arr,intm,intn); //household
分解函数
voidGivens(float*arr,intm,intn); //givens分解
函数
voiddivide(); //分解
矩阵总函数
intmain()
{
divide();return0;
}
voiddivide()
{
intm,n,method,flag=1;charc='0';
while(c!
='#')
{
printf("请输入mXn矩阵维度m:
");scanf("%d",&m);
printf("请输入mXn矩阵维度n:
");scanf("%d",&n);
float*arr=(float*)malloc(sizeof(float)*m*n);ScanMatrix(arr,m,n);
while(flag)
{
printf("请输入矩阵的分解方式:
\n1-PA=LU,\n2-Gram-Schmit正交化A=QR,\n3-Household
分解A=QU,\n4-Givens分解A=QU。
\n5-结束分解.\n");scanf("%d",&method);
switch(method)
{
case1:
PALU(arr,m,n);break;
case2:
Gram_Schmint(arr,m,n);break;
case3:
Household(arr,m,n);break;
case4:
Givens(arr,m,n);break;
case5:
flag=0;break;
}
}
flag=1;
printf("输入#结束,否则继续重新输入矩阵开始分解:
");getchar();
c=getchar();free(arr);
}
}
voidPALU(float*arr,intm,intn)
{
inti,j,k,p,q,exsit=1,change1=0,change2=0;floatt,s;
float*arrt=(float*)malloc(sizeof(float)*m*n);float*P=(float*)malloc(sizeof(float)*m*m);float*L=(float*)malloc(sizeof(float)*m*m);for(i=0;i for(j=0;j { *(arrt+i*n+j)=*(arr+i*n+j); } for(i=0;i for(j=0;j { *(P+i*m+j)=*(L+i*m+j)=0; } for(i=0;i { *(P+i*m+i)=*(L+i*m+i)=1; } for(i=0;i { if(*(arrt+i*n+i)==0)//判断主元是否为0; { for(j=i+1;j { if(*(arrt+j*n+i)>0&&j { for(k=0;k { t=*(arrt+j*n+k); //arr的i行与j行逐个进行值交换; *(arrt+j*n+k)=*(arrt+i*n+k); *(arrt+i*n+k)=t;change1=i;change2=j; } for(k=0;k { } break; } if(j==n) { t=*(P+j*m+k); //P的i行与j行逐个进行值交换; *(P+j*m+k)=*(P+i*m+k); *(P+i*m+k)=t; exsit=0; printf("矩阵不满秩,不可分解! \n");//不存在非零主元;break; } } } for(p=i+1;p { s=*(arrt+p*n+i); *(arrt+p*n+i)=0; *(L+p*m+i)=s/(*(arrt+i*n+i)); for(q=i+1;q { *(arrt+p*n+q)=*(arrt+p*n+q)-(s)*(*(arrt+i*n+q))/(*(arrt+i*n+i));//消元过程 } } for(k=0;k { if(change1! =0&&k { t=*(L+change2*m+k); //L的i行与j行逐个进行值交换; *(L+change2*m+k)=*(L+change1*m+k); *(L+change1*m+k)=t; } } change1=0;change2=0; } if(exsit) { printf("PA=LU的分解矩阵如下: \n");printm('P',P,m,m); printm('L',L,m,m); printm('U',arrt,m,n); } free(L); free(P);free(arrt); } voidGram_Schmint(float*arr,intm,intn) { inti,j,l; floatmod,ip; float*arrt=(float*)malloc(sizeof(float)*m*n);float*Q=(float*)malloc(sizeof(float)*m*n);float*R=(float*)malloc(sizeof(float)*n*n);for(i=0;i for(j=0;j { *(arrt+i*n+j)=*(arr+i*n+j); *(Q+i*n+j)=0; } for(i=0;i for(j=0;j { *(R+i*n+j)=0; } for(j=0;j { mod=0;for(l=0;l { ip=0; for(i=0;i { ip=ip+(*(arr+i*n+j))*(*(Q+i*n+l)); } *(R+l*n+j)=ip; //i行l行的Q值; for(i=0;i { *(arrt+i*n+j)=*(arrt+i*n+j)-ip*(*(Q+i*n+l)); } } for(i=0;i { mod=sqrt((*(arrt+i*n+j))*(*(arrt+i*n+j))+mod*mod); } if(! mod) { printf("A矩阵并不是行间线性无关! ! 无法QR分解! ! \n");break; } *(R+j*n+j)=mod;for(i=0;i { *(Q+i*n+j)=(*(arrt+i*n+j))/mod; } } if(mod) { printf("A=QR的分解矩阵如下: \n");printm('Q',Q,m,n); printm('R',R,n,n); } free(Q); free(R);free(arrt); } voidHousehold(float*arr,intm,intn) { inti,j,l,flag=1;floatmod,modu; float*arrt=(float*)malloc(sizeof(float)*m*n);float*H=(float*)malloc(sizeof(float)*m*m);float*R=(float*)malloc(sizeof(float)*m*m);float*RT=(float*)malloc(sizeof(float)*m*m);float*Rsum=(float*)malloc(sizeof(float)*m*m); for(i=0;i for(j=0;j { *(arrt+i*n+j)=*(arr+i*n+j); } for(i=0;i for(j=0;j { *(H+i*m+j)=*(R+i*m+j)=0; } for(i=0;i { *(H+i*m+i)=*(R+i*m+i)=1; } for(j=0;j { mod=0;modu=0; for(i=j;i { mod=sqrt((*(arrt+i*n+j))*(*(arrt+i*n+j))+mod*mod); } if(*(arrt+j*n+j)-mod) *(arrt+j*n+j)=*(arrt+j*n+j)-mod; else if((*(arrt+j*n+j)+mod)! =0) *(arrt+j*n+j)=*(arrt+j*n+j)+mod;else { printf("矩阵分解有误,存在全为0元素列! ");flag=0; break; } for(i=j;i { modu=sqrt((*(arrt+i*n+j))*(*(arrt+i*n+j))+modu*modu); } for(i=j;i { *(R+i*m+l)=-2*(*(arrt+i*n+j))*(*(arrt+l*n+j))/modu/modu; } for(i=j;i { *(R+i*m+i)=1+(*(R+i*m+i)); } for(i=0;i for(l=0;l { if(i==l) { } else { } } *(R+i*m+l)=1; *(R+i*m+l)=*(R+l*m+i)=0; matrixmul(R,H,Rsum,m,m,m);for(i=0;i for(l=0;l { *(H+i*m+l)=*(Rsum+i*m+l); } matrixmul(Rsum,arr,arrt,m
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 分解 PALU QR Household Givens 程序