compute_jacobian_standard Subroutine

    subroutine compute_jacobian_standard(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 (size
                                                !! `num_nonzero_elements`)

    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

    if (me%exception_raised) return ! check for exceptions

    ! initialize:
    jac = zero

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

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

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

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

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

            ! put result into the output vector:
            jac(pack(me%sparsity%indices,mask=me%sparsity%icol==i)) = &
                df(nonzero_elements_in_col)

        end if

    end do

    end subroutine compute_jacobian_standard