33#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
44#:set RC_KINDS_TYPES = R_KINDS_TYPES + C_KINDS_TYPES
55
6- #:def cnjg(type,expression)
7- #:if 'complex' in type
8- conjg(${expression}$)
9- #:else
10- ${expression}$
11- #:endif
12- #:enddef
13-
146! This module is based on https://github.com/jalvesz/fast_math
157module stdlib_intrinsics
168 !!Replacement for certain Fortran intrinsic functions offering either faster and/or more accurate implementations.
@@ -21,30 +13,52 @@ module stdlib_intrinsics
2113
2214 interface fsum
2315 #: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
2625 #:endfor
2726 end interface
2827 public :: fsum
2928
3029 interface fsum_kahan
3130 #: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
3440 #:endfor
3541 end interface
3642 public :: fsum_kahan
3743
3844 interface fprod
3945 #: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
4151 #:endfor
4252 end interface
4353 public :: fprod
4454
4555 interface fprod_kahan
4656 #: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
4862 #:endfor
4963 end interface
5064 public :: fprod_kahan
@@ -55,174 +69,27 @@ module stdlib_intrinsics
5569 module procedure :: vkahan_m_${rs}$
5670 #:endfor
5771 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
6673
6774contains
6875
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
21582 y = a - c
21683 t = s + y
21784 c = (t - s) - y
21885 s = t
21986end 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
22491 logical, intent(in) :: m
225- ${t1 }$ :: t, y
92+ ${rt }$ :: t, y
22693 y = a - c
22794 t = s + y
22895 c = (t - s) - y
0 commit comments