11* > \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+ *
211* > \author Weslley S. Pereira, University of Colorado Denver, U.S.
312*
413* > \verbatim
4251* >
4352* > \endverbatim
4453*
54+ * > \ingroup auxOTHERauxiliary
55+ *
56+ * =====================================================================
4557 program zdiv
58+ *
59+ * -- LAPACK test routine --
60+ * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
4661
62+ * ..
63+ * .. Local parameters ..
4764 logical debug
4865 parameter ( debug = .false. )
49-
50- integer N, i, nNaN, nInf, min, Max, m
66+ integer N, nNaN, nInf
5167 parameter ( N = 4 , nNaN = 3 , nInf = 5 )
52-
53- double precision X( N ), threeFourth, fiveFourth, aInf, aNaN, b,
54- $ eps, blueMin, blueMax, OV, Xj, stepX(N), limX(N)
68+ double precision threeFourth, fiveFourth
5569 parameter ( threeFourth = 3.0d0 / 4 ,
5670 $ fiveFourth = 5.0d0 / 4 )
57-
58- double complex Y, Y2, R, cInf( nInf ), cNaN( nNaN ), czero,
59- $ cone
71+ double complex czero, cone
6072 parameter ( czero = DCMPLX( 0.0d0 , 0.0d0 ),
6173 $ 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 )
6282*
83+ * .. Intrinsic Functions ..
6384 intrinsic DCONJG, DBLE, RADIX, CEILING, TINY, DIGITS,
6485 $ MAXEXPONENT, MINEXPONENT, FLOOR, HUGE, DCMPLX,
6586 $ EPSILON
6687
67- integer subnormalTreatedAs0, caseAFails, caseBFails,
68- $ caseCFails, caseDFails
6988*
89+ * .. Initialize error counts ..
7090 subnormalTreatedAs0 = 0
7191 caseAFails = 0
7292 caseBFails = 0
@@ -75,6 +95,7 @@ program zdiv
7595 caseEFails = 0
7696 caseFFails = 0
7797*
98+ * .. Initialize machine constants ..
7899 min = MINEXPONENT (0.0d0 )
79100 Max = MAXEXPONENT (0.0d0 )
80101 m = DIGITS (0.0d0 )
@@ -84,20 +105,40 @@ program zdiv
84105 blueMax = b** FLOOR( (Max - m + 1 ) * 0.5d0 )
85106 OV = HUGE (0.0d0 )
86107*
108+ * .. Vector X ..
87109 X(1 ) = TINY (0.0d0 ) * b** ( DBLE (1 - m) )
88110 X(2 ) = TINY (0.0d0 )
89111 X(3 ) = OV
90112 X(4 ) = b** ( DBLE (Max-1 ) )
91113*
114+ * .. Then modify X using the step ..
92115 stepX(1 ) = 2.0
93116 stepX(2 ) = 2.0
94117 stepX(3 ) = 0.0
95118 stepX(4 ) = 0.5
96119*
120+ * .. Up to the value ..
97121 limX(1 ) = X(2 )
98122 limX(2 ) = 1.0
99123 limX(3 ) = 0.0
100124 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 ..
101142*
102143 if ( debug ) then
103144 print * , ' # X :=' , X
@@ -107,26 +148,21 @@ program zdiv
107148*
108149 Xj = X(1 )
109150 if ( Xj .eq. 0.0d0 ) then
110- print * , " # Subnormal numbers treated as 0"
151+ subnormalTreatedAs0 = subnormalTreatedAs0 + 1
152+ if ( debug .or. subnormalTreatedAs0 .eq. 1 ) then
153+ print * , " !! fl( subnormal ) may be 0"
154+ endif
111155 else
112156 do 100 i = 1 , N
113157 Xj = X(i)
114- if ( Xj .eq. 0.0d0 ) print * ,
115- $ " # Subnormal numbers may be treated as 0"
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
116164 100 continue
117165 endif
118- *
119- aInf = OV * 2
120- cInf(1 ) = DCMPLX( aInf, 0.0d0 )
121- cInf(2 ) = DCMPLX(- aInf, 0.0d0 )
122- cInf(3 ) = DCMPLX( 0.0d0 , aInf )
123- cInf(4 ) = DCMPLX( 0.0d0 ,- aInf )
124- cInf(5 ) = DCMPLX( aInf, aInf )
125- *
126- aNaN = aInf / aInf
127- cNaN(1 ) = DCMPLX( aNaN, 0.0d0 )
128- cNaN(2 ) = DCMPLX( 0.0d0 , aNaN )
129- cNaN(3 ) = DCMPLX( aNaN, aNaN )
130166*
131167* Test (a) y = x + 0 * I, y/y = 1
132168 do 10 i = 1 , N
@@ -296,11 +332,6 @@ program zdiv
296332 WRITE ( * , FMT = 9998 ) ' ic' ,i, Y, Y, R, ' NaN'
297333 endif
298334 70 continue
299- *
300- if ( (caseAFails .gt. 0 ) .or. (caseBFails .gt. 0 ) .or.
301- $ (caseCFails .gt. 0 ) .or. (caseDFails .gt. 0 ) .or.
302- $ (caseEFails .gt. 0 ) .or. (caseFFails .gt. 0 ) )
303- $ print * , " # Please check the failed divisions in [stderr]"
304335*
305336* Test (h) NaNs
306337 do 80 i = 1 , nNaN
@@ -319,6 +350,13 @@ program zdiv
319350 endif
320351 80 continue
321352*
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 ..
322360 9998 FORMAT ( ' [' ,A2,I1, ' ] X = ' , ES24.16E3 , ' : ' , A15, ' = ' ,
323361 $ (ES24.16E3 ,SP,ES24.16E3 ," *I" ), ' differs from ' , A10 )
324362*
0 commit comments