module active_module
use w2f__types
implicit none
private
public :: deriv, active, saxpy, sax, setderiv, zero_deriv, &
& getderiv, init_adjoint, set_adjoint, get_adjoint, &
& adjoint_interpreter, dump_tape, dump_adjoint, &
& init_fm_rm, finish
integer, parameter :: tape_chunk_size=100000
integer :: tape_counter=1
integer, parameter :: addresses_chunk_size=100
integer :: addresses_counter=0
integer :: nr_adjoints
! need a static here if active common blocks present like in box_timestep
integer, parameter :: nr_tangents=6
type deriv
sequence
double precision, dimension(nr_tangents) :: d
end type deriv
type adj
sequence
double precision, dimension(:), pointer :: d=>NULL()
end type adj
!
! active needs to be a sequence type
! with no initialization
!
type active
sequence
double precision :: v
type(deriv) :: d
end type active
!
! x=x+a*y
!
type adjoint_saxpy
integer :: type ! 0=zero_deriv, 1=setderiv, 2=saxpy, 3=sax
integer :: addr_x
integer :: addr_y
double precision :: a
end type adjoint_saxpy
!
! can make this allocatable etc.
!
type(adjoint_saxpy), dimension(:), allocatable :: tape
integer, dimension(:), allocatable :: addresses
type(adj), dimension(:), allocatable :: adjoint
interface init_fm_rm
module procedure init_fm_rm
end interface
interface finish
module procedure finish
end interface
interface saxpy
module procedure saxpy_a_a
end interface
interface setderiv
module procedure setderiv_a_a, setderiv_a_c
end interface
interface zero_deriv
module procedure zero_deriv_a
end interface
interface sax
module procedure sax_d_a_a, sax_i_a_a
end interface
interface getderiv
module procedure getderiv
end interface
interface set_adjoint
module procedure set_adjoint_scal, set_adjoint_vec
end interface
interface get_adjoint
module procedure get_adjoint
end interface
interface adjoint_interpreter
module procedure adjoint_interpreter
end interface
interface dump_tape
module procedure dump_tape
end interface
interface dump_adjoint
module procedure dump_adjoint
end interface
contains
subroutine init_fm_rm(nr_tan, nr_adj)
integer, intent(in) :: nr_tan, nr_adj
! nr_tangents=nr_tan
nr_adjoints=nr_adj
end subroutine init_fm_rm
subroutine finish
deallocate(tape)
deallocate(addresses)
deallocate(adjoint)
end subroutine finish
subroutine check_addresses_size
implicit none
integer, dimension(:), allocatable :: addresses_tmp
integer s
if (.not.allocated(addresses)) then
allocate(addresses(addresses_chunk_size))
endif
s=size(addresses)
if (addresses_counter>=s) then
allocate(addresses_tmp(s))
addresses_tmp=addresses
deallocate(addresses)
allocate(addresses(s+addresses_chunk_size))
addresses=addresses_tmp
deallocate(addresses_tmp)
endif
end subroutine check_addresses_size
integer function hash(x)
integer, intent(in) :: x
integer i
call check_addresses_size
do i=1,addresses_counter
if (addresses(i)==x) then
hash=i
return
endif
end do
addresses_counter=addresses_counter+1
addresses(addresses_counter)=x
hash=addresses_counter
return
end function hash
!
! Local routine to handle dynamic allocation of tape
!
subroutine check_tape_size
type(adjoint_saxpy), dimension(:), allocatable :: tape_tmp
integer s
if (.not.allocated(tape)) then
allocate(tape(tape_chunk_size))
endif
s=size(tape)
if (tape_counter>s) then
allocate(tape_tmp(s))
tape_tmp=tape
deallocate(tape)
allocate(tape(s+tape_chunk_size))
tape(1:s)=tape_tmp
deallocate(tape_tmp)
endif
end subroutine check_tape_size
!
! chain rule saxpy to be used in forward and reverse modes
!
subroutine saxpy_a_a(a,x,y)
double precision, intent(in) :: a
type(active), intent(in) :: x
type(active), intent(inout) :: y
type(adjoint_saxpy) :: adj_saxpy
integer iaddr
call check_tape_size
adj_saxpy%type=2
adj_saxpy%addr_x=hash(iaddr(x))
adj_saxpy%addr_y=hash(iaddr(y))
adj_saxpy%a=a
tape(tape_counter)=adj_saxpy
tape_counter=tape_counter+1
y%d%d=y%d%d+x%d%d*a
end subroutine saxpy_a_a
!
! chain rule saxpy to be used in forward and reverse modes
! derivative component of y is equal to zero initially
! note: y needs to be inout as otherwise value component gets
! zeroed out
!
subroutine sax_d_a_a(a,x,y)
double precision, intent(in) :: a
type(active), intent(in) :: x
type(active), intent(inout) :: y
type(adjoint_saxpy) :: adj_saxpy
integer iaddr
call check_tape_size
adj_saxpy%type=3
adj_saxpy%addr_x=hash(iaddr(x))
adj_saxpy%addr_y=hash(iaddr(y))
adj_saxpy%a=a
tape(tape_counter)=adj_saxpy
tape_counter=tape_counter+1
y%d%d=x%d%d*a
end subroutine sax_d_a_a
subroutine sax_i_a_a(a,x,y)
integer(kind=w2f__i8), intent(in) :: a
type(active), intent(in) :: x
type(active), intent(inout) :: y
type(adjoint_saxpy) :: adj_saxpy
integer iaddr
call check_tape_size
adj_saxpy%type=3
adj_saxpy%addr_x=hash(iaddr(x))
adj_saxpy%addr_y=hash(iaddr(y))
adj_saxpy%a=a
tape(tape_counter)=adj_saxpy
tape_counter=tape_counter+1
y%d%d=x%d%d*a
end subroutine sax_i_a_a
!
! set derivative of y to be equal to derivative of x
! note: making y inout allows for already existing active
! variables to become the target of a derivative assignment
!
subroutine setderiv_a_a(y,x)
type(active), intent(inout) :: y
type(active), intent(in) :: x
type(adjoint_saxpy) :: adj_saxpy
integer iaddr
call check_tape_size
adj_saxpy%type=1
adj_saxpy%addr_x=hash(iaddr(x))
adj_saxpy%addr_y=hash(iaddr(y))
adj_saxpy%a=1.D0
tape(tape_counter)=adj_saxpy
tape_counter=tape_counter+1
y%d%d=x%d%d
end subroutine setderiv_a_a
!
! initialization of derivative component with a constant
!
subroutine setderiv_a_c(y,x)
type(active), intent(inout) :: y
double precision, dimension(:), intent(in) :: x
y%d%d=x
end subroutine setderiv_a_c
!
! set derivative components to 0.0
!
subroutine zero_deriv_a(x)
type(active) :: x
type(adjoint_saxpy) :: adj_saxpy
integer iaddr
call check_tape_size
adj_saxpy%type=0
adj_saxpy%addr_x=hash(iaddr(x))
adj_saxpy%addr_y=0
adj_saxpy%a=0.D0
tape(tape_counter)=adj_saxpy
tape_counter=tape_counter+1
x%d%d=0.0d0
end subroutine zero_deriv_a
!
! return double value of deriv
!
function getderiv(x) result(res)
type(active), intent(in) :: x
double precision, dimension(nr_tangents) :: res
res=x%d%d
return
end function getderiv
subroutine init_adjoint
integer i
allocate(adjoint(addresses_counter))
do i=1,addresses_counter
allocate(adjoint(i)%d(nr_adjoints))
adjoint(i)%d=0.D0
end do
end subroutine init_adjoint
!
! initialization of adjoint component with a constant
!
subroutine set_adjoint_vec(y,x)
type(active), intent(inout) :: y
double precision, dimension(nr_adjoints), &
& intent(in) :: x
integer iaddr
adjoint(hash(iaddr(y)))%d=x
end subroutine set_adjoint_vec
subroutine set_adjoint_scal(y,x)
type(active), intent(inout) :: y
double precision, intent(in) :: x
integer iaddr
adjoint(hash(iaddr(y)))%d=x
end subroutine set_adjoint_scal
function get_adjoint(y) result(res)
type(active), intent(in) :: y
double precision, dimension(nr_adjoints) :: res
integer iaddr
res=adjoint(hash(iaddr(y)))%d
return
end function get_adjoint
subroutine adjoint_interpreter
integer i
do i=tape_counter-1,1,-1
select case (tape(i)%type)
case (0)
adjoint(tape(i)%addr_x)%d=0.D0
case (1)
adjoint(tape(i)%addr_x)%d=adjoint(tape(i)%addr_x)%d &
& +adjoint(tape(i)%addr_y)%d
adjoint(tape(i)%addr_y)%d=0.D0
case (2)
adjoint(tape(i)%addr_x)%d= adjoint(tape(i)%addr_x)%d &
& +adjoint(tape(i)%addr_y)%d* &
& tape(i)%a
case (3)
adjoint(tape(i)%addr_x)%d= adjoint(tape(i)%addr_x)%d &
& +adjoint(tape(i)%addr_y)%d &
& *tape(i)%a
adjoint(tape(i)%addr_y)%d=0.D0
end select
end do
end subroutine adjoint_interpreter
subroutine dump_tape
integer i
print*, "tape"
do i=1,tape_counter-1
print*, tape(i)
end do
end subroutine dump_tape
subroutine dump_adjoint
integer i
print*, "adjoint"
do i=1,addresses_counter-1
print*, adjoint(i)%d
end do
end subroutine dump_adjoint
end module