Forwarddiff - A Fortran library for forward mode automatic differentiation

Announcements
1

Recently, I put together the Forwarddiff package: GitHub - Nicholaswogan/forwarddiff

It enables the calculation of derivatives, gradients and Jacobians of fortran subroutines/functions using forward mode automatic differentiation. To create the package, I borrowed code, syntax and inspiration from DNAD, ForwardDiff.jl, and a lecture series by Chris Rackauckas.

The package can efficiently compute dense and sparse (banded and block-diagonal) Jacobians.

I built the package to compute the Jacobian of a reaction diffusion system (~10,000 equations and a Jacobian bandwidth of ~200). It seems to perform well. Also, I did some micro benchmarks comparing it to ForwardDiff.jl, and the performance appears comparable.

I probably have not included compatibility with all Fortran intrinsics. So, I’m happy to accept pull requests, if you want to add support for a new function!

16 Likes
2

Interesting. I’ve tried the same thing using both parameterized types and allocatable arrays to set the derivative sizes. My implementation for using allocatable arrays keeps returning a SEGFAULT when I try add two dual values. PDTs have a weird problem with oneAPI ifort and ifx (2024.2). For the following PDT

Type :: ad_t(nd)

    Integer, LEN :: nd
    Real(WP)     :: f
    Real(WP)     :: df(nd)

  End Type

The following test function

function gfun(x)

 implicit none

 type(ad_t(*)),    dimension(3), intent(in) :: x
 type(ad_t(x%nd)), dimension(3)             :: gfun

!Compute the function g.

#ifdef ASSIGN_TO_FUNCTION
! this doesn't work
 gfun(1) = x(1)**2 + x(2)*5.0_WP + 2.0_WP**x(3) + SQRT(x(2))
 gfun(2) = x(2)**2 + x(3)*5.0_WP + MAX(x(1),x(2)) + SIGN(1.0_WP,x(2)) + ABS(x(3))
 gfun(3) = x(3)**2 + x(1)*5.0_WP + MIN(x(1),x(3)) + x(1)/x(3) + MIN(x(1),x(3))
#else
! this works
 type(ad_t(x%nd)), dimension(3) :: g
 g(1) = x(1)**2 + x(2)*5.0_WP + 2.0_WP**x(3) + SQRT(x(2))
 g(2) = x(2)**2 + x(3)*5.0_WP + MAX(x(1),x(2)) + SIGN(1.0_WP,x(2)) + ABS(x(3))
 g(3) = x(3)**2 + x(1)*5.0_WP + MIN(x(1),x(3)) + x(1)/x(3) + MIN(x(1),x(3))
 gfun(1)%f = g(1)%f
 gfun(2)%f = g(2)%f
 gfun(3)%f = g(3)%f
 gfun(1)%df = g(1)%df
 gfun(2)%df = g(2)%df
 gfun(3)%df = g(3)%df
#endif
end function gfun

gives the following error if I try to assign directly to the function result (gfun).

ad32.x: malloc.c:2379: sysmalloc: Assertion `(old_top == initial_top (av) && old_size == 0) || ((unsigned long) (old_size) >= MINSIZE && prev_inuse (old_top) && ((unsigned long) old_end & (pagesize - 1)) == 0)' failed.
Aborted (core dumped)

It works correctly if I create a local temporary (g) to hold the results of the expression and them assign each g to the corresponding gfun.

I’m perfectly willing to accept that I’m doing something stupid due to my ignorance of just how PDTs are suppose to work but I find this behavior very strange

Here are the results when it works with the x values initialized to
x(1)%f = 2.0
x(2)%f =5.0
x(3)%f = 7.0

 g(1)%f=    1.592360679774998E+02 g(1)%df(1), g(1)%df(2), g(1)%df(3) =     4.000000000000000E+00    5.223606797749979E+00    8.872283911167300E+01
 g(2)%f=    7.300000000000000E+01 g(2)%df(1), g(2)%df(2), g(2)%df(3) =     0.000000000000000E+00    1.100000000000000E+01    6.000000000000000E+00
 g(3)%f=    6.328571428571428E+01 g(3)%df(1), g(3)%df(2), g(3)%df(3) =     7.142857142857143E+00    0.000000000000000E+00    1.395918367346939E+01
 
 Exact dg/dx:
 dg1/dx1, dg1/dx2, dg1/dx3 =     4.000000000000000E+00    5.223606797749979E+00    8.872283911167300E+01
 dg2/dx1, dg2/dx2, dg2/dx3 =     0.000000000000000E+00    1.100000000000000E+01    6.000000000000000E+00
 dg3/dx1, dg3/dx2, dg3/dx3 =     7.142857142857143E+00    0.000000000000000E+00    1.395918367346939E+01

I’ll look at your code to see if I can figure out why my allocatable array version SEGFAULTS. I think it has something to do with using ELEMENTAL functions and relying on allocate on assignment

3

Yes, I rely mostly upon allocation by assignment. But, to deal with some intrinsics in forwarddiff, I had to invoke the allocate command.