dlinsd Subroutine

subroutine dlinsd( &
   neq, y, t, ydot, cj, tscale, p, pnorm, wt, &
   lsoff, stptol, iret, res, ires, rwm, iwm, &
   fnorm, icopt, idy, r, ynew, ydotnew, icnflg, &
   icnstr, rlx, rpar, ipar)
!! This routine uses a linesearch algorithm to calculate a new \( (y,\dot{y}) \) pair,
!! denoted \( (y_{new},\dot{y}_{new}) \), such that:
!!
!!   $$ f(y_{new},\dot{y}_{new}) \leq (1 - 2 \alpha \lambda) f(y,\dot{y})$$
!!
!! where \(0 < \lambda \le 1\), and \( f(y,\dot{y}) \) is defined as:
!!
!! $$ f(y,\dot{y}) = (1/2) \| (J^{-1} G(t,y,\dot{y}) \|^2 $$
!!
!! where the norm is the weighted root mean square vector norm, \(G\) is the DAE system
!! residual function, and \(J\) is the system iteration matrix (Jacobian).

   use daskr_kinds, only: rk, zero, one
   use daskr
   implicit none
   external :: xerrwd

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(inout) :: y(neq)
      !! On entry, the current solution vector. On return, the new solution vector.
   real(rk), intent(in) :: t
      !! Independent variable.
   real(rk), intent(inout) :: ydot(neq)
      !! On entry, the current derivative vector. On return, the new derivative vector.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: tscale ! @todo: what is "t"?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(inout) :: p(neq)
      !! Approximate Newton step used in backtracking.
   real(rk), intent(inout) :: pnorm
      !! Weighted root mean square norm of `p`.
   real(rk), intent(inout) :: wt(neq)
      !! Scaling factors.
   integer, intent(in) :: lsoff ! @todo: convert to logical
      !! Flag showing whether the linesearch algorithm is to be invoked.
      !! `0`: do the linesearch.
      !! `1`: turn off linesearch.
   real(rk), intent(in) :: stptol
      !! Tolerance used in calculating the minimum lambda (`rl`) value allowed.
   integer, intent(out) :: iret
      !! Return flag.
      !! `0`: a satisfactory (y,y') was found.
      !! `1`: the routine failed to find a new (y,y') pair that was sufficiently distinct
      !! from the current one.
      !! `2`: a failure in `res`.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   real(rk), intent(inout) :: rwm(*)
      !! Real workspace for the linear system solver.
   integer, intent(inout) :: iwm(*)
      !! Integer workspace for the linear system solver.
   real(rk), intent(inout) :: fnorm
      !! Value of \( \sqrt{2f} \) for the current \( (y,\dot{y}) \) on input and output.
   integer, intent(in) :: icopt
      !! Initial condition flag.
   integer, intent(in) :: idy(neq)
      !! Array indicating which variables are differential and which are algebraic.
   real(rk), intent(out) :: r(neq)
      !! Scaled residual vector \(J^{-1} G(t,y,\dot{y})\).
   real(rk), intent(out) :: ynew(neq)
      !! New solution vector.
   real(rk), intent(out) :: ydotnew(neq)
      !! New derivative of solution vector.
   integer, intent(in) :: icnflg
      !! Constraint flag. If nonzero, then constraint violations in the proposed new
      !! approximate solution will be checked for, and the maximum step length will be
      !! adjusted accordingly.
   integer, intent(out) :: icnstr(neq)
      !! Flags for checking constraints.
   real(rk), intent(in) :: rlx
      !! Factor for restricting update size in [[dcnstr]].
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   real(rk), parameter :: alpha = 1e-4_rk
   integer :: ivar, kprint
   real(rk) :: f1norm, f1normp, fnormp, ratio, ratio1, rl, rlmin, slpi, tau
   character(len=80) :: msg

   kprint = iwm(iwloc_kprint)
   f1norm = (fnorm**2)/2
   ratio = one

   if (kprint .ge. 2) then
      msg = '------ IN ROUTINE DLINSD-- PNRM = (R1)'
      call xerrwd(msg, 38, 901, 0, 0, 0, 0, 1, pnorm, zero)
   end if
   tau = pnorm
   rl = one

   ! Check for violations of the constraints, if any are imposed.
   ! If any violations are found, the step vector P is rescaled, and the
   ! constraint check is repeated, until no violations are found.
   if (icnflg .ne. 0) then
      do
         call dyypnw(neq, y, ydot, cj, rl, p, icopt, idy, ynew, ydotnew)
         call dcnstr(neq, y, ynew, icnstr, tau, rlx, iret, ivar)
         if (iret /= 1) exit
         ratio1 = tau/pnorm
         ratio = ratio*ratio1
         p = p*ratio1
         pnorm = tau
         if (kprint .ge. 2) then
            msg = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)'
            call xerrwd(msg, 50, 902, 0, 1, ivar, 0, 1, pnorm, zero)
         end if
         if (pnorm .le. stptol) then
            iret = 1
            return
         end if
      end do
   end if

   slpi = -2*f1norm*ratio
   rlmin = stptol/pnorm
   if ((lsoff .eq. 0) .and. (kprint .ge. 2)) then
      msg = '------ MIN. LAMBDA = (R1)'
      call xerrwd(msg, 25, 903, 0, 0, 0, 0, 1, rlmin, zero)
   end if

   ! Begin iteration to find RL value satisfying alpha-condition.
   ! If RL becomes less than RLMIN, then terminate with IRET = 1.
   do
      call dyypnw(neq, y, ydot, cj, rl, p, icopt, idy, ynew, ydotnew)
      call dfnrmd(neq, ynew, t, ydotnew, r, cj, tscale, wt, res, ires, &
                  fnormp, rwm, iwm, rpar, ipar)
      iwm(iwloc_nre) = iwm(iwloc_nre) + 1
      if (ires .ne. 0) then
         iret = 2
         return
      end if

      if (lsoff .eq. 1) goto 150

      f1normp = (fnormp**2)/2
      if (kprint .ge. 2) then
         msg = '------ LAMBDA = (R1)'
         call xerrwd(msg, 20, 904, 0, 0, 0, 0, 1, rl, zero)
         msg = '------ NORM(F1) = (R1),  NORM(F1NEW) = (R2)'
         call xerrwd(msg, 43, 905, 0, 0, 0, 0, 2, f1norm, f1normp)
      end if

      if (f1normp .gt. f1norm + alpha*slpi*rl) goto 200

      ! Alpha-condition is satisfied, or linesearch is turned off.
      ! Copy YNEW,YDOTNEW to Y,YDOT and return.
150   continue
      iret = 0
      y = ynew
      ydot = ydotnew
      fnorm = fnormp
      if (kprint .ge. 1) then
         msg = '------ LEAVING ROUTINE DLINSD, FNRM = (R1)'
         call xerrwd(msg, 42, 906, 0, 0, 0, 0, 1, fnorm, zero)
      end if
      return

      ! Alpha-condition not satisfied. Perform backtrack to compute new RL
      ! value. If no satisfactory YNEW,YDOTNEW can be found sufficiently
      ! distinct from Y,YDOT, then return IRET = 1.
200   continue
      if (rl .lt. rlmin) then
         iret = 1
         return
      end if

      rl = rl/2

   end do

end subroutine dlinsd