compute_jacobian_partitioned Subroutine

private subroutine compute_jacobian_partitioned(me, x, dx, jac, fx)

Compute the Jacobian using finite differences, (using the partitioned sparsity pattern to compute multiple columns at a time).

Arguments

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

vector of variables (size n)

real(kind=wp), intent(in), dimension(me%n) :: dx

absolute perturbation (>0) for each variable

real(kind=wp), intent(out), dimension(:) :: jac

sparse jacobian vector

real(kind=wp), intent(in), optional, dimension(me%m) :: fx

function value at x, if already known


Calls

proc~~compute_jacobian_partitioned~~CallsGraph proc~compute_jacobian_partitioned compute_jacobian_partitioned proc~columns_in_partition_group sparsity_pattern%columns_in_partition_group proc~compute_jacobian_partitioned->proc~columns_in_partition_group proc~init_xwork numdiff_type%init_xwork proc~compute_jacobian_partitioned->proc~init_xwork proc~perturb_x_and_compute_f_partitioned numdiff_type%perturb_x_and_compute_f_partitioned proc~compute_jacobian_partitioned->proc~perturb_x_and_compute_f_partitioned proc~raise_exception numdiff_type%raise_exception proc~compute_jacobian_partitioned->proc~raise_exception proc~select_finite_diff_method_for_partition_group numdiff_type%select_finite_diff_method_for_partition_group proc~compute_jacobian_partitioned->proc~select_finite_diff_method_for_partition_group proc~compute_nominal_function numdiff_type%compute_nominal_function proc~perturb_x_and_compute_f_partitioned->proc~compute_nominal_function interface~unique unique proc~compute_nominal_function->interface~unique 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

Source Code

    subroutine compute_jacobian_partitioned(me,x,dx,jac,fx)

    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
    real(wp),dimension(me%m),intent(in),optional :: fx !! function value at `x`, if already known

    integer                          :: i             !! column counter
    integer                          :: j             !! function evaluation counter
    integer                          :: igroup        !! group number counter
    integer                          :: n_cols        !! number of columns in a group
    integer,dimension(:),allocatable :: cols          !! array of column indices in a group
    integer,dimension(:),allocatable :: nonzero_rows  !! the indices of the nonzero Jacobian
                                                      !! elementes (row numbers) in a group
    integer,dimension(:),allocatable :: indices       !! nonzero indices in `jac` for a group
    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
    real(wp),dimension(me%m)         :: f0           !! function value at the nominal `x`
    logical                          :: f0_computed   !! if `f0` has been computed

    if (me%exception_raised) return ! check for exceptions

    ! initialize:
    jac = zero
    f0_computed = present(fx)
    if (f0_computed) f0 = fx

    ! df is only cleared once. after each group, only the rows
    ! used by that group are reset (to avoid an O(m) reset per group)
    df = zero
    call me%init_xwork(x)

    ! compute by group:
    do igroup = 1, me%sparsity%maxgrp

        ! get the columns in this group:
        call me%sparsity%columns_in_partition_group(igroup,n_cols,cols,&
                                                    nonzero_rows,indices,status_ok)
        if (.not. status_ok) then
            call me%raise_exception(25,'compute_jacobian_partitioned',&
                                       'the partition has not been computed.')
            return
        end if

        if (n_cols>0) then

            if (allocated(nonzero_rows)) then

                select case (me%mode)
                case(1) ! use the specified methods

                    ! note: all the methods must be the same within a group

                    ! compute the columns of the Jacobian in this group:
                    do j = 1, size(me%meth(1)%dx_factors)
                         if (associated(me%info_function)) call me%info_function(cols,j,x)
                         call me%perturb_x_and_compute_f_partitioned(x,me%meth(1)%dx_factors(j),&
                                                         dx,me%meth(1)%df_factors(j),&
                                                         cols,nonzero_rows,df,&
                                                         f0,f0_computed)
                         if (me%exception_raised) return ! check for exceptions
                    end do
                    ! divide by the denominator, which can be different for each column:
                    do i = 1, n_cols
                        associate (rows => me%sparsity%irow(me%sparsity%col_idx(&
                                            me%sparsity%col_ptr(cols(i)):me%sparsity%col_ptr(cols(i)+1)-1)))
                            df(rows) = df(rows) / (me%meth(1)%df_den_factor*dx(cols(i)))
                        end associate
                    end do

                case(2) ! select the method from the class so as not to violate
                        ! the bounds on *any* of the variables in the group

                    ! note: all the classes must be the same within a group

                    call me%select_finite_diff_method_for_partition_group( &
                                    x(cols),me%xlow(cols),me%xhigh(cols),&
                                    dx(cols),me%class_meths(1),fd,status_ok)

                    if (.not. status_ok) then
                        if (me%print_messages) then
                            ! will not consider this a fatal error for now:
                            write(error_unit,'(A,1X,I5,1X,A,1X,*(I5,1X))') &
                                'Error in compute_jacobian_partitioned: '//&
                                'variable bounds violated for group: ',&
                                igroup,'. columns: ',cols
                        end if
                    end if

                    ! compute the columns of the Jacobian in this group:
                    do j = 1, size(fd%dx_factors)
                        if (associated(me%info_function)) call me%info_function(cols,j,x)
                        call me%perturb_x_and_compute_f_partitioned(x,fd%dx_factors(j),&
                                                        dx,fd%df_factors(j),&
                                                        cols,nonzero_rows,df,&
                                                        f0,f0_computed)
                        if (me%exception_raised) return ! check for exceptions
                    end do
                    ! divide by the denominator, which can be different for each column:
                    do i = 1, n_cols
                        associate (rows => me%sparsity%irow(me%sparsity%col_idx(&
                                            me%sparsity%col_ptr(cols(i)):me%sparsity%col_ptr(cols(i)+1)-1)))
                            df(rows) = df(rows) / (fd%df_den_factor*dx(cols(i)))
                        end associate
                    end do

                case default
                    call me%raise_exception(26,'compute_jacobian_partitioned','invalid mode')
                    return
                end select

                ! put result into the output vector:
                jac(indices) = df(nonzero_rows)
                df(nonzero_rows) = zero ! reset for the next group

            end if

        end if

    end do

    end subroutine compute_jacobian_partitioned