qiuni.txt

设计一个矩阵相乘的程序

#include <stdio.h> 
#include <stdlib.h> 
#include <malloc.h> 
#include <math.h> 
#define MAX 255 

void MatrixMul(a,b,m,n,k,c) /*实矩阵相乘*/ 
int m,n,k; /*m:矩阵A的行数, n:矩阵B的行数, k:矩阵B的列数*/ 
double a[],b[],c[]; /*a为A矩阵, b为B矩阵, c为结果，即C = AB */ 
{ 
int i,j,l,u; 
/*逐行逐列计算乘积*/ 
for (i=0; i<=m-1; i++) 
for (j=0; j<=k-1; j++) 
{ 
u=i*k+j; c[u]=0.0; 
for (l=0; l<=n-1; l++) 
c[u]=c[u]+a[i*n+l]*b[l*k+j]; 
} 
return; 
} 
int brinv(a,n) /*求矩阵的逆矩阵*/ 
int n; /*矩阵的阶数*/ 
double a[]; /*矩阵C*/ 
{ 
int *is,*js,i,j,k,l,u,v; 
double d,p; 
is=malloc(n*sizeof(int)); 
js=malloc(n*sizeof(int)); 
for (k=0; k<=n-1; k++) 
{ 
d=0.0; 
for (i=k; i<=n-1; i++) 
/*全选主元，即选取绝对值最大的元素*/ 
for (j=k; j<=n-1; j++) 
{ 
l=i*n+j; p=fabs(a[l]); 
if (p>d) { d=p; is[k]=i; js[k]=j;} 
} 
/*全部为0，此时为奇异矩阵*/ 
if (d+1.0==1.0) 
{ 
free(is); free(js); printf(" >> This is a singular matrix, can't be inversed!\n"); 
return(0); 
} 
/*行交换*/ 
if (is[k]!=k) 
for (j=0; j<=n-1; j++) 
{ 
u=k*n+j; v=is[k]*n+j; 
p=a[u]; a[u]=a[v]; a[v]=p; 
} 
/*列交换*/ 
if (js[k]!=k) 
for (i=0; i<=n-1; i++) 
{ 
u=i*n+k; v=i*n+js[k]; 
p=a[u]; a[u]=a[v]; a[v]=p; 
} 
l=k*n+k; 
a[l]=1.0/a[l]; /*求主元的倒数*/ 
/* a[kj]a[kk] -> a[kj] */ 
for (j=0; j<=n-1; j++) 
if (j!=k) 
{ 
u=k*n+j; a[u]=a[u]*a[l]; 
} 
/* a[ij] - a[ik]a[kj] -> a[ij] */ 
for (i=0; i<=n-1; i++) 
if (i!=k) 
for (j=0; j<=n-1; j++) 
if (j!=k) 
{ 
u=i*n+j; 
a[u]=a[u]-a[i*n+k]*a[k*n+j]; 
} 
/* -a[ik]a[kk] -> a[ik] */ 
for (i=0; i<=n-1; i++) 
if (i!=k) 
{ 
u=i*n+k; a[u]=-a[u]*a[l]; 
} 
} 
for (k=n-1; k>=0; k--) 
{ 
/*恢复列*/ 
if (js[k]!=k) 
for (j=0; j<=n-1; j++) 
{ 
u=k*n+j; v=js[k]*n+j; 
p=a[u]; a[u]=a[v]; a[v]=p; 
} 
/*恢复行*/ 
if (is[k]!=k) 
for (i=0; i<=n-1; i++) 
{ 
u=i*n+k; v=i*n+is[k]; 
p=a[u]; a[u]=a[v]; a[v]=p; 
} 
} 
free(is); free(js); 
return(1); 
} 

print_matrix(a,n)/*打印方阵a的元素*/ 
int n; /*矩阵的阶数*/ 
double a[]; /*矩阵C或D*/ 
{ 
int i,j; 
for (i=0; i<n; i++) 
{ 
for (j=0; j<n; j++) 
printf("%13.7f\t",a[i*n+j]); 
printf("\n"); 
} 
} 

int main() 
{ 
int i,j,k,n=0,m; 
double A[MAX],B[MAX],C[MAX],D[MAX],E[MAX]; /* 题中没有说明，这里做实数看待 */ 
system("cls"); 
puts("**********************************************************"); 
puts("* This program is to inverse a square matrix(nxn). *"); 
puts("**********************************************************"); 
while(n<=0) 
{ /* 设A是n*m矩阵则B必为m*n矩阵，这样AB才会是方阵存在逆矩阵 */ 
printf(" >> Please input the number of rows in A, n= "); 
scanf("%d",&n); 
printf(" >> Please input the number of cols in A, m= "); 
scanf("%d",&m); 
} 
printf(" >> Please input the elements of the matrix A one by one:\n >> "); 
for(i=0;i<n*m;i++) 
scanf("%lf",&A[i]); 

printf(" >> Please input the elements of the matrix B one by one:\n >> "); 
for(i=0;i<m*n;i++) 
scanf("%lf",&B[i]); 

MatrixMul(A,B,n,m,n,C);/* A*B的结果存在C里面 */ 
for(i=0;i<n*n;i++) /* D[n*n]是C[n*n](也就是AB)的副本，当C[n*n]变换为本身的逆矩阵后再和D[n*n]做乘积运算 */ 
D[i]=C[i]; /* 如果乘积结果为单位矩阵，则说明计算结果正确 */ 
i=brinv(C,n); 
if (i!=0) 
{ 
printf(" Matrix AB:\n"); 
print_matrix(D,n); 
printf("\n"); 
printf(" AB's Inverse Matrix AB-:\n"); 
print_matrix(C,n); 
printf("\n"); 
printf(" Product of AB and AB- :\n"); 
MatrixMul(D,C,n,n,n,E); 
print_matrix(E,n); 
} 
system("pause"); 
return 0; 
}