The Discussion

The subroutine is implemented in the form of a large DO loop that moves down the diagonal while performing the operations of the Gauss Jordan elimination on the resulting submatrix, i.e., on all that between the current point and the bottom right corner of the matrix.

The first step is to locate the pivot element. This is done by calling Fortran intrinsic MAXLOC:

     irc = MAXLOC(ABS(a), outerand(lpiv, lpiv))

The location of the pivot is placed in the array irc, whose entries are pointed to by irow and icol.

Next we swap rows pointed to by irow and icol, so that the pivot ends up on the diagonal in the (icol, icol) location. The original placement of the pivot term is memorised on indxr and indxc for unscrambling at the end of the procedure:

     indxr(i) = irow     ! We are now ready to divid the pivot row by the
     indxc(i) = icol     ! pivot element, located at irow and icol.

If that term, i.e., the pivot term happens to be 0 then the matrix, obviously, is singular, so we must abort raising an error flag:

     IF (a(icol, icol) == 0.0) &
          CALL nrerror('gaussj: singular matrix(2)')

Now we perform the first division of the Gauss-Jordan elimination procedure, i.e., we devide by A_kk:

     pivinv=1.0_sp/a(icol, icol)
     a(icol, icol)=1.0
     a(icol, :)=a(icol, :)*pivinv
     b(icol, :)=b(icol, :)*pivinv

and then reduce all other rows of the matrix. Observe that the whole column a(:, icol) is memorized first on dumc. We will need it, because, in the meantime we're setting that column of A to

$$\left( \begin{array}{c} 0 \\ 0 \\ \vdots \\ 1/A_{kl} \\ \vdots \\ 0 \end{array} \right)$$

where 1/A_kl is the inverse of the pivot.

     dumc=a(:, icol)     ! Next we reduce the rows, except for the pivot one.
     a(:, icol)=0.0
     a(icol, icol)=pivinv

Why do we do that? That's because we no longer need that column in its previous form, and instead we're now setting it to what the corresponding column of an identity matrix would have become if it was subject to the Gauss-Jordan operations performed so far. This way, as we keep using A we're simultaneously building its inverse in the same space!

Finally we perform the reductions on b and at the same time we build the inverse on a thusly:

     a(1:icol-1,:)=a(1:icol-1,:) - outerprod(dumc(1:icol-1), a(icol,:))
     b(1:icol-1,:)=b(1:icol-1,:) - outerprod(dumc(1:icol-1), b(icol,:))
     a(icol+1:,:)=a(icol+1:,:) - outerprod(dumc(icol+1:), a(icol,:))
     b(icol+1:,:)=b(icol+1:,:) - outerprod(dumc(icol+1:), b(icol,:))

The pivoting would have scrambled the matrix, which now needs to be set back in the right order. But we have memorized the sequence of scrambling on indxr and indxc, so now, we can undo it:

  DO l = n, 1, -1
     CALL swap(a(:, indxr(l)), a(:, indxc(l)))
  END DO

The safety check right at the beginning:

     ipiv(icol) = ipiv(icol) + 1
     IF (ipiv(icol) > 1) CALL nrerror('gaussj: singular matrix(1)')

ensures that we do not hit the same pivot row twice. If that is to happen the matrix must be singular and in that case we abort flagging an error message.