SUBROUTINE MINV(VARX,VINV,NGDO,DET,LINDX)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C MINV CALCULATES THE INVERSE OF THE MATRIX VARX AND PLACES IT IN THE
C  MATRIX VINV. IT ALSO CALCULATES THE DETERMINANT, DET, OF VARX.
C
C  SOFTWARE WRITTEN BY Benjamin T. Marshall.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C      (c) Copyright 1987-1999 by GATS, Inc.
C          11864 Canon Blvd., Suite 101, Newport News, VA 23606
C          Phone: (757) 873-5920
C
C  All Rights Reserved. No part of this software or publication may be
C  reproduced, stored in a retrieval system, or transmitted, in any form
C  or by any means, electronic, mechanical, photocopying, recording, or
C  otherwise without the prior written permission of GATS, Inc.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT NONE
      INTEGER NGDO
      LOGICAL FLAG
      REAL*8 VARX(NGDO,NGDO),VINV(NGDO,NGDO),DET
      INTEGER LINDX(NGDO),IG,IG2
      CALL LUDCMP(VARX,NGDO,NGDO,LINDX,DET,VINV)
      DO IG=1,NGDO
        DO IG2=1,NGDO
          VINV(IG,IG2)=0.
        ENDDO
        VINV(IG,IG)=1.
      ENDDO
      DO IG=1,NGDO
        DET=DET*VARX(IG,IG)
        CALL LUBKSB(VARX,NGDO,NGDO,LINDX,VINV(1,IG))
      ENDDO
      RETURN
      END
      SUBROUTINE LUDCMP(A,N,NP,INDX,D,VV)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C LUDCMP REPLACES THE MATRIX A WITH ITS LU DECOMPOSITION. INDX IS AN
C  OUTPUT VECTOR WHICH RECORDS THE ROW PERMUTATION EFFECTED BY THE
C  PARTIAL PIVOTING, D IS OUTPUT AS +OR- 1 DEPENDING ON WHETHER THE
C  NUMBER OF ROW INTERCHANGES WAS EVEN OR ODD. MATRIX A IS NPxNP IN SIZE
C  OF WHICH NxN ELEMENTS ARE USED. THE DECOMPOSED A AND THE OUTPUT INDX
C  CAN THEN BE USED WITH LUBKSB TO SOLVE LINEAR EQUATIONS OR TO INVERT A
C  MATRIX.
C
C     THIS ROUTINE IS DERIVED FROM
C     "NUMERICAL RECIPES THE ART OF SCIENTIFIC COMPUTING" by
C     W.H. PRESS, B.P. FLANNERY, S.A. TEUKOLSKY AND W.T. VETTERLING;
C     CAMBRIDGE UNIVERSITY PRESS, 510 NORTH AVE. NEW ROCHELLE, NY 10801.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (TINY=1.0D-99)
      DIMENSION A(NP,NP),INDX(N),VV(N)
      D=1.
      DO 12 I=1,N
        AAMAX=0.
        DO 11 J=1,N
          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11      CONTINUE
        IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
        VV(I)=1./AAMAX
12    CONTINUE
      DO 19 J=1,N
        IF (J.GT.1) THEN
          DO 14 I=1,J-1
            SUM=A(I,J)
            IF (I.GT.1)THEN
              DO 13 K=1,I-1
                SUM=SUM-A(I,K)*A(K,J)
13            CONTINUE
              A(I,J)=SUM
            ENDIF
14        CONTINUE
        ENDIF
        AAMAX=0.
        DO 16 I=J,N
          SUM=A(I,J)
          IF (J.GT.1)THEN
            DO 15 K=1,J-1
              SUM=SUM-A(I,K)*A(K,J)
15          CONTINUE
            A(I,J)=SUM
          ENDIF
          DUM=VV(I)*ABS(SUM)
          IF (DUM.GE.AAMAX) THEN
            IMAX=I
            AAMAX=DUM
          ENDIF
16      CONTINUE
        IF (J.NE.IMAX)THEN
          DO 17 K=1,N
            DUM=A(IMAX,K)
            A(IMAX,K)=A(J,K)
            A(J,K)=DUM
17        CONTINUE
          D=-D
          VV(IMAX)=VV(J)
        ENDIF
        INDX(J)=IMAX
        IF(J.NE.N)THEN
          IF(A(J,J).EQ.0.)A(J,J)=TINY
          DUM=1./A(J,J)
          DO 18 I=J+1,N
            A(I,J)=A(I,J)*DUM
18        CONTINUE
        ENDIF
19    CONTINUE
      IF(A(N,N).EQ.0.)A(N,N)=TINY
      RETURN
      END
      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C LUBKSB SOLVES A SET OF LINEAR EQUATIONS C*X=B, WHERE THE INPUT A IS
C  THE LU DECOMPOSITION OF C AS DETERMINED BY LUDCMP. INDX IS INPUT AS 
C  THE PERMUTATION VECTOR RETURNED BY LUDCMP.
C
C     THIS ROUTINE IS DERIVED FROM
C     "NUMERICAL RECIPES THE ART OF SCIENTIFIC COMPUTING" by
C     W.H. PRESS, B.P. FLANNERY, S.A. TEUKOLSKY AND W.T. VETTERLING;
C     CAMBRIDGE UNIVERSITY PRESS, 510 NORTH AVE. NEW ROCHELLE, NY 10801.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(NP,NP),INDX(N),B(N)
      II=0
      DO 12 I=1,N
        LL=INDX(I)
        SUM=B(LL)
        B(LL)=B(I)
        IF (II.NE.0)THEN
          DO 11 J=II,I-1
            SUM=SUM-A(I,J)*B(J)
11        CONTINUE
        ELSE IF (SUM.NE.0.) THEN
          II=I
        ENDIF
        B(I)=SUM
12    CONTINUE
      DO 14 I=N,1,-1
        SUM=B(I)
        IF(I.LT.N)THEN
          DO 13 J=I+1,N
            SUM=SUM-A(I,J)*B(J)
13        CONTINUE
        ENDIF
        B(I)=SUM/A(I,I)
14    CONTINUE
      RETURN
      END