initialize_numdiff_for_diff Subroutine

    subroutine initialize_numdiff_for_diff(me,n,m,xlow,xhigh,&
                                    problem_func,sparsity_mode,info,&
                                    chunk_size,eps,acc,cache_size,&
                                    xlow_for_sparsity,xhigh_for_sparsity,&
                                    dpert_for_sparsity,sparsity_perturb_mode,&
                                    linear_sparsity_tol,function_precision_tol,&
                                    num_sparsity_points,print_messages)

    implicit none

    class(numdiff_type),intent(inout):: me
    integer,intent(in)               :: n             !! number of `x` variables
    integer,intent(in)               :: m             !! number of `f` functions
    real(wp),dimension(n),intent(in) :: xlow          !! lower bounds on `x`
    real(wp),dimension(n),intent(in) :: xhigh         !! upper bounds on `x`
    procedure(func)                  :: problem_func  !! the user function that defines the problem
                                                      !! (returns `m` functions)
    integer,intent(in)               :: sparsity_mode !! the sparsity computation method:
                                                      !!
                                                      !! **1** - assume dense,
                                                      !! **2** - three-point simple method,
                                                      !! **3** - will be specified by the user in
                                                      !! a subsequent call to [[set_sparsity_pattern]].
                                                      !! **4** - computes a two-point jacobian
                                                      !! at `num_sparsity_points` points.
    procedure(info_f),optional       :: info          !! a function the user can define
                                                      !! which is called when each column
                                                      !! of the jacobian is computed.
                                                      !! It can be used to perform any
                                                      !! setup operations.
    integer,intent(in),optional      :: chunk_size    !! chunk size for allocating the arrays
                                                      !! (must be >0) [default is 100]
    real(wp),intent(in),optional     :: eps           !! tolerance parameter for [[diff]]
                                                      !! if not present, default is `1.0e-9_wp`
    real(wp),intent(in),optional     :: acc           !! tolerance parameter for [[diff]]
                                                      !! if not present, default is `0.0_wp`
    integer,intent(in),optional      :: cache_size    !! if present, this is the cache size
                                                      !! for the function cache
                                                      !! (default is not to enable cache)
    real(wp),dimension(n),intent(in),optional :: xlow_for_sparsity   !! lower bounds on `x` used for
                                                                     !! sparsity computation (when
                                                                     !! `sparsity_mode` is 2). If not
                                                                     !! present, then `xlow` is used.
    real(wp),dimension(n),intent(in),optional :: xhigh_for_sparsity  !! upper bounds on `x` used for
                                                                     !! sparsity computation (when
                                                                     !! `sparsity_mode` is 2). If not
                                                                     !! present, then `xhigh` is used.
    real(wp),dimension(n),intent(in),optional :: dpert_for_sparsity  !! required if `sparsity_mode=4`
    integer,intent(in),optional :: sparsity_perturb_mode   !! perturbation mode (required if `sparsity_mode=4`):
                                                           !!
                                                           !! **1** - perturbation is `dx=dpert`,
                                                           !! **2** - perturbation is `dx=dpert*x`,
                                                           !! **3** - perturbation is `dx=dpert*(1+x)`
    integer,intent(in),optional :: num_sparsity_points  !! for `sparsity_mode=4`, the number of jacobian
                                                        !! evaluations used to estimate the sparsity pattern.
    real(wp),intent(in),optional :: linear_sparsity_tol !! the equality tolerance for derivatives to
                                                        !! indicate a constant jacobian element (linear sparsity)
    real(wp),intent(in),optional :: function_precision_tol  !! the function precision. two functions values
                                                            !! that are the within this tolerance are
                                                            !! considered the same value. This is used
                                                            !! when estimating the sparsity pattern when
                                                            !! `sparsity_mode=2` in [[compute_sparsity_random]]
    logical,intent(in),optional :: print_messages !! if true, print error messages to `error_unit`.
                                                  !! default is True.

    logical :: cache  !! if the cache is to be used

    me%use_diff = .true.

    ! cache:
    cache = present(cache_size)
    if (cache) cache = cache_size>0
    if (cache) then
        me%compute_function  => compute_function_with_cache
        call me%cache%initialize(cache_size,n,m)
    else
        me%compute_function  => problem_func
        call me%cache%destroy()
    end if

    ! functions:
    me%problem_func      => problem_func
    me%jacobian_function => compute_jacobian_with_diff

    ! size of the problem:
    me%n = n
    me%m = m

    ! input variable bounds:
    call me%set_numdiff_bounds(xlow,xhigh)

    ! set the sparsity function
    call me%set_sparsity_mode(sparsity_mode,xlow_for_sparsity,xhigh_for_sparsity)

    if (sparsity_mode==4) then
        ! these must be present, since we don't have a dpert and perturb mode for the gradients:
        if (present(dpert_for_sparsity) .and. present(sparsity_perturb_mode)) then
            me%dpert_for_sparsity = abs(dpert_for_sparsity)
            ! perturbation options:
            select case (sparsity_perturb_mode)
            case(1:3)
                me%sparsity_perturb_mode = sparsity_perturb_mode
            case default
                call me%raise_exception(1,'initialize_numdiff_for_diff',&
                                          'sparsity_perturb_mode must be 1, 2, or 3.')
                return
            end select
        else
            call me%raise_exception(2,'initialize_numdiff_for_diff',&
                                      'missing required inputs for sparsity mode 4')
            return
        end if
    end if
    ! if these aren't present, they will just keep the defaults:
    if (present(linear_sparsity_tol))    me%linear_sparsity_tol    = linear_sparsity_tol
    if (present(function_precision_tol)) me%function_precision_tol = function_precision_tol
    if (present(num_sparsity_points))    me%num_sparsity_points    = num_sparsity_points

    ! optional:
    if (present(chunk_size))     me%chunk_size = abs(chunk_size)
    if (present(eps))            me%eps = eps
    if (present(acc))            me%acc = acc
    if (present(info))           me%info_function => info
    if (present(print_messages)) me%print_messages = print_messages

    end subroutine initialize_numdiff_for_diff