compute_sparsity_random Subroutine

private subroutine compute_sparsity_random(me, x)

Compute the sparsity pattern by computing the function at three "random" points in the [xlow_for_sparsity, xhigh_for_sparsity] interval and checking if the function values are the same.

Arguments

Type IntentOptional Attributes Name
class(numdiff_type), intent(inout) :: me
real(kind=wp), intent(in), dimension(:) :: x

vector of variables (size n) (not used here)

Calls

proc~~compute_sparsity_random~~CallsGraph proc~compute_sparsity_random numerical_differentiation_module::compute_sparsity_random interface~expand_vector numdiff_utilities_module::expand_vector proc~compute_sparsity_random->interface~expand_vector proc~compute_indices numerical_differentiation_module::sparsity_pattern%compute_indices proc~compute_sparsity_random->proc~compute_indices proc~destroy_sparsity_pattern numerical_differentiation_module::numdiff_type%destroy_sparsity_pattern proc~compute_sparsity_random->proc~destroy_sparsity_pattern proc~dsm_wrapper numerical_differentiation_module::sparsity_pattern%dsm_wrapper proc~compute_sparsity_random->proc~dsm_wrapper proc~equal_within_tol numdiff_utilities_module::equal_within_tol proc~compute_sparsity_random->proc~equal_within_tol proc~raise_exception numerical_differentiation_module::numdiff_type%raise_exception proc~compute_sparsity_random->proc~raise_exception proc~resize_sparsity_vectors numerical_differentiation_module::numdiff_type%resize_sparsity_vectors proc~compute_sparsity_random->proc~resize_sparsity_vectors proc~expand_vector_int numdiff_utilities_module::expand_vector_int interface~expand_vector->proc~expand_vector_int proc~expand_vector_real numdiff_utilities_module::expand_vector_real interface~expand_vector->proc~expand_vector_real proc~destroy_sparsity numerical_differentiation_module::sparsity_pattern%destroy_sparsity proc~destroy_sparsity_pattern->proc~destroy_sparsity proc~dsm dsm_module::dsm proc~dsm_wrapper->proc~dsm proc~resize_sparsity_vectors->interface~expand_vector proc~degr dsm_module::degr proc~dsm->proc~degr proc~ido dsm_module::ido proc~dsm->proc~ido proc~numsrt dsm_module::numsrt proc~dsm->proc~numsrt proc~seq dsm_module::seq proc~dsm->proc~seq proc~setr dsm_module::setr proc~dsm->proc~setr proc~slo dsm_module::slo proc~dsm->proc~slo proc~srtdat dsm_module::srtdat proc~dsm->proc~srtdat proc~ido->proc~numsrt

Source Code

    subroutine compute_sparsity_random(me,x)

    implicit none

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

    integer :: i !! column counter
    integer :: j !! row counter
    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
    integer,dimension(me%m) :: idx !! indices to compute `[1,2,...,m]`
    real(wp),dimension(me%n) :: x1 !! perturbed variable vector
    real(wp),dimension(me%n) :: x2 !! perturbed variable vector
    real(wp),dimension(me%n) :: x3 !! perturbed variable vector
    real(wp),dimension(me%n) :: x4 !! perturbed variable vector
    real(wp),dimension(me%m) :: f1 !! function evaluation
    real(wp),dimension(me%m) :: f2 !! function evaluation
    real(wp),dimension(me%m) :: f3 !! function evaluation
    real(wp),dimension(me%m) :: f4 !! function evaluation
    real(wp) :: dfdx1 !! for linear sparsity estimation
    real(wp) :: dfdx2 !! for linear sparsity estimation
    real(wp) :: dfdx3 !! for linear sparsity estimation
    real(wp) :: dfdx  !! for linear sparsity estimation
    integer :: info !! status output form [[dsm]]

    real(wp),dimension(4),parameter :: coeffs = [0.20123456787654321_wp,&
                                                 0.40123456787654321_wp,&
                                                 0.60123456787654321_wp,&
                                                 0.80123456787654321_wp]
        !! Pick three pseudo-random points roughly equally spaced.
        !! (add some noise in attempt to avoid freak zeros)
        !!````
        !! xlow---|----|--x--|---xhigh
        !!        1    2     3
        !!````
        !!
        !! Also using an extra point to estimate the
        !! constant elements of the Jacobian.

    ! 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(18,'compute_sparsity_random',&
                                   'the x bounds have not been set.')
        return
    end if

    associate (xlow => me%xlow_for_sparsity, xhigh => me%xhigh_for_sparsity)

        ! we will compute all the functions:
        idx = [(i,i=1,me%m)]

        n_icol = 0  ! initialize vector size counters
        n_irow = 0
        n_linear_icol = 0
        n_linear_irow = 0
        n_linear_vals = 0

        ! define a nominal point roughly in the middle:
        x2 = xlow + (xhigh-xlow)*coeffs(2)
        call me%compute_function(x2,f2,idx)
        if (me%exception_raised) return ! check for exceptions

        if (me%compute_linear_sparsity_pattern) then
            ! we need another point where we perturb all the variables
            ! to check to make sure it is linear in only one variable.
            x4 = xlow + (xhigh-xlow)*coeffs(4)
            call me%compute_function(x4,f4,idx)
            if (me%exception_raised) return ! check for exceptions
        end if

        do i = 1, me%n  ! columns of Jacobian

            ! restore nominal:
            x1 = x2
            x3 = x2

            x1(i) = xlow(i) + (xhigh(i)-xlow(i))*coeffs(1)
            x3(i) = xlow(i) + (xhigh(i)-xlow(i))*coeffs(3)

            call me%compute_function(x1,f1,idx)
            if (me%exception_raised) return ! check for exceptions
            call me%compute_function(x3,f3,idx)
            if (me%exception_raised) return ! check for exceptions

            do j = 1, me%m ! each function (rows of Jacobian)
                if (equal_within_tol([f1(j),f2(j),f3(j)], me%function_precision_tol)) then
                    ! no change in the function, so no sparsity element here.
                    cycle
                else
                    if (me%compute_linear_sparsity_pattern) then
                        ! if computing the linear pattern separately.
                        ! do the three points lie on a line? (x1,f1) -> (x2, f2), -> (x3,f3)
                        ! [check that the slopes are equal within a tolerance]
                        dfdx1 = (f1(j)-f2(j)) / (x1(i)-x2(i)) ! slope of line from 1->2
                        dfdx2 = (f1(j)-f3(j)) / (x1(i)-x3(i)) ! slope of line from 1->3
                        dfdx3 = (f1(j)-f4(j)) / (x1(i)-x4(i)) ! slope of line from 1->4
                        if (equal_within_tol([dfdx1,dfdx2,dfdx3],me%linear_sparsity_tol)) then
                            ! this is a linear element (constant value)
                            dfdx = (dfdx1 + dfdx2 + dfdx3) / 3.0_wp ! just take the average and use that
                            call expand_vector(me%sparsity%linear_icol,n_linear_icol,me%chunk_size,val=i)
                            call expand_vector(me%sparsity%linear_irow,n_linear_irow,me%chunk_size,val=j)
                            call expand_vector(me%sparsity%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(me%sparsity%icol,n_icol,me%chunk_size,val=i)
                    call expand_vector(me%sparsity%irow,n_irow,me%chunk_size,val=j)
                end if
            end do

        end do

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

    end associate

    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(19,'compute_sparsity_random',&
                                       'error partitioning sparsity pattern.')
            return
        end if
    end if

    ! finished:
    me%sparsity%sparsity_computed = .true.

    end subroutine compute_sparsity_random