Skip to content

Commit 2bc7af9

Browse files
jalveszjalvesz
committed
add fast dot_product and start tests
1 parent 08ec0aa commit 2bc7af9

File tree

5 files changed

+190
-13
lines changed

5 files changed

+190
-13
lines changed

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ set(fppFiles
1616
stdlib_hash_64bit_fnv.fypp
1717
stdlib_hash_64bit_pengy.fypp
1818
stdlib_hash_64bit_spookyv2.fypp
19+
stdlib_intrinsics.fypp
1920
stdlib_io.fypp
2021
stdlib_io_npy.fypp
2122
stdlib_io_npy_load.fypp

‎src/stdlib_intrinsics.fypp

Lines changed: 87 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,32 +1,58 @@
11
#:include "common.fypp"
22
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
33
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4-
#:set RANKS = range(2, 2)
4+
#:set RC_KINDS_TYPES = R_KINDS_TYPES + C_KINDS_TYPES
5+
6+
#:def cnjg(type,expression)
7+
#:if 'complex' in type
8+
conjg(${expression}$)
9+
#:else
10+
${expression}$
11+
#:endif
12+
#:enddef
13+
14+
! This module is based on https://github.com/jalvesz/fast_math
515
module stdlib_intrinsics
16+
!!Replacement for certain Fortran intrinsic functions offering either faster and/or more accurate implementations.
17+
!! ([Specification](../page/specs/stdlib_intrinsics.html))
618
use stdlib_kinds
719
implicit none
820
private
921

1022
interface fsum
11-
#:for rk, rt in REAL_KINDS_TYPES
12-
module procedure :: fsum_1d_${rk}$
13-
module procedure :: fsum_1d_${rk}$_mask
23+
#:for rk, rt, rs in RC_KINDS_TYPES
24+
module procedure :: fsum_1d_${rs}$
25+
module procedure :: fsum_1d_${rs}$_mask
1426
#:endfor
1527
end interface
1628
public :: fsum
1729

1830
interface fsum_kahan
19-
#:for rk, rt in REAL_KINDS_TYPES
20-
module procedure :: fsum_kahan_1d_${rk}$
21-
module procedure :: fsum_kahan_1d_${rk}$_mask
31+
#:for rk, rt, rs in RC_KINDS_TYPES
32+
module procedure :: fsum_kahan_1d_${rs}$
33+
module procedure :: fsum_kahan_1d_${rs}$_mask
2234
#:endfor
2335
end interface
2436
public :: fsum_kahan
2537

38+
interface fprod
39+
#:for rk, rt, rs in RC_KINDS_TYPES
40+
module procedure :: fprod_${rs}$
41+
#:endfor
42+
end interface
43+
public :: fprod
44+
45+
interface fprod_kahan
46+
#:for rk, rt, rs in RC_KINDS_TYPES
47+
module procedure :: fprod_kahan_${rs}$
48+
#:endfor
49+
end interface
50+
public :: fprod_kahan
51+
2652
interface vkahan
27-
#:for rk, rt in REAL_KINDS_TYPES
28-
module procedure :: vkahan_${rk}$
29-
module procedure :: vkahan_m_${rk}$
53+
#:for rk, rt, rs in RC_KINDS_TYPES
54+
module procedure :: vkahan_${rs}$
55+
module procedure :: vkahan_m_${rs}$
3056
#:endfor
3157
end interface
3258

@@ -40,7 +66,7 @@ module stdlib_intrinsics
4066

4167
contains
4268

43-
#:for k1, t1, s1 in R_KINDS_TYPES
69+
#:for k1, t1, s1 in RC_KINDS_TYPES
4470
pure function fsum_1d_${s1}$(a) result(s)
4571
${t1}$, intent(in) :: a(:)
4672
${t1}$ :: s
@@ -85,7 +111,7 @@ pure function fsum_1d_${s1}$_mask(a,mask) result(s)
85111
end function
86112
#:endfor
87113

88-
#:for k1, t1, s1 in R_KINDS_TYPES
114+
#:for k1, t1, s1 in RC_KINDS_TYPES
89115
pure function fsum_kahan_1d_${s1}$(a) result(s)
90116
${t1}$, intent(in) :: a(:)
91117
${t1}$ :: s
@@ -132,7 +158,55 @@ pure function fsum_kahan_1d_${s1}$_mask(a,mask) result(s)
132158
end function
133159
#:endfor
134160

135-
#:for k1, t1, s1 in R_KINDS_TYPES
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
136210
elemental subroutine vkahan_${s1}$(a,s,c)
137211
${t1}$, intent(in) :: a
138212
${t1}$, intent(inout) :: s

‎test/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ add_subdirectory(constants)
1717
add_subdirectory(hash_functions)
1818
add_subdirectory(hash_functions_perf)
1919
add_subdirectory(hashmaps)
20+
add_subdirectory(intrinsics)
2021
add_subdirectory(io)
2122
add_subdirectory(linalg)
2223
add_subdirectory(logger)

‎test/intrinsics/CMakeLists.txt

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
set(
2+
fppFiles
3+
"test_intrinsics.fypp"
4+
)
5+
fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
6+
7+
ADDTEST(intrinsics)

‎test/intrinsics/test_intrinsics.fypp

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
#:include "common.fypp"
2+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
3+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4+
#:set RC_KINDS_TYPES = R_KINDS_TYPES + C_KINDS_TYPES
5+
6+
module test_intrinsics
7+
use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test
8+
use stdlib_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64
9+
use stdlib_intrinsics
10+
use stdlib_math, only: swap
11+
implicit none
12+
13+
contains
14+
15+
!> Collect all exported unit tests
16+
subroutine collect_suite(testsuite)
17+
!> Collection of tests
18+
type(unittest_type), allocatable, intent(out) :: testsuite(:)
19+
20+
testsuite = [ &
21+
new_unittest('sum', test_sum) &
22+
!new_unittest('dot_product', test_dot_product) &
23+
]
24+
end subroutine
25+
26+
subroutine test_sum(error)
27+
!> Error handling
28+
type(error_type), allocatable, intent(out) :: error
29+
30+
!> Internal parameters and variables
31+
integer, parameter :: n = 1e3, ncalc = 3, niter = 100
32+
real(sp) :: u
33+
integer :: iter, i, j
34+
!====================================================================================
35+
#:for k1, t1, s1 in R_KINDS_TYPES
36+
block
37+
${t1}$, allocatable :: x(:)
38+
${t1}$, parameter :: total_sum = 4*atan(1._${k1}$), tolerance = epsilon(1._${k1}$)*100
39+
${t1}$ :: xsum(ncalc), meanval(ncalc), err(ncalc)
40+
41+
allocate(x(n))
42+
do i = 1, n
43+
x(i) = 8*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$)/real(n,kind=${k1}$)**2
44+
end do
45+
! scramble array
46+
do i = 1, n
47+
call random_number(u)
48+
j = 1 + FLOOR(n*u)
49+
call swap( x(i), x(j) )
50+
end do
51+
52+
err(:) = 0._${k1}$
53+
do iter = 1, niter
54+
xsum(1) = sum(x)
55+
xsum(2) = fsum_kahan(x)
56+
xsum(3) = fsum(x)
57+
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-xsum(1:ncalc)/total_sum)
58+
end do
59+
err(1:ncalc) = err(1:ncalc) / niter
60+
print *, err
61+
call check(error, all(err(:)<tolerance) )
62+
if (allocated(error)) return
63+
end block
64+
#:endfor
65+
66+
end subroutine
67+
68+
end module test_intrinsics
69+
70+
program tester
71+
use, intrinsic :: iso_fortran_env, only : error_unit
72+
use testdrive, only : run_testsuite, new_testsuite, testsuite_type
73+
use test_intrinsics, only : collect_suite
74+
implicit none
75+
integer :: stat, is
76+
type(testsuite_type), allocatable :: testsuites(:)
77+
character(len=*), parameter :: fmt = '("#", *(1x, a))'
78+
79+
stat = 0
80+
81+
testsuites = [ &
82+
new_testsuite("sparse", collect_suite) &
83+
]
84+
85+
do is = 1, size(testsuites)
86+
write(error_unit, fmt) "Testing:", testsuites(is)%name
87+
call run_testsuite(testsuites(is)%collect, error_unit, stat)
88+
end do
89+
90+
if (stat > 0) then
91+
write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
92+
error stop
93+
end if
94+
end program

0 commit comments

Comments
 (0)