#include <math.h> #include "..CommonFiles rutil.h" #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} /* 完全主元法--高斯消去法*/ /* gdb 查看数组*(a+1)[2]@4*/ 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; indxc=ivector(1,n); /* 主元所在的列 */ indxr=ivector(1,n); /* 主元所在的行 */ ipiv=ivector(1,n); /* ipiv用于记录主元*/ for (j=1; j<=n; j++) ipiv[j]=0; for (i=1; i<=n; i++) { // 约化列的主循环 big=0.0; for (j=1; j<=n; j++) // 用于寻求主元素的外层循环 if (ipiv[j] != 1) // 判断是否被选过主元 for (k=1; k<=n; k++) { if (ipiv[k] == 0) { if (fabs(a[j][k]) >= big) { big=fabs(a[j][k]); irow=j; // 寻求主元核心成果就是求出irow和icol icol=k; } } else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1"); } ++(ipiv[icol]); /* ** 至此已经求得主元素,如果需要的话,进行行交换把主元素放到对角线位置上. ** 也就是把主元所在的行位置换到与主元所在列相等的行位置上. */ if (irow != icol) { for (l=1; l<=n; l++) SWAP(a[irow][l],a[icol][l]) for (l=1; l<=m; l++) SWAP(b[irow][l],b[icol][l]) } /* ** indxc[i] 为第一个主元所在的列 ** indxr[r] 为第一个主元所在的行 */ indxr[i]=irow; indxc[i]=icol; if (a[icol][icol] == 0.0) nrerror("gaussj: 奇异矩阵Singular Matrix-2"); pivinv=1.0/a[icol][icol]; /* 对第一个主元求倒数 */ a[icol][icol]=1.0; for (l=1; l<=n; l++) a[icol][l] *= pivinv; for (l=1; l<=m; l++) b[icol][l] *= pivinv; for (ll=1; ll<=n; ll++) if (ll != icol) { dum=a[ll][icol]; a[ll][icol]=0.0; for (l=1; l<=n; l++) a[ll][l] -= a[icol][l]*dum; for (l=1; l<=m; l++) b[ll][l] -= b[icol][l]*dum; } } /* 这是列约化的结尾为了使解向量保持原来的顺序,在根据交换的相反顺序交换个列*/ for (l=n; l>=1; l--) { if (indxr[l] != indxc[l])/* 第一个主元所在的行和列不相等*/ for (k=1; 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); } #undef SWAP