@(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]]);
}
}