C =============================================================================
C MATRIX INVERSION SUBROUTINES FROM "NUMERICAL RECIPES"

      SUBROUTINE LUDCMP(A,N,INDX,D)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NMAX=1000,TINY=1.0D-20)
      DIMENSION A(N,N),INDX(N),VV(NMAX)
      D = 1.0D0
      DO 12 I = 1,N
        AAMAX = 0.0D0
        DO 11 J = 1,N
          IF (ABS(A(I,J)).GT.AAMAX) AAMAX = ABS(A(I,J))
 11     CONTINUE
        IF (AAMAX.EQ.0.0D0) PAUSE 'Singular matrix.'
        VV(I) = 1.0D0/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.0D0
        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.0D0) A(J,J) = TINY
          DUM = 1.0D0/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.0D0) A(N,N) = TINY
      RETURN
      END

C -----------------------------------------------------------------------------

      SUBROUTINE LUBKSB(A,N,INDX,B)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N,N),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.0D0) 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

C -----------------------------------------------------------------------------

      SUBROUTINE MPROVE(A,ALUD,N,INDX,B,X)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NMAX=1000)
      DIMENSION A(N,N),ALUD(N,N),INDX(N),B(N),X(N),R(NMAX)
      DO 12 I = 1,N
        SDP = -B(I)
        DO 11 J = 1,N
          SDP = SDP + A(I,J)*X(J)
 11     CONTINUE
        R(I) = SDP
 12   CONTINUE
      CALL LUBKSB(ALUD,N,INDX,R)
      DO 13 I = 1,N
        X(I) = X(I) - R(I)
 13   CONTINUE
      RETURN
      END

C =============================================================================