compute_jacobian_for_sparsity Subroutine

private subroutine compute_jacobian_for_sparsity(me, i, class_meths, x, jac)

A separate version of compute_jacobian to be used only when computing the sparsity pattern in compute_sparsity_random_2. It uses 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
type(meth_array), intent(in) :: class_meths

set of finite diff methods to use
real(kind=wp), intent(in), dimension(:) :: x

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

sparse jacobian vector

Calls

proc~~compute_jacobian_for_sparsity~~CallsGraph proc~compute_jacobian_for_sparsity numerical_differentiation_module::numdiff_type%compute_jacobian_for_sparsity proc~compute_sparsity_perturbation_vector numerical_differentiation_module::numdiff_type%compute_sparsity_perturbation_vector proc~compute_jacobian_for_sparsity->proc~compute_sparsity_perturbation_vector proc~perturb_x_and_compute_f numerical_differentiation_module::numdiff_type%perturb_x_and_compute_f proc~compute_jacobian_for_sparsity->proc~perturb_x_and_compute_f proc~select_finite_diff_method numerical_differentiation_module::numdiff_type%select_finite_diff_method proc~compute_jacobian_for_sparsity->proc~select_finite_diff_method proc~compute_perturb_vector numerical_differentiation_module::numdiff_type%compute_perturb_vector proc~compute_sparsity_perturbation_vector->proc~compute_perturb_vector proc~raise_exception numerical_differentiation_module::numdiff_type%raise_exception proc~compute_perturb_vector->proc~raise_exception

Called by

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

Source Code

    subroutine compute_jacobian_for_sparsity(me,i,class_meths,x,jac)

    implicit none

    class(numdiff_type),intent(inout)             :: me
    integer,intent(in)                            :: i           !! the column being computed
    type(meth_array),intent(in)                   :: class_meths !! set of finite diff methods to use
    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%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

    ! 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),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:
        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)
            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