C     $Id: minv.f 3304 2011-01-17 15:25:59Z brideout $
C
      SUBROUTINE MINV(A,N,D,L,M)
C
C     Inverts general matrix A (overwrites A) using standard
C     gauss-jordan method. The determinant is also calculated. a
C     determinant of zero indicates that the matrix is singular.
C
C      Input:
C        n - order of matrix a
C        l - work vector of length n
C        m - work vector of length n
C
C      Input, Output:
C        a - input matrix, destroyed in computation and replaced by
C            resultant inverse.
C
C      Output:
C        d - resultant determinant
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION D
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(*)
      INTEGER L(*),M(*)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION BIGA,HOLD
      INTEGER I,IJ,IK,IZ,J,JI,JK,JP,JQ,JR,K,KI,KJ,KK,NK
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS
C     ..
C     search for largest element
      D = 1.0D0
      NK = -N
      DO 90 K = 1,N
         NK = NK + N
         L(K) = K
         M(K) = K
         KK = NK + K
         BIGA = A(KK)
         DO 20 J = K,N
            IZ = N*(J-1)
            DO 10 I = K,N
               IJ = IZ + I
               IF (ABS(BIGA).LT.ABS(A(IJ))) THEN
                  BIGA = A(IJ)
                  L(K) = I
                  M(K) = J
               END IF
   10       CONTINUE
   20    CONTINUE
C
C        interchange rows
C
         J = L(K)
         IF (J.GT.K) THEN
            KI = K - N
            DO 30 I = 1,N
               KI = KI + N
               HOLD = -A(KI)
               JI = KI - K + J
               A(KI) = A(JI)
               A(JI) = HOLD
   30       CONTINUE
         END IF
C
C        interchange columns
C
         I = M(K)
         IF (I.GT.K) THEN
            JP = N*(I-1)
            DO 40 J = 1,N
               JK = NK + J
               JI = JP + J
               HOLD = -A(JK)
               A(JK) = A(JI)
               A(JI) = HOLD
   40       CONTINUE
         END IF
C
C        divide column by minus pivot (value of pivot element is
C        contained in biga)
C
         IF (BIGA.EQ.0) GO TO 130
         DO 50 I = 1,N
            IF (I.NE.K) THEN
               IK = NK + I
               A(IK) = A(IK)/(-BIGA)
            END IF
   50    CONTINUE
C
C        reduce matrix
C
         DO 70 I = 1,N
            IK = NK + I
            HOLD = A(IK)
            IJ = I - N
            DO 60 J = 1,N
               IJ = IJ + N
               IF (I.NE.K) THEN
                  IF (J.NE.K) THEN
                     KJ = IJ - I + K
                     A(IJ) = HOLD*A(KJ) + A(IJ)
                  END IF
               END IF
   60       CONTINUE
   70    CONTINUE
C
C        divide row by pivot
C
         KJ = K - N
         DO 80 J = 1,N
            KJ = KJ + N
            IF (J.NE.K) A(KJ) = A(KJ)/BIGA
   80    CONTINUE
C
C        product of pivots
C
         D = D*BIGA
C
C        replace pivot by reciprocal
C
         A(KK) = 1.0D0/BIGA
   90 CONTINUE
C
C        final row and column interchange
C
      K = N
  100 CONTINUE
      K = (K-1)
      IF (K.LE.0) GO TO 140
      I = L(K)
      IF (I.GT.K) THEN
         JQ = N*(K-1)
         JR = N*(I-1)
         DO 110 J = 1,N
            JK = JQ + J
            HOLD = A(JK)
            JI = JR + J
            A(JK) = -A(JI)
            A(JI) = HOLD
  110    CONTINUE
      END IF
      J = M(K)
      IF (J.GT.K) THEN
         KI = K - N
         DO 120 I = 1,N
            KI = KI + N
            HOLD = A(KI)
            JI = KI - K + J
            A(KI) = -A(JI)
            A(JI) = HOLD
  120    CONTINUE
      END IF
      GO TO 100
  130 CONTINUE
      D = 0.0D0
  140 CONTINUE
      RETURN
      END