Skip to content

Commit 19fafd4

Browse files
weslleyspereiraweslleyspereira
committed
Test the Fortran intrinsic ABS and complex divisions and report during compile time. Currently only working with the CMake build!
1 parent 79aa0f2 commit 19fafd4

File tree

6 files changed

+641
-0
lines changed

6 files changed

+641
-0
lines changed

‎BLAS/SRC/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -105,3 +105,7 @@ set_target_properties(
105105
SOVERSION ${LAPACK_MAJOR_VERSION}
106106
)
107107
lapack_install_library(${BLASLIB})
108+
109+
if( TEST_FORTRAN_COMPILER )
110+
add_dependencies( ${BLASLIB} run_test_zcomplexabs run_test_zcomplexdiv )
111+
endif()

‎CMakeLists.txt‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,26 @@ if(_is_coverage_build)
4141
find_package(codecov)
4242
endif()
4343

44+
# By default test Fortran compiler complex abs and complex division
45+
option(TEST_FORTRAN_COMPILER "Test Fortran compiler complex abs and complex division" ON)
46+
if( TEST_FORTRAN_COMPILER )
47+
48+
add_executable( test_zcomplexabs ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexabs.f )
49+
add_custom_target( run_test_zcomplexabs
50+
COMMAND test_zcomplexabs 2> test_zcomplexabs.err
51+
WORKING_DIRECTORY ${LAPACK_BINARY_DIR}/INSTALL
52+
COMMENT "Running test_zcomplexabs in ${LAPACK_BINARY_DIR}/INSTALL with stderr: test_zcomplexabs.err"
53+
SOURCES ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexabs.f )
54+
55+
add_executable( test_zcomplexdiv ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexdiv.f )
56+
add_custom_target( run_test_zcomplexdiv
57+
COMMAND test_zcomplexdiv 2> test_zcomplexdiv.err
58+
WORKING_DIRECTORY ${LAPACK_BINARY_DIR}/INSTALL
59+
COMMENT "Running test_zcomplexdiv in ${LAPACK_BINARY_DIR}/INSTALL with stderr: test_zcomplexdiv.err"
60+
SOURCES ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexdiv.f )
61+
62+
endif()
63+
4464
# By default static library
4565
option(BUILD_SHARED_LIBS "Build shared libraries" OFF)
4666

‎INSTALL/CMakeLists.txt‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,3 +7,10 @@ add_executable(secondtst_INT_ETIME second_INT_ETIME.f secondtst.f)
77
add_executable(secondtst_INT_CPU_TIME second_INT_CPU_TIME.f secondtst.f)
88
add_executable(testieee tstiee.f)
99
add_executable(testversion ilaver.f LAPACK_version.f)
10+
11+
if( NOT TARGET test_zcomplexabs )
12+
add_executable( test_zcomplexabs test_zcomplexabs.f )
13+
endif()
14+
if( NOT TARGET test_zcomplexdiv )
15+
add_executable( test_zcomplexdiv test_zcomplexdiv.f )
16+
endif()

‎INSTALL/test_zcomplexabs.f‎

Lines changed: 275 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,275 @@
1+
*> \brief zabs tests the robustness and precision of the intrinsic ABS for double complex
2+
*> \author Weslley S. Pereira, University of Colorado Denver, U.S.
3+
*
4+
*> \verbatim
5+
*>
6+
*> Real values for test:
7+
*> (1) x = 2**m, where m = MINEXPONENT-DIGITS, ..., MINEXPONENT-1. Stop on the first success.
8+
*> Mind that not all platforms might implement subnormal numbers.
9+
*> (2) x = 2**m, where m = MINEXPONENT, ..., 0. Stop on the first success.
10+
*> (3) x = OV, where OV is the overflow threshold. OV^2 overflows but the norm is OV.
11+
*> (4) x = 2**m, where m = MAXEXPONENT-1, ..., 1. Stop on the first success.
12+
*>
13+
*> Tests:
14+
*> (a) y = x + 0 * I, |y| = x
15+
*> (b) y = 0 + x * I, |y| = x
16+
*> (c) y = (3/4)*x + x * I, |y| = (5/4)*x whenever (3/4)*x and (5/4)*x can be exactly stored
17+
*> (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2) whenever (1/2)*x can be exactly stored
18+
*>
19+
*> Special cases:
20+
*>
21+
*> (i) Inf propagation
22+
*> (1) y = Inf + 0 * I, |y| is Inf.
23+
*> (2) y =-Inf + 0 * I, |y| is Inf.
24+
*> (3) y = 0 + Inf * I, |y| is Inf.
25+
*> (4) y = 0 - Inf * I, |y| is Inf.
26+
*> (5) y = Inf + Inf * I, |y| is Inf.
27+
*>
28+
*> (n) NaN propagation
29+
*> (1) y = NaN + 0 * I, |y| is NaN.
30+
*> (2) y = 0 + NaN * I, |y| is NaN.
31+
*> (3) y = NaN + NaN * I, |y| is NaN.
32+
*>
33+
*> \endverbatim
34+
*
35+
program zabs
36+
37+
logical debug
38+
parameter ( debug = .false. )
39+
40+
integer N, i, nNaN, nInf, min, Max, m
41+
parameter ( N = 4, nNaN = 3, nInf = 5 )
42+
43+
double precision X( N ), R, threeFourth, fiveFourth, answerC,
44+
$ answerD, oneHalf, aInf, aNaN, relDiff, b,
45+
$ eps, blueMin, blueMax, Xj, stepX(N), limX(N)
46+
parameter ( threeFourth = 3.0d0 / 4,
47+
$ fiveFourth = 5.0d0 / 4,
48+
$ oneHalf = 1.0d0 / 2 )
49+
50+
double complex Y, cInf( nInf ), cNaN( nNaN )
51+
intrinsic ABS, DBLE, RADIX, CEILING, TINY, DIGITS, SQRT,
52+
$ MAXEXPONENT, MINEXPONENT, FLOOR, HUGE, DCMPLX,
53+
$ EPSILON
54+
55+
integer subnormalTreatedAs0, caseAFails, caseBFails,
56+
$ caseCFails, caseDFails
57+
*
58+
subnormalTreatedAs0 = 0
59+
caseAFails = 0
60+
caseBFails = 0
61+
caseCFails = 0
62+
caseDFails = 0
63+
*
64+
min = MINEXPONENT(0.0d0)
65+
Max = MAXEXPONENT(0.0d0)
66+
m = DIGITS(0.0d0)
67+
b = DBLE(RADIX(0.0d0))
68+
eps = EPSILON(0.0d0)
69+
blueMin = b**CEILING( (min - 1) * 0.5d0 )
70+
blueMax = b**FLOOR( (Max - m + 1) * 0.5d0 )
71+
*
72+
X(1) = TINY(0.0d0) * b**( DBLE(1-m) )
73+
X(2) = TINY(0.0d0)
74+
X(3) = HUGE(0.0d0)
75+
X(4) = b**( DBLE(Max-1) )
76+
*
77+
stepX(1) = 2.0
78+
stepX(2) = 2.0
79+
stepX(3) = 0.0
80+
stepX(4) = 0.5
81+
*
82+
limX(1) = X(2)
83+
limX(2) = 1.0
84+
limX(3) = 0.0
85+
limX(4) = 2.0
86+
*
87+
if( debug ) then
88+
print *, '# X :=', X
89+
print *, '# Blue min constant :=', blueMin
90+
print *, '# Blue max constant :=', blueMax
91+
endif
92+
*
93+
Xj = X(1)
94+
if( Xj .eq. 0.0d0 ) then
95+
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
96+
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
97+
print *, "!! fl( subnormal ) may be 0"
98+
endif
99+
else
100+
do 100 i = 1, N
101+
Xj = X(i)
102+
if( Xj .eq. 0.0d0 ) then
103+
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
104+
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
105+
print *, "!! fl( subnormal ) may be 0"
106+
endif
107+
endif
108+
100 continue
109+
endif
110+
*
111+
aInf = X(3) * 2
112+
cInf(1) = DCMPLX( aInf, 0.0d0 )
113+
cInf(2) = DCMPLX(-aInf, 0.0d0 )
114+
cInf(3) = DCMPLX( 0.0d0, aInf )
115+
cInf(4) = DCMPLX( 0.0d0,-aInf )
116+
cInf(5) = DCMPLX( aInf, aInf )
117+
*
118+
aNaN = aInf / aInf
119+
cNaN(1) = DCMPLX( aNaN, 0.0d0 )
120+
cNaN(2) = DCMPLX( 0.0d0, aNaN )
121+
cNaN(3) = DCMPLX( aNaN, aNaN )
122+
*
123+
* Test (a) y = x + 0 * I, |y| = x
124+
do 10 i = 1, N
125+
Xj = X(i)
126+
if( Xj .eq. 0.0d0 ) then
127+
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
128+
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
129+
print *, "!! [a] fl( subnormal ) may be 0"
130+
endif
131+
else
132+
do while( Xj .ne. limX(i) )
133+
Y = DCMPLX( Xj, 0.0d0 )
134+
R = ABS( Y )
135+
if( R .ne. Xj ) then
136+
caseAFails = caseAFails + 1
137+
if( caseAFails .eq. 1 ) then
138+
print *, "!! Some ABS(x+0*I) differ from ABS(x)"
139+
endif
140+
WRITE( 0, FMT = 9999 ) 'a',i, Xj, '(1+0*I)', R, Xj
141+
endif
142+
Xj = Xj * stepX(i)
143+
end do
144+
endif
145+
10 continue
146+
*
147+
* Test (b) y = 0 + x * I, |y| = x
148+
do 20 i = 1, N
149+
Xj = X(i)
150+
if( Xj .eq. 0.0d0 ) then
151+
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
152+
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
153+
print *, "!! [b] fl( subnormal ) may be 0"
154+
endif
155+
else
156+
do while( Xj .ne. limX(i) )
157+
Y = DCMPLX( 0.0d0, Xj )
158+
R = ABS( Y )
159+
if( R .ne. Xj ) then
160+
caseBFails = caseBFails + 1
161+
if( caseBFails .eq. 1 ) then
162+
print *, "!! Some ABS(0+x*I) differ from ABS(x)"
163+
endif
164+
WRITE( 0, FMT = 9999 ) 'b',i, Xj, '(0+1*I)', R, Xj
165+
endif
166+
Xj = Xj * stepX(i)
167+
end do
168+
endif
169+
20 continue
170+
*
171+
* Test (c) y = (3/4)*x + x * I, |y| = (5/4)*x
172+
do 30 i = 1, N
173+
if( i .eq. 3 ) go to 30
174+
if( i .eq. 1 ) then
175+
Xj = 4*X(i)
176+
else
177+
Xj = X(i)
178+
endif
179+
if( Xj .eq. 0.0d0 ) then
180+
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
181+
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
182+
print *, "!! [c] fl( subnormal ) may be 0"
183+
endif
184+
else
185+
do while( Xj .ne. limX(i) )
186+
answerC = fiveFourth * Xj
187+
Y = DCMPLX( threeFourth * Xj, Xj )
188+
R = ABS( Y )
189+
if( R .ne. answerC ) then
190+
caseCFails = caseCFails + 1
191+
if( caseCFails .eq. 1 ) then
192+
print *,
193+
$ "!! Some ABS(x*(3/4+I)) differ from (5/4)*ABS(x)"
194+
endif
195+
WRITE( 0, FMT = 9999 ) 'c',i, Xj, '(3/4+I)', R,
196+
$ answerC
197+
endif
198+
Xj = Xj * stepX(i)
199+
end do
200+
endif
201+
30 continue
202+
*
203+
* Test (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2)
204+
do 40 i = 1, N
205+
if( i .eq. 1 ) then
206+
Xj = 2*X(i)
207+
else
208+
Xj = X(i)
209+
endif
210+
if( Xj .eq. 0.0d0 ) then
211+
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
212+
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
213+
print *, "!! [d] fl( subnormal ) may be 0"
214+
endif
215+
else
216+
do while( Xj .ne. limX(i) )
217+
answerD = (oneHalf * Xj) * SQRT(2.0d0)
218+
if( answerD .eq. 0.0d0 ) then
219+
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
220+
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
221+
print *, "!! [d] fl( subnormal ) may be 0"
222+
endif
223+
else
224+
Y = DCMPLX( oneHalf * Xj, oneHalf * Xj )
225+
R = ABS( Y )
226+
relDiff = ABS(R-answerD)/answerD
227+
if( relDiff .ge. (0.5*eps) ) then
228+
caseDFails = caseDFails + 1
229+
if( caseDFails .eq. 1 ) then
230+
print *,
231+
$ "!! Some ABS(x*(1+I)) differ from sqrt(2)*ABS(x)"
232+
endif
233+
WRITE( 0, FMT = 9999 ) 'd',i, (oneHalf*Xj),
234+
$ '(1+1*I)', R, answerD
235+
endif
236+
endif
237+
Xj = Xj * stepX(i)
238+
end do
239+
endif
240+
40 continue
241+
*
242+
* Test (e) Infs
243+
do 50 i = 1, nInf
244+
Y = cInf(i)
245+
R = ABS( Y )
246+
if( .not.(R .gt. HUGE(0.0d0)) ) then
247+
WRITE( *, FMT = 9997 ) 'i',i, Y, R
248+
endif
249+
50 continue
250+
*
251+
* Test (f) NaNs
252+
do 60 i = 1, nNaN
253+
Y = cNaN(i)
254+
R = ABS( Y )
255+
if( R .eq. R ) then
256+
WRITE( *, FMT = 9998 ) 'n',i, Y, R
257+
endif
258+
60 continue
259+
*
260+
if( (caseAFails .gt. 0) .or. (caseBFails .gt. 0) .or.
261+
$ (caseCFails .gt. 0) .or. (caseDFails .gt. 0) )
262+
$ print *, "# Please check the failed ABS(a+b*I) in [stderr]"
263+
*
264+
9997 FORMAT( '[',A1,I1, '] ABS(', (ES8.1,SP,ES8.1,"*I"), ' ) = ',
265+
$ ES8.1, ' differs from Inf' )
266+
*
267+
9998 FORMAT( '[',A1,I1, '] ABS(', (ES8.1,SP,ES8.1,"*I"), ' ) = ',
268+
$ ES8.1, ' differs from NaN' )
269+
*
270+
9999 FORMAT( '[',A1,I1, '] ABS(', ES24.16E3, ' * ', A7, ' ) = ',
271+
$ ES24.16E3, ' differs from ', ES24.16E3 )
272+
*
273+
* End of zabs
274+
*
275+
END

0 commit comments

Comments
 (0)