1 |
|
2 |
module active_module |
3 |
use w2f__types |
4 |
implicit none |
5 |
private |
6 |
public :: active, saxpy, sax, setderiv, zero_deriv, & |
7 |
&convert_p2a_scalar, convert_a2p_scalar, & |
8 |
&convert_p2a_vector, convert_a2p_vector, & |
9 |
&convert_p2a_matrix, convert_a2p_matrix |
10 |
|
11 |
|
12 |
! |
13 |
! active needs to be a sequence type |
14 |
! with no initialization |
15 |
! |
16 |
type active |
17 |
sequence |
18 |
double precision :: v |
19 |
! initialization does not work for active variables |
20 |
! inside of common block, such as in boxmodel |
21 |
! initialization is required for correct adjoint |
22 |
double precision :: d=0.0 |
23 |
! double precision :: d |
24 |
end type active |
25 |
|
26 |
interface saxpy |
27 |
module procedure saxpy_a_a |
28 |
end interface |
29 |
|
30 |
interface setderiv |
31 |
module procedure setderiv_a_a |
32 |
end interface |
33 |
|
34 |
interface zero_deriv |
35 |
module procedure zero_deriv_a |
36 |
end interface |
37 |
|
38 |
interface sax |
39 |
module procedure sax_d_a_a, sax_i_a_a |
40 |
end interface |
41 |
|
42 |
interface convert_p2a_scalar |
43 |
module procedure convert_p2a_scalar_impl |
44 |
end interface |
45 |
interface convert_a2p_scalar |
46 |
module procedure convert_a2p_scalar_impl |
47 |
end interface |
48 |
interface convert_p2a_vector |
49 |
module procedure convert_sp2a_vector_impl |
50 |
end interface |
51 |
interface convert_a2p_vector |
52 |
module procedure convert_a2sp_vector_impl |
53 |
end interface |
54 |
interface convert_p2a_matrix |
55 |
module procedure convert_sp2a_matrix_impl |
56 |
end interface |
57 |
interface convert_a2p_matrix |
58 |
module procedure convert_a2sp_matrix_impl |
59 |
end interface |
60 |
|
61 |
contains |
62 |
|
63 |
! |
64 |
! chain rule saxpy to be used in forward and reverse modes |
65 |
! |
66 |
|
67 |
subroutine saxpy_a_a(a,x,y) |
68 |
double precision, intent(in) :: a |
69 |
type(active), intent(in) :: x |
70 |
type(active), intent(inout) :: y |
71 |
|
72 |
y%d=y%d+x%d*a |
73 |
end subroutine saxpy_a_a |
74 |
|
75 |
! |
76 |
! chain rule saxpy to be used in forward and reverse modes |
77 |
! derivative component of y is equal to zero initially |
78 |
! note: y needs to be inout as otherwise value component gets |
79 |
! zeroed out |
80 |
! |
81 |
|
82 |
subroutine sax_d_a_a(a,x,y) |
83 |
double precision, intent(in) :: a |
84 |
type(active), intent(in) :: x |
85 |
type(active), intent(inout) :: y |
86 |
|
87 |
y%d=x%d*a |
88 |
end subroutine sax_d_a_a |
89 |
|
90 |
subroutine sax_i_a_a(a,x,y) |
91 |
integer(kind=w2f__i8), intent(in) :: a |
92 |
type(active), intent(in) :: x |
93 |
type(active), intent(inout) :: y |
94 |
|
95 |
y%d=x%d*a |
96 |
end subroutine sax_i_a_a |
97 |
|
98 |
! |
99 |
! set derivative of y to be equal to derivative of x |
100 |
! note: making y inout allows for already existing active |
101 |
! variables to become the target of a derivative assignment |
102 |
! |
103 |
|
104 |
subroutine setderiv_a_a(y,x) |
105 |
type(active), intent(inout) :: y |
106 |
type(active), intent(in) :: x |
107 |
|
108 |
y%d=x%d |
109 |
end subroutine setderiv_a_a |
110 |
|
111 |
! |
112 |
! set derivative components to 0.0 |
113 |
! |
114 |
subroutine zero_deriv_a(x) |
115 |
type(active), intent(out) :: x |
116 |
|
117 |
x%d=0.0d0 |
118 |
end subroutine zero_deriv_a |
119 |
|
120 |
subroutine convert_a2p_scalar_impl(convertTo, convertFrom) |
121 |
double precision, intent(out) :: convertTo |
122 |
type(active), intent(in) :: convertFrom |
123 |
convertTo=convertFrom%v |
124 |
end subroutine |
125 |
|
126 |
subroutine convert_p2a_scalar_impl(convertTo, convertFrom) |
127 |
double precision, intent(in) :: convertFrom |
128 |
type(active), intent(out) :: convertTo |
129 |
convertTo%v=convertFrom |
130 |
end subroutine |
131 |
|
132 |
subroutine convert_a2sp_vector_impl(convertTo, convertFrom) |
133 |
type(active), dimension(:), intent(in) :: convertFrom |
134 |
real, dimension(:), intent(out) :: convertTo |
135 |
integer i |
136 |
do i=lbound(convertFrom,1),ubound(convertFrom,1) |
137 |
convertTo(i)=convertFrom(i)%v |
138 |
end do |
139 |
end subroutine |
140 |
|
141 |
subroutine convert_sp2a_vector_impl(convertTo, convertFrom) |
142 |
real, dimension(:), intent(in) :: convertFrom |
143 |
type(active), dimension(:), intent(out) :: convertTo |
144 |
integer i |
145 |
do i=lbound(convertFrom,1),ubound(convertFrom,1) |
146 |
convertTo(i)%v=convertFrom(i) |
147 |
end do |
148 |
end subroutine |
149 |
|
150 |
subroutine convert_a2sp_matrix_impl(convertTo, convertFrom) |
151 |
type(active), dimension(:,:), intent(in) :: convertFrom |
152 |
real, dimension(:,:), intent(out) :: convertTo |
153 |
integer i,j |
154 |
do i=lbound(convertFrom,1),ubound(convertFrom,1) |
155 |
do j=lbound(convertFrom,2),ubound(convertFrom,2) |
156 |
convertTo(i,j)=convertFrom(i,j)%v |
157 |
end do |
158 |
end do |
159 |
end subroutine |
160 |
|
161 |
subroutine convert_sp2a_matrix_impl(convertTo, convertFrom) |
162 |
real, dimension(:,:), intent(in) :: convertFrom |
163 |
type(active), dimension(:,:), intent(out) :: convertTo |
164 |
integer i,j |
165 |
do i=lbound(convertFrom,1),ubound(convertFrom,1) |
166 |
do j=lbound(convertFrom,2),ubound(convertFrom,2) |
167 |
convertTo(i,j)%v=convertFrom(i,j) |
168 |
end do |
169 |
end do |
170 |
end subroutine |
171 |
|
172 |
end module |
173 |
|
174 |
|