intdy.f

Go to the documentation of this file.

      SUBROUTINE intdy (T, K, YH, NYH, DKY, IFLAG)
CLLL. OPTIMIZE
      INTEGER K, NYH, IFLAG
      INTEGER ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
     1   MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, CNYH,
     2   IALTH, IPUP, LMAX, MEO, NQNYH, NSLP
      INTEGER ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
     2   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
      INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
      DOUBLE PRECISION T, YH, DKY
      DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
     1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
      DOUBLE PRECISION C, R, S, TP
      dimension yh(nyh,*), dky(*)
      COMMON /ls0001/ conit, crate, el(13), elco(13,12),
     1   hold, rmax, tesco(3,12),
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     2   illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
     3   mxstep, mxhnil, nhnil, ntrep, nslast, cnyh,
     3   ialth, ipup, lmax, meo, nqnyh, nslp,
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
C-----------------------------------------------------------------------
C INTDY COMPUTES INTERPOLATED VALUES OF THE K-TH DERIVATIVE OF THE
C DEPENDENT VARIABLE VECTOR Y, AND STORES IT IN DKY.  THIS ROUTINE
C IS CALLED WITHIN THE PACKAGE WITH K = 0 AND T = TOUT, BUT MAY
C ALSO BE CALLED BY THE USER FOR ANY K UP TO THE CURRENT ORDER.
C (SEE DETAILED INSTRUCTIONS IN THE USAGE DOCUMENTATION.)
C-----------------------------------------------------------------------
C THE COMPUTED VALUES IN DKY ARE GOTTEN BY INTERPOLATION USING THE
C NORDSIECK HISTORY ARRAY YH.  THIS ARRAY CORRESPONDS UNIQUELY TO A
C VECTOR-VALUED POLYNOMIAL OF DEGREE NQCUR OR LESS, AND DKY IS SET
C TO THE K-TH DERIVATIVE OF THIS POLYNOMIAL AT T.
C THE FORMULA FOR DKY IS..
C              Q
C  DKY(I)  =  SUM  C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1)
C             J=K
C WHERE  C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR.
C THE QUANTITIES  NQ = NQCUR, L = NQ+1, N = NEQ, TN, AND H ARE
C COMMUNICATED BY COMMON.  THE ABOVE SUM IS DONE IN REVERSE ORDER.
C IFLAG IS RETURNED NEGATIVE IF EITHER K OR T IS OUT OF BOUNDS.
C-----------------------------------------------------------------------
      iflag = 0
      IF (k .LT. 0 .OR. k .GT. nq) GO TO 80
      tp = tn - hu -  100.0d0*uround*(tn + hu)
      IF ((t-tp)*(t-tn) .GT. 0.0d0) GO TO 90
C
      s = (t - tn)/h
      ic = 1
      IF (k .EQ. 0) GO TO 15
      jj1 = l - k
      DO 10 jj = jj1,nq
 10     ic = ic*jj
 15   c = dble(ic)
      DO 20 i = 1,n
 20     dky(i) = c*yh(i,l)
      IF (k .EQ. nq) GO TO 55
      jb2 = nq - k
      DO 50 jb = 1,jb2
        j = nq - jb
        jp1 = j + 1
        ic = 1
        IF (k .EQ. 0) GO TO 35
        jj1 = jp1 - k
        DO 30 jj = jj1,j
 30       ic = ic*jj
 35     c = dble(ic)
        DO 40 i = 1,n
 40       dky(i) = c*yh(i,jp1) + s*dky(i)
 50     CONTINUE
      IF (k .EQ. 0) RETURN
 55   r = h**(-k)
      DO 60 i = 1,n
 60     dky(i) = r*dky(i)
      RETURN
C
 80   CALL xerrwd('INTDY--  K (=I1) ILLEGAL      ',
     1   30, 51, 0, 1, k, 0, 0, 0.0d0, 0.0d0)
      iflag = -1
      RETURN
 90   CALL xerrwd('INTDY--  T (=R1) ILLEGAL      ',
     1   30, 52, 0, 0, 0, 0, 1, t, 0.0d0)
      CALL xerrwd(
     1  '      T NOT IN INTERVAL TCUR - HU (= R1) TO TCUR (=R2)      ',
     1   60, 52, 0, 0, 0, 0, 2, tp, tn)
      iflag = -2
      RETURN
C----------------------- END OF SUBROUTINE INTDY -----------------------
      END