LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
test_zcomplexabs.f
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  double precision x( n ), r, answerc,
64  $ answerd, ainf, anan, reldiff, b,
65  $ eps, bluemin, bluemax, xj, stepx(n), limx(n)
66  double complex y, cinf( ninf ), cnan( nnan )
67 *
68 * .. Intrinsic Functions ..
69  intrinsic abs, dble, radix, ceiling, tiny, digits, sqrt,
70  $ maxexponent, minexponent, floor, huge, dcmplx,
71  $ epsilon
72 
73 *
74 * .. Initialize error counts ..
75  subnormaltreatedas0 = 0
76  caseafails = 0
77  casebfails = 0
78  casecfails = 0
79  casedfails = 0
80 *
81 * .. Initialize machine constants ..
82  min = minexponent(0.0d0)
83  max = maxexponent(0.0d0)
84  m = digits(0.0d0)
85  b = dble(radix(0.0d0))
86  eps = epsilon(0.0d0)
87  bluemin = b**ceiling( (min - 1) * 0.5d0 )
88  bluemax = b**floor( (max - m + 1) * 0.5d0 )
89 *
90 * .. Vector X ..
91  x(1) = tiny(0.0d0) * b**( dble(1-m) )
92  x(2) = tiny(0.0d0)
93  x(3) = huge(0.0d0)
94  x(4) = b**( dble(max-1) )
95 *
96 * .. Then modify X using the step ..
97  stepx(1) = 2.0
98  stepx(2) = 2.0
99  stepx(3) = 0.0
100  stepx(4) = 0.5
101 *
102 * .. Up to the value ..
103  limx(1) = x(2)
104  limx(2) = 1.0
105  limx(3) = 0.0
106  limx(4) = 2.0
107 *
108 * .. Inf entries ..
109  ainf = x(3) * 2
110  cinf(1) = dcmplx( ainf, 0.0d0 )
111  cinf(2) = dcmplx(-ainf, 0.0d0 )
112  cinf(3) = dcmplx( 0.0d0, ainf )
113  cinf(4) = dcmplx( 0.0d0,-ainf )
114  cinf(5) = dcmplx( ainf, ainf )
115 *
116 * .. NaN entries ..
117  anan = ainf / ainf
118  cnan(1) = dcmplx( anan, 0.0d0 )
119  cnan(2) = dcmplx( 0.0d0, anan )
120  cnan(3) = dcmplx( anan, anan )
121 
122 *
123 * .. Tests ..
124 *
125  if( debug ) then
126  print *, '# X :=', x
127  print *, '# Blue min constant :=', bluemin
128  print *, '# Blue max constant :=', bluemax
129  endif
130 *
131  xj = x(1)
132  if( xj .eq. 0.0d0 ) then
133  subnormaltreatedas0 = subnormaltreatedas0 + 1
134  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
135  print *, "!! fl( subnormal ) may be 0"
136  endif
137  else
138  do 100 i = 1, n
139  xj = x(i)
140  if( xj .eq. 0.0d0 ) then
141  subnormaltreatedas0 = subnormaltreatedas0 + 1
142  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
143  print *, "!! fl( subnormal ) may be 0"
144  endif
145  endif
146  100 continue
147  endif
148 *
149 * Test (a) y = x + 0 * I, |y| = x
150  do 10 i = 1, n
151  xj = x(i)
152  if( xj .eq. 0.0d0 ) then
153  subnormaltreatedas0 = subnormaltreatedas0 + 1
154  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
155  print *, "!! [a] fl( subnormal ) may be 0"
156  endif
157  else
158  do while( xj .ne. limx(i) )
159  y = dcmplx( xj, 0.0d0 )
160  r = abs( y )
161  if( r .ne. xj ) then
162  caseafails = caseafails + 1
163  if( caseafails .eq. 1 ) then
164  print *, "!! Some ABS(x+0*I) differ from ABS(x)"
165  endif
166  WRITE( 0, fmt = 9999 ) 'a',i, xj, '(1+0*I)', r, xj
167  endif
168  xj = xj * stepx(i)
169  end do
170  endif
171  10 continue
172 *
173 * Test (b) y = 0 + x * I, |y| = x
174  do 20 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 *, "!! [b] fl( subnormal ) may be 0"
180  endif
181  else
182  do while( xj .ne. limx(i) )
183  y = dcmplx( 0.0d0, xj )
184  r = abs( y )
185  if( r .ne. xj ) then
186  casebfails = casebfails + 1
187  if( casebfails .eq. 1 ) then
188  print *, "!! Some ABS(0+x*I) differ from ABS(x)"
189  endif
190  WRITE( 0, fmt = 9999 ) 'b',i, xj, '(0+1*I)', r, xj
191  endif
192  xj = xj * stepx(i)
193  end do
194  endif
195  20 continue
196 *
197 * Test (c) y = (3/4)*x + x * I, |y| = (5/4)*x
198  do 30 i = 1, n
199  if( i .eq. 3 ) go to 30
200  if( i .eq. 1 ) then
201  xj = 4*x(i)
202  else
203  xj = x(i)
204  endif
205  if( xj .eq. 0.0d0 ) then
206  subnormaltreatedas0 = subnormaltreatedas0 + 1
207  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
208  print *, "!! [c] fl( subnormal ) may be 0"
209  endif
210  else
211  do while( xj .ne. limx(i) )
212  answerc = fivefourth * xj
213  y = dcmplx( threefourth * xj, xj )
214  r = abs( y )
215  if( r .ne. answerc ) then
216  casecfails = casecfails + 1
217  if( casecfails .eq. 1 ) then
218  print *,
219  $ "!! Some ABS(x*(3/4+I)) differ from (5/4)*ABS(x)"
220  endif
221  WRITE( 0, fmt = 9999 ) 'c',i, xj, '(3/4+I)', r,
222  $ answerc
223  endif
224  xj = xj * stepx(i)
225  end do
226  endif
227  30 continue
228 *
229 * Test (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2)
230  do 40 i = 1, n
231  if( i .eq. 1 ) then
232  xj = 2*x(i)
233  else
234  xj = x(i)
235  endif
236  if( xj .eq. 0.0d0 ) then
237  subnormaltreatedas0 = subnormaltreatedas0 + 1
238  if( debug .or. subnormaltreatedas0 .eq. 1 ) then
239  print *, "!! [d] fl( subnormal ) may be 0"
240  endif
241  else
242  do while( xj .ne. limx(i) )
243  answerd = (onehalf * xj) * sqrt(2.0d0)
244  if( answerd .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  y = dcmplx( onehalf * xj, onehalf * xj )
251  r = abs( y )
252  reldiff = abs(r-answerd)/answerd
253  if( reldiff .ge. (0.5*eps) ) then
254  casedfails = casedfails + 1
255  if( casedfails .eq. 1 ) then
256  print *,
257  $ "!! Some ABS(x*(1+I)) differ from sqrt(2)*ABS(x)"
258  endif
259  WRITE( 0, fmt = 9999 ) 'd',i, (onehalf*xj),
260  $ '(1+1*I)', r, answerd
261  endif
262  endif
263  xj = xj * stepx(i)
264  end do
265  endif
266  40 continue
267 *
268 * Test (e) Infs
269  do 50 i = 1, ninf
270  y = cinf(i)
271  r = abs( y )
272  if( .not.(r .gt. huge(0.0d0)) ) then
273  WRITE( *, fmt = 9997 ) 'i',i, y, r
274  endif
275  50 continue
276 *
277 * Test (f) NaNs
278  do 60 i = 1, nnan
279  y = cnan(i)
280  r = abs( y )
281  if( r .eq. r ) then
282  WRITE( *, fmt = 9998 ) 'n',i, y, r
283  endif
284  60 continue
285 *
286 * If anything was written to stderr, print the message
287  if( (caseafails .gt. 0) .or. (casebfails .gt. 0) .or.
288  $ (casecfails .gt. 0) .or. (casedfails .gt. 0) )
289  $ print *, "# Please check the failed ABS(a+b*I) in [stderr]"
290 *
291 * .. Formats ..
292  9997 FORMAT( '[',a1,i1, '] ABS(', (es8.1,sp,es8.1,"*I"), ' ) = ',
293  $ es8.1, ' differs from Inf' )
294 *
295  9998 FORMAT( '[',a1,i1, '] ABS(', (es8.1,sp,es8.1,"*I"), ' ) = ',
296  $ es8.1, ' differs from NaN' )
297 *
298  9999 FORMAT( '[',a1,i1, '] ABS(', es24.16e3, ' * ', a7, ' ) = ',
299  $ es24.16e3, ' differs from ', es24.16e3 )
300 *
301 * End of zabs
302 *
303  END