c     this is source code that simulates the old NAG routine 
c     f01abf by calling matrix inversion routine DINV supplied
c     by Peter Martinez

      SUBROUTINE DINV(N,A,IDIM,IR,IFAIL)
      REAL T1,T2,T3
      DIMENSION IR(N)
      DOUBLE PRECISION A(IDIM,N),DET,TEMP,S,
     $                 C11,C12,C13,C21,C22,C23,C31,C32,C33
      CHARACTER*6 NAME
      DATA NAME/'DINV'/,KPRNT/0/
***   PATCHED CERN LIBRARY ROUTINE F010
***   FOR ASTRONOMY, UCT, INTERNAL USE ONLY
***   SINCE WE DID NOT WRITE THIS ROUTINE WE DO NOT HAVE THE RIGHT TO 
***   GIVE IT AWAY   P.M.

C     ******************************************************************

C     REPLACES A BY ITS INVERSE.

C     (PARAMETERS AS FOR DEQINV.)

C     CALLS ... DFACT, DFINV, F010PR, ABEND.

C     ******************************************************************


C     TEST FOR PARAMETER ERRORS.

      IF((N.LT.1).OR.(N.GT.IDIM)) GO TO 7

C     TEST FOR N.LE.3.

      IF(N.GT.3) GO TO 6
      IFAIL=0
      IF(N.LT.3) GO TO 4

C     N=3 CASE.
C
C     COMPUTE COFACTORS.
      C11=A(2,2)*A(3,3)-A(2,3)*A(3,2)
      C12=A(2,3)*A(3,1)-A(2,1)*A(3,3)
      C13=A(2,1)*A(3,2)-A(2,2)*A(3,1)
      C21=A(3,2)*A(1,3)-A(3,3)*A(1,2)
      C22=A(3,3)*A(1,1)-A(3,1)*A(1,3)
      C23=A(3,1)*A(1,2)-A(3,2)*A(1,1)
      C31=A(1,2)*A(2,3)-A(1,3)*A(2,2)
      C32=A(1,3)*A(2,1)-A(1,1)*A(2,3)
      C33=A(1,1)*A(2,2)-A(1,2)*A(2,1)
      T1=ABS(SNGL(A(1,1)))
      T2=ABS(SNGL(A(2,1)))
      T3=ABS(SNGL(A(3,1)))

C     (SET TEMP=PIVOT AND DET=PIVOT*DET.)
      IF(T1.GE.T2) GO TO 1
      IF(T3.GE.T2) GO TO 2
C     (PIVOT IS A21)
      TEMP=A(2,1) 
      DET=C13*C32-C12*C33
      GO TO 3

    1 IF(T3.GE.T1) GO TO 2
C     (PIVOT IS A11)
      TEMP=A(1,1)
      DET=C22*C33-C23*C32
       GO TO 3
C     (PIVOT IS A31)
 2    TEMP=A(3,1)
      DET=C23*C12-C22*C13

C     SET ELEMENTS OF INVERSE IN A.
    3 IF(DET.EQ.0D0) GO TO 8
      S=TEMP/DET
      A(1,1)=S*C11
      A(1,2)=S*C21
      A(1,3)=S*C31
      A(2,1)=S*C12
      A(2,2)=S*C22
      A(2,3)=S*C32
      A(3,1)=S*C13
      A(3,2)=S*C23
      A(3,3)=S*C33
      RETURN

    4 IF(N.LT.2) GO TO 5

C     N=2 CASE BY CRAMERS RULE.

      DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)
      IF(DET.EQ.0D0) GO TO 8
      S=1D0/DET
      C11   =S*A(2,2)
      A(1,2)=-S*A(1,2)
      A(2,1)=-S*A(2,1)
      A(2,2)=S*A(1,1)
      A(1,1)=C11
      RETURN

C     N=1 CASE.

    5 IF(A(1,1).EQ.0D0) GO TO 8
      A(1,1)=1D0/A(1,1)
      RETURN

C     N.GT.3 CASES.  FACTORIZE MATRIX AND INVERT.

    6 CONTINUE
      CALL DFACT(N,A,IDIM,IR,IFAIL,DET,JFAIL)
      IF(IFAIL.NE.0) RETURN
      CALL DFINV(N,A,IDIM,IR)
      RETURN

C     ERROR EXITS.

    7 IFAIL=+1
*      CALL F010PR(NAME,N,IDIM,K,KPRNT)
      PRINT*,'*** ERROR ERROR *** CERN ERR ROUTINE  F010PR REQD '
      RETURN

    8 IFAIL=-1
      RETURN

      END     



C----------------------------------------------------------------------
         SUBROUTINE          DFACT(N,A,IDIM,IR,IFAIL,DET,JFAIL)
*        INTEGER             IR(9),    IPAIRF
*        DOUBLE PRECISION    A(IDIM,9),DET,      ZERO,     ONE,X,Y,TF
         INTEGER             IR(N),    IPAIRF
         DOUBLE PRECISION    A(IDIM,N),DET,      ZERO,     ONE,X,Y,TF
         REAL                G1,       G2
         REAL                PIVOTF,   P,        Q,        SIZEF,  T
         DOUBLE PRECISION    S11, S12, DOTF
         CHARACTER*6         HNAME
         IPAIRF(J,K)  =  J*2**12 + K
         PIVOTF(X)    =  ABS(SNGL(X))
         SIZEF(X)     =  ABS(SNGL(X))
         DOTF(X,Y,S11)  =  X * Y + S11
         DATA      G1, G2              /  1.E-19,  1.E19  /
         DATA      HNAME               /  ' DFACT'  /
         DATA      ZERO, ONE           /  0.D0, 1.D0  /
         DATA      NORMAL, IMPOSS      /  0, -1  /
         DATA      JRANGE, JOVER, JUNDER  /  0, +1, -1  /
         

         IF(IDIM .GE. N  .AND.  N .GT. 0)  GOTO 110
*        CALL TMPRNT(HNAME,N,IDIM,0)
         PRINT*,'*** ERROR ERROR *** CERN ERR ROUTINE TMPRNT REQD'
         RETURN
         
 110     IFAIL  =  NORMAL
         JFAIL  =  JRANGE
         NXCH   =  0
         DET    =  ONE
         DO 144    J  =  1, N
 120        K  =  J
            P  =  PIVOTF(A(J,J))
            IF(J .EQ. N)  GOTO 122
            JP1  =  J+1
            DO 121    I  =  JP1, N
               Q  =  PIVOTF(A(I,J))
               IF(Q .LE. P)  GOTO 121
               K  =  I
               P  =  Q
121         CONTINUE
            IF(K .NE. J)  GOTO 123
122         IF(P .GT. 0.)  GOTO 130
            DET    =  ZERO
            IFAIL  =  IMPOSS
            JFAIL  =  JRANGE
            PRINT 888,'DFACT IN MATRIX INVERSION - DET ZERO'
            PRINT 999,N,J,P,Q,T,DET
888         FORMAT(1X,A)
999         FORMAT(1X,'N ',I3,'   ROW ',I3,'   P= ',E15.7/
     *             1X,'Q= ',E15.7,'   T= ',E15.7,'   DET= ',E15.7)
            RETURN
            
123         DO 124    L  =  1, N
               TF      =  A(J,L)
               A(J,L)  =  A(K,L)
               A(K,L)  =  TF
124         CONTINUE
            NXCH      =  NXCH + 1
            IR(NXCH)  =  IPAIRF(J,K)
130         DET     =  DET * A(J,J)
            A(J,J)  =  ONE / A(J,J)
            T  =  SIZEF(DET)
            IF(T .LT. G1)  GOTO 132
            IF(T .GT. G2)  GOTO 133
131         IF(J .EQ. N)  GOTO 144
            GOTO 140
            
132         DET    =  ZERO
            IF(JFAIL .EQ. JRANGE)  JFAIL  =  JUNDER
C           PRINT 888,'DFACT IN MATRIX INVERSION UNDER-FLOWING'
C           PRINT 999,N,J,P,Q,T,DET
            GOTO 131
            
133         DET    =  ONE
            IF(JFAIL .EQ. JRANGE)  JFAIL  =  JOVER
C           PRINT 888,'DFACT IN MATRIX INVERSION OVER-FLOWING'
C           PRINT 999,N,J,P,Q,T,DET
            GOTO 131
            
140         JM1  =  J-1
            JP1  =  J+1
            DO 143   K  =  JP1, N
               S11  =  -A(J,K)
               S12  =  -A(K,J+1)
               IF(J .EQ. 1)  GOTO 142
               DO 141  I  =  1, JM1
                  S11  =  DOTF(A(I,K),A(J,I),S11)
                  S12  =  DOTF(A(I,J+1),A(K,I),S12)
141            CONTINUE
142            A(J,K)    =  -S11 * A(J,J)
               A(K,J+1)  =  -DOTF(A(J,J+1),A(K,J),S12)
143         CONTINUE
144      CONTINUE
150      IF(MOD(NXCH,2) .NE. 0)  DET  =  -DET
         IF(JFAIL .NE. JRANGE)   DET  =  ZERO
         IR(N)  =  NXCH
         RETURN
         END
C----------------------------------------------------------------------
         SUBROUTINE          DFINV(N,A,IDIM,IR)
*        INTEGER             IR(9)
*        DOUBLE PRECISION    A(IDIM,9),ZERO,     X, Y, TI
         INTEGER             IR(N)
         DOUBLE PRECISION    A(IDIM,N),ZERO,     X, Y, TI
         DOUBLE PRECISION    S31, S32, S33, S34, DOTF
         CHARACTER*6         HNAME
         DATA      HNAME               /  ' DFINV'  /
         DOTF(X,Y,S31)  =  X*Y + S31
         DATA      ZERO      /  0.D0  /
         IF(IDIM .GE. N  .AND.  N .GT. 0)  GOTO 310
*        CALL TMPRNT(HNAME,N,IDIM,0)
         PRINT*,'*** ERROR ERROR *** CERN ERR ROUTINE TMPRNT REQUIRED '
         RETURN                        
         
310      IF(N .EQ. 1)  RETURN
         A(2,1)  =  -A(2,2) * DOTF(A(1,1),A(2,1),ZERO)
         A(1,2)  =  -A(1,2)
         IF(N .EQ. 2)  GOTO 330
         DO 314    I  =  3, N
            IM2  =  I-2
            DO 312 J  =  1, IM2
               S31  =  ZERO
               S32  =  A(J,I)
               DO 311  K  =  J, IM2
                  S31  =  DOTF(A(K,J),A(I,K),S31)
                  S32  =  DOTF(A(J,K+1),A(K+1,I),S32)
311            CONTINUE
               A(I,J)  =  -A(I,I) * DOTF(A(I-1,J),A(I,I-1),S31)
               A(J,I)  =  -S32
312         CONTINUE
            A(I,I-1)  =  -A(I,I) * DOTF(A(I-1,I-1),A(I,I-1),ZERO)
            A(I-1,I)  =  -A(I-1,I)
314      CONTINUE
330      NM1  =  N-1
         DO 335   I  =  1, NM1
            NMI  =  N-I
            DO 332   J  =  1, I
               S33  =  A(I,J)
               DO 331   K  =  1, NMI
                  S33  =  DOTF(A(I+K,J),A(I,I+K),S33)
331            CONTINUE
               A(I,J)  =  S33
332         CONTINUE
            DO 334   J  =  1, NMI
               S34  =  ZERO
               DO 333   K  =  J, NMI
                  S34  =  DOTF(A(I+K,I+J),A(I,I+K),S34)
333            CONTINUE
               A(I,I+J)  =  S34
334         CONTINUE
335      CONTINUE
         NXCH  =  IR(N)
         IF(NXCH .EQ. 0)  RETURN
         DO 342 M  =  1, NXCH
            K   =  NXCH - M+1
            IJ  =  IR(K)
            I   =  IJ / 4096
            J   =  MOD(IJ,4096)
            DO 341  K  =  1, N
               TI      =  A(K,I)
               A(K,I)  =  A(K,J)
               A(K,J)  =  TI
341         CONTINUE
342      CONTINUE
         RETURN
         END
C----------------------------------------------------------------------
*** ERROR HANDLING SUBROUTINE **
*          SUBROUTINE          TMPRNT(NAME,N,IDIM,K)
*          CHARACTER*6         NAME
*          LOGICAL             MFLAG,    RFLAG
*          IF(NAME(2:2) .EQ. 'S') THEN
*             CALL KERMTR('F012.1',LGFILE,MFLAG,RFLAG)
*          ELSE
*             CALL KERMTR('F011.1',LGFILE,MFLAG,RFLAG)
*          ENDIF
*          IF(NAME(3:6) .EQ. 'FEQN') ASSIGN 1002 TO IFMT
*          IF(NAME(3:6) .NE. 'FEQN') ASSIGN 1001 TO IFMT
*          IF(MFLAG) THEN
*             IF(LGFILE .EQ. 0) THEN
*                IF(NAME(3:6) .EQ. 'FEQN') THEN
C                   WRITE(*,IFMT) NAME, N, IDIM, K
*                ELSE
C                   WRITE(*,IFMT) NAME, N, IDIM
*                ENDIF
*             ELSE
*              IF(NAME(3:6) .EQ. 'FEQN') THEN
*                  WRITE(LGFILE,IFMT) NAME, N, IDIM, K
*                ELSE
*                   WRITE(LGFILE,IFMT) NAME, N, IDIM
*                ENDIF
*             ENDIF
*          ENDIF
*          IF(.NOT. RFLAG) CALL ABEND
*          RETURN
*1001      FORMAT(7X, 31H PARAMETER ERROR IN SUBROUTINE , A6,
*     +             27H ... (N.LT.1 OR IDIM.LT.N).,
*     +             5X, 3HN =, I4, 5X, 6HIDIM =, I4, 1H. )
*1002      FORMAT(7X, 31H PARAMETER ERROR IN SUBROUTINE , A6,
*     +             37H ... (N.LT.1 OR IDIM.LT.N OR K.LT.1).,
*     +             5X, 3HN =, I4, 5X, 6HIDIM =, I4, 5X, 3HK =, I4,1H.)
*          END
*----------------------------------------------------------------------
