compute_jacobian_for_sparsity Subroutine

private subroutine compute_jacobian_for_sparsity(me, i, x, jac, f0, f0_computed, dx_col)

A separate version of compute_jacobian to be used only when computing the sparsity pattern in compute_sparsity_random_2. It uses sparsity_class_meths and the sparsity dperts and bounds.

Note

Based on compute_jacobian. The index manipulation here could be greatly simplified, since we realdy know we are computed all the elements in one column.

Type Bound

numdiff_type

Arguments

Type IntentOptional Attributes Name
class(numdiff_type), intent(inout) :: me
integer, intent(in) :: i

the column being computed
real(kind=wp), intent(in), dimension(:) :: x

vector of variables (size n)
real(kind=wp), intent(out), dimension(:), allocatable :: jac

sparse jacobian vector
real(kind=wp), intent(inout), dimension(me%m) :: f0

function value at x (computed on first use and reused for subsequent columns)
logical, intent(inout) :: f0_computed

if f0 has been computed
real(kind=wp), intent(out) :: dx_col

the perturbation used for column i

Calls

proc~~compute_jacobian_for_sparsity~~CallsGraph proc~compute_jacobian_for_sparsity numdiff_type%compute_jacobian_for_sparsity proc~compute_sparsity_perturbation_vector numdiff_type%compute_sparsity_perturbation_vector proc~compute_jacobian_for_sparsity->proc~compute_sparsity_perturbation_vector proc~init_xwork numdiff_type%init_xwork proc~compute_jacobian_for_sparsity->proc~init_xwork proc~perturb_x_and_compute_f numdiff_type%perturb_x_and_compute_f proc~compute_jacobian_for_sparsity->proc~perturb_x_and_compute_f proc~select_finite_diff_method numdiff_type%select_finite_diff_method proc~compute_jacobian_for_sparsity->proc~select_finite_diff_method proc~compute_perturb_vector numdiff_type%compute_perturb_vector proc~compute_sparsity_perturbation_vector->proc~compute_perturb_vector proc~compute_nominal_function numdiff_type%compute_nominal_function proc~perturb_x_and_compute_f->proc~compute_nominal_function interface~unique unique proc~compute_nominal_function->interface~unique proc~raise_exception numdiff_type%raise_exception proc~compute_perturb_vector->proc~raise_exception proc~unique_int unique_int interface~unique->proc~unique_int proc~unique_real unique_real interface~unique->proc~unique_real interface~expand_vector expand_vector proc~unique_int->interface~expand_vector interface~sort_ascending sort_ascending proc~unique_int->interface~sort_ascending proc~unique_real->interface~expand_vector proc~unique_real->interface~sort_ascending proc~expand_vector_int expand_vector_int interface~expand_vector->proc~expand_vector_int proc~expand_vector_real expand_vector_real interface~expand_vector->proc~expand_vector_real proc~sort_ascending_int sort_ascending_int interface~sort_ascending->proc~sort_ascending_int proc~sort_ascending_real sort_ascending_real interface~sort_ascending->proc~sort_ascending_real interface~swap swap proc~sort_ascending_int->interface~swap proc~sort_ascending_real->interface~swap proc~swap_int swap_int interface~swap->proc~swap_int proc~swap_real swap_real interface~swap->proc~swap_real

Called by

proc~~compute_jacobian_for_sparsity~~CalledByGraph proc~compute_jacobian_for_sparsity numdiff_type%compute_jacobian_for_sparsity proc~compute_sparsity_random_2 compute_sparsity_random_2 proc~compute_sparsity_random_2->proc~compute_jacobian_for_sparsity

Source Code

    subroutine compute_jacobian_for_sparsity(me,i,x,jac,f0,f0_computed,dx_col)

    implicit none

    class(numdiff_type),intent(inout)             :: me
    integer,intent(in)                            :: i           !! the column being computed
    real(wp),dimension(:),intent(in)              :: x           !! vector of variables (size `n`)
    real(wp),dimension(:),allocatable,intent(out) :: jac         !! sparse jacobian vector
    real(wp),dimension(me%m),intent(inout)        :: f0          !! function value at `x`
                                                                 !! (computed on first use and
                                                                 !! reused for subsequent columns)
    logical,intent(inout)                         :: f0_computed !! if `f0` has been computed
    real(wp),intent(out)                          :: dx_col      !! the perturbation used for column `i`

    real(wp),dimension(me%n) :: dx  !! absolute perturbation (>0) for each variable
    integer,dimension(:),allocatable :: nonzero_elements_in_col  !! the indices of the
                                                                 !! nonzero Jacobian
                                                                 !! elements in a column
    integer                  :: j   !! function evaluation counter
    real(wp),dimension(me%m) :: df  !! accumulated function
    type(finite_diff_method) :: fd  !! a finite different method (when
                                    !! specifying class rather than the method)
    logical :: status_ok                    !! error flag
    integer :: num_nonzero_elements_in_col  !! number of nonzero elements in a column

    if (me%exception_raised) return ! check for exceptions

    ! Note that a sparsity pattern has already been set

    ! size the jacobian vector:
    allocate(jac(me%sparsity%num_nonzero_elements))

    ! compute the perturbation vector (really we only need dx(i)):
    call me%compute_sparsity_perturbation_vector(x,dx)
    if (me%exception_raised) return ! check for exceptions
    dx_col = dx(i)

    ! initialize:
    jac = zero

    ! determine functions to compute for this column:
    num_nonzero_elements_in_col = count(me%sparsity%icol==i)
    if (num_nonzero_elements_in_col/=0) then ! there are functions to compute

        nonzero_elements_in_col = pack(me%sparsity%irow,mask=me%sparsity%icol==i)

        call me%select_finite_diff_method(x(i),me%xlow_for_sparsity(i),me%xhigh_for_sparsity(i),&
                                            dx(i),me%sparsity_class_meths,fd,status_ok)
        if (.not. status_ok) then
            if (me%print_messages) then
                write(error_unit,'(A,1X,I5)') &
                'Error in compute_jacobian_for_sparsity: variable bounds violated for column: ',i
            end if
        end if

        ! compute this column of the Jacobian:
        call me%init_xwork(x)
        df = zero
        do j = 1, size(fd%dx_factors)
            if (associated(me%info_function)) call me%info_function([i],j,x)
            call me%perturb_x_and_compute_f(x,fd%dx_factors(j),&
                                            dx,fd%df_factors(j),&
                                            i,nonzero_elements_in_col,df,&
                                            f0,f0_computed)
            if (me%exception_raised) return ! check for exceptions
        end do
        df(nonzero_elements_in_col) = df(nonzero_elements_in_col) / &
                                        (fd%df_den_factor*dx(i))

        ! put result into the output vector:
        jac(pack(me%sparsity%indices,mask=me%sparsity%icol==i)) = &
            df(nonzero_elements_in_col)

    end if

    end subroutine compute_jacobian_for_sparsity