C**********************************************************************
C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/exinv.for
C  ROUTINE:   EXINV1 FORTRAN
C  LANGUAGE LEVEL: Fortran IV or Fortran 77.
C  PURPOSE:   To Illustrate robust Matrix decomposition and
C             Solution of AX=B or Inverse if B=vectors of I where
C             A=Coefficient Matrix, B is RHS, X=Solution Vector,
C             I=An Identity Matrix Vector. Computes Matrix Condition
C             and uses Machine tolerance to return reliable results
C             or an error (in most cases). Machine tolerance isn't
C             actually needed since its definition is the smallest
C             floating-point number Eps such that on that machine:
C             Eps+1.0 .LT. 1.0  and  Eps-1.0 .LT. 1.0  (see line 50
C             of this file where that's done machine independently).
C  PACKAGE:   MOLER
C  CALLS:     DECOMP, SOLVE
C  CALLED BY: on PSUVM issue: WATFOR77 EXINV1
C               or on the PC: WATFOR77 EXINV1.FOR
C  AUTHORS:  George E. Forsythe and Cleve B. Moler.
C  TEXT:     COMPUTER SOLUTION OF LINEAR ALGEBRAIC SYSTEMS
C   Englewood Cliffs, N.J., Prentice-Hall, "1967". xi, 148 p.
C   Series: Prentice-Hall series in automatic computation.
C   Bibliography: p. 138-141.
C   1. Numerical analysis -- Data processing. 2. Matrices.
C   PSU Pattee Library Call#: QA297.F57
C**********************************************************************
      DOUBLE PRECISION A(10,10), B(10), WORK(10), COND, CONDP1
      INTEGER IPVT(10), I, J, N, NDIM
      NDIM = 10
      N = 3
      A(1,1) = 10
      A(2,1) = -3
      A(3,1) = 5
      A(1,2) = -7
      A(2,2) = 2
      A(3,2) = -1
      A(1,3) = 0
      A(2,3) = 6
      A(3,3) = 5
      WRITE(6,*) 'Coefficient matrix A...'
      DO 1 I = 1, N
         WRITE(6,2) (A(I,J), J=1,N)
    1 CONTINUE
    2 FORMAT(1X ,10F5.0)
      WRITE(6,7)
      WRITE(6,*) '** Solving AX=B ** '
      CALL DECOMP(NDIM, N, A, COND, IPVT, WORK)
      WRITE(6,3) COND
    3 FORMAT(' Condition Number =', E16.8 )
      WRITE(6,7)
      CONDP1 = COND + 1.
      IF (CONDP1 .EQ. COND) WRITE(6,4)
    4 FORMAT(' $$Matrix A is singular to Machine Precision.')
      IF (CONDP1 .EQ. COND) STOP 99
      WRITE(6,*) 'Right=Hand Side B...'
      B(1) = 7
      B(2) = 4
      B(3) = 6
      DO 5 I = 1, N
         WRITE(6,2) B(I)
    5 CONTINUE
      WRITE(6,7)
      CALL SOLVE(NDIM, N, A, B, IPVT)
      WRITE(6,*) 'Solution Vector X...'
      DO 6 I = 1, N
         WRITE(6,7) B(I)
    6 CONTINUE
    7 FORMAT(1X ,F10.5)
      STOP
      END
      SUBROUTINE DECOMP(NDIM,N,A,COND,IPVT,WORK)
C
      INTEGER NDIM,N
      DOUBLE PRECISION A(NDIM,N),COND,WORK(N)
      INTEGER IPVT(N)
C
C     DECOMPOSES A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION
C     AND ESTIMATES THE CONDITION OF THE MATRIX.
C
C     USE SOLVE TO COMPUTE SOLUTIONS TO LINEAR SYSTEMS.
C
C     INPUT..
C
C        NDIM = DECLARED ROW DIMENSION OF THE ARRAY CONTAINING  A.
C
C        N = ORDER OF THE MATRIX.
C
C        A = MATRIX TO BE TRIANGULARIZED.
C
C     OUTPUT..
C
C        A  CONTAINS AN UPPER TRIANGULAR MATRIX  U  AND A PERMUTED
C          VERSION OF A LOWER TRIANGULAR MATRIX  I-L  SO THAT
C          (PERMUTATION MATRIX)*A = L*U .
C
C        COND = AN ESTIMATE OF THE CONDITION OF  A .
C           FOR THE LINEAR SYSTEM  A*X = B, CHANGES IN  A  AND  B
C           MAY CAUSE CHANGES  COND  TIMES AS LARGE IN  X .
C           IF  COND+1.0 .EQ. COND , A IS SINGULAR TO WORKING
C           PRECISION.  COND  IS SET TO  1.0D+32  IF EXACT
C           SINGULARITY IS DETECTED.
C
C        IPVT = THE PIVOT VECTOR.
C           IPVT(K) = THE INDEX OF THE K-TH PIVOT ROW
C           IPVT(N) = (-1)**(NUMBER OF INTERCHANGES)
C
C     WORK SPACE..  THE VECTOR  WORK  MUST BE DECLARED AND INCLUDED
C                   IN THE CALL.  ITS INPUT CONTENTS ARE IGNORED.
C                   ITS OUTPUT CONTENTS ARE USUALLY UNIMPORTANT.
C
C     THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY
C        DET(A) = IPVT(N) * A(1,1) * A(2,2) * ... * A(N,N).
C
      DOUBLE PRECISION EK, T, ANORM, YNORM, ZNORM
      INTEGER NM1, I, J, K, KP1, KB, KM1, M
      DOUBLE PRECISION DABS
C
      IPVT(N) = 1
      IF (N .EQ. 1) GO TO 80
      NM1 = N - 1
C
C     COMPUTE 1-NORM OF A
C
      ANORM = 0.0D0
      DO 10 J = 1, N
         T = 0.0D0
         DO 5 I = 1, N
            T = T + DABS(A(I,J))
    5    CONTINUE
         IF (T .GT. ANORM) ANORM = T
   10 CONTINUE
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
      DO 35 K = 1,NM1
         KP1= K+1
C
C        FIND PIVOT
C
         M = K
         DO 15 I = KP1,N
            IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
   15    CONTINUE
         IPVT(K) = M
         IF (M .NE. K) IPVT(N) = -IPVT(N)
         T = A(M,K)
         A(M,K) = A(K,K)
         A(K,K) = T
C
C        SKIP STEP IF PIVOT IS ZERO
C
         IF (T .EQ. 0.0D0) GO TO 35
C
C        COMPUTE MULTIPLIERS
C
         DO 20 I = KP1,N
             A(I,K) = -A(I,K)/T
   20    CONTINUE
C
C        INTERCHANGE AND ELIMINATE BY COLUMNS
C
         DO 30 J = KP1,N
             T = A(M,J)
             A(M,J) = A(K,J)
             A(K,J) = T
             IF (T .EQ. 0.0D0) GO TO 30
             DO 25 I = KP1,N
                A(I,J) = A(I,J) + A(I,K)*T
   25        CONTINUE
   30    CONTINUE
   35 CONTINUE
C
C     COND = (1-NORM OF A)*(AN ESTIMATE OF 1-NORM OF A-INVERSE)
C     ESTIMATE OBTAINED BY ONE STEP OF INVERSE ITERATION FOR THE
C     SMALL SINGULAR VECTOR.  THIS INVOLVES SOLVING TWO SYSTEMS
C     OF EQUATIONS, (A-TRANSPOSE)*Y = E  AND  A*Z = Y  WHERE  E
C     IS A VECTOR OF +1 OR -1 CHOSEN TO CAUSE GROWTH IN Y.
C     ESTIMATE = (1-NORM OF Z)/(1-NORM OF Y)
C
C     SOLVE (A-TRANSPOSE)*Y = E
C
      DO 50 K = 1, N
         T = 0.0D0
         IF (K .EQ. 1) GO TO 45
         KM1 = K-1
         DO 40 I = 1, KM1
            T = T + A(I,K)*WORK(I)
   40    CONTINUE
   45    EK = 1.0D0
         IF (T .LT. 0.0D0) EK = -1.0D0
         IF (A(K,K) .EQ. 0.0D0) GO TO 90
         WORK(K) = -(EK + T)/A(K,K)
   50 CONTINUE
      DO 60 KB = 1, NM1
         K = N - KB
         T = 0.0D0
         KP1 = K+1
         DO 55 I = KP1, N
            T = T + A(I,K)*WORK(K)
   55    CONTINUE
         WORK(K) = T
         M = IPVT(K)
         IF (M .EQ. K) GO TO 60
         T = WORK(M)
         WORK(M) = WORK(K)
         WORK(K) = T
   60 CONTINUE
C
      YNORM = 0.0D0
      DO 65 I = 1, N
         YNORM = YNORM + DABS(WORK(I))
   65 CONTINUE
C
C     SOLVE A*Z = Y
C
      CALL SOLVE(NDIM, N, A, WORK, IPVT)
C
      ZNORM = 0.0D0
      DO 70 I = 1, N
         ZNORM = ZNORM + DABS(WORK(I))
   70 CONTINUE
C
C     ESTIMATE CONDITION
C
      COND = ANORM*ZNORM/YNORM
      IF (COND .LT. 1.0D0) COND = 1.0D0
      RETURN
C
C     1-BY-1
C
   80 COND = 1.0D0
      IF (A(1,1) .NE. 0.0D0) RETURN
C
C     EXACT SINGULARITY
C
   90 COND = 1.0D+32
      RETURN
      END
      SUBROUTINE SOLVE(NDIM, N, A, B, IPVT)
C
      INTEGER NDIM, N, IPVT(N)
      DOUBLE PRECISION A(NDIM,N),B(N)
C
C   SOLUTION OF LINEAR SYSTEM, A*X = B .
C   DO NOT USE IF DECOMP HAS DETECTED SINGULARITY.
C
C   INPUT..
C
C     NDIM = DECLARED ROW DIMENSION OF ARRAY CONTAINING A .
C
C     N = ORDER OF MATRIX.
C
C     A = TRIANGULARIZED MATRIX OBTAINED FROM DECOMP .
C
C     B = RIGHT HAND SIDE VECTOR.
C
C     IPVT = PIVOT VECTOR OBTAINED FROM DECOMP .
C
C   OUTPUT..
C
C     B = SOLUTION VECTOR, X .
C
      INTEGER KB, KM1, NM1, KP1, I, K, M
      DOUBLE PRECISION T
C
C     FORWARD ELIMINATION
C
      IF (N .EQ. 1) GO TO 50
      NM1 = N-1
      DO 20 K = 1, NM1
         KP1 = K+1
         M = IPVT(K)
         T = B(M)
         B(M) = B(K)
         B(K) = T
         DO 10 I = KP1, N
             B(I) = B(I) + A(I,K)*T
   10    CONTINUE
   20 CONTINUE
C
C     BACK SUBSTITUTION
C
      DO 40 KB = 1,NM1
         KM1 = N-KB
         K = KM1+1
         B(K) = B(K)/A(K,K)
         T = -B(K)
         DO 30 I = 1, KM1
             B(I) = B(I) + A(I,K)*T
   30    CONTINUE
   40 CONTINUE
   50 B(1) = B(1)/A(1,1)
      RETURN
      END