SUBROUTINE gaussj(a, n, b, m)

!  Gauss-Jordan solution of linear equations in quadruple precision.
!  Adapted from the routine GAUSSJ in the Numerical Recipes package.
!  A is an NxN quadruple-precision matrix, declared as quadruple-precision
!  A(:,:).
!  B contains M right-hand side vectors, where M = 1 in many cases,
!  and is declared as quadruple-precision B(:,:).
!  On output, A has been overwritten by its inverse, and B has been
!  replaced with the solutions.

! Programmer : Alan Miller, CSIRO Mathematical Information Sciences
! e-mail: alan @ mel.dms.csiro.au   WWW: http://www.mel.dms.csiro.au/~alan
! Fax: (+61) 3-9545-8080

! Latest revision - 6 December 1996

USE quadruple_precision
IMPLICIT NONE

INTEGER, INTENT(IN)                         :: n, m
TYPE (quad), INTENT(IN OUT), DIMENSION(:,:) :: a, b

!       Local variables

INTEGER, PARAMETER :: nmax = 50
INTEGER            :: ipiv(nmax), indxr(nmax), indxc(nmax), j, i, irow,  &
                      icol, k, l, ll
TYPE (quad)        :: pivinv, dum, one, zero
REAL (dp)          :: big

one%hi = 1._dp
one%lo = 0._dp
zero%hi = 0._dp
zero%lo = 0._dp

ipiv(1:n) = 0

!       Start of the main loop

DO i = 1, n
  big = 0.
  DO j = 1, n
    IF(ipiv(j) /= 1) THEN
      DO k = 1, n
        IF (ipiv(k) == 0) THEN
          IF (ABS(a(j,k)%hi) >= big) THEN
            big = ABS(a(j,k)%hi)
            irow = j
            icol = k
          END IF
        ELSE IF (ipiv(k) > 1) THEN
          WRITE(*, *) 'Singular matrix'
          RETURN
        END IF
      END DO
    END IF
  END DO
  ipiv(icol) = ipiv(icol)+1

!       Pivot has been selected.   Interchange rows, if necessary, to
!       put the element on the diagonal.

  IF (irow /= icol) THEN
    DO l = 1, n
      dum = a(irow,l)
      a(irow,l) = a(icol,l)
      a(icol,l) = dum
    END DO
    DO l = 1, m
      dum = b(irow,l)
      b(irow,l) = b(icol,l)
      b(icol,l) = dum
    END DO
  END IF
  indxr(i) = irow
  indxc(i) = icol

!       Now do the pivoting

  IF (a(icol,icol)%hi == 0.) THEN
    WRITE(*, *) 'Singular matrix.'
    RETURN
  END IF
  pivinv = one / a(icol,icol)
  a(icol,icol) = one
  DO l = 1, n
    a(icol, l) = a(icol, l) * pivinv
  END DO
  DO l = 1, m
    b(icol, l) = b(icol, l) * pivinv
  END DO
  DO ll = 1, n
    IF(ll /= icol)THEN
    dum = a(ll,icol)
    a(ll,icol) = zero
    DO l = 1, n
      a(ll, l) = a(ll, l) - a(icol, l) * dum
    END DO
    DO l = 1, m
      b(ll, l) = b(ll, l) - b(icol, l) * dum
    END DO
    END IF
  END DO
END DO

!       Finally unscramble A by re-arranging the columns

DO l = n, 1, -1
  IF(indxr(l) /= indxc(l))THEN
    DO k = 1, n
      dum = a(k,indxr(l))
      a(k,indxr(l)) = a(k,indxc(l))
      a(k,indxc(l)) = dum
    END DO
  END IF
END DO

RETURN
END SUBROUTINE gaussj