compute_jacobian_partitioned Subroutine

    subroutine compute_jacobian_partitioned(me,x,dx,jac)

    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

    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
    integer,dimension(:),allocatable :: col_indices   !! nonzero indices in `jac` for a column
    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

    if (me%exception_raised) return ! check for exceptions

    ! initialize:
    jac = zero

    ! 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:
                    df = zero
                    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)
                         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
                        num_nonzero_elements_in_col = count(me%sparsity%icol==cols(i))
                        if (allocated(col_indices)) deallocate(col_indices)
                        allocate(col_indices(num_nonzero_elements_in_col))
                        ! col_indices = pack(me%sparsity%indices,mask=me%sparsity%icol==cols(i))
                        block
                            integer :: j,n
                            n = 0
                            do j = 1, size(me%sparsity%icol)
                                if (me%sparsity%icol(j)==cols(i)) then
                                    n = n + 1
                                    col_indices(n) = j
                                end if
                            end do
                        end block

                        df(me%sparsity%irow(col_indices)) = df(me%sparsity%irow(col_indices)) / &
                                                            (me%meth(1)%df_den_factor*dx(cols(i)))
                    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:
                    df = zero
                    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)
                        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
                        num_nonzero_elements_in_col = count(me%sparsity%icol==cols(i))
                        if (allocated(col_indices)) deallocate(col_indices)
                        allocate(col_indices(num_nonzero_elements_in_col))
                        col_indices = pack(me%sparsity%indices,mask=me%sparsity%icol==cols(i))
                        df(me%sparsity%irow(col_indices)) = df(me%sparsity%irow(col_indices)) / &
                                                            (fd%df_den_factor*dx(cols(i)))
                    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)

            end if

        end if

    end do

    end subroutine compute_jacobian_partitioned