share/odepack/fortran/dlhin.f

   1 *DECK DLHIN
   2       SUBROUTINE DLHIN (NEQ, N, T0, Y0, YDOT, F, TOUT, UROUND,
   3      1   EWT, ITOL, ATOL, Y, TEMP, H0, NITER, IER)
   4       EXTERNAL F
   5       DOUBLE PRECISION T0, Y0, YDOT, TOUT, UROUND, EWT, ATOL, Y,
   6      1   TEMP, H0
   7       INTEGER NEQ, N, ITOL, NITER, IER
   8       DIMENSION NEQ(*), Y0(*), YDOT(*), EWT(*), ATOL(*), Y(*), TEMP(*)
   9 C-----------------------------------------------------------------------
  10 C Call sequence input -- NEQ, N, T0, Y0, YDOT, F, TOUT, UROUND,
  11 C                        EWT, ITOL, ATOL, Y, TEMP
  12 C Call sequence output -- H0, NITER, IER
  13 C Common block variables accessed -- None
  14 C
  15 C Subroutines called by DLHIN: F, DCOPY
  16 C Function routines called by DLHIN: DVNORM
  17 C-----------------------------------------------------------------------
  18 C This routine computes the step size, H0, to be attempted on the
  19 C first step, when the user has not supplied a value for this.
  20 C
  21 C First we check that TOUT - T0 differs significantly from zero.  Then
  22 C an iteration is done to approximate the initial second derivative
  23 C and this is used to define H from WRMS-norm(H**2 * yddot / 2) = 1.
  24 C A bias factor of 1/2 is applied to the resulting h.
  25 C The sign of H0 is inferred from the initial values of TOUT and T0.
  26 C
  27 C Communication with DLHIN is done with the following variables:
  28 C
  29 C NEQ    = NEQ array of solver, passed to F.
  30 C N      = size of ODE system, input.
  31 C T0     = initial value of independent variable, input.
  32 C Y0     = vector of initial conditions, input.
  33 C YDOT   = vector of initial first derivatives, input.
  34 C F      = name of subroutine for right-hand side f(t,y), input.
  35 C TOUT   = first output value of independent variable
  36 C UROUND = machine unit roundoff
  37 C EWT, ITOL, ATOL = error weights and tolerance parameters
  38 C                   as described in the driver routine, input.
  39 C Y, TEMP = work arrays of length N.
  40 C H0     = step size to be attempted, output.
  41 C NITER  = number of iterations (and of f evaluations) to compute H0,
  42 C          output.
  43 C IER    = the error flag, returned with the value
  44 C          IER = 0  if no trouble occurred, or
  45 C          IER = -1 if TOUT and t0 are considered too close to proceed.
  46 C-----------------------------------------------------------------------
  47 C
  48 C Type declarations for local variables --------------------------------
  49 C
  50       DOUBLE PRECISION AFI, ATOLI, DELYI, HALF, HG, HLB, HNEW, HRAT,
  51      1     HUB, HUN, PT1, T1, TDIST, TROUND, TWO, DVNORM, YDDNRM
  52       INTEGER I, ITER
  53 C-----------------------------------------------------------------------
  54 C The following Fortran-77 declaration is to cause the values of the
  55 C listed (local) variables to be saved between calls to this integrator.
  56 C-----------------------------------------------------------------------
  57       SAVE HALF, HUN, PT1, TWO
  58       DATA HALF /0.5D0/, HUN /100.0D0/, PT1 /0.1D0/, TWO /2.0D0/
  59 C
  60       NITER = 0
  61       TDIST = ABS(TOUT - T0)
  62       TROUND = UROUND*MAX(ABS(T0),ABS(TOUT))
  63       IF (TDIST .LT. TWO*TROUND) GO TO 100
  64 C
  65 C Set a lower bound on H based on the roundoff level in T0 and TOUT. ---
  66       HLB = HUN*TROUND
  67 C Set an upper bound on H based on TOUT-T0 and the initial Y and YDOT. -
  68       HUB = PT1*TDIST
  69       ATOLI = ATOL(1)
  70       DO 10 I = 1,N
  71         IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
  72         DELYI = PT1*ABS(Y0(I)) + ATOLI
  73         AFI = ABS(YDOT(I))
  74         IF (AFI*HUB .GT. DELYI) HUB = DELYI/AFI
  75  10     CONTINUE
  76 C
  77 C Set initial guess for H as geometric mean of upper and lower bounds. -
  78       ITER = 0
  79       HG = SQRT(HLB*HUB)
  80 C If the bounds have crossed, exit with the mean value. ----------------
  81       IF (HUB .LT. HLB) THEN
  82         H0 = HG
  83         GO TO 90
  84       ENDIF
  85 C
  86 C Looping point for iteration. -----------------------------------------
  87  50   CONTINUE
  88 C Estimate the second derivative as a difference quotient in f. --------
  89       T1 = T0 + HG
  90       DO 60 I = 1,N
  91  60     Y(I) = Y0(I) + HG*YDOT(I)
  92       CALL F (NEQ, T1, Y, TEMP)
  93       DO 70 I = 1,N
  94  70     TEMP(I) = (TEMP(I) - YDOT(I))/HG
  95       YDDNRM = DVNORM (N, TEMP, EWT)
  96 C Get the corresponding new value of H. --------------------------------
  97       IF (YDDNRM*HUB*HUB .GT. TWO) THEN
  98         HNEW = SQRT(TWO/YDDNRM)
  99       ELSE
 100         HNEW = SQRT(HG*HUB)
 101       ENDIF
 102       ITER = ITER + 1
 103 C-----------------------------------------------------------------------
 104 C Test the stopping conditions.
 105 C Stop if the new and previous H values differ by a factor of .lt. 2.
 106 C Stop if four iterations have been done.  Also, stop with previous H
 107 C if hnew/hg .gt. 2 after first iteration, as this probably means that
 108 C the second derivative value is bad because of cancellation error.
 109 C-----------------------------------------------------------------------
 110       IF (ITER .GE. 4) GO TO 80
 111       HRAT = HNEW/HG
 112       IF ( (HRAT .GT. HALF) .AND. (HRAT .LT. TWO) ) GO TO 80
 113       IF ( (ITER .GE. 2) .AND. (HNEW .GT. TWO*HG) ) THEN
 114         HNEW = HG
 115         GO TO 80
 116       ENDIF
 117       HG = HNEW
 118       GO TO 50
 119 C
 120 C Iteration done.  Apply bounds, bias factor, and sign. ----------------
 121  80   H0 = HNEW*HALF
 122       IF (H0 .LT. HLB) H0 = HLB
 123       IF (H0 .GT. HUB) H0 = HUB
 124  90   H0 = SIGN(H0, TOUT - T0)
 125 C Restore Y array from Y0, then exit. ----------------------------------
 126       CALL DCOPY (N, Y0, 1, Y, 1)
 127       NITER = ITER
 128       IER = 0
 129       RETURN
 130 C Error return for TOUT - T0 too small. --------------------------------
 131  100  IER = -1
 132       RETURN
 133 C----------------------- End of Subroutine DLHIN -----------------------
 134       END