ddeabm_initialize Subroutine

   subroutine ddeabm_initialize(me, neq, maxnum, df, rtol, atol, report, initial_step_mode, initial_step_size, &
                                max_step_size, max_order)

      implicit none

      class(ddeabm_class), intent(inout)   :: me
      integer, intent(in)                  :: neq     !! number of equations to be integrated
      integer, intent(in)                  :: maxnum  !! maximum number of integration steps allowed
      procedure(deriv_func)               :: df      !! derivative function dx/dt
      real(wp), dimension(:), intent(in)    :: rtol    !! relative tolerance for integration (see [[ddeabm]])
      real(wp), dimension(:), intent(in)    :: atol    !! absolution tolerance for integration (see [[ddeabm]])
      procedure(report_func), optional     :: report  !! reporting function
      integer, intent(in), optional         :: initial_step_mode !! how to choose the initial step `h`:
                                                             !!
                                                             !! 1. Use [[dhstrt]].
                                                             !! 2. Use the older (quicker) algorithm.
                                                             !! 3. Use the user-specified value `initial_step_size` (>0).
      real(wp), intent(in), optional        :: initial_step_size !! for `initial_step_mode=3`
      real(wp), intent(in), optional        :: max_step_size !! maximum step size magnitude (>0).
                                                             !! If not present, the step size is not limited.
      integer, intent(in), optional         :: max_order !! maximum order of the method
                                                         !! (`1` to `ddeabm_max_order_limit`). If not present,
                                                         !! `ddeabm_default_max_order` (12) is used.
                                                         !! A lower value can improve stability on some problems.
                                                         !! A higher value can reduce the number of steps for very
                                                         !! tight tolerances when using higher precision (`REAL128`).

      logical :: vector_tols
      integer :: i, j, k
      integer :: kmax !! maximum order of the method
      character(len=100) :: msg !! error message

      !initialize the class:
      call me%destroy()

      !maximum order of the method:
      kmax = ddeabm_default_max_order
      if (present(max_order)) then
         if (max_order < 1 .or. max_order > ddeabm_max_order_limit) then
            write (msg, '(A,I0,A,I0)') 'Error in ddeabm_initialize: max_order must be in the range 1 to ', &
               ddeabm_max_order_limit, '. It was ', max_order
            error stop trim(msg)
         end if
         kmax = max_order
      end if
      me%max_order = kmax

      !number of equations:
      me%neq = neq

      !maximum number of steps:
      me%maxnum = maxnum

      !set the derivative routine pointer:
      me%df => df

      !set the report function (optional):
      if (present(report)) then
         me%report => report
      else
         nullify (me%report)
      end if

      !allocate the work arrays:
      allocate (me%ypout(neq))
      allocate (me%yp(neq))
      allocate (me%yy(neq))
      allocate (me%wt(neq))
      allocate (me%p(neq))
      allocate (me%phi(neq, max(kmax + 4, 6)))  ! (dhstrt uses columns 3:6 as work arrays)
      allocate (me%alpha(kmax), me%beta(kmax), me%psi(kmax), me%v(kmax), me%w(kmax), me%gi(kmax))
      allocate (me%sig(kmax + 1), me%g(kmax + 1))
      allocate (me%iv(kmax))
      me%alpha = 0.0_wp
      me%beta = 0.0_wp
      me%psi = 0.0_wp
      me%v = 0.0_wp
      me%w = 0.0_wp
      me%gi = 0.0_wp
      me%sig = 0.0_wp
      me%g = 0.0_wp
      me%iv = 0
      ! so that interpolation before the first step is well defined:
      me%ypout = 0.0_wp
      me%yp = 0.0_wp
      me%yy = 0.0_wp
      me%wt = 0.0_wp
      me%p = 0.0_wp
      me%phi = 0.0_wp

      !tolerances:
      ! [for now, we are considering these unchangeable, although they do not have to be]
      ! each of rtol and atol can be size 1 (used for all equations) or size neq.

      if (size(rtol) /= 1 .and. size(rtol) /= neq) &
         error stop 'Error in ddeabm_initialize: rtol must be size 1 or neq'
      if (size(atol) /= 1 .and. size(atol) /= neq) &
         error stop 'Error in ddeabm_initialize: atol must be size 1 or neq'

      vector_tols = (size(rtol) == neq .or. size(atol) == neq) .and. neq > 1
      me%scalar_tols = .not. vector_tols

      if (me%scalar_tols) then
         allocate (me%rtol(1)); allocate (me%rtol_tmp(1))
         allocate (me%atol(1)); allocate (me%atol_tmp(1))
         me%rtol = rtol(1)
         me%atol = atol(1)
      else
         ! a size-1 tolerance is used for all the equations
         allocate (me%rtol(neq)); allocate (me%rtol_tmp(neq))
         allocate (me%atol(neq)); allocate (me%atol_tmp(neq))
         if (size(rtol) == 1) then
            me%rtol = rtol(1)
         else
            me%rtol = rtol
         end if
         if (size(atol) == 1) then
            me%atol = atol(1)
         else
            me%atol = atol
         end if
      end if

      !Compute the two and gstr arrays:
      !  Note: this is a modification of the original dsteps code.
      !  The full-precision coefficients are used here, instead
      !  of the less precise ones in the original.
      !  These are computed recursively using the equation on p. 159 of
      !  Shampine/Gordon, "Computer Solution of Ordinary Differential Equations", 1975.
      allocate (me%two(kmax + 1))
      allocate (me%gstr(0:kmax + 1))
      do i = 1, kmax + 1
         me%two(i) = 2.0_wp**i
      end do
      me%gstr = 0.0_wp
      me%gstr(0) = 1.0_wp
      do i = 1, kmax + 1
         k = 1
         do j = i - 1, 0, -1
            k = k + 1
            me%gstr(i) = me%gstr(i) + 1.0_wp/real(k, wp)*me%gstr(j)
         end do
         me%gstr(i) = -me%gstr(i)
      end do
      me%gstr = abs(me%gstr)

      ! initial step size (can have initial_step_mode 1, 2, or 3):
      me%initial_step_size = 0.0_wp ! dummy value (only used if initial_step_mode is 3)
      if (present(initial_step_mode)) then
         select case (initial_step_mode)
         case (1, 2, 3)
            me%initial_step_mode = initial_step_mode
            if (initial_step_mode == 3) then
               if (present(initial_step_size)) then
                  if (initial_step_size > 0.0_wp) then
                     me%initial_step_size = initial_step_size
                  else
                     error stop 'Error in ddeabm_initialize: initial_step_size must be > 0'
                  end if
               else
                  error stop 'Error in ddeabm_initialize: initial_step_size must be input when initial_step_mode=3'
               end if
            end if
         case default
            error stop 'Error in ddeabm_initialize: invalid input for initial_step_mode'
         end select
      else
         me%initial_step_mode = 1 ! original method (dhstrt) if none is specified
      end if

      ! maximum step size (optional):
      me%max_step_size = 0.0_wp ! no limit
      if (present(max_step_size)) then
         if (max_step_size > 0.0_wp) then
            me%max_step_size = max_step_size
         else
            error stop 'Error in ddeabm_initialize: max_step_size must be > 0'
         end if
      end if

   end subroutine ddeabm_initialize