compute_jacobian_with_diff Subroutine

    subroutine compute_jacobian_with_diff(me,x,dx,jac)

    implicit none

    class(numdiff_type),intent(inout)   :: me
    real(wp),dimension(:),intent(in)    :: x    !! vector of variables (size `n`)
    real(wp),dimension(me%n),intent(in) :: dx   !! absolute perturbation (>0)
                                                !! for each variable
    real(wp),dimension(:),intent(out)   :: jac  !! sparse jacobian vector (size
                                                !! `num_nonzero_elements`)

    integer,parameter  :: iord  = 1  !! tells [[diff]] to compute first derivative

    integer                  :: i      !! counter for nonzero elements of jacobian
    type(diff_func)          :: d      !! the [[diff]] class to use
    real(wp)                 :: x0     !! value of ith variable for [[diff]]
    real(wp)                 :: xmin   !! ith variable lower bound
    real(wp)                 :: xmax   !! ith variable upper bound
    real(wp)                 :: deriv  !! derivative value df/dx from [[diff]]
    real(wp)                 :: error  !! estimated error from [[diff]]
    integer                  :: ifail  !! [[diff]] error flag
    real(wp),dimension(me%n) :: xp     !! used by [[dfunc]].
    real(wp),dimension(me%m) :: fvec   !! used by [[dfunc]].
    integer                  :: ir     !! row index of next non-zero element
    integer                  :: ic     !! column index of next non-zero element

    integer :: ic_prev  !! previous column perturbed
    integer :: icount   !! count of number of times a column has been perturbed
    logical :: use_info !! if we are reporting to the user

    if (me%exception_raised) return ! check for exceptions

    ! initialize:
    jac = zero
    use_info = associated(me%info_function)
    if (use_info) ic_prev = -1

    ! set the function for diff:
    call d%set_function(dfunc)

    ! each element is computed one by one
    do i = 1, me%sparsity%num_nonzero_elements

        ir   = me%sparsity%irow(i)
        ic   = me%sparsity%icol(i)
        x0   = x(ic)
        xmin = me%xlow(ic)
        xmax = me%xhigh(ic)

        ! if reporting to the user, have to keep track
        ! of which column is being perturbed, and how
        ! many times:
        if (use_info) then
            if (ic/=ic_prev) icount = 0
        end if

        call d%compute_derivative(iord,x0,xmin,xmax,me%eps,me%acc,deriv,error,ifail)

        if (ifail==0 .or. ifail==1) then
            jac(i) = deriv
        else if (ifail == -1) then
            ! this indicates an exception was
            ! already raised, so just return.
            return
        else
            call me%raise_exception(24, 'compute_jacobian_with_diff',&
                                        'error computing derivative with DIFF.')
            return
        end if

        if (use_info) ic_prev = ic

    end do

    contains

        function dfunc(this,xval) result(fx)

        !! interface to user function for [[diff]] routine.

        implicit none

        class(diff_func),intent(inout)  :: this
        real(wp),intent(in)             :: xval  !! input variable (`ic` variable)
        real(wp)                        :: fx    !! derivative of `ir` function
                                                 !! w.r.t. `xval` variable

        if (use_info) then
            icount = icount + 1
            call me%info_function([ic],icount,x)
        end if

        xp = x
        xp(ic) = xval
        call me%compute_function(xp,fvec,funcs_to_compute=[ir])

        if (me%exception_raised) then ! check for exceptions
            fx = 0.0_wp
            call me%terminate()
        else
            fx = fvec(ir)
        end if

        end function dfunc

    end subroutine compute_jacobian_with_diff