高斯—约当消去法求解矩阵的逆和解向量,求解的逆矩阵保存在A中,解向量保存在B中,其中AX=B,不知为什么老是出错,函数传递有问题,特向高手求救!!谢谢了
#include <math.h>
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#define NR_END 1
#define FREE_ARG char*
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
int *ivector(long nl, long nh)
{ int *v;
v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if(!v) printf("allocation failure in ivector()");
return v-nl+NR_END;
}
void free_ivector(int *v, long nl, long nh)
{ free((FREE_ARG) (v+nl-NR_END));
}
void gaussj(float **a, int n, float **b, int m)
{
int *indxc,*indxr,*ipiv;
int i,icol,irow,j,k,l,ll;
float big,dum,pivinv,temp;
for(j=0;j<n;j++)
{for(i=0;i<m;i++)
printf("a[%d][%d]=%f",j,i,b[j][i]);
printf("\n");}
indxc=ivector(1,n);
indxr=ivector(1,n);
ipiv=ivector(1,n);
for(j=0;j<n;j++) ipiv[j]=0;
for(i=0;i<n;i++) {
big=0;
for(j=0;j<n;j++)
if(ipiv[j]!=1)
for(k=0;k<n;k++) {
if(ipiv[k]==0) {
if(fabs(a[j][k])>=big) {
big=fabs(a[j][k]);
irow=j;
icol=k;
}
}
}
++(ipiv[icol]);
if(irow!=icol) {
for(l=0;l<n;l++) SWAP(a[irow][l],a[icol][l])
for(l=0;l<m;l++) SWAP(b[irow][l],b[icol][l])
}
indxr[i]=irow;
indxc[i]=icol;
if(a[icol][icol]==0.0) printf("gaussj: Singular Matrix");
pivinv=1.0/a[icol][icol];
a[icol][icol]=1.0;
for(l=0;l<n;l++) a[icol][l] *=pivinv;
for(l=0;l<m;l++) b[icol][l] *=pivinv;
for(ll=0;ll<n;ll++)
if(ll!=icol) {
dum=a[ll][icol];
a[ll][icol]=0.0;
for(l=0;l<n;l++) a[ll][l] -=a[icol][l]*dum;
for(l=0;l<m;l++) b[ll][l] -=b[icol][l]*dum;
}
}
for(l=n-1;l>=0;l--) {
if(indxr[l]!=indxc[l])
for(k=0;k<n;k++)
SWAP(a[k][indxr[l]],a[k][indxc[l]]);
}
free_ivector(ipiv,1,n);
free_ivector(indxr,1,n);
free_ivector(indxc,1,n);
}
main()
{ int i;
static float a[4][4]={{1.0,3.0,2.0,13.0},{7.0,2.0,1.0,-2.0},{9.0,15.0,3.0,-2.0},{-2.0,-2.0,11.0,5.0}};
static float b[4][2]={{9.0,0.0},{6.0,4.0},{11.0,7.0},{-2.0,-1.0}};
gaussj(a,4,b,2);
for(i=0;i<=3;i++)
printf("x(%d)=%13.7e, %13.7e\n", i, b[i][0], b[i][1]);
}