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