So I'm working on an automatic differentiation toolbox in Fortran using operator overloading. I have previously implemented this in C++ but really need to get it to work in Fortran.
I have the following module defined in Fortran:
module adopov
integer :: indexcount
integer, parameter :: tape_size = 1000
!
!....... ADtype
public :: ADtype
type ADtype
integer :: index = -1
real :: v = 0.0
!
contains
procedure :: oo_asg
generic, public :: assignment(=) => oo_asg
end type ADtype
!
!....... class tape
public :: ADtape
type ADtape
real :: v = 0.0
end type ADtape
!
!....... interface(s)
interface assignment(=)
module procedure oo_asg
end interface
!
type (ADtape), dimension(tape_size) :: tape
!
!....... definitions
contains
!
!....... assignment
subroutine oo_asg (x,y)
implicit none
class(ADtype), intent(out) :: x
class(ADtype), intent(in) :: y
!
tape(indexcount)%v = y%v
indexcount = indexcount + 1
x%v = y%v
x%index = indexcount
end subroutine oo_asg
!
end module adopov
In C++, I have a similar user-defined type as
class ADType {
public:
int index;
double v;
ADType() : index(-1), v(0) {};
ADType(const double&);
ADType& operator=(const ADType&);
};
where the constructor sets the initial values for the index and value parts. Next, I have a constructor for passive values or constants (of type double) so that I can define a new variable of class (ADType) whenever I have a double variable. For example, when I have:
ADType x;
x = 2.0;
initially a new variable of type ADType is created with value set to 2.0, let's say var1 = 2.0 and next (according to the assignment operator (=) defined in the class ADType) I will assign that variable to x, i.e. x = var1. This entire process is being recorded in a tape that counts operations and records the values and indices.
Now, you may say "why do you have to do this?". Well, during the adjoint method of automatic differentiation using operator overloading, this is a necessary step.
The way I do it in C++ is that I simply have the following two constructors:
ADType:: ADType(const double& x): v(x) {
tape[indexcounter].v = x;
indexcounter++;
};
ADType& ADType::operator=(const ADType& x) {
if (this==&x) return *this;
tape[indexcounter].v = v = x.v;
indexcounter++;
return *this;
}
but I don't know how to implement the constructor for passive values and constants in Fortran.
You have two options, that you might use in combination.
Overload the defined assignment with a procedure that corresponds to having a REAL object on the right hand side.
TYPE ADTYPE
...
CONTAINS
PROCEDURE :: OO_ASG
PROCEDURE :: ASSIGN_FROM_REAL
GENERIC, PUBLIC :: ASSIGNMENT(=) => OO_ASG, ASSIGN_FROM_REAL
END TYPE AREAL
! Users of the module shoudn't be calling this procedure
! (or the specific binding, really) directly.
PRIVATE :: ASSIGN_FROM_REAL
...
SUBROUTINE ASSIGN_FROM_REAL(x,y)
CLASS(ADTYPE), INTENT(OUT) :: x
REAL, INTENT(IN) :: y
! Do what you have to do...
...
x%V = y
END SUBROUTINE ASSIGN_FROM_REAL
! Example of use...
TYPE(ADTYPE) :: z
z = 2.0
Use a structure constructor, or the overloaded procedural equivalent.
INTERFACE ADTYPE
MODULE PROCEDURE ADTYPE_construct_from_REAL
END INTERFACE ADTYPE
PRIVATE :: ADTYPE_construct_from_REAL
...
FUNCTION ADTYPE_construct_from_REAL(y) RESULT(x)
REAL, INTENT(IN) :: y
TYPE(ADTYPE) :: x
! Do what you have to do.
...
x%V = y
END FUNCTION ADTYPE_construct_from_REAL
! Example of use:
TYPE(ADTYPE) :: z
z = ADTYPE(3.0)
If you made the invocation of the constructor that takes a double explicit in the source code in your C++ example (i.e ADType x; x = ADType(2.0);, then you have the equivalent to the second of these two approaches - Fortran does not have implicit conversions between objects of derived type.
(Your example code shows both a type bound assignment and a Fortran 90 style stand-alone interface. It does not make sense to have both - that module should not even compile.)
Here's a complete working suggestion to your problem. Please note that every variable inside a module automatically inherits the save attribute. If you're eventually interested in concurrency you may have to enclose indexcounter inside ADtape with an appropriate type-bound procedure for bookkeeping.
module adopov
use, intrinsic :: ISO_C_binding, only: &
ip => C_INT, &
wp => C_DOUBLE
! Explicit typing only
implicit none
! Everything is private unless stated otherwise
private
public :: Adtype, wp
! Declare derived data types
type ADtape
real(wp) :: v = 0.0_wp
end type ADtape
type, public :: ADtype
integer(ip) :: index = -1
real(wp) :: v = 0.0_wp
contains
procedure, private :: asgn_from_type, asgn_from_real, asgn_from_int
generic, public :: assignment(=) => asgn_from_type, asgn_from_real, asgn_from_int
end type ADtype
! Set user-defined constructor
interface ADtype
module procedure :: ADtype_constructor
end interface ADtype
! Variables confined to the module.
! Please note that every variable
! in a module implicitly inherits the save attribute.
integer :: indexcount = 0 ! Your original code left this uninitialized
integer, parameter :: TAPE_SIZE = 1000
type (ADtape) :: tape(TAPE_SIZE)
contains
pure function ADtype_constructor(x, indx) result (return_value)
real (wp), intent(in) :: x
integer (ip), intent (in), optional :: indx
type (ADtype) :: return_value
return_value%v = x
if (present(indx)) return_value%index = indx
end function ADtype_constructor
subroutine update_tape(float)
real (wp), intent (in) :: float
tape(indexcount)%v = float
indexcount = indexcount + 1
end subroutine update_tape
subroutine asgn_from_type(this, y_type)
class(ADtype), intent(out) :: this
class(ADtype), intent(in) :: y_type
associate( &
v => this%v, &
indx => this%index &
)
call update_tape(y_type%v)
v = y_type%v
indx = indexcount
end associate
end subroutine asgn_from_type
subroutine asgn_from_real(this, y_real)
class(ADtype), intent(out) :: this
real(wp), intent(in) :: y_real
associate( &
v => this%v, &
indx => this%index &
)
call update_tape(y_real)
v = y_real
indx = indexcount
end associate
end subroutine asgn_from_real
subroutine asgn_from_int(this, y_int)
class(ADtype), intent(out) :: this
integer(ip), intent(in) :: y_int
associate( &
v => this%v, &
indx => this%index, &
float => real(y_int, kind=wp) &
)
call update_tape(float)
v = float
indx = indexcount
end associate
end subroutine asgn_from_int
end module adopov
program main
use, intrinsic :: ISO_Fortran_env, only: &
stdout => OUTPUT_UNIT, &
compiler_version, &
compiler_options
use adopov, only: &
ADtype, wp
! Explicit typing only
implicit none
type(ADtype) :: foo, bar, woo
! Invoke the user-defined constructor
foo = ADtype(42.0_wp)
bar = ADtype(42.0_wp, -6)
woo = foo
print *, foo
print *, bar
print *, woo
write( stdout, '(/4a/)') &
' This file was compiled using ', compiler_version(), &
' using the options ', compiler_options()
end program main
This yields
gfortran -Wall -o main.exe adopov.f90 main.f90
./main.exe
1 42.000000000000000
2 42.000000000000000
3 42.000000000000000
This file was compiled using GCC version 6.1.1 20160802 using the options -mtune=generic -march=x86-64 -Wall
