LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
test_zcomplexdiv.f
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  double precision x( n ), ainf, anan, b,
80  $ eps, bluemin, bluemax, ov, xj, stepx(n), limx(n)
81  double complex y, y2, r, cinf( ninf ), cnan( nnan )
82 *
83 * .. Intrinsic Functions ..
84  intrinsic dconjg, dble, radix, ceiling, tiny, digits,
85  $ maxexponent, minexponent, floor, huge, dcmplx,
86  $ epsilon
87 
88 *
89 * .. Initialize error counts ..
90  subnormaltreatedas0 = 0
91  caseafails = 0
92  casebfails = 0
93  casecfails = 0
94  casedfails = 0
95  caseefails = 0
96  caseffails = 0
97 *
98 * .. Initialize machine constants ..
99  min = minexponent(0.0d0)
100  max = maxexponent(0.0d0)
101  m = digits(0.0d0)
102  b = dble(radix(0.0d0))
103  eps = epsilon(0.0d0)
104  bluemin = b**ceiling( (min - 1) * 0.5d0 )
105  bluemax = b**floor( (max - m + 1) * 0.5d0 )
106  ov = huge(0.0d0)
107 *
108 * .. Vector X ..
109  x(1) = tiny(0.0d0) * b**( dble(1-m) )
110  x(2) = tiny(0.0d0)
111  x(3) = ov
112  x(4) = b**( dble(max-1) )
113 *
114 * .. Then modify X using the step ..
115  stepx(1) = 2.0
116  stepx(2) = 2.0
117  stepx(3) = 0.0
118  stepx(4) = 0.5
119 *
120 * .. Up to the value ..
121  limx(1) = x(2)
122  limx(2) = 1.0
123  limx(3) = 0.0
124  limx(4) = 2.0
125 *
126 * .. Inf entries ..
127  ainf = ov * 2
128  cinf(1) = dcmplx( ainf, 0.0d0 )
129  cinf(2) = dcmplx(-ainf, 0.0d0 )
130  cinf(3) = dcmplx( 0.0d0, ainf )
131  cinf(4) = dcmplx( 0.0d0,-ainf )
132  cinf(5) = dcmplx( ainf, ainf )
133 *
134 * .. NaN entries ..
135  anan = ainf / ainf
136  cnan(1) = dcmplx( anan, 0.0d0 )
137  cnan(2) = dcmplx( 0.0d0, anan )
138  cnan(3) = dcmplx( anan, anan )
139 
140 *
141 * .. Tests ..
142 *
143  if( debug ) then
144  print *, '# X :=', x
145  print *, '# Blue min constant :=', bluemin
146  print *, '# Blue max constant :=', bluemax
147  endif
148 *
149  xj = x(1)
150  if( xj .eq. 0.0d0 ) then
151  subnormaltreatedas0 = subnormaltreatedas0 + 1
152  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
153  print *, "!! fl( subnormal ) may be 0"
154  endif
155  else
156  do 100 i = 1, n
157  xj = x(i)
158  if( xj .eq. 0.0d0 ) then
159  subnormaltreatedas0 = subnormaltreatedas0 + 1
160  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
161  print *, "!! fl( subnormal ) may be 0"
162  endif
163  endif
164  100 continue
165  endif
166 *
167 * Test (a) y = x + 0 * I, y/y = 1
168  do 10 i = 1, n
169  xj = x(i)
170  if( xj .eq. 0.0d0 ) then
171  subnormaltreatedas0 = subnormaltreatedas0 + 1
172  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
173  print *, "!! [a] fl( subnormal ) may be 0"
174  endif
175  else
176  do while( xj .ne. limx(i) )
177  y = dcmplx( xj, 0.0d0 )
178  r = y / y
179  if( r .ne. 1.0d0 ) then
180  caseafails = caseafails + 1
181  if( caseafails .eq. 1 ) then
182  print *, "!! Some (x+0*I)/(x+0*I) differ from 1"
183  endif
184  WRITE( 0, fmt = 9999 ) 'a',i, xj,
185  $ '(x+0*I)/(x+0*I)', r, 1.0d0
186  endif
187  xj = xj * stepx(i)
188  end do
189  endif
190  10 continue
191 *
192 * Test (b) y = 0 + x * I, y/y = 1
193  do 20 i = 1, n
194  xj = x(i)
195  if( xj .eq. 0.0d0 ) then
196  subnormaltreatedas0 = subnormaltreatedas0 + 1
197  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
198  print *, "!! [b] fl( subnormal ) may be 0"
199  endif
200  else
201  do while( xj .ne. limx(i) )
202  y = dcmplx( 0.0d0, xj )
203  r = y / y
204  if( r .ne. 1.0d0 ) then
205  casebfails = casebfails + 1
206  if( casebfails .eq. 1 ) then
207  print *, "!! Some (0+x*I)/(0+x*I) differ from 1"
208  endif
209  WRITE( 0, fmt = 9999 ) 'b',i, xj,
210  $ '(0+x*I)/(0+x*I)', r, 1.0d0
211  endif
212  xj = xj * stepx(i)
213  end do
214  endif
215  20 continue
216 *
217 * Test (c) y = x + x * I, y/y = 1
218  do 30 i = 1, n
219  xj = x(i)
220  if( xj .eq. 0.0d0 ) then
221  subnormaltreatedas0 = subnormaltreatedas0 + 1
222  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
223  print *, "!! [c] fl( subnormal ) may be 0"
224  endif
225  else
226  do while( xj .ne. limx(i) )
227  y = dcmplx( xj, xj )
228  r = y / y
229  if( r .ne. 1.0d0 ) then
230  casecfails = casecfails + 1
231  if( casecfails .eq. 1 ) then
232  print *, "!! Some (x+x*I)/(x+x*I) differ from 1"
233  endif
234  WRITE( 0, fmt = 9999 ) 'c',i, xj,
235  $ '(x+x*I)/(x+x*I)', r, 1.0d0
236  endif
237  xj = xj * stepx(i)
238  end do
239  endif
240  30 continue
241 *
242 * Test (d) y1 = 0 + x * I, y2 = x + 0 * I, y1/y2 = I
243  do 40 i = 1, n
244  xj = x(i)
245  if( xj .eq. 0.0d0 ) then
246  subnormaltreatedas0 = subnormaltreatedas0 + 1
247  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
248  print *, "!! [d] fl( subnormal ) may be 0"
249  endif
250  else
251  do while( xj .ne. limx(i) )
252  y = dcmplx( 0.0d0, xj )
253  y2 = dcmplx( xj, 0.0d0 )
254  r = y / y2
255  if( r .ne. dcmplx(0.0d0,1.0d0) ) then
256  casedfails = casedfails + 1
257  if( casedfails .eq. 1 ) then
258  print *, "!! Some (0+x*I)/(x+0*I) differ from I"
259  endif
260  WRITE( 0, fmt = 9999 ) 'd',i, xj, '(0+x*I)/(x+0*I)',
261  $ r, dcmplx(0.0d0,1.0d0)
262  endif
263  xj = xj * stepx(i)
264  end do
265  endif
266  40 continue
267 *
268 * Test (e) y1 = 0 + x * I, y2 = x + 0 * I, y2/y1 = -I
269  do 50 i = 1, n
270  xj = x(i)
271  if( xj .eq. 0.0d0 ) then
272  subnormaltreatedas0 = subnormaltreatedas0 + 1
273  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
274  print *, "!! [e] fl( subnormal ) may be 0"
275  endif
276  else
277  do while( xj .ne. limx(i) )
278  y = dcmplx( 0.0d0, xj )
279  y2 = dcmplx( xj, 0.0d0 )
280  r = y2 / y
281  if( r .ne. dcmplx(0.0d0,-1.0d0) ) then
282  caseefails = caseefails + 1
283  if( caseefails .eq. 1 ) then
284  print *,"!! Some (x+0*I)/(0+x*I) differ from -I"
285  endif
286  WRITE( 0, fmt = 9999 ) 'e',i, xj, '(x+0*I)/(0+x*I)',
287  $ r, dcmplx(0.0d0,-1.0d0)
288  endif
289  xj = xj * stepx(i)
290  end do
291  endif
292  50 continue
293 *
294 * Test (f) y = x + x * I, y/conj(y) = I
295  do 60 i = 1, n
296  xj = x(i)
297  if( xj .eq. 0.0d0 ) then
298  subnormaltreatedas0 = subnormaltreatedas0 + 1
299  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
300  print *, "!! [f] fl( subnormal ) may be 0"
301  endif
302  else
303  do while( xj .ne. limx(i) )
304  y = dcmplx( xj, xj )
305  r = y / dconjg( y )
306  if( r .ne. dcmplx(0.0d0,1.0d0) ) then
307  caseffails = caseffails + 1
308  if( caseffails .eq. 1 ) then
309  print *, "!! Some (x+x*I)/(x-x*I) differ from I"
310  endif
311  WRITE( 0, fmt = 9999 ) 'f',i, xj, '(x+x*I)/(x-x*I)',
312  $ r, dcmplx(0.0d0,1.0d0)
313  endif
314  xj = xj * stepx(i)
315  end do
316  endif
317  60 continue
318 *
319 * Test (g) Infs
320  do 70 i = 1, ninf
321  y = cinf(i)
322  r = czero / y
323  if( (r .ne. czero) .and. (r .eq. r) ) then
324  WRITE( *, fmt = 9998 ) 'ia',i, czero, y, r, 'NaN and 0'
325  endif
326  r = cone / y
327  if( (r .ne. czero) .and. (r .eq. r) ) then
328  WRITE( *, fmt = 9998 ) 'ib',i, cone, y, r, 'NaN and 0'
329  endif
330  r = y / y
331  if( r .eq. r ) then
332  WRITE( *, fmt = 9998 ) 'ic',i, y, y, r, 'NaN'
333  endif
334  70 continue
335 *
336 * Test (h) NaNs
337  do 80 i = 1, nnan
338  y = cnan(i)
339  r = czero / y
340  if( r .eq. r ) then
341  WRITE( *, fmt = 9998 ) 'na',i, czero, y, r, 'NaN'
342  endif
343  r = cone / y
344  if( r .eq. r ) then
345  WRITE( *, fmt = 9998 ) 'nb',i, cone, y, r, 'NaN'
346  endif
347  r = y / y
348  if( r .eq. r ) then
349  WRITE( *, fmt = 9998 ) 'nc',i, y, y, r, 'NaN'
350  endif
351  80 continue
352 *
353 * If anything was written to stderr, print the message
354  if( (caseafails .gt. 0) .or. (casebfails .gt. 0) .or.
355  $ (casecfails .gt. 0) .or. (casedfails .gt. 0) .or.
356  $ (caseefails .gt. 0) .or. (caseffails .gt. 0) )
357  $ print *, "# Please check the failed divisions in [stderr]"
358 *
359 * .. Formats ..
360  9998 FORMAT( '[',a2,i1, '] ', (es24.16e3,sp,es24.16e3,"*I"), ' * ',
361  $ (es24.16e3,sp,es24.16e3,"*I"), ' = ',
362  $ (es24.16e3,sp,es24.16e3,"*I"), ' differs from ', a10 )
363 *
364  9999 FORMAT( '[',a2,i1, '] X = ', es24.16e3, ' : ', a15, ' = ',
365  $ (es24.16e3,sp,es24.16e3,"*I"), ' differs from ',
366  $ (es24.16e3,sp,es24.16e3,"*I") )
367 *
368 * End of zdiv
369 *
370  END