Example Program

Here is an example program:

SUBROUTINE gaussj(a, b)

  USE nrtype; USE nrutil, ONLY : assert_eq, nrerror, outerand, outerprod, swap
  IMPLICIT NONE
  REAL(sp), DIMENSION(:,:), INTENT(inout) :: a, b

  ! Linear equation solution by Gauss-Jordan elimination.
  ! a is an NxN input coefficient matrix. b is an NxM input matrix
  ! containing M right-hand-side vectors. On output, a is replaced
  ! by its matrix inverse, and b is replaced by the corresponding
  ! set of solution vectors.

  INTEGER(i4b), DIMENSION(SIZE(a, 1)) :: ipiv, indxr, indxc

  ! These arrays are used for bookkeeping on the pivoting

  LOGICAL(lgt), DIMENSION(SIZE(a, 1)) :: lpiv
  REAL(sp) :: pivinv
  REAL(sp), DIMENSION(SIZE(a, 1)) :: dumc
  INTEGER(i4b), TARGET :: irc(2)
  INTEGER(i4b) :: i, l, n
  INTEGER(i4b), POINTER :: irow, icol

  n = assert_eq(SIZE(a, 1), SIZE(a, 2), SIZE(b, 1), 'gaussj')
  irow => irc(1)
  icol => irc(2)
  ipiv = 0
  DO i = 1, n
     lpiv = (ipiv == 0)
     irc = MAXLOC(ABS(a), outerand(lpiv, lpiv))
     ipiv(icol) = ipiv(icol) + 1
     IF (ipiv(icol) > 1) CALL nrerror('gaussj: singular matrix(1)')

     ! We now have the pivot element, so we interchange rows, if needed,
     ! to put the pivot element on the diagonal. The columns are not 
     ! physically interchanged, only relabeled: indxc(i), the column of
     ! the ith pivot element, is the ith column that is reduced, while
     ! indxr(i) is the row in which that pivot element was originally 
     ! located. If indxr(i) \= indxc(i) there is an implied column
     ! interchange. With this form of bookkeeping, the solution b's
     ! will end up in the correct order, and the inverse matrix will be
     ! scrambled by columns.

     IF (irow /= icol) THEN
        CALL swap(a(irow, :), a(icol, :))
        CALL swap(b(irow, :), b(icol, :))
     END IF
     indxr(i) = irow     ! We are now ready to divide the pivot row by the
     indxc(i) = icol     ! pivot element, located at irow and icol.
     IF (a(icol, icol) == 0.0) &
          CALL nrerror('gaussj: singular matrix(2)')
     pivinv=1.0_sp/a(icol, icol)
     a(icol, icol)=1.0
     a(icol, :)=a(icol, :)*pivinv
     b(icol, :)=b(icol, :)*pivinv
     dumc=a(:, icol)     ! Next we reduce the rows, except for the pivot one.
     a(:, icol)=0.0
     a(icol, icol)=pivinv
     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,:))
  END DO

  ! It only remains to unscramble the solution in view of the column
  ! interchanges. We do this by interchanging pairs of columns in the
  ! revers order that the permutation was built up.

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

END SUBROUTINE gaussj