compute_sparsity_random_2 Subroutine

    subroutine compute_sparsity_random_2(me,x)

    implicit none

    class(numdiff_type),intent(inout) :: me
    real(wp),dimension(:),intent(in)  :: x !! vector of variables (size `n`)
                                           !! (not used here)

    type :: jac_type
        !! so we can define an array of jacobian columns
        real(wp),dimension(:),allocatable :: jac !! a column of the jacobian
    end type jac_type

    real(wp),dimension(:),allocatable :: coeffs           !! coefficients for computing `xp`
    real(wp),dimension(:,:),allocatable :: xp             !! perturbed `x` array (size `n,num_sparsity_points`)
    type(jac_type),dimension(:),allocatable :: jac_array  !! array of jacobian columns
    type(sparsity_pattern) :: tmp_sparsity_pattern        !! will accumulate the pattern here
                                                          !! and copy it into the class when
                                                          !! finished.
    integer :: i              !! counter for the number of points
    integer :: j              !! counter
    integer :: icol           !! column counter
    integer :: irow           !! row counter
    real(wp) :: dfdx          !! value of a jacobian element
    integer :: n_icol         !! `icol` size counter
    integer :: n_irow         !! `irow` size counter
    integer :: n_linear_icol  !! `linear_icol` size counter
    integer :: n_linear_irow  !! `linear_irow` size counter
    integer :: n_linear_vals  !! `linear_vals` size counter
    type(meth_array) :: class_meths  !! set of finite diff methods to use
    real(wp),dimension(:),allocatable :: jac !! array of jacobian element values
    integer :: info !! status output form [[dsm]]

    ! initialize:
    call me%destroy_sparsity_pattern()

    if (me%exception_raised) return ! check for exceptions

    ! error check:
    if (.not. allocated(me%xlow_for_sparsity) .or. .not. allocated(me%xhigh_for_sparsity)) then
        call me%raise_exception(20,'compute_sparsity_random_2',&
                                   'the x bounds have not been set.')
        return
    end if

    ! compute the coefficients:
    coeffs = divide_interval(me%num_sparsity_points)

    ! save the initial settings:
    tmp_sparsity_pattern = me%sparsity

    ! we create a provisional sparsity pattern here,
    ! to compute one row at a time:
    me%sparsity%irow = [(irow,irow=1,me%m)]
    me%sparsity%sparsity_computed = .true.
    me%sparsity%num_nonzero_elements = me%m
    call me%sparsity%compute_indices()
    if (me%partition_sparsity_pattern) call me%generate_dense_sparsity_partition()
    n_icol = 0
    n_irow = 0
    n_linear_icol = 0
    n_linear_irow = 0
    n_linear_vals = 0
    allocate(jac(me%num_sparsity_points))

    ! the idea here is to compute the (dense) jacobian
    ! one column at at time, and keep a running track of
    ! the ones that have changed.

    ! first precompute the perturbation arrays for each point:
    allocate(xp(me%n,me%num_sparsity_points))
    do i = 1, me%num_sparsity_points
        xp(:,i) = me%xlow_for_sparsity + (me%xhigh_for_sparsity-me%xlow_for_sparsity)*coeffs(i)
    end do

    ! we will use 2-point methods (simple differences):
    class_meths = get_all_methods_in_class(2)

    do icol = 1, me%n  ! column loop

        ! size the array:
        if (allocated(jac_array)) deallocate(jac_array)
        allocate(jac_array(me%num_sparsity_points))

        ! compute the ith column of the jacobian:
        me%sparsity%icol = [(icol, j=1,me%m)]
        do i = 1, me%num_sparsity_points
            call me%compute_jacobian_for_sparsity( icol, class_meths, xp(:,i), jac_array(i)%jac )
            if (me%exception_raised) return ! check for exceptions
        end do

        ! check each row:
        do irow = 1, me%m

            ! get the jacobian values for this row,col for all the points:
            do j = 1, me%num_sparsity_points
                jac(j) = jac_array(j)%jac(irow)
            end do

            ! put the results into the tmp_sparsity_pattern
            if (equal_within_tol([0.0_wp,jac],me%linear_sparsity_tol)) then
                ! they are all zero
                cycle
            else
                if (me%compute_linear_sparsity_pattern) then
                    if (equal_within_tol(jac,me%linear_sparsity_tol)) then
                        ! this is a linear element (constant value)
                        dfdx = sum(jac) / me%num_sparsity_points ! just take the average and use that
                        call expand_vector(tmp_sparsity_pattern%linear_icol,n_linear_icol,me%chunk_size,val=icol)
                        call expand_vector(tmp_sparsity_pattern%linear_irow,n_linear_irow,me%chunk_size,val=irow)
                        call expand_vector(tmp_sparsity_pattern%linear_vals,n_linear_vals,me%chunk_size,val=dfdx)
                        cycle ! this element will not be added to the nonlinear pattern
                    end if
                end if
                ! otherwise, add it to the nonlinear pattern:
                call expand_vector(tmp_sparsity_pattern%icol,n_icol,me%chunk_size,val=icol)
                call expand_vector(tmp_sparsity_pattern%irow,n_irow,me%chunk_size,val=irow)
            end if
        end do

    end do

    ! copy over the computed pattern:
    me%sparsity = tmp_sparsity_pattern

    ! resize to correct size:
    call me%resize_sparsity_vectors(n_icol,n_irow,n_linear_icol,&
                                        n_linear_irow,n_linear_vals)

    call me%sparsity%compute_indices()
    if (me%partition_sparsity_pattern) then
        call me%sparsity%dsm_wrapper(me%n,me%m,info)
        if (info/=1) then
            call me%raise_exception(21,'compute_sparsity_random_2',&
                                       'error partitioning sparsity pattern.')
            return
        end if
    end if
    me%sparsity%sparsity_computed = .true.

    end subroutine compute_sparsity_random_2