矩阵四种分解：PALU,QR,Household,Givens的c程序

modu=sqrt((*(arrt+i*n+j))*(*(arrt+i*n+j))+modu*modu); } for(i=j;i<m;i++) //求单次的R矩阵 for(l=j;l<m;l++) { *(R+i*m+l)=-2*(*(arrt+i*n+j))*(*(arrt+l*n+j))/modu/modu; } for(i=j;i<m;i++) { *(R+i*m+i)=1+(*(R+i*m+i)); } for(i=0;i<j;i++) for(l=0;l<m;l++) { if(i==l) { *(R+i*m+l)=1; } else { *(R+i*m+l)=*(R+l*m+i)=0; } } matrixmul(R,H,Rsum,m,m,m); for(i=0;i<m;i++) for(l=0;l<m;l++) { *(H+i*m+l)=*(Rsum+i*m+l); } matrixmul(Rsum,arr,arrt,m,m,n); } for(i=0;i<m;i++) for(l=0;l<m;l++) { *(RT+i*m+l)=*(Rsum+l*m+i); }
wk.baidu.com
{ case 1: PALU(arr,m,n); break; case 2: Gram_Schmint(arr,m,n); break; case 3: Household(arr,m,n); break; case 4: Givens(arr,m,n); break; case 5: flag=0; break; } } flag=1; printf("输入#结束，否则继续重新输入矩阵开始分解:"); getchar(); c=getchar(); free(arr); } } void PALU(float *arr,int m,int n) { int i,j,k,p,q,exsit=1,change1=0,change2=0; float t,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<m;i++) for(j=0;j<n;j++) { *(arrt+i*n+j)=*(arr+i*n+j);
if(n&&flag) { printf("A=QT的分解矩阵如下（其中Q为PA=T中P的逆矩阵）：\n"); printm('Q',RT,m,m); printm('T',arrt,m,n); printm('P',Rsum,m,m); } free(arrt); free(H); free(RT); free(Rsum); free(R); } void Givens(float *arr,int m,int n) { int i,j,l,k,flag=1; float mod,modu,c,s; float *arrt= (float *)malloc(sizeof (float)* m* n); float *G= (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<m;i++) for(j=0;j<n;j++) { *(arrt+i*n+j)=*(arr+i*n+j); } for(i=0;i<m;i++) for(j=0;j<m;j++) { *(G+i*m+j)=*(R+i*m+j)=0; } for(i=0;i<m;i++) { *(G+i*m+i)=*(R+i*m+i)=1; }
exsit=0; printf("矩阵不满秩,不可分解!\n");//不存在非零主元; break; } } } for(p=i+1;p<m;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<n;q++)//对所选列所有的元素开始消元； { *(arrt+p*n+q)=*(arrt+p*n+q)-(s)*(*(arrt+i*n+q))/(*(arrt+i*n+i));// 消元过程 } } for(k=0;k<m;k++) { if(change1!=0&&k<change1) { 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);
float *Rsum= (float *)malloc(sizeof (float) * m*m); for(i=0;i<m;i++) //矩阵初始化 for(j=0;j<n;j++) { *(arrt+i*n+j)=*(arr+i*n+j); } for(i=0;i<m;i++) //矩阵初始化 for(j=0;j<m;j++) { *(H+i*m+j)=*(R+i*m+j)=0; } for(i=0;i<m;i++) { *(H+i*m+i)=*(R+i*m+i)=1; } for(j=0;j<n;j++) //对第A（n-j)阶子矩阵操作 { mod=0; modu=0; for(i=j;i<m;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<m;i++) //求列模值 {
#include "stdio.h" #include "malloc.h" #include "math.h" void ScanMatrix(float *arr,int m,int n); //输入矩阵 void matrixmul(float*arr1,float*arr2,float*arr3,int m,int n,int k); //矩阵 数乘函数 void printm(char c,float*arr,int m,int n); //矩阵打印函数 void PALU(float *arr,int m,int n); //PA=LU分解函数 void Gram_Schmint(float *arr,int m,int n); //QR分解函数 void Household(float *arr,int m,int n); //household分解 函数 void Givens(float *arr,int m,int n); //givens分解函数 void divide(); //分解矩阵总函数 int main() { divide(); return 0; } void divide() { int m,n,method,flag=1; char c='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)
*(arrt+i*n+j)=*(arrt+i*n+j)-ip*(*(Q+i*n+l)); } } for(i=0;i<m;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<m;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); } void Household(float *arr,int m,int n) { int i,j,l,flag=1; float mod,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);
} for(i=0;i<m;i++) for(j=0;j<m;j++) { *(P+i*m+j)=*(L+i*m+j)=0; } for(i=0;i<m;i++) { *(P+i*m+i)=*(L+i*m+i)=1; } for(i=0;i<m&&i<n&&exsit==1;i++) { if(*(arrt+i*n+i)==0) //判断主元是否为0； { for(j=i+1;j<m+1;j++)//对主元下面的元素进行考察； { if(*(arrt+j*n+i)>0&&j<m)//存在大于0的主元； { for(k=0;k<n;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<m;k++) { t=*(P+j*m+k); //P的i行与j行逐个进行值交换； *(P+j*m+k)=*(P+i*m+k); *(P+i*m+k)=t; } break; } if(j==n) {
} free(L); free(P); free(arrt); } void Gram_Schmint(float *arr,int m,int n) { int i,j,l; float mod,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<m;i++) for(j=0;j<n;j++) { *(arrt+i*n+j)=*(arr+i*n+j); *(Q+i*n+j)=0; } for(i=0;i<n;i++) for(j=0;j<n;j++) { *(R+i*n+j)=0; } for(j=0;j<n;j++) { mod=0; for(l=0;l<j;l++) { ip=0; for(i=0;i<m;i++) //A[*][j]与Q[*][l]做一次内积； { ip=ip+(*(arr+i*n+j))*(*(Q+i*n+l)); } *(R+l*n+j)=ip; //i行l行的Q值； for(i=0;i<m;i++) //A[*][l]存放消除前面基底的分量； {