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

pivinv=1.0_sp/a(icol, icol) a(icol, icol)=1.0 a(icol, :)=a(icol, :)*pivinv b(icol, :)=b(icol, :)*pivinvand 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
where 1/

dumc=a(:, icol) ! Next we reduce the rows, except for the pivot one. a(:, icol)=0.0 a(icol, icol)=pivinvWhy 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

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.