      SUBROUTINE PADE8M (MAX,N,A,EA,MDIG,IDIG,WK,IERR)
C
C ****
C
C  FUNCTION - COMPUTES THE EXPONENTIAL OF THE MATRIX A USING THE EIGHTH
C             DIAGONAL PADE TABLE APPROXIMATION AND RETURNS AN ESTIMATE
C             OF THE ACCURACY.
C
C ****
C  NOTE:  PADE8M IS A SLIGHT MODIFICATION OF WARD'S PADE8.
C         NOTE THAT NEITHER PADE8 NOR PADE8M ARE PORTABLE
C         BUT ARE SUITABLE FOR IBM 360/370 OR AMDAHL 470.
C
C  INPUT PARAMETERS
C
C    MAX - THE ROW (OR FIRST) DIMENSION OF THE ARRAYS A AND EA AS IT
C          APPEARS IN THE PROGRAM CALLING PADE8. THIS PARAMETER IS USED
C          TO SPECIFY THE VARIABLE DIMENSION OF THESE ARRAYS. MAX MUST
C          BE GREATER THAN OR EQUAL TO N.
C
C    N - THE ORDER OF THE INPUT MATRIX A. THE OUTPUT EXPONENTIAL MATRIX
C        EA WILL ALSO BE OF ORDER N.
C
C    A - A TWO-DIMENSIONAL ARRAY WHICH SPECIFIES THE (I,J) COMPONENTS
C        OF THE INPUT MATRIX WHOSE EXPONENTIAL IS DESIRED. THE ROW
C        DIMENSION OF A MUST BE MAX AND THE COLUMN DIMENSION MUST BE
C        GREATER THAN OR EQUAL TO N. THE LAST DIGIT OF THE ELEMENTS
C        OF A MAY BE ALTERED ON RETURN.
C
C    EA - A TWO-DIMENSIONAL ARRAY PROVIDED BY THE USER FOR RETURNING
C         THE EXPONENTIAL OF THE MATRIX A. THE ROW DIMENSION OF EA MUST
C         BE MAX AND THE COLUMN DIMENSION MUST BE GREATER THAN OR EQUAL
C         TO N.
C
C    WK - AN ARRAY OF ANY DIMENSION THAT MUST BE PROVIDED BY THE USER
C         FOR WORK SPACE. WK MUST CONTAIN AT LEAST N**2 + 8*N STORAGE
C         LOCATIONS.
C
C  OUTPUT PARAMETERS
C
C    EA - CONTAINS THE (I,J) COMPONENT OF THE EXPONENTIAL OF THE MATRIX
C         A IN THE LOCATION EA(I,J).
C
C    MDIG - AN ESTIMATE FOR THE MINIMUM NUMBER OF DIGITS ACCURATE IN
C           THE NORM OF EA.
C
C    IDIG - AN ESTIMATE FOR THE NUMBER OF DIGITS ACCURATE IN THE NORM
C           OF EA AT THE 95 PER CENT CONFIDENCE LEVEL.
C
C    IERR - AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS OR
C           FAILURE TO APPROXIMATE THE EXPONENTIAL ACCURATELY.
C
C           = 0   NO ERROR  .
C           = 1   MAX .LT. N  .
C           = 2   NORM OF A .GT. 2**40  .
C           = 3   MDIG .EQ. 0  AND  IDIG .GT. 0  .
C           = 4   MDIG .EQ. 0  AND  IDIG .EQ. 0  .
C
C           FOR ERROR NUMBERS 1 AND 2, AN APPROXIMATION IS NOT
C           ATTEMPTED. ERROR NUMBER 3 SHOULD ACT AS A WARNING OF
C           POSSIBLE INACCURACY. ERROR NUMBER 4 SHOULD ACT AS A SEVERE
C           WARNING OF POSSIBLE INACCURACY.
C
C  REQUIRED SUBPROGRAMS - BALANC(EISPACK),BALINV,GAUSEL
C
C  REQUIRED FUNCTIONS - DABS,DEXP,DBLE,DLOG10,IFIX,MOD,SNGL
C
C  AUTHOR/IMPLEMENTER/MODIFIER -  R.C.WARD / R.C.WARD / A.J.LAUB
C
C  LANGUAGE - FORTRAN (IBM DOUBLE PRECISION)
C
C  DATE RELEASED - NOV. 14, 1975
C
C  LATEST VERSION - NOV. 19, 1981
C
C ****
C
      DOUBLE PRECISION A(MAX,N),EA(MAX,N),WK(N,1),C(8)
      DOUBLE PRECISION TR,AVGEV,FACTOR,S,U,P,EAVGEV,EPS,SS,RERL,ANORM,
     1       EMNORM,BIG,SMAL,ENORM,XN,TEMP,VAR,VAREPS,FN,GN,SIZE,SIZEJ,
     2       BD,SD2,SUM2D,EPSP1,ABS,RERR
      DOUBLE PRECISION DBLE,DABS,DEXP,DLOG10
C   *********************************************
      DATA C / .5D0,.11666666666666667D0,1.6666666666666667D-2,
     1         1.6025641025641026D-3,1.0683760683760684D-4,
     2         4.8562548562548563D-6,1.3875013875013875D-7,
     3         1.9270852604185938D-9 /
C
C     DETERMINE MACHINE PRECISION
C
      EPS = 1.0D0
 3    EPS = EPS/2.0D0
      EPSP1 = EPS+1.0D0
      IF (EPSP1 .GT. 1.0D0) GO TO 3
      EPS = 2.0D0*EPS
C
      IERR = 0
      IF (N .EQ. 1) GO TO 450
      IF (MAX .GE. N) GO TO 5
      IERR = 1
      GO TO 500
 5    XN = DBLE(N)
      FN = (4.D0*XN) / ((XN+2.D0)*(XN+1.D0))
      GN = (2.D0*(XN**2)+10.D0*XN-4.D0) / (((XN+2.D0)**2) *
     1                                     ((XN+1.D0)**2))
      VAREPS = (255.D0*(16.D0**(-26))) / (24.D0*DLOG(16.D0))
      ISCA = N + 8
C ****
C     SHIFT ORIGIN BY EIGENVALUE AVERAGE AND BALANCE
C ****
      TR = 0.D0
      DO 10 I=1,N
         TR = TR + A(I,I)
 10   CONTINUE
      AVGEV = TR/XN
      DO 20 I=1,N
         A(I,I) = A(I,I) - AVGEV
 20   CONTINUE
      CALL BALANC (MAX,N,A,LOW,IGH,WK(1,ISCA))
C ****
C     CALCULATE NORM OF BALANCED,SHIFTED A
C ****
      ANORM = 0.D0
      DO 40 J=1,N
         SS = 0.D0
         DO 30 I=1,N
            SS = SS + DABS(A(I,J))
 30      CONTINUE
         IF (SS .GT. ANORM) ANORM = SS
 40   CONTINUE
C ****
C     DETERMINE POWER-OF-TWO NORMALIZATION FACTOR
C ****
      M = 0
      IF (ANORM .LE. 1.D0) GO TO 80
      FACTOR = 2.D0
      DO 50 M=1,40
         IF (ANORM .LE. FACTOR) GO TO 60
         FACTOR = FACTOR * 2.D0
 50   CONTINUE
      IERR = 2
      GO TO 500
C ****
C     NORMALIZE A
C ****
 60   DO 70 I=1,N
         DO 70 J=1,N
            A(I,J) = A(I,J) / FACTOR
 70   CONTINUE
 80   IFIR = N+1
      DO 160 J=1,N
C    ****
C        COMPUTE J-TH COLUMN OF FIRST EIGHT POWERS OF A
C    ****
         DO 100 I=1,N
            S = 0.D0
            DO 90 L=1,N
               S = S + A(I,L)*A(L,J)
 90         CONTINUE
            WK(I,IFIR) = S
 100     CONTINUE
         DO 120 K=1,6
            KP1 = K+1
            IPOW = N+K
            ISUB = IPOW+1
            DO 120 I=1,N
               S = 0.D0
               DO 110 L=1,N
                 S = S + A(I,L)*WK(L,IPOW)
 110           CONTINUE
               WK(I,ISUB) = S
 120     CONTINUE
C    ****
C        COLLECT TERMS FOR J-TH COLUMN OF PADE NUMERATOR AND DENOMINATOR
C    ****
         DO 160 I=1,N
            S = 0.D0
            U = 0.D0
            DO 140 L=1,7
               K = 8-L
               KP1 = K+1
               IPOW = N+K
               P = C(KP1)*WK(I,IPOW)
               S = S + P
               KEO = MOD(KP1,2)
               IF (KEO .EQ. 0) GO TO 130
               U = U - P
               GO TO 140
 130           U = U + P
 140        CONTINUE
            P = C(1)*A(I,J)
            S = S + P
            U = U - P
            IF (I .NE. J) GO TO 150
            S = S + 1.D0
            U = U + 1.D0
 150        EA(I,J) = S
            WK(I,J) = U
 160  CONTINUE
C ****
C     CALCULATE NORMALIZED EXP(A) BY WK * EXP(A) = EA
C ****
      CALL GAUSEL (MAX,N,WK,N,EA,IERR)
      ENORM = 0.D0
      DO 200 J=1,N
         SS = 0.D0
          DO 190 I=1,N
            SS = SS + DABS(EA(I,J))
 190      CONTINUE
          IF (SS .GT. ENORM) ENORM = SS
 200  CONTINUE
      ABS = (19.D0 * XN + 47.D0) * EPS * ENORM
      VAR = (12.D0*XN)*VAREPS
      IF (M .EQ. 0) GO TO 285
C ****
C     REMOVE EFFECT OF NORMALIZATION ON EXP(A) - TESTING ROUNDOFF ERROR
C ****
      VAR = XN*VAREPS
      DO 230 K=1,M
         DO 220 I=1,N
            DO 220 J=1,N
               S = 0.D0
               DO 210 L=1,N
                 S = S + EA(I,L)*EA(L,J)
 210           CONTINUE
               WK(I,J) = S
 220     CONTINUE
         TEMP = ENORM**2
         ABS = 2.D0*ABS*ENORM + XN*EPS*TEMP + ABS**2
         VAR = (FN*VAR + GN*TEMP*VAREPS) * TEMP
         ENORM = 0.D0
         DO 230 J=1,N
            SS = 0.D0
            DO 225 I=1,N
               EA(I,J) = WK(I,J)
               SS = SS + DABS(EA(I,J))
 225        CONTINUE
            IF (SS .GT. ENORM) ENORM = SS
 230  CONTINUE
C ****
C     RETURN A TO ORIGINAL INPUT
C ****
      DO 280 I=1,N
         DO 280 J=1,N
            A(I,J) = A(I,J)*FACTOR
 280  CONTINUE
 285  CALL BALINV (MAX,N,LOW,IGH,WK(1,ISCA),A)
      DO 290 I=1,N
         A(I,I) = A(I,I) + AVGEV
 290  CONTINUE
C ****
C     CALCULATE 2 STANDARD DEVIATION BOUND ON PROBABILISTIC ERROR
C ****
      BIG = 1.D0
      SMAL = 1.D0
      SUM2D = 0.D0
      IF (LOW .EQ. IGH) GO TO 320
      SUM2D = -1.D0
      DO 310 I=LOW,IGH
         U = WK(I,ISCA)
         SUM2D = SUM2D + U**2
         IF (BIG .GE. U) GO TO 300
         BIG = U
 300     IF (SMAL .LE. U) GO TO 310
         SMAL = U
 310  CONTINUE
 320  SUM2D = SUM2D + XN + DBLE(LOW - IGH)
      SD2 = 2.D0*DSQRT(SUM2D*VAR + (8.D0*DBLE(IGH-LOW)*(ENORM**2)
     1                                  *SUM2D*VAREPS) / (XN**2))
C ****
C     REMOVE EFFECT OF BALANCE AND ORIGIN SHIFT ON EXP(A) AND DETERMINE
C       MAXIMUM PROBABILISTIC ERROR IN ITS NORM
C ****
      CALL BALINV (MAX,N,LOW,IGH,WK(1,ISCA),EA)
      EAVGEV = DEXP(AVGEV)
      EMNORM = 0.D0
      SIZE = 0.D0
      DO 340 J=1,N
         SS = 0.D0
         DO 330 I=1,N
            SS = SS + DABS(EA(I,J))
            EA(I,J) = EA(I,J)*EAVGEV
 330     CONTINUE
         IF (SS .GT. EMNORM) EMNORM = SS
         IF (LOW .EQ. IGH) GO TO 340
         BD = WK(J,ISCA)
         IF (J .LT. LOW .OR. IGH .LT. J) BD = 1.D0
         SIZEJ = SS + SD2/BD
         IF (SIZEJ .GT. SIZE) SIZE = SIZEJ
 340  CONTINUE
C ****
C     COMPUTE RETURNED INTEGRAL ACCURACY ESTIMATES
C ****
      RERR = (BIG*ABS) / (SMAL*EMNORM*EPS)
      IF (LOW .EQ. IGH) SIZE = EMNORM + SD2
      RERL = (SIZE-EMNORM) / (EMNORM*EPS)
      MDIG = 17 - IFIX(SNGL(DLOG10(RERR)) + .5)
      IDIG = 17 - IFIX(SNGL(DLOG10(RERL)) + .5)
      IF (IDIG .GT. 15) IDIG = 15
      IF (MDIG .GT. 0) GO TO 500
      MDIG = 0
      IERR = 3
      IF (IDIG .GT. 0) GO TO 500
      IDIG = 0
      IERR = 4
      GO TO 500
C ****
C     A IS A SCALAR (1X1 MATRIX)
C ****
 450  EA(1,1) = DEXP(A(1,1))
      MDIG = 16
      IDIG = 16
C ****
C     EXIT ROUTINE
C ****
 500  CONTINUE
      RETURN
      END
