initialize_numdiff Subroutine

    subroutine initialize_numdiff(me,n,m,xlow,xhigh,perturb_mode,dpert,&
                        problem_func,sparsity_mode,jacobian_method,jacobian_methods,&
                        class,classes,info,chunk_size,partition_sparsity_pattern,&
                        cache_size,xlow_for_sparsity,xhigh_for_sparsity,&
                        dpert_for_sparsity,sparsity_perturb_mode,&
                        linear_sparsity_tol,function_precision_tol,&
                        num_sparsity_points)

    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`
    integer,intent(in)                  :: perturb_mode    !! perturbation mode:
                                                           !! **1** - perturbation is `dx=dpert`,
                                                           !! **2** - perturbation is `dx=dpert*x`,
                                                           !! **3** - perturbation is `dx=dpert*(1+x)`
    real(wp),dimension(n),intent(in)    :: dpert           !! perturbation vector for `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.
    integer,intent(in),optional         :: jacobian_method !! `id` code for the finite difference method
                                                           !! to use for all `n` variables.
                                                           !! see [[get_finite_difference_method]]
    integer,dimension(n),intent(in),optional :: jacobian_methods !! `id` codes for the finite difference method
                                                                 !! to use for each variable.
                                                                 !! see [[get_finite_difference_method]]
    integer,intent(in),optional :: class  !! method class for the finite difference method
                                          !! to use for all `n` variables.
                                          !! see [[get_finite_difference_method]]
    integer,dimension(n),intent(in),optional :: classes   !! method class for the finite difference methods
                                                          !! to use for each variable.
                                                          !! see [[get_finite_difference_method]]
    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]
    logical,intent(in),optional :: partition_sparsity_pattern  !! if the sparisty pattern is to
                                                               !! be partitioned using [[DSM]]
                                                               !! [default is False]
    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  !! for `sparsity_mode=4`, the perturbation
    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)`
    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]]
    integer,intent(in),optional :: num_sparsity_points  !! for `sparsity_mode=4`, the number of jacobian
                                                        !! evaluations used to estimate the sparsity pattern.

    integer :: i      !! counter
    logical :: found  !! flag for [[get_finite_difference_method]]
    logical :: cache  !! if the cache is to be used

    me%use_diff = .false.

    ! 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

    if (present(partition_sparsity_pattern)) then
        me%partition_sparsity_pattern = partition_sparsity_pattern
    else
        me%partition_sparsity_pattern = .false.
    end if

    ! method:
    if (allocated(me%meth)) deallocate(me%meth)
    if (allocated(me%class)) deallocate(me%class)
    if (allocated(me%class_meths)) deallocate(me%class_meths)

    if (      present(jacobian_method) .and. .not. present(jacobian_methods) .and. &
        .not. present(class) .and. .not. present(classes)) then
        ! use the same for all variable
        me%mode = 1
        allocate(me%meth(n))
        do i=1,n
            call get_finite_difference_method(jacobian_method,me%meth(i),found)
            if (.not. found) then
                call me%raise_exception(8,'initialize_numdiff',&
                                          'invalid jacobian_method')
                return
            end if
        end do
    elseif (.not. present(jacobian_method) .and. present(jacobian_methods) .and. &
            .not. present(class) .and. .not. present(classes)) then
        ! specify a separate method for each variable
        me%mode = 1
        allocate(me%meth(n))
        do i=1,n
            call get_finite_difference_method(jacobian_methods(i),me%meth(i),found)
            if (.not. found) then
                call me%raise_exception(9,'initialize_numdiff',&
                                          'invalid jacobian_methods')
                return
            end if
        end do
        if (me%partition_sparsity_pattern) then
            call me%raise_exception(10,'initialize_numdiff',&
                                        'when using partitioned sparsity pattern, '//&
                                        'all columns must use the same '//&
                                        'finite diff method.')
            return
        end if
    elseif (.not. present(jacobian_method) .and. .not. present(jacobian_methods) .and. &
                  present(class) .and. .not. present(classes)) then
        ! use the class for all variables
        me%mode = 2
        allocate(me%class(n))
        me%class = class
        allocate(me%class_meths(n))
        me%class_meths(1) = get_all_methods_in_class(class)
        if (n>1) me%class_meths(2:n) = me%class_meths(1)  ! just copy them over
    elseif (.not. present(jacobian_method) .and. .not. present(jacobian_methods) .and. &
            .not. present(class) .and. present(classes)) then
        ! specify a separate class for each variable
        me%mode = 2
        me%class = classes
        allocate(me%class_meths(n))
        do i=1,n
            me%class_meths(i) = get_all_methods_in_class(me%class(i))
        end do
        if (me%partition_sparsity_pattern) then
            call me%raise_exception(11,'initialize_numdiff',&
                                        'when using partitioned sparsity pattern, '//&
                                        'all columns must use the same '//&
                                        'finite diff method.')
        end if
    else
        call me%raise_exception(12,'initialize_numdiff',&
                                    'must specify one of either '//&
                                    'jacobian_method, jacobian_methods, class, or classes.')
    end if

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

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

    ! perturbation options:
    select case (perturb_mode)
    case(1:3)
        me%perturb_mode = perturb_mode
    case default
        call me%raise_exception(13,'initialize_numdiff',&
                                   'perturb_mode must be 1, 2, or 3.')
        return
    end select
    if (allocated(me%dpert)) deallocate(me%dpert)
    allocate(me%dpert(n))
    me%dpert = abs(dpert)

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

    ! set the jacobian function, depending on the options:
    if (me%partition_sparsity_pattern) then
        me%jacobian_function => compute_jacobian_partitioned
    else
        me%jacobian_function => compute_jacobian_standard
    end if

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

    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

    if (present(dpert_for_sparsity)) then
        me%dpert_for_sparsity = abs(dpert_for_sparsity)
    else
        me%dpert_for_sparsity = dpert ! use the same dpert as the jacobian
    end if

    if (present(sparsity_perturb_mode)) then
        select case (sparsity_perturb_mode)
        case(1:3)
            me%sparsity_perturb_mode = sparsity_perturb_mode
        case default
            call me%raise_exception(14,'initialize_numdiff',&
                                       'sparsity_perturb_mode must be 1, 2, or 3.')
            return
        end select
    else
        me%sparsity_perturb_mode = perturb_mode ! use the same perturb mode as the jacobian
    end if

    end subroutine initialize_numdiff