dstoka Subroutine

subroutine dstoka(Neq,Y,Yh,Nyh,Yh1,Ewt,Savf,Savx,Acor,Wm,Iwm,f,jac,psol)

integer, dimension(*) :: Neq
real(kind=dp), dimension(*) :: Y
integer, intent(in) :: Nyh
real(kind=dp), intent(inout), dimension(Nyh,*) :: Yh
real(kind=dp), intent(inout), dimension(*) :: Yh1
real(kind=dp), dimension(*) :: Ewt
real(kind=dp), intent(inout), dimension(*) :: Savf
real(kind=dp), intent(inout), dimension(*) :: Savx
real(kind=dp), intent(inout), dimension(*) :: Acor
real(kind=dp), dimension(*) :: Wm
integer, dimension(*) :: Iwm
external f
external jac
external psol

real(kind=dp) :: dcon, ddn, del, delp, dfnrm, drc, dsm, dup, exdn, exsm, exup, r, rh, rhdn, rhsm, rhup, roc, stiff, told
integer :: i, i1, iredo, iret, j, jb, jok, m, ncf, newq, nslow

dls1%kflag = 0
told = dls1%tn
ncf = 0
dls1%ierpj = 0
dls1%iersl = 0
dls1%jcur = 0
dls1%icf = 0
delp = 0.0D0
if ( dls1%jstart>0 ) goto 400
if ( dls1%jstart==-1 ) then
!-----------------------------------------------------------------------
!  The following block handles preliminaries needed when JSTART = -1.
!  IPUP is set to MITER to force a matrix update.
!  If an order increase is about to be considered (IALTH = 1),
!  IALTH is reset to 2 to postpone consideration one more step.
!  If the caller has changed METH, DCFODE is called to reset
!  the coefficients of the method.
!  If the caller has changed MAXORD to a value less than the current
!  order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
!  If H is to be changed, YH must be rescaled.
!  If H or METH is being changed, IALTH is reset to L = NQ + 1
!  to prevent further changes in H for that many steps.
!-----------------------------------------------------------------------
   dls1%ipup = dls1%miter
   dls1%lmax = dls1%maxord + 1
   if ( dls1%ialth==1 ) dls1%ialth = 2
   if ( dls1%meth/=dls1%meo ) then
      call dcfode(dls1%meth,dls1%elco,dls1%tesco)
      dls1%meo = dls1%meth
      if ( dls1%nq<=dls1%maxord ) then
         dls1%ialth = dls1%l
         iret = 1
         goto 100
      endif
   elseif ( dls1%nq<=dls1%maxord ) then
      goto 200
   endif
   dls1%nq = dls1%maxord
   dls1%l = dls1%lmax
   do i = 1, dls1%l
      dls1%el(i) = dls1%elco(i,dls1%nq)
   enddo
   dls1%nqnyh = dls1%nq*Nyh
   dls1%rc = dls1%rc*dls1%el(1)/dls1%el0
   dls1%el0 = dls1%el(1)
   dls1%conit = 0.5D0/(dls1%nq+2)
   dlpk%epcon = dls1%conit*dls1%tesco(2,dls1%nq)
   ddn = dvnorm(dls1%n,Savf,Ewt)/dls1%tesco(1,dls1%l)
   exdn = 1.0D0/dls1%l
   rhdn = 1.0D0/(1.3D0*ddn**exdn+0.0000013D0)
   rh = min(rhdn,1.0D0)
   iredo = 3
   if ( dls1%h==dls1%hold ) then
      rh = max(rh,dls1%hmin/abs(dls1%h))
   else
      rh = min(rh,abs(dls1%h/dls1%hold))
      dls1%h = dls1%hold
   endif
   goto 300
else
   if ( dls1%jstart==-2 ) goto 200
!-----------------------------------------------------------------------
!  On the first call, the order is set to 1, and other variables are
!  initialized.  RMAX is the maximum ratio by which H can be increased
!  in a single step.  It is initially 1.E4 to compensate for the small
!  initial H, but then is normally equal to 10.  If a failure
!  occurs (in corrector convergence or error test), RMAX is set at 2
!  for the next increase.
!-----------------------------------------------------------------------
   dls1%lmax = dls1%maxord + 1
   dls1%nq = 1
   dls1%l = 2
   dls1%ialth = 2
   dls1%rmax = 10000.0D0
   dls1%rc = 0.0D0
   dls1%el0 = 1.0D0
   dls1%crate = 0.7D0
   dls1%hold = dls1%h
   dls1%meo = dls1%meth
   dls1%nslp = 0
   dls%nslj = 0
   dls1%ipup = 0
   iret = 3
   dls%newt = 0
   dls%stifr = 0.0D0
!-----------------------------------------------------------------------
!  DCFODE is called to get all the integration coefficients for the
!  current METH.  Then the EL vector and related constants are reset
!  whenever the order NQ is changed, or at the start of the problem.
!-----------------------------------------------------------------------
   call dcfode(dls1%meth,dls1%elco,dls1%tesco)
endif
 100  continue
do i = 1, dls1%l
   dls1%el(i) = dls1%elco(i,dls1%nq)
enddo
dls1%nqnyh = dls1%nq*Nyh
dls1%rc = dls1%rc*dls1%el(1)/dls1%el0
dls1%el0 = dls1%el(1)
dls1%conit = 0.5D0/(dls1%nq+2)
dlpk%epcon = dls1%conit*dls1%tesco(2,dls1%nq)
select case (iret)
case (2)
   rh = max(rh,dls1%hmin/abs(dls1%h))
   goto 300
case (3)
   goto 400
case default
endselect
!-----------------------------------------------------------------------
!  If H is being changed, the H ratio RH is checked against
!  RMAX, HMIN, and HMXI, and the YH array rescaled.  IALTH is set to
!  L = NQ + 1 to prevent a change of H for that many steps, unless
!  forced by a convergence or error test failure.
!-----------------------------------------------------------------------
 200  continue
if ( dls1%h==dls1%hold ) goto 400
rh = dls1%h/dls1%hold
dls1%h = dls1%hold
iredo = 3
 300  continue
rh = min(rh,dls1%rmax)
rh = rh/max(1.0D0,abs(dls1%h)*dls1%hmxi*rh)
r = 1.0D0
do j = 2, dls1%l
   r = r*rh
   do i = 1, dls1%n
      Yh(i,j) = Yh(i,j)*r
   enddo
enddo
dls1%h = dls1%h*rh
dls1%rc = dls1%rc*rh
dls1%ialth = dls1%l
if ( iredo==0 ) then
   dls1%rmax = 10.0D0
   goto 1200
endif
!-----------------------------------------------------------------------
!  This section computes the predicted values by effectively
!  multiplying the YH array by the Pascal triangle matrix.
!  The flag IPUP is set according to whether matrix data is involved
!  (NEWT .gt. 0 .and. JACFLG .ne. 0) or not, to trigger a call to DSETPK.
!  IPUP is set to MITER when RC differs from 1 by more than CCMAX,
!  and at least every MSBP steps, when JACFLG = 1.
!  RC is the ratio of new to old values of the coefficient  H*EL(1).
!-----------------------------------------------------------------------
 400  continue
if ( dls%newt==0 .or. dlpk%jacflg==0 ) then
   drc = 0.0D0
   dls1%ipup = 0
   dls1%crate = 0.7D0
else
   drc = abs(dls1%rc-1.0D0)
   if ( drc>dls1%ccmax ) dls1%ipup = dls1%miter
   if ( dls1%nst>=dls1%nslp+dls1%msbp ) dls1%ipup = dls1%miter
endif
dls1%tn = dls1%tn + dls1%h
i1 = dls1%nqnyh + 1
do jb = 1, dls1%nq
   i1 = i1 - Nyh
! DIR$ IVDEP
   do i = i1, dls1%nqnyh
      Yh1(i) = Yh1(i) + Yh1(i+Nyh)
   enddo
enddo
!-----------------------------------------------------------------------
!  Up to MAXCOR corrector iterations are taken.  A convergence test is
!  made on the RMS-norm of each correction, weighted by the error
!  weight vector EWT.  The sum of the corrections is accumulated in the
!  vector ACOR(i).  The YH array is not altered in the corrector loop.
!  Within the corrector loop, an estimated rate of convergence (ROC)
!  and a stiffness ratio estimate (STIFF) are kept.  Corresponding
!  global estimates are kept as CRATE and dls%stifr.
!-----------------------------------------------------------------------
 500  continue
m = 0
dlpk%mnewt = 0
stiff = 0.0D0
roc = 0.05D0
nslow = 0
do i = 1, dls1%n
   Y(i) = Yh(i,1)
enddo
call f(Neq,dls1%tn,Y,Savf)
dls1%nfe = dls1%nfe + 1
if ( dls%newt/=0 .and. dls1%ipup>0 ) then
!-----------------------------------------------------------------------
!  If indicated, DSETPK is called to update any matrix data needed,
!  before starting the corrector iteration.
!  JOK is set to indicate if the matrix data need not be recomputed.
!  IPUP is set to 0 as an indicator that the matrix data is up to date.
!-----------------------------------------------------------------------
   jok = 1
   if ( dls1%nst==0 .or. dls1%nst>dls%nslj+50 ) jok = -1
   if ( dls1%icf==1 .and. drc<0.2D0 ) jok = -1
   if ( dls1%icf==2 ) jok = -1
   if ( jok==-1 ) then
      dls%nslj = dls1%nst
      dls%njev = dls%njev + 1
   endif
   call dsetpk(Neq,Y,Yh1,Ewt,Acor,Savf,jok,Wm,Iwm,f,jac)
   dls1%ipup = 0
   dls1%rc = 1.0D0
   drc = 0.0D0
   dls1%nslp = dls1%nst
   dls1%crate = 0.7D0
   if ( dls1%ierpj/=0 ) goto 800
endif
do i = 1, dls1%n
   Acor(i) = 0.0D0
enddo
 600  continue
if ( dls%newt/=0 ) then
!-----------------------------------------------------------------------
!  In the case of the chord method, compute the corrector error,
!  and solve the linear system with that as right-hand side and
!  P as coefficient matrix.  STIFF is set to the ratio of the norms
!  of the residual and the correction vector.
!-----------------------------------------------------------------------
   do i = 1, dls1%n
      Savx(i) = dls1%h*Savf(i) - (Yh(i,2)+Acor(i))
   enddo
   dfnrm = dvnorm(dls1%n,Savx,Ewt)
   call dsolpk(Neq,Y,Savf,Savx,Ewt,Wm,Iwm,f,psol)
   if ( dls1%iersl<0 ) goto 800
   if ( dls1%iersl>0 ) goto 700
   del = dvnorm(dls1%n,Savx,Ewt)
   if ( del>1.0D-8 ) stiff = max(stiff,dfnrm/del)
   do i = 1, dls1%n
      Acor(i) = Acor(i) + Savx(i)
      Y(i) = Yh(i,1) + dls1%el(1)*Acor(i)
   enddo
else
!-----------------------------------------------------------------------
!  In the case of functional iteration, update Y directly from
!  the result of the last function evaluation, and STIFF is set to 1.0.
!-----------------------------------------------------------------------
   do i = 1, dls1%n
      Savf(i) = dls1%h*Savf(i) - Yh(i,2)
      Y(i) = Savf(i) - Acor(i)
   enddo
   del = dvnorm(dls1%n,Y,Ewt)
   do i = 1, dls1%n
      Y(i) = Yh(i,1) + dls1%el(1)*Savf(i)
      Acor(i) = Savf(i)
   enddo
   stiff = 1.0D0
endif
!-----------------------------------------------------------------------
!  Test for convergence.  If M .gt. 0, an estimate of the convergence
!  rate constant is made for the iteration switch, and is also used
!  in the convergence test.   If the iteration seems to be diverging or
!  converging at a slow rate (.gt. 0.8 more than once), it is stopped.
!-----------------------------------------------------------------------
if ( m/=0 ) then
   roc = max(0.05D0,del/delp)
   dls1%crate = max(0.2D0*dls1%crate,roc)
endif
dcon = del*min(1.0D0,1.5D0*dls1%crate)/dlpk%epcon
if ( dcon<=1.0D0 ) then
!-----------------------------------------------------------------------
!  The corrector has converged.  JCUR is set to 0 to signal that the
!  preconditioner involved may need updating later.
!  The stiffness ratio STIFR is updated using the latest STIFF value.
!  The local error test is made and control passes to statement 500
!  if it fails.
!-----------------------------------------------------------------------
   dls1%jcur = 0
   if ( dls%newt>0 ) dls%stifr = 0.5D0*(dls%stifr+stiff)
   if ( m==0 ) dsm = del/dls1%tesco(2,dls1%nq)
   if ( m>0 ) dsm = dvnorm(dls1%n,Acor,Ewt)/dls1%tesco(2,dls1%nq)
   if ( dsm>1.0D0 ) then
!-----------------------------------------------------------------------
!  The error test failed.  KFLAG keeps track of multiple failures.
!  Restore TN and the YH array to their previous values, and prepare
!  to try the step again.  Compute the optimum step size for this or
!  one lower order.  After 2 or more failures, H is forced to decrease
!  by a factor of 0.2 or less.
!-----------------------------------------------------------------------
      dls1%kflag = dls1%kflag - 1
      dls1%tn = told
      i1 = dls1%nqnyh + 1
      do jb = 1, dls1%nq
         i1 = i1 - Nyh
! DIR$ IVDEP
         do i = i1, dls1%nqnyh
            Yh1(i) = Yh1(i) - Yh1(i+Nyh)
         enddo
      enddo
      dls1%rmax = 2.0D0
      if ( abs(dls1%h)<=dls1%hmin*1.00001D0 ) then
!-----------------------------------------------------------------------
!  All returns are made through this section.  H is saved in HOLD
!  to allow the caller to change H on the next step.
!-----------------------------------------------------------------------
         dls1%kflag = -1
         goto 1300
      elseif ( dls1%kflag<=-3 ) then
!-----------------------------------------------------------------------
!  Control reaches this section if 3 or more failures have occured.
!  If 10 failures have occurred, exit with KFLAG = -1.
!  It is assumed that the derivatives that have accumulated in the
!  YH array have errors of the wrong order.  Hence the first
!  derivative is recomputed, and the order is set to 1.  Then
!  H is reduced by a factor of 10, and the step is retried,
!  until it succeeds or H reaches HMIN.
!-----------------------------------------------------------------------
         if ( dls1%kflag==-10 ) then
            dls1%kflag = -1
            goto 1300
         else
            rh = 0.1D0
            rh = max(dls1%hmin/abs(dls1%h),rh)
            dls1%h = dls1%h*rh
            do i = 1, dls1%n
               Y(i) = Yh(i,1)
            enddo
            call f(Neq,dls1%tn,Y,Savf)
            dls1%nfe = dls1%nfe + 1
            do i = 1, dls1%n
               Yh(i,2) = dls1%h*Savf(i)
            enddo
            dls1%ipup = dls1%miter
            dls1%ialth = 5
            if ( dls1%nq==1 ) goto 400
            dls1%nq = 1
            dls1%l = 2
            iret = 3
            goto 100
         endif
      else
         iredo = 2
         rhup = 0.0D0
         goto 900
      endif
   else
!-----------------------------------------------------------------------
!  After a successful step, update the YH array.
!  If Newton iteration is being done and STIFR is less than 1.5,
!  then switch to functional iteration.
!  Consider changing H if IALTH = 1.  Otherwise decrease IALTH by 1.
!  If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
!  use in a possible order increase on the next step.
!  If a change in H is considered, an increase or decrease in order
!  by one is considered also.  A change in H is made only if it is by a
!  factor of at least 1.1.  If not, IALTH is set to 3 to prevent
!  testing for that many steps.
!-----------------------------------------------------------------------
      dls1%kflag = 0
      iredo = 0
      dls1%nst = dls1%nst + 1
      if ( dls%newt==0 ) dls%nsfi = dls%nsfi + 1
      if ( dls%newt>0 .and. dls%stifr<1.5D0 ) dls%newt = 0
      dls1%hu = dls1%h
      dls1%nqu = dls1%nq
      do j = 1, dls1%l
         do i = 1, dls1%n
            Yh(i,j) = Yh(i,j) + dls1%el(j)*Acor(i)
         enddo
      enddo
      dls1%ialth = dls1%ialth - 1
      if ( dls1%ialth==0 ) then
!-----------------------------------------------------------------------
!  Regardless of the success or failure of the step, factors
!  RHDN, RHSM, and RHUP are computed, by which H could be multiplied
!  at order NQ - 1, order NQ, or order NQ + 1, respectively.
!  in the case of failure, RHUP = 0.0 to avoid an order increase.
!  the largest of these is determined and the new order chosen
!  accordingly.  If the order is to be increased, we compute one
!  additional scaled derivative.
!-----------------------------------------------------------------------
         rhup = 0.0D0
         if ( dls1%l/=dls1%lmax ) then
            do i = 1, dls1%n
               Savf(i) = Acor(i) - Yh(i,dls1%lmax)
            enddo
            dup = dvnorm(dls1%n,Savf,Ewt)/dls1%tesco(3,dls1%nq)
            exup = 1.0D0/(dls1%l+1)
            rhup = 1.0D0/(1.4D0*dup**exup+0.0000014D0)
         endif
         goto 900
      else
         if ( dls1%ialth<=1 ) then
            if ( dls1%l/=dls1%lmax ) then
               do i = 1, dls1%n
                  Yh(i,dls1%lmax) = Acor(i)
               enddo
            endif
         endif
         goto 1200
      endif
   endif
else
   m = m + 1
   if ( m/=dls1%maxcor ) then
      if ( m<2 .or. del<=2.0D0*delp ) then
         if ( roc<=10.0D0 ) then
            if ( roc>0.8D0 ) nslow = nslow + 1
            if ( nslow<2 ) then
               dlpk%mnewt = m
               delp = del
               call f(Neq,dls1%tn,Y,Savf)
               dls1%nfe = dls1%nfe + 1
               goto 600
            endif
         endif
      endif
   endif
endif
!-----------------------------------------------------------------------
!  The corrector iteration failed to converge.
!  If functional iteration is being done (NEWT = 0) and MITER .gt. 0
!  (and this is not the first step), then switch to Newton
!  (NEWT = MITER), and retry the step.  (Setting STIFR = 1023 insures
!  that a switch back will not occur for 10 step attempts.)
!  If Newton iteration is being done, but using a preconditioner that
!  is out of date (JACFLG .ne. 0 .and. JCUR = 0), then signal for a
!  re-evalutation of the preconditioner, and retry the step.
!  In all other cases, the YH array is retracted to its values
!  before prediction, and H is reduced, if possible.  If H cannot be
!  reduced or MXNCF failures have occurred, exit with KFLAG = -2.
!-----------------------------------------------------------------------
 700  continue
dls1%icf = 1
if ( dls%newt==0 ) then
   if ( dls1%nst==0 ) goto 800
   if ( dls1%miter==0 ) goto 800
   dls%newt = dls1%miter
   dls%stifr = 1023.0D0
   dls1%ipup = dls1%miter
   goto 500
endif
if ( dls1%jcur/=1 .and. dlpk%jacflg/=0 ) then
   dls1%ipup = dls1%miter
   goto 500
endif
 800  continue
dls1%icf = 2
ncf = ncf + 1
dlpk%ncfn = dlpk%ncfn + 1
dls1%rmax = 2.0D0
dls1%tn = told
i1 = dls1%nqnyh + 1
do jb = 1, dls1%nq
   i1 = i1 - Nyh
! DIR$ IVDEP
   do i = i1, dls1%nqnyh
      Yh1(i) = Yh1(i) - Yh1(i+Nyh)
   enddo
enddo
if ( dls1%ierpj<0 .or. dls1%iersl<0 ) then
   dls1%kflag = -3
   goto 1300
elseif ( abs(dls1%h)<=dls1%hmin*1.00001D0 ) then
   dls1%kflag = -2
   goto 1300
elseif ( ncf==dls1%mxncf ) then
   dls1%kflag = -2
   goto 1300
else
   rh = 0.5D0
   dls1%ipup = dls1%miter
   iredo = 1
   rh = max(rh,dls1%hmin/abs(dls1%h))
   goto 300
endif
 900  continue
exsm = 1.0D0/dls1%l
rhsm = 1.0D0/(1.2D0*dsm**exsm+0.0000012D0)
rhdn = 0.0D0
if ( dls1%nq/=1 ) then
   ddn = dvnorm(dls1%n,Yh(1,dls1%l),Ewt)/dls1%tesco(1,dls1%nq)
   exdn = 1.0D0/dls1%nq
   rhdn = 1.0D0/(1.3D0*ddn**exdn+0.0000013D0)
endif
if ( rhsm>=rhup ) then
   if ( rhsm>=rhdn ) then
      newq = dls1%nq
      rh = rhsm
      goto 1000
   endif
elseif ( rhup>rhdn ) then
   newq = dls1%l
   rh = rhup
   if ( rh<1.1D0 ) then
      dls1%ialth = 3
      goto 1200
   else
      r = dls1%el(dls1%l)/dls1%l
      do i = 1, dls1%n
         Yh(i,newq+1) = Acor(i)*r
      enddo
      goto 1100
   endif
endif
newq = dls1%nq - 1
rh = rhdn
if ( dls1%kflag<0 .and. rh>1.0D0 ) rh = 1.0D0
 1000 continue
if ( (dls1%kflag==0) .and. (rh<1.1D0) ) then
   dls1%ialth = 3
   goto 1200
else
   if ( dls1%kflag<=-2 ) rh = min(rh,0.2D0)
!-----------------------------------------------------------------------
!  If there is a change of order, reset NQ, L, and the coefficients.
!  In any case H is reset according to RH and the YH array is rescaled.
!  Then exit from 690 if the step was OK, or redo the step otherwise.
!-----------------------------------------------------------------------
   if ( newq==dls1%nq ) then
      rh = max(rh,dls1%hmin/abs(dls1%h))
      goto 300
   endif
endif
 1100 continue
dls1%nq = newq
dls1%l = dls1%nq + 1
iret = 2
goto 100
 1200 continue
r = 1.0D0/dls1%tesco(2,dls1%nqu)
do i = 1, dls1%n
   Acor(i) = Acor(i)*r
enddo
 1300 continue
dls1%hold = dls1%h
dls1%jstart = 1
end subroutine dstoka