Next: The Discussion Up: Gauss-Jordan Elimination Previous: Gauss-Jordan Elimination

## 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
```

Zdzislaw Meglicki
2001-02-26