Compute the Jacobian using finite differences, (using the partitioned sparsity pattern to compute multiple columns at a time).
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(numdiff_type), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(in), | dimension(:) | :: | x |
vector of variables (size |
|
| 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 |
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