compute_jacobian_standard Subroutine

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

    integer,dimension(:),allocatable :: nonzero_elements_in_col  !! the indices of the
                                                                 !! nonzero Jacobian
                                                                 !! elements in a column
    integer                  :: i   !! column counter
    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
    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

    if (.not. allocated(me%sparsity%col_ptr)) then
        call me%raise_exception(31,'compute_jacobian_standard',&
                                   'the sparsity column index has not been computed.')
        return
    end if

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

    ! compute Jacobian matrix column-by-column:
    do i=1,me%n

        ! determine functions to compute for this column:
        associate (col_indices => me%sparsity%col_idx(me%sparsity%col_ptr(i):me%sparsity%col_ptr(i+1)-1))

        num_nonzero_elements_in_col = size(col_indices)
        if (num_nonzero_elements_in_col/=0) then ! there are functions to compute

            nonzero_elements_in_col = me%sparsity%irow(col_indices)

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

                ! compute this column of the Jacobian:
                do j = 1, size(me%meth(i)%dx_factors)
                    if (associated(me%info_function)) call me%info_function([i],j,x)
                    call me%perturb_x_and_compute_f(x,me%meth(i)%dx_factors(j),&
                                                    dx,me%meth(i)%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) / &
                                              (me%meth(i)%df_den_factor*dx(i))

            case(2) ! select the method from the class so as not to violate the bounds

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

                ! compute this column of the Jacobian:
                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))

            case default
                call me%raise_exception(23,'compute_jacobian_standard',&
                                           'invalid mode')
                return
            end select

            ! put result into the output vector:
            jac(col_indices) = df(nonzero_elements_in_col)
            df(nonzero_elements_in_col) = zero ! reset for the next column

        end if

        end associate

    end do

    end subroutine compute_jacobian_standard