@(Cabinet)[cpp_num, aca_book, published_gitbook]

date: 2016-05-24

Numerical Recipes, Ch2

Digest of Numerical Recipes

Gauss-Jordan Elimination

void gaussj(MatDoub_IO &a, MatDoub_IO &b)
{
    Int i,icol,irow,j,k,l,ll,n=a.nrows(),m=b.ncols();
    Doub big,dum,pivinv;
    //bookkeeping for pivoting
    //indxc and indxr store pivot positions for every column
    VecInt indxc(n),indxr(n),ipiv(n);
    //initialize ipiv[]
    for (j=0;j<n;j++) ipiv[j]=0;
    //main loop over the columns to be reduced
    for (i=0;i<n;i++) {
        big=0.0;
        //outer loop searching for a pivot element
        for (j=0;j<n;j++)
            //pass, only if ipiv[j] is a pivot, i.e. ==1
            if (ipiv[j] != 1)
            //search over all columns, find the biggest one
            //as the pivot
                for (k=0;k<n;k++) {
                //check if ipiv[k] has been reduced before, per column
                    if (ipiv[k] == 0) {
                        if (abs(a[j][k]) >= big) {
                            big=abs(a[j][k]);
                            irow=j;
                            icol=k;
                        }
                    }
                }
        //found the pivot
        //set the flag=1
        ++(ipiv[icol]);
        //interchange rows
        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) throw("gaussj: Singular Matrix");
        pivinv=1.0/a[icol][icol];
        a[icol][icol]=1.0;
        //divide the pivot row by the pivot element
        for (l=0;l<n;l++) a[icol][l] *= pivinv;
        for (l=0;l<m;l++) b[icol][l] *= pivinv;
        //reduce the rows, except for the pivot one
        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;
            }
    }
    //unscramble the solution via interchanging columns
    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]]);
    }
}

results matching ""

    No results matching ""