SUBROUTINE PADE(NM,NB,NQ,N,M,A,B,Q,DELTA,TOL,ICODE,MAXQ,MDEG,
     1                ISCALE,RHO,E,N33,D33,N23,D23,D12,D13,X,R,Z,
     2                N34,D34,N14,D14,N24,D24,D,V,Y,WORK,IPVT)
C
C     *****PARAMETERS:
      INTEGER NM,NB,NQ,N,M,ICODE,MAXQ,ISCALE,IPVT(N)
      DOUBLE PRECISION A(NM,N),B(NB,M),Q(NQ,N),DELTA,TOL,RHO,E(5),
     1                 N33(NM,N),D33(NM,N),N23(NM,N),D23(NM,N),
     2                 D12(NM,N),D13(NM,N),X(NM,N),R(NM,N),Z(NM,N),
     3                 N34(NM,M),D34(NM,M),N14(NM,M),D14(NM,M),
     4                 N24(NM,M),D24(NM,M),D(NM,M),V(NM,M),Y(NM,M),
     5                 WORK(N)
C
C     *****LOCAL VARIABLES:
      INTEGER I,J,K,L,MDEG
      DOUBLE PRECISION ALFA,ANORM,BETA,BNORM,CK,CKS,COND,CNORM,
     +                 EPS,QNORM,S,SIGN,SUM,T
C
C     *****FUNCTIONS:
      INTEGER IDINT
      DOUBLE PRECISION DEXP,DBLE,DSQRT,DMAX1,DLOG10
C
C     *****SUBROUTINES CALLED:
C     FNORM,SAVE,TRNATB,MLINEQ,MMUL,MADD,MULWOB
C
C     ------------------------------------------------------------------
C
C     *****PURPOSE:
C     THIS SUBROUTINE COMPUTES APPROXIMATES F,H,QD,M,W TO
C
C        F(DELTA)  = EXP(A*DELTA)
C
C        H(DELTA)  = INTEGRAL FROM 0 TO DELTA OF   F(S)B
C                                                      T
C        QD(DELTA) = INTEGRAL FROM 0 TO DELTA OF   F(S)  Q F(S)
C                                                      T
C        M(DELTA)  = INTEGRAL FROM 0 TO DELTA OF   F(S)  Q H(S)
C                                                      T
C        W(DELTA)  = INTEGRAL FROM 0 TO DELTA OF   H(S)  Q H(S)
C
C     *****PARAMETER DESCRIPTION:
C
C     ON INPUT:
C        NM               ROW DIMENSION OF THE ARRAYS CONTAINING
C                         A (AND N33,D33,N23,D23,D12,D13,X,R,Z,N34,
C                         D34,N14,D14,N24,D24,D,V,Y) AS DECLARED
C                         IN THE CALLING PROGRAM DIMENSION STATEMENT;
C
C        NB               ROW DIMENSION OF THE ARRAY CONTAINING B
C                         AS DECLARED IN THE CALLING PROGRAM DIMENSION
C                         STATEMENT;
C
C        NQ               ROW DIMENSION OF THE ARRAY CONTAINING Q
C                         AS DECLARED IN THE CALLING PROGRAM DIMENSION
C                         STATEMENT;
C
C        N                ORDER OF THE SQUARE MATRIX A;
C
C        M                COLUMN DIMENSION OF B;
C
C        A                N X N MATRIX;
C
C        B                N X M MATRIX;
C
C        Q                N X N MATRIX;
C
C        DELTA            INTERVAL LENGTH;
C
C        TOL              APPROXIMATE RELATIVE ERROR BOUND FOR
C                         COMPUTED MATRICES;
C
C        ICODE            INTEGER VARIABLE WHICH DETERMINES WHICH
C                         OF F,H,QD,M,W ARE COMPUTED.
C
C                         IF ICODE EQUALS 1, F IS COMPUTED, AND THE
C                         ARRAYS B,Q,D23,N34,D34,N14,D14,N24,D24,
C                         D12,D13,V,Y,Z,D ARE NOT USED.
C
C                         IF ICODE EQUALS 2, F,H ARE COMPUTED,
C                         AND THE ARRAYS Q,N23,D23,N14,D14,N24,D24,
C                         D12,D13,V,Z ARE NOT USED.
C
C                         IF ICODE EQUALS 3, F,QD ARE COMPUTED,
C                         AND THE ARRAYS N34,D34,N14,D14,N24,D24,
C                         D12,D23,D13,V,Y,Z,D ARE NOT USED.
C
C                         IF ICODE EQUALS 4, F,H,QD,M ARE COMPUTED,
C                         AND THE ARRAYS N14,D14,D13,D23,D12,Z ARE NOT
C                         USED.
C
C                         IF ICODE EQUALS 5, F,H,QD,M,W ARE COMPUTED;
C
C        MAXQ             MAXIMUM DEGREE OF PADE APPROXIMANT.
C                         (MAXQ = 10 FOR IBM LONG ARITHMETIC).
C
C     ON OUTPUT:
C
C        MDEG             DEGREE OF PADE APPROXIMANT;
C
C        ISCALE           SCALE PARAMETER;
C
C        RHO              REAL QUANTITY USED IN COMPUTING APPROXIMATE
C                         ERROR BOUNDS. RHO IS USUALLY A SMALL NUMBER
C                         AND THEREFORE KNOWLEDGE OF ITS EXACT VALUE IS
C                         NOT TOO CRITICAL.  FURTHERMORE, ONE IS OFTEN
C                         INTERESTED IN RELATIVE ERROR BOUNDS, AND IN
C                         THIS SENSE TOL IS A FAIR MEASURE;
C
C        E                REAL VECTOR OF LENGTH 5 CONTAINING APPROXIMATE
C                         RELATIVE ERROR BOUNDS WHICH SHOULD BE LESS
C                         THAN TOL. IF THEY ARE NOT, THIS MEANS THAT
C                         THE ACCURACY COULD NOT BE ACHIEVED FOR THE
C                         GIVEN MAXQ;
C
C        N33              N X N ARRAY CONTAINING F.
C                         (APPROXIMATE ERROR BOUND = E(1)*RHO);
C
C        N23              N X N ARRAY CONTAINING QD.
C                         (APPROXIMATE ERROR BOUND = E(2)*RHO);
C
C        N34              N X M ARRAY CONTAINING H.
C                         (APPROXIMATE ERROR BOUND = E(3)*RHO*RHO);
C
C        N14              N X M ARRAY CONTAINING W IN ITS UPPER LEFT
C                         M X M BLOCK.
C                         (APPROXIMATE ERROR BOUND = E(5)*RHO*RHO);
C
C        N24              N X M ARRAY CONTAINING M.
C                         (APPROXIMATE ERROR BOUND = E(4)*RHO*RHO);
C
C        D33,D23,D13,     N X N REAL SCRATCH ARRAYS;
C        D12,X,R,Z
C
C        D34,D14,D24,     N X M REAL SCRATCH ARRAYS;
C        D,V,Y
C
C        WORK             REAL SCRATCH VECTOR OF LENGTH N;
C
C        IPVT             INTEGER SCRATCH VECTOR OF LENGTH N.
C
C     *****REFERENCES:
C        (1)  VAN LOAN, CHARLES F., COMPUTING INTEGRALS INVOLVING THE
C             MATRIX EXPONENTIAL, IEEE TRANSACTIONS ON AUTOMATIC
C             CONTROL, VOL. AC-23, NO.3, JUNE 1978, PAGES 395-404.
C
C     *****HISTORY:
C     WRITTEN BY C.F. VANLOAN,
C     MODIFIED BY J.A.K. CARRIG, LIDS, MIT, NOV. 1978,
C     MODIFIED BY ALAN J. LAUB, UCSB, FEB. 1982; FEB. 1984.
C     MOST RECENT VERSION: FEB. 21, 1984.
C
C     ------------------------------------------------------------------
C
C     DETERMINE PADE PARAMETER MDEG AND SCALE POWER ISCALE
C
      BNORM=0.0D0
      QNORM=0.0D0
      CALL FNORM (NM,N,N,A,ANORM)
      IF(ICODE.EQ.2 .OR. ICODE.GE.4)   CALL FNORM(NB,N,M,B,BNORM)
      IF(ICODE.GE.3)   CALL FNORM(NQ,N,N,Q,QNORM)
      ALFA = DMAX1(BNORM,QNORM)
      S = DMAX1(ANORM,ALFA)
      IF (S .EQ. 0.0D0) GO TO 5
      ANORM = ANORM/S
      BNORM = BNORM/S
      QNORM = QNORM/S
    5 CONTINUE
      CNORM = DSQRT(ANORM**2 + BNORM**2 + QNORM**2)
      IF(ICODE.GE.3)   CNORM = DSQRT(CNORM**2 + ANORM**2)
      IF(ICODE.EQ.5)   CNORM = DSQRT(CNORM**2 + ANORM**2 +
     +                               (DBLE(N)/S)**2)
      CNORM = S*CNORM
      ISCALE = IDINT(2.0D0 + DLOG10(CNORM*DELTA)/DLOG10(2.0D0))
      IF(ISCALE.LT.0)   ISCALE = 1
C
      DO 10 I=1,5
         E(I) = 0.0D0
   10 CONTINUE
      EPS = 8.0D0*CNORM
      DO 20 K=1,MAXQ
         MDEG = K
         EPS = EPS/DBLE(16*(4*K*K-1))
         T = EPS*DELTA*DEXP(EPS*DELTA)
         E(1) = T
         IF(ICODE.EQ.2.OR.ICODE.GE.4) E(2) = T*(1.0D0+ALFA*DELTA/2.0D0)
         T = T*DEXP(EPS*DELTA)
         IF(ICODE.GE.3) E(3) = T*(1.0D0+ALFA*DELTA)
         IF(ICODE.GE.4) E(4) = T*(1.0D0+EPS+ALFA*DELTA)**2
         IF(ICODE.EQ.5) E(5) = 4.0D0*T*(1.0D0+1.5D0*(ALFA+EPS)*DELTA)**3
         IF(E(ICODE).LE.TOL)   GO TO 30
   20 CONTINUE
   30 CONTINUE
      BETA=DELTA/(2.0D0**ISCALE)
      SIGN=-1.0D0
      CK=0.5D0*BETA
      CKS=-CK
C
C     INITIALIZE NUMERATOR AND DENOMINATOR MATRICES
C
C     IF ICODE .GE. 1, INITIALIZE N33,D33,X
C
      CALL SAVE(NM,NM,N,N,A,X)
      DO 50 J=1,N
         DO 40 I=1,N
            N33(I,J)=CK*X(I,J)
            D33(I,J)=CKS*X(I,J)
   40    CONTINUE
         N33(J,J)=N33(J,J)+1.0D0
         D33(J,J)=D33(J,J)+1.0D0
   50 CONTINUE
C
C     IF ICODE .EQ. 2 OR ICODE .GE.4, INITIALIZE N34,D34
C
      IF(ICODE.EQ.1.OR.ICODE.EQ.3) GO TO 80
C
      CALL SAVE(NB,NM,N,M,B,D)
      DO 70 J=1,M
         DO 70 I=1,N
            N34(I,J)=CK*D(I,J)
            D34(I,J)=CKS*D(I,J)
   70 CONTINUE
C
C     IF ICODE .GE. 3, INITIALIZE N23,D23,R
C
   80 IF(ICODE.LT.3) GO TO 90
C
      CALL SAVE(NQ,NM,N,N,Q,R)
      CALL SAVE(NM,NM,N,N,R,N23)
      CALL SAVE(NM,NM,N,N,R,D23)
      CALL MSCALE(NM,N,N,CK,N23)
      CALL MSCALE(NM,N,N,CKS,D23)
C
C     IF ICODE .GE. 4, INITIALIZE N24,D24
C
   90 IF(ICODE.LT.4) GO TO 110
C
      DO 100 J=1,M
         DO 100 I=1,N
            N24(I,J)=0.0D0
            D24(I,J)=0.0D0
  100 CONTINUE
C
C     IF ICODE .EQ. 5, INITIALIZE N14,D14,D12
C
  110 IF(ICODE.NE.5) GO TO 150
C
      DO 140 J=1,N
         DO 120 I=1,N
            D12(I,J)=0.0D0
            D13(I,J)=0.0D0
  120    CONTINUE
         DO 130 I=1,M
            D14(J,I)=0.0D0
            N14(J,I)=0.0D0
  130    CONTINUE
         D12(J,J)=-CK
  140 CONTINUE
  150 CONTINUE
C
C     BUILD UP PADE NUMERATOR AND DENOMINATOR
C
      IF(MAXQ .EQ. 1) GO TO 560
      DO 550 K=2,MDEG
         SIGN = -SIGN
         CK = CK*BETA*DBLE(MDEG-K+1)/DBLE(K*(2*MDEG-K+1))
         CKS=SIGN*CK
C
         IF(ICODE .LT. 5)   GO TO 360
C
C        SKIP IF DO NOT WANT W
C
         IF(K.EQ.2) GO TO 290
         DO 270 J=1,M
            DO 230 I=1,N
               WORK(I)=Y(I,J)
  230       CONTINUE
            DO 260 I=1,N
              SUM=0.0D0
              IF(K.EQ.3) GO TO 250
              DO 240 L=1,N
                 SUM=SUM-A(L,I)*WORK(L)
  240         CONTINUE
  250         Y(I,J)=SUM+V(I,J)
  260       CONTINUE
  270    CONTINUE
         DO 280 J=1,M
            DO 280 I=1,N
               N14(I,J)=N14(I,J)+CK*Y(I,J)
               D14(I,J)=D14(I,J)+CKS*Y(I,J)
  280    CONTINUE
  290    CONTINUE
         DO 340 J=1,N
            DO 300 I=1,N
               WORK(I)=Z(I,J)
  300       CONTINUE
            DO 330 I=1,N
               SUM=0.0D0
               IF(K.EQ.2) GO TO 320
               DO 310 L=1,N
                  SUM=SUM-A(L,I)*WORK(L)
  310          CONTINUE
  320          Z(I,J)=SUM+R(I,J)
  330       CONTINUE
  340    CONTINUE
C
         T = -SIGN*DBLE(K)
         DO 350 J=1,N
            DO 350 I=1,N
               D12(I,J)=D12(I,J)+CKS*T*X(J,I)
               D13(I,J)=D13(I,J)+CKS*Z(I,J)
  350    CONTINUE
C
  360    IF(ICODE.LT.4)   GO TO 440
C
C        SKIP IF DO NOT WANT M OR BOTH H AND Q
C
         IF(K.EQ.2) CALL MMUL(NQ,NB,NM,M,N,N,Q,B,V)
         IF(K.EQ.2) GO TO 420
         DO 410 J=1,M
            DO 370 I=1,N
               WORK(I)=V(I,J)
  370       CONTINUE
            DO 400 I=1,N
               SUM=0.0D0
               DO 380 L=1,N
                  SUM=SUM-A(L,I)*WORK(L)
  380          CONTINUE
               DO 390 L=1,N
                  SUM=SUM+Q(I,L)*D(L,J)
  390          CONTINUE
               V(I,J)=SUM
  400       CONTINUE
  410    CONTINUE
  420    CONTINUE
         DO 430 J=1,M
            DO 430 I=1,N
               N24(I,J)=N24(I,J)+CK*V(I,J)
               D24(I,J)=D24(I,J)+CKS*V(I,J)
  430    CONTINUE
C
  440    CONTINUE
         IF(ICODE.LT.3)   GO TO 510
C
C        SKIP IF DO NOT WANT Q
C
         DO 490 J=1,N
            DO 450 I=1,N
               WORK(I)=R(I,J)
  450       CONTINUE
            DO 480 I=1,N
               SUM=0.0D0
               DO 460 L=1,N
                  SUM=SUM-A(L,I)*WORK(L)
  460          CONTINUE
               DO 470 L=1,N
                  SUM=SUM+Q(I,L)*X(L,J)
  470          CONTINUE
               R(I,J)=SUM
  480       CONTINUE
  490    CONTINUE
         DO 500 J=1,N
            DO 500 I=1,N
               N23(I,J)=N23(I,J)+CK*R(I,J)
               D23(I,J)=D23(I,J)+CKS*R(I,J)
  500    CONTINUE
C
  510    IF(ICODE.EQ.3 .OR. ICODE.EQ.1)   GO TO 530
C
C        SKIP IF DO NOT WANT H
C
         CALL MMUL(NM,NB,NM,M,N,N,X,B,D)
         DO 520 J=1,M
            DO 520 I=1,N
               N34(I,J)=N34(I,J)+CK*D(I,J)
               D34(I,J)=D34(I,J)+CKS*D(I,J)
  520    CONTINUE
  525    CONTINUE
C
  530    CALL MULWOB(NM,NM,N,A,X,WORK)
         DO 540 J=1,N
            DO 540 I=1,N
               N33(I,J)=N33(I,J)+CK*X(I,J)
               D33(I,J)=D33(I,J)+CKS*X(I,J)
  540    CONTINUE
  550 CONTINUE
C
C     SOLVE THE LINEAR SYSTEM
C
  560 IF(ICODE.GE.3) CALL TRNATB(NM,NM,N,N,N33,R)
C
      CALL SAVE(NM,NM,N,N,D33,X)
      CALL MLINEQ(NM,NM,N,N,X,N33,COND,IPVT,WORK)
C
      IF(ICODE.EQ.1.OR.ICODE.EQ.3) GO TO 600
C
      CALL MSUB(NM,NM,NM,N,M,N34,D34,N34)
      CALL SAVE(NM,NM,N,N,D33,X)
      CALL MLINEQ(NM,NM,N,M,X,N34,COND,IPVT,WORK)
C
  600 IF(ICODE.LT.3) GO TO 610
C
      CALL MMUL(NM,NM,NM,N,N,N,D23,N33,D33)
      CALL MSUB(NM,NM,NM,N,N,N23,D33,N23)
      CALL SAVE(NM,NM,N,N,R,X)
      CALL MLINEQ(NM,NM,N,N,X,N23,COND,IPVT,WORK)
C
  610 IF(ICODE.LT.4) GO TO 620
C
      CALL MMUL(NM,NM,NM,M,N,N,D23,N34,D33)
      CALL MSUB(NM,NM,NM,N,M,N24,D33,N24)
      CALL MSUB(NM,NM,NM,N,M,N24,D24,N24)
      CALL SAVE(NM,NM,N,N,R,X)
      CALL MLINEQ(NM,NM,N,M,X,N24,COND,IPVT,WORK)
C
  620 IF(ICODE.LT.5) GO TO 630
C
      CALL MMUL(NM,NM,NM,M,N,N,D12,N24,D33)
      CALL MSUB(NM,NM,NM,N,M,N14,D33,N14)
      CALL MMUL(NM,NM,NM,M,N,N,D13,N34,D33)
      CALL MSUB(NM,NM,NM,N,M,N14,D33,N14)
      CALL MSUB(NM,NM,NM,N,M,N14,D14,N14)
      CALL SAVE(NM,NM,N,N,R,X)
      CALL MLINEQ(NM,NM,N,M,X,N14,COND,IPVT,WORK)
C
C     FORM SCALED INTEGRALS
C
  630 IF(ICODE.LE.2) GO TO 700
C
      CALL TRNATB(NM,NM,N,N,N33,X)
      CALL MULWOB(NM,NM,N,X,N23,WORK)
C
      IF(ICODE.LE.3) GO TO 700
C
      CALL TRNATB(NM,NM,N,N,N33,X)
      CALL MMUL(NM,NM,NM,M,N,N,X,N24,V)
      CALL SAVE(NM,NM,N,M,V,N24)
C
      IF(ICODE.LE.4) GO TO 700
C
      CALL TRNATB(NM,NM,N,N,N33,X)
      CALL MMUL(NM,NM,NM,M,N,N,X,N14,Y)
      CALL TRNATB(NB,NM,N,M,B,X)
      CALL MMUL(NM,NM,NM,M,M,N,X,Y,N14)
      CALL TRNATB(NM,NM,M,M,N14,D14)
      CALL MADD(NM,NM,NM,M,M,N14,D14,N14)
C
  700 CONTINUE
C
C     APPLY DOUBLING FORMULAE
C
      CALL FNORM(NM,N,N,N33,RHO)
C
      IF(ISCALE.EQ.0)   RETURN
C
      DO 850 K=1,ISCALE
C
         IF(ICODE.LT.5)   GO TO 810
C
         CALL MADD(NM,NM,NM,M,M,N14,N14,N14)
         CALL TRNATB(NM,NM,N,M,N34,X)
         CALL MMUL(NM,NM,NM,M,M,N,X,N24,Y)
         CALL MADD(NM,NM,NM,M,M,N14,Y,N14)
         CALL TRNATB(NM,NM,M,M,Y,X)
         CALL MADD(NM,NM,NM,M,M,X,N14,N14)
         CALL MMUL(NM,NM,NM,M,N,N,N23,N34,Y)
         CALL TRNATB(NM,NM,N,M,N34,R)
         CALL MMUL(NM,NM,NM,M,M,N,R,Y,X)
         CALL MADD(NM,NM,NM,M,M,X,N14,N14)
C
  810    IF(ICODE.LT.4)   GO TO 820
C
         CALL MMUL(NM,NM,NM,M,N,N,N23,N34,D24)
         CALL MADD(NM,NM,NM,N,M,D24,N24,D24)
         CALL TRNATB(NM,NM,N,N,N33,X)
         CALL MMUL(NM,NM,NM,M,N,N,X,D24,V)
         CALL MADD(NM,NM,NM,N,M,V,N24,N24)
C
  820    IF(ICODE.LT.3)   GO TO 830
C
         CALL MQF(NM,NM,NM,N,N,N23,N33,X,WORK)
         CALL MADD(NM,NM,NM,N,N,X,N23,N23)
C
  830    IF(ICODE.EQ.1 .OR. ICODE.EQ.3)   GO TO 840
C
         CALL MMUL(NM,NM,NM,M,N,N,N33,N34,D34)
         CALL MADD(NM,NM,NM,N,M,N34,D34,N34)
C
  840    CALL SAVE(NM,NM,N,N,N33,X)
         CALL MMUL(NM,NM,NM,N,N,N,N33,X,D33)
         CALL SAVE(NM,NM,N,N,D33,N33)
         CALL FNORM(NM,N,N,N33,T)
         IF(RHO.LT.T)   RHO=T
  850 CONTINUE
C
      RETURN
C
C     LAST LINE OF PADE
C
      END