LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
test_zcomplexabs.f
Go to the documentation of this file.
1*> \brief zabs tests the robustness and precision of the intrinsic ABS for double complex
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \author Weslley S. Pereira, University of Colorado Denver, U.S.
9*
10*> \verbatim
11*>
12*> Real values for test:
13*> (1) x = 2**m, where m = MINEXPONENT-DIGITS, ..., MINEXPONENT-1. Stop on the first success.
14*> Mind that not all platforms might implement subnormal numbers.
15*> (2) x = 2**m, where m = MINEXPONENT, ..., 0. Stop on the first success.
16*> (3) x = OV, where OV is the overflow threshold. OV^2 overflows but the norm is OV.
17*> (4) x = 2**m, where m = MAXEXPONENT-1, ..., 1. Stop on the first success.
18*>
19*> Tests:
20*> (a) y = x + 0 * I, |y| = x
21*> (b) y = 0 + x * I, |y| = x
22*> (c) y = (3/4)*x + x * I, |y| = (5/4)*x whenever (3/4)*x and (5/4)*x can be exactly stored
23*> (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2) whenever (1/2)*x can be exactly stored
24*>
25*> Special cases:
26*>
27*> (i) Inf propagation
28*> (1) y = Inf + 0 * I, |y| is Inf.
29*> (2) y =-Inf + 0 * I, |y| is Inf.
30*> (3) y = 0 + Inf * I, |y| is Inf.
31*> (4) y = 0 - Inf * I, |y| is Inf.
32*> (5) y = Inf + Inf * I, |y| is Inf.
33*>
34*> (n) NaN propagation
35*> (1) y = NaN + 0 * I, |y| is NaN.
36*> (2) y = 0 + NaN * I, |y| is NaN.
37*> (3) y = NaN + NaN * I, |y| is NaN.
38*>
39*> \endverbatim
40*
41*> \ingroup auxOTHERauxiliary
42*
43* =====================================================================
44 program zabs
45*
46* -- LAPACK test routine --
47* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
48
49* ..
50* .. Local parameters ..
51 logical debug
52 parameter( debug = .false. )
53 integer n, nnan, ninf
54 parameter( n = 4, nnan = 3, ninf = 5 )
55 double precision threefourth, fivefourth, onehalf
56 parameter( threefourth = 3.0d0 / 4,
57 $ fivefourth = 5.0d0 / 4,
58 $ onehalf = 1.0d0 / 2 )
59* ..
60* .. Local Variables ..
61 integer i, min, max, m, subnormaltreatedas0,
62 $ caseafails, casebfails, casecfails, casedfails,
63 $ caseefails, caseffails, nfailingtests, ntests
64 double precision x( n ), r, answerc,
65 $ answerd, ainf, anan, reldiff, b,
66 $ eps, bluemin, bluemax, xj, stepx(n), limx(n)
67 double complex y, cinf( ninf ), cnan( nnan )
68*
69* .. Intrinsic Functions ..
70 intrinsic abs, dble, radix, ceiling, tiny, digits, sqrt,
71 $ maxexponent, minexponent, floor, huge, dcmplx,
72 $ epsilon
73
74*
75* .. Initialize error counts ..
76 subnormaltreatedas0 = 0
77 caseafails = 0
78 casebfails = 0
79 casecfails = 0
80 casedfails = 0
81 caseefails = 0
82 caseffails = 0
83 nfailingtests = 0
84 ntests = 0
85*
86* .. Initialize machine constants ..
87 min = minexponent(0.0d0)
88 max = maxexponent(0.0d0)
89 m = digits(0.0d0)
90 b = dble(radix(0.0d0))
91 eps = epsilon(0.0d0)
92 bluemin = b**ceiling( (min - 1) * 0.5d0 )
93 bluemax = b**floor( (max - m + 1) * 0.5d0 )
94*
95* .. Vector X ..
96 x(1) = tiny(0.0d0) * b**( dble(1-m) )
97 x(2) = tiny(0.0d0)
98 x(3) = huge(0.0d0)
99 x(4) = b**( dble(max-1) )
100*
101* .. Then modify X using the step ..
102 stepx(1) = 2.0
103 stepx(2) = 2.0
104 stepx(3) = 0.0
105 stepx(4) = 0.5
106*
107* .. Up to the value ..
108 limx(1) = x(2)
109 limx(2) = 1.0
110 limx(3) = 0.0
111 limx(4) = 2.0
112*
113* .. Inf entries ..
114 ainf = x(3) * 2
115 cinf(1) = dcmplx( ainf, 0.0d0 )
116 cinf(2) = dcmplx(-ainf, 0.0d0 )
117 cinf(3) = dcmplx( 0.0d0, ainf )
118 cinf(4) = dcmplx( 0.0d0,-ainf )
119 cinf(5) = dcmplx( ainf, ainf )
120*
121* .. NaN entries ..
122 anan = ainf / ainf
123 cnan(1) = dcmplx( anan, 0.0d0 )
124 cnan(2) = dcmplx( 0.0d0, anan )
125 cnan(3) = dcmplx( anan, anan )
126
127*
128* .. Tests ..
129*
130 if( debug ) then
131 print *, '# X :=', x
132 print *, '# Blue min constant :=', bluemin
133 print *, '# Blue max constant :=', bluemax
134 endif
135*
136 xj = x(1)
137 if( xj .eq. 0.0d0 ) then
138 subnormaltreatedas0 = subnormaltreatedas0 + 1
139 if( debug .or. subnormaltreatedas0 .eq. 1 ) then
140 print *, "!! fl( subnormal ) may be 0"
141 endif
142 else
143 do 100 i = 1, n
144 xj = x(i)
145 if( xj .eq. 0.0d0 ) then
146 subnormaltreatedas0 = subnormaltreatedas0 + 1
147 if( debug .or. subnormaltreatedas0 .eq. 1 ) then
148 print *, "!! fl( subnormal ) may be 0"
149 endif
150 endif
151 100 continue
152 endif
153*
154* Test (a) y = x + 0 * I, |y| = x
155 do 10 i = 1, n
156 xj = x(i)
157 if( xj .eq. 0.0d0 ) then
158 subnormaltreatedas0 = subnormaltreatedas0 + 1
159 if( debug .or. subnormaltreatedas0 .eq. 1 ) then
160 print *, "!! [a] fl( subnormal ) may be 0"
161 endif
162 else
163 do while( xj .ne. limx(i) )
164 ntests = ntests + 1
165 y = dcmplx( xj, 0.0d0 )
166 r = abs( y )
167 if( r .ne. xj ) then
168 caseafails = caseafails + 1
169 if( caseafails .eq. 1 ) then
170 print *, "!! Some ABS(x+0*I) differ from ABS(x)"
171 endif
172 WRITE( 0, fmt = 9999 ) 'a',i, xj, '(1+0*I)', r, xj
173 endif
174 xj = xj * stepx(i)
175 end do
176 endif
177 10 continue
178*
179* Test (b) y = 0 + x * I, |y| = x
180 do 20 i = 1, n
181 xj = x(i)
182 if( xj .eq. 0.0d0 ) then
183 subnormaltreatedas0 = subnormaltreatedas0 + 1
184 if( debug .or. subnormaltreatedas0 .eq. 1 ) then
185 print *, "!! [b] fl( subnormal ) may be 0"
186 endif
187 else
188 do while( xj .ne. limx(i) )
189 ntests = ntests + 1
190 y = dcmplx( 0.0d0, xj )
191 r = abs( y )
192 if( r .ne. xj ) then
193 casebfails = casebfails + 1
194 if( casebfails .eq. 1 ) then
195 print *, "!! Some ABS(0+x*I) differ from ABS(x)"
196 endif
197 WRITE( 0, fmt = 9999 ) 'b',i, xj, '(0+1*I)', r, xj
198 endif
199 xj = xj * stepx(i)
200 end do
201 endif
202 20 continue
203*
204* Test (c) y = (3/4)*x + x * I, |y| = (5/4)*x
205 do 30 i = 1, n
206 if( i .eq. 3 ) go to 30
207 if( i .eq. 1 ) then
208 xj = 4*x(i)
209 else
210 xj = x(i)
211 endif
212 if( xj .eq. 0.0d0 ) then
213 subnormaltreatedas0 = subnormaltreatedas0 + 1
214 if( debug .or. subnormaltreatedas0 .eq. 1 ) then
215 print *, "!! [c] fl( subnormal ) may be 0"
216 endif
217 else
218 do while( xj .ne. limx(i) )
219 ntests = ntests + 1
220 answerc = fivefourth * xj
221 y = dcmplx( threefourth * xj, xj )
222 r = abs( y )
223 if( r .ne. answerc ) then
224 casecfails = casecfails + 1
225 if( casecfails .eq. 1 ) then
226 print *,
227 $ "!! Some ABS(x*(3/4+I)) differ from (5/4)*ABS(x)"
228 endif
229 WRITE( 0, fmt = 9999 ) 'c',i, xj, '(3/4+I)', r,
230 $ answerc
231 endif
232 xj = xj * stepx(i)
233 end do
234 endif
235 30 continue
236*
237* Test (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2)
238 do 40 i = 1, n
239 if( i .eq. 1 ) then
240 xj = 2*x(i)
241 else
242 xj = x(i)
243 endif
244 if( xj .eq. 0.0d0 ) then
245 subnormaltreatedas0 = subnormaltreatedas0 + 1
246 if( debug .or. subnormaltreatedas0 .eq. 1 ) then
247 print *, "!! [d] fl( subnormal ) may be 0"
248 endif
249 else
250 do while( xj .ne. limx(i) )
251 answerd = (onehalf * xj) * sqrt(2.0d0)
252 if( answerd .eq. 0.0d0 ) then
253 subnormaltreatedas0 = subnormaltreatedas0 + 1
254 if( debug .or. subnormaltreatedas0 .eq. 1 ) then
255 print *, "!! [d] fl( subnormal ) may be 0"
256 endif
257 else
258 ntests = ntests + 1
259 y = dcmplx( onehalf * xj, onehalf * xj )
260 r = abs( y )
261 reldiff = abs(r-answerd)/answerd
262 if( reldiff .ge. (0.5*eps) ) then
263 casedfails = casedfails + 1
264 if( casedfails .eq. 1 ) then
265 print *,
266 $ "!! Some ABS(x*(1+I)) differ from sqrt(2)*ABS(x)"
267 endif
268 WRITE( 0, fmt = 9999 ) 'd',i, (onehalf*xj),
269 $ '(1+1*I)', r, answerd
270 endif
271 endif
272 xj = xj * stepx(i)
273 end do
274 endif
275 40 continue
276*
277* Test (e) Infs
278 do 50 i = 1, ninf
279 ntests = ntests + 1
280 y = cinf(i)
281 r = abs( y )
282 if( .not.(r .gt. huge(0.0d0)) ) then
283 caseefails = caseefails + 1
284 WRITE( *, fmt = 9997 ) 'i',i, y, r
285 endif
286 50 continue
287*
288* Test (f) NaNs
289 do 60 i = 1, nnan
290 ntests = ntests + 1
291 y = cnan(i)
292 r = abs( y )
293 if( r .eq. r ) then
294 caseffails = caseffails + 1
295 WRITE( *, fmt = 9998 ) 'n',i, y, r
296 endif
297 60 continue
298*
299* If any test fails, displays a message
300 nfailingtests = caseafails + casebfails + casecfails + casedfails
301 $ + caseefails + caseffails
302 if( nfailingtests .gt. 0 ) then
303 print *, "# ", ntests-nfailingtests, " tests out of ", ntests,
304 $ " pass for ABS(a+b*I),", nfailingtests, " tests fail."
305 else
306 print *, "# All tests pass for ABS(a+b*I)"
307 endif
308*
309* If anything was written to stderr, print the message
310 if( (caseafails .gt. 0) .or. (casebfails .gt. 0) .or.
311 $ (casecfails .gt. 0) .or. (casedfails .gt. 0) ) then
312 print *, "# Please check the failed ABS(a+b*I) in [stderr]"
313 endif
314*
315* .. Formats ..
316 9997 FORMAT( '[',a1,i1, '] ABS(', (es8.1,sp,es8.1,"*I"), ' ) = ',
317 $ es8.1, ' differs from Inf' )
318*
319 9998 FORMAT( '[',a1,i1, '] ABS(', (es8.1,sp,es8.1,"*I"), ' ) = ',
320 $ es8.1, ' differs from NaN' )
321*
322 9999 FORMAT( '[',a1,i1, '] ABS(', es24.16e3, ' * ', a7, ' ) = ',
323 $ es24.16e3, ' differs from ', es24.16e3 )
324*
325* End of zabs
326*
327 END
zabs
program zabs
zabs tests the robustness and precision of the intrinsic ABS for double complex
Definition test_zcomplexabs.f:44