diff --git a/CMakeLists.txt b/CMakeLists.txt index 68bf7a1..86dd52a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ cmake_minimum_required(VERSION "3.14" ) -project(FORWARDDIFF LANGUAGES Fortran VERSION "0.1.1") +project(FORWARDDIFF LANGUAGES Fortran VERSION "0.1.2") set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_BINARY_DIR}/modules") add_subdirectory(src) diff --git a/src/forwarddiff_derivative.f90 b/src/forwarddiff_derivative.f90 index cff38b5..ca60b9b 100644 --- a/src/forwarddiff_derivative.f90 +++ b/src/forwarddiff_derivative.f90 @@ -9,43 +9,57 @@ module forwarddiff_derivative public :: jacobian, jacobian_sig abstract interface + !> Interface for the function input to the `derivative` + !> routine for computing the derivative of scalar + !> functions. function derivative_sig(x) result(res) import :: dual - type(dual), intent(in) :: x - type(dual) :: res + type(dual), intent(in) :: x !! Input scalar + type(dual) :: res !! Resulting scalar end function + !> Interface for the function input to the `gradient` + !> subroutine for computing the gradient of a function + !> mapping a vector to a scalar. function gradient_sig(x) result(res) import :: dual - type(dual), target, intent(in) :: x(:) - type(dual) :: res + type(dual), target, intent(in) :: x(:) !! Input vector + type(dual) :: res !! Resulting scalar end function + !> Interface for the function input to the `jacobian` + !> subroutine for computing the jacobian of a function + !> mapping a vector to a vector. subroutine jacobian_sig(x, f) import :: dual - type(dual), target, intent(in) :: x(:) - type(dual), target, intent(out) :: f(:) + type(dual), target, intent(in) :: x(:) !! Input vector + type(dual), target, intent(out) :: f(:) !! Resulting vector end subroutine end interface contains + !> Computes the derivative of input scalar function + !> `fcn` at the input `x`. subroutine derivative(fcn, x, f, dfdx) - procedure(derivative_sig) :: fcn + procedure(derivative_sig) :: fcn !! Input scalar function real(wp), intent(in) :: x - real(wp), intent(out) :: f - real(wp), intent(out) :: dfdx + real(wp), intent(out) :: f !! `fcn` evaluated at `x` + real(wp), intent(out) :: dfdx !! The derivative of `fcn` at `x` type(dual) :: ff ff = fcn(dual(x, [1.0_wp])) f = ff%val dfdx = ff%der(1) end subroutine + !> Computes the gradient of the input function + !> `fcn` at the input `x`. subroutine gradient(fcn, x, f, dfdx, err) - procedure(gradient_sig) :: fcn + procedure(gradient_sig) :: fcn !! Input function mapping a vector to a scalar real(wp), intent(in) :: x(:) - real(wp), intent(out) :: f - real(wp), intent(out) :: dfdx(:) + real(wp), intent(out) :: f !! `fcn` evaluated at `x` + real(wp), intent(out) :: dfdx(:) !! The gradient of `fcn` at `x` + !> If an error occurs, `err` will be allocated with an error message. character(:), allocatable, intent(out) :: err type(dual) :: xx(size(x)) @@ -78,16 +92,55 @@ subroutine gradient(fcn, x, f, dfdx, err) end subroutine + !> Computes the Jacobian of the input function `fcn + !> at the input `x`, with support for some sparse Jacobians. Only can + !> consider square Jacobians. subroutine jacobian(fcn, x, f, dfdx, jt, bandwidth, blocksize, err) use forwarddiff_const, only: DenseJacobian, BandedJacobian, BlockDiagonalJacobian - procedure(jacobian_sig) :: fcn + procedure(jacobian_sig) :: fcn !! Input function mapping a vector to a vector real(wp), intent(in) :: x(:) - real(wp), intent(out) :: f(:) - real(wp), intent(out) :: dfdx(:,:) + real(wp), intent(out) :: f(:) !! `fcn` evaluated at `x` + real(wp), intent(out) :: dfdx(:,:) + !! The Jacobian of `fcn` evaluated at `x`. For dense Jacobians (`jt == DenseJacobian`), + !! `dfdx` has shape (n,n) where n is `size(x)`. The first dimension indexes rows + !! of the Jacobian, while the second dimension indexes columns. For example, + !! + !! | df(1,1) df(1,2) df(1,3) | + !! J = | df(2,1) df(2,2) df(2,3) | + !! | df(3,1) df(3,2) df(3,3) | + !! + !! If the Jacobian is banded (`jt == BandedJacobian`), then `dfdx` has shape + !! (bandwidth,n), where bandwidth is the Jacobian bandwidth. In this case, the + !! The diagonals of the Jacobian are loaded into the rows of `dfdx` + !! + !! | df(1,1) df(1,2) 0 0 0 | + !! | df(2,1) df(2,2) df(2,3) 0 0 | -> | 0 df(1,2) df(2,3) df(3,4) df(4,5) | + !! J = | 0 df(3,2) df(3,3) df(3,4) 0 | -> | df(1,1) df(2,2) df(3,3) df(4,4) df(5,5) | + !! | 0 0 df(4,3) df(4,4) df(4,5) | -> | df(2,1) df(3,2) df(4,3) df(5,4) 0 | + !! | 0 0 0 df(5,4) df(5,5) | + !! + !! If the Jacobian is block-diagonal (`jt == BlockDiagonalJacobian`), then `dfdx` has + !! shape (blocksize,n) where blocksize is the blocksize of the Jacobian. In this case, + !! `dfdx` stores the Jacobian entries in the following way: + !! + !! | df(1,1) df(1,2) 0 0 | + !! J = | df(2,1) df(2,2) 0 0 | -> | df(1,1) df(1,2) df(3,3) df(3,4) | + !! | 0 0 df(3,3) df(3,4) | -> | df(2,1) df(2,2) df(4,3) df(4,4) | + !! | 0 0 df(4,3) df(4,4) | integer, optional, intent(in) :: jt + !! Jacobian sparsity indicator. The default is `jt = DenseJacobian`. + !! - If `jt == DenseJacobian`, then the algorithm assumes the Jacobian is dense. + !! - If `jt == BandedJacobian`, then the algorithm assumes the Jacobian is banded + !! with `bandwidth`. + !! - If `jt == BlockDiagonalJacobian`, then the algorithm assumes the Jacobian is + !! block diagonal with `blocksize`. Note that size(x) must be an integer multiple + !! of `blocksize`. integer, optional, intent(in) :: bandwidth + !! The Jacobian bandwidth for `jt == BandedJacobian` integer, optional, intent(in) :: blocksize + !! The Jacobian blocksize for `jt == BlockDiagonalJacobian` character(:), allocatable, intent(out) :: err + !! If an error occurs, `err` will be allocated with an error message. integer :: jt_ @@ -121,6 +174,9 @@ subroutine jacobian(fcn, x, f, dfdx, jt, bandwidth, blocksize, err) endif call jacobian_blockdiagonal(fcn, x, f, dfdx, blocksize, err) if (allocated(err)) return + else + err = 'Invalid value for the Jacobian type indicator `jt`.' + return endif end subroutine