2-1 Gaussj

#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