C语言实现矩阵的LU分解、施密特正交化、Givens分解、Householder分解

C语言实现矩阵的LU分解、施密特正交化、Givens分解、Householder分解

By Kim.Wang，UCAS

#include

#include

#include

#define HS 10

#define LS 10

int 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;

else

q[i][j]=0;

}

for(j=0;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;

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 < 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;

else

dw[i][j]=q[i][j]=0;

}

for(k=0;k

{

for (i = 0; i < m; i++)

{

if(i

u[i]=0;

else if(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 < 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

{

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("\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;

else

p[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

{

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;

else

q[i][j]=0;

}

for(j=0;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;

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 < 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];

else

jdz=r[n-1][n-1];

if(m

{

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);

}

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)

{

case 0:return 0; break;

case 1:schmidt(); break;

case 2:householder(); break;

case 3:givens(); break;

}

}

while(zl!=0);

}