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