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