3
3
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4
4
#:set RC_KINDS_TYPES = R_KINDS_TYPES + C_KINDS_TYPES
5
5
6
- #:def cnjg(type,expression)
7
- #:if 'complex' in type
8
- conjg(${expression}$)
9
- #:else
10
- ${expression}$
11
- #:endif
12
- #:enddef
13
-
14
6
! This module is based on https://github.com/jalvesz/fast_math
15
7
module stdlib_intrinsics
16
8
!!Replacement for certain Fortran intrinsic functions offering either faster and/or more accurate implementations.
@@ -21,30 +13,52 @@ module stdlib_intrinsics
21
13
22
14
interface fsum
23
15
#:for rk, rt, rs in RC_KINDS_TYPES
24
- module procedure :: fsum_1d_${rs}$
25
- module procedure :: fsum_1d_${rs}$_mask
16
+ pure module function fsum_1d_${rs}$(a) result(s)
17
+ ${rt}$, intent(in) :: a(:)
18
+ ${rt}$ :: s
19
+ end function
20
+ pure module function fsum_1d_${rs}$_mask(a,mask) result(s)
21
+ ${rt}$, intent(in) :: a(:)
22
+ logical, intent(in) :: mask(:)
23
+ ${rt}$ :: s
24
+ end function
26
25
#:endfor
27
26
end interface
28
27
public :: fsum
29
28
30
29
interface fsum_kahan
31
30
#:for rk, rt, rs in RC_KINDS_TYPES
32
- module procedure :: fsum_kahan_1d_${rs}$
33
- module procedure :: fsum_kahan_1d_${rs}$_mask
31
+ pure module function fsum_kahan_1d_${rs}$(a) result(s)
32
+ ${rt}$, intent(in) :: a(:)
33
+ ${rt}$ :: s
34
+ end function
35
+ pure module function fsum_kahan_1d_${rs}$_mask(a,mask) result(s)
36
+ ${rt}$, intent(in) :: a(:)
37
+ logical, intent(in) :: mask(:)
38
+ ${rt}$ :: s
39
+ end function
34
40
#:endfor
35
41
end interface
36
42
public :: fsum_kahan
37
43
38
44
interface fprod
39
45
#:for rk, rt, rs in RC_KINDS_TYPES
40
- module procedure :: fprod_${rs}$
46
+ pure module function fprod_${rs}$(a,b) result(p)
47
+ ${rt}$, intent(in) :: a(:)
48
+ ${rt}$, intent(in) :: b(:)
49
+ ${rt}$ :: p
50
+ end function
41
51
#:endfor
42
52
end interface
43
53
public :: fprod
44
54
45
55
interface fprod_kahan
46
56
#:for rk, rt, rs in RC_KINDS_TYPES
47
- module procedure :: fprod_kahan_${rs}$
57
+ pure module function fprod_kahan_${rs}$(a,b) result(p)
58
+ ${rt}$, intent(in) :: a(:)
59
+ ${rt}$, intent(in) :: b(:)
60
+ ${rt}$ :: p
61
+ end function
48
62
#:endfor
49
63
end interface
50
64
public :: fprod_kahan
@@ -55,174 +69,27 @@ module stdlib_intrinsics
55
69
module procedure :: vkahan_m_${rs}$
56
70
#:endfor
57
71
end interface
58
-
59
- #:for k1, t1, s1 in (R_KINDS_TYPES)
60
- ${t1}$, parameter :: zero_${s1}$ = 0._${k1}$
61
- #:endfor
62
- #:for k1, t1, s1 in (C_KINDS_TYPES)
63
- ${t1}$, parameter :: zero_${s1}$ = (0._${k1}$,0._${k1}$)
64
- #:endfor
65
- integer, parameter :: chunk = 64
72
+ public :: vkahan
66
73
67
74
contains
68
75
69
- #:for k1, t1, s1 in RC_KINDS_TYPES
70
- pure function fsum_1d_${s1}$(a) result(s)
71
- ${t1}$, intent(in) :: a(:)
72
- ${t1}$ :: s
73
- ${t1}$ :: abatch(chunk)
74
- integer :: i, dr, rr
75
- ! -----------------------------
76
- dr = size(a)/chunk
77
- rr = size(a) - dr*chunk
78
-
79
- abatch = zero_${s1}$
80
- do i = 1, dr
81
- abatch(1:chunk) = abatch(1:chunk) + a(chunk*i-chunk+1:chunk*i)
82
- end do
83
- abatch(1:rr) = abatch(1:rr) + a(size(a)-rr+1:size(a))
84
-
85
- s = zero_${s1}$
86
- do i = 1, chunk/2
87
- s = s + abatch(i)+abatch(chunk/2+i)
88
- end do
89
- end function
90
-
91
- pure function fsum_1d_${s1}$_mask(a,mask) result(s)
92
- ${t1}$, intent(in) :: a(:)
93
- logical, intent(in) :: mask(:)
94
- ${t1}$ :: s
95
- ${t1}$ :: abatch(chunk)
96
- integer :: i, dr, rr
97
- ! -----------------------------
98
- dr = size(a)/chunk
99
- rr = size(a) - dr*chunk
100
-
101
- abatch = zero_${s1}$
102
- do i = 1, dr
103
- abatch(1:chunk) = abatch(1:chunk) + merge( zero_${s1}$ , a(chunk*i-chunk+1:chunk*i) , mask(chunk*i-chunk+1:chunk*i) )
104
- end do
105
- abatch(1:rr) = abatch(1:rr) + merge( zero_${s1}$ , a(size(a)-rr+1:size(a)) , mask(size(a)-rr+1:size(a)) )
106
-
107
- s = zero_${s1}$
108
- do i = 1, chunk/2
109
- s = s + abatch(i)+abatch(chunk/2+i)
110
- end do
111
- end function
112
- #:endfor
113
-
114
- #:for k1, t1, s1 in RC_KINDS_TYPES
115
- pure function fsum_kahan_1d_${s1}$(a) result(s)
116
- ${t1}$, intent(in) :: a(:)
117
- ${t1}$ :: s
118
- ${t1}$ :: sbatch(chunk)
119
- ${t1}$ :: cbatch(chunk)
120
- integer :: i, dr, rr
121
- ! -----------------------------
122
- dr = size(a)/(chunk)
123
- rr = size(a) - dr*chunk
124
- sbatch = zero_${s1}$
125
- cbatch = zero_${s1}$
126
- do i = 1, dr
127
- call vkahan( a(chunk*i-chunk+1:chunk*i) , sbatch(1:chunk) , cbatch(1:chunk) )
128
- end do
129
- call vkahan( a(size(a)-rr+1:size(a)) , sbatch(1:rr) , cbatch(1:rr) )
130
-
131
- s = zero_${s1}$
132
- do i = 1,chunk
133
- call vkahan( sbatch(i) , s , cbatch(i) )
134
- end do
135
- end function
136
-
137
- pure function fsum_kahan_1d_${s1}$_mask(a,mask) result(s)
138
- ${t1}$, intent(in) :: a(:)
139
- logical, intent(in) :: mask(:)
140
- ${t1}$ :: s
141
- ${t1}$ :: sbatch(chunk)
142
- ${t1}$ :: cbatch(chunk)
143
- integer :: i, dr, rr
144
- ! -----------------------------
145
- dr = size(a)/(chunk)
146
- rr = size(a) - dr*chunk
147
- sbatch = zero_${s1}$
148
- cbatch = zero_${s1}$
149
- do i = 1, dr
150
- call vkahan( a(chunk*i-chunk+1:chunk*i) , sbatch(1:chunk) , cbatch(1:chunk) , mask(chunk*i-chunk+1:chunk*i) )
151
- end do
152
- call vkahan( a(size(a)-rr+1:size(a)) , sbatch(1:rr) , cbatch(1:rr) , mask(size(a)-rr+1:size(a)) )
153
-
154
- s = zero_${s1}$
155
- do i = 1,chunk
156
- call vkahan( sbatch(i) , s , cbatch(i) )
157
- end do
158
- end function
159
- #:endfor
160
-
161
- #:for k1, t1, s1 in RC_KINDS_TYPES
162
- pure function fprod_${s1}$(a,b) result(p)
163
- ${t1}$, intent(in) :: a(:)
164
- ${t1}$, intent(in) :: b(:)
165
- ${t1}$ :: p
166
- ${t1}$ :: abatch(chunk)
167
- integer :: i, dr, rr
168
- ! -----------------------------
169
- dr = size(a)/chunk
170
- rr = size(a) - dr*chunk
171
-
172
- abatch = zero_${s1}$
173
- do i = 1, dr
174
- abatch(1:chunk) = abatch(1:chunk) + a(chunk*i-chunk+1:chunk*i)*${cnjg(t1,'b(chunk*i-chunk+1:chunk*i)')}$
175
- end do
176
- abatch(1:rr) = abatch(1:rr) + a(size(a)-rr+1:size(a))*${cnjg(t1,'b(size(a)-rr+1:size(a))')}$
177
-
178
- p = zero_${s1}$
179
- do i = 1, chunk/2
180
- p = p + abatch(i)+abatch(chunk/2+i)
181
- end do
182
- end function
183
-
184
- pure function fprod_kahan_${s1}$(a,b) result(p)
185
- ${t1}$, intent(in) :: a(:)
186
- ${t1}$, intent(in) :: b(:)
187
- ${t1}$ :: p
188
- ${t1}$ :: pbatch(chunk)
189
- ${t1}$ :: cbatch(chunk)
190
- integer :: i, dr, rr
191
- ! -----------------------------
192
- dr = size(a)/(chunk)
193
- rr = size(a) - dr*chunk
194
- pbatch = zero_${s1}$
195
- cbatch = zero_${s1}$
196
- do i = 1, dr
197
- call vkahan( a(chunk*i-chunk+1:chunk*i)*${cnjg(t1,'b(chunk*i-chunk+1:chunk*i)')}$ , pbatch(1:chunk) , cbatch(1:chunk) )
198
- end do
199
- call vkahan( a(size(a)-rr+1:size(a))*${cnjg(t1,'b(size(a)-rr+1:size(a))')}$ , pbatch(1:rr) , cbatch(1:rr) )
200
-
201
- p = zero_${s1}$
202
- do i = 1,chunk
203
- call vkahan( pbatch(i) , p , cbatch(i) )
204
- end do
205
- end function
206
-
207
- #:endfor
208
-
209
- #:for k1, t1, s1 in RC_KINDS_TYPES
210
- elemental subroutine vkahan_${s1}$(a,s,c)
211
- ${t1}$, intent(in) :: a
212
- ${t1}$, intent(inout) :: s
213
- ${t1}$, intent(inout) :: c
214
- ${t1}$ :: t, y
76
+ #:for rk, rt, rs in RC_KINDS_TYPES
77
+ elemental subroutine vkahan_${rs}$(a,s,c)
78
+ ${rt}$, intent(in) :: a
79
+ ${rt}$, intent(inout) :: s
80
+ ${rt}$, intent(inout) :: c
81
+ ${rt}$ :: t, y
215
82
y = a - c
216
83
t = s + y
217
84
c = (t - s) - y
218
85
s = t
219
86
end subroutine
220
- elemental subroutine vkahan_m_${s1 }$(a,s,c,m)
221
- ${t1 }$, intent(in) :: a
222
- ${t1 }$, intent(inout) :: s
223
- ${t1 }$, intent(inout) :: c
87
+ elemental subroutine vkahan_m_${rs }$(a,s,c,m)
88
+ ${rt }$, intent(in) :: a
89
+ ${rt }$, intent(inout) :: s
90
+ ${rt }$, intent(inout) :: c
224
91
logical, intent(in) :: m
225
- ${t1 }$ :: t, y
92
+ ${rt }$ :: t, y
226
93
y = a - c
227
94
t = s + y
228
95
c = (t - s) - y
0 commit comments