You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
401 lines
13 KiB
401 lines
13 KiB
*> \brief zdiv tests the robustness and precision of the double complex division
|
|
*
|
|
* =========== DOCUMENTATION ===========
|
|
*
|
|
* Online html documentation available at
|
|
* http://www.netlib.org/lapack/explore-html/
|
|
*
|
|
* Authors:
|
|
* ========
|
|
*
|
|
*> \author Weslley S. Pereira, University of Colorado Denver, U.S.
|
|
*
|
|
*> \verbatim
|
|
*>
|
|
*> Real values for test:
|
|
*> (1) x = 2**m, where m = MINEXPONENT-DIGITS, ..., MINEXPONENT-1.
|
|
*> Mind that not all platforms might implement subnormal numbers.
|
|
*> (2) x = 2**m, where m = MINEXPONENT, ..., 0.
|
|
*> (3) x = OV, where OV is the overflow threshold. OV^2 overflows but the norm is OV.
|
|
*> (4) x = 2**m, where m = MAXEXPONENT-1, ..., 1.
|
|
*>
|
|
*> Tests:
|
|
*> (a) y = x + 0 * I, y/y = 1
|
|
*> (b) y = 0 + x * I, y/y = 1
|
|
*> (c) y = x + x * I, y/y = 1
|
|
*> (d) y1 = 0 + x * I, y2 = x + 0 * I, y1/y2 = I
|
|
*> (e) y1 = 0 + x * I, y2 = x + 0 * I, y2/y1 = -I
|
|
*> (f) y = x + x * I, y/conj(y) = I
|
|
*>
|
|
*> Special cases:
|
|
*>
|
|
*> (i) Inf inputs:
|
|
*> (1) y = ( Inf + 0 * I)
|
|
*> (2) y = ( 0 + Inf * I)
|
|
*> (3) y = (-Inf + 0 * I)
|
|
*> (4) y = ( 0 - Inf * I)
|
|
*> (5) y = ( Inf + Inf * I)
|
|
*> Tests:
|
|
*> (a) 0 / y is either 0 or NaN.
|
|
*> (b) 1 / y is either 0 or NaN.
|
|
*> (c) y / y is NaN.
|
|
*>
|
|
*> (n) NaN inputs:
|
|
*> (1) y = (NaN + 0 * I)
|
|
*> (2) y = (0 + NaN * I)
|
|
*> (3) y = (NaN + NaN * I)
|
|
*> Tests:
|
|
*> (a) 0 / y is NaN.
|
|
*> (b) 1 / y is NaN.
|
|
*> (c) y / y is NaN.
|
|
*>
|
|
*> \endverbatim
|
|
*
|
|
*> \ingroup auxOTHERauxiliary
|
|
*
|
|
* =====================================================================
|
|
program zdiv
|
|
*
|
|
* -- LAPACK test routine --
|
|
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|
|
|
* ..
|
|
* .. Local parameters ..
|
|
logical debug
|
|
parameter ( debug = .false. )
|
|
integer N, nNaN, nInf
|
|
parameter ( N = 4, nNaN = 3, nInf = 5 )
|
|
double precision threeFourth, fiveFourth
|
|
parameter ( threeFourth = 3.0d0 / 4,
|
|
$ fiveFourth = 5.0d0 / 4 )
|
|
double complex czero, cone
|
|
parameter ( czero = DCMPLX( 0.0d0, 0.0d0 ),
|
|
$ cone = DCMPLX( 1.0d0, 0.0d0 ) )
|
|
* ..
|
|
* .. Local Variables ..
|
|
integer i, min, Max, m,
|
|
$ subnormalTreatedAs0, caseAFails, caseBFails,
|
|
$ caseCFails, caseDFails, caseEFails, caseFFails,
|
|
$ caseInfFails, caseNaNFails, nFailingTests,
|
|
$ nTests
|
|
double precision X( N ), aInf, aNaN, b,
|
|
$ eps, blueMin, blueMax, OV, Xj, stepX(N), limX(N)
|
|
double complex Y, Y2, R, cInf( nInf ), cNaN( nNaN )
|
|
*
|
|
* .. Intrinsic Functions ..
|
|
intrinsic DCONJG, DBLE, RADIX, CEILING, TINY, DIGITS,
|
|
$ MAXEXPONENT, MINEXPONENT, FLOOR, HUGE, DCMPLX,
|
|
$ EPSILON
|
|
|
|
*
|
|
* .. Initialize error counts ..
|
|
subnormalTreatedAs0 = 0
|
|
caseAFails = 0
|
|
caseBFails = 0
|
|
caseCFails = 0
|
|
caseDFails = 0
|
|
caseEFails = 0
|
|
caseFFails = 0
|
|
caseInfFails = 0
|
|
caseNaNFails = 0
|
|
nFailingTests = 0
|
|
nTests = 0
|
|
*
|
|
* .. Initialize machine constants ..
|
|
min = MINEXPONENT(0.0d0)
|
|
Max = MAXEXPONENT(0.0d0)
|
|
m = DIGITS(0.0d0)
|
|
b = DBLE(RADIX(0.0d0))
|
|
eps = EPSILON(0.0d0)
|
|
blueMin = b**CEILING( (min - 1) * 0.5d0 )
|
|
blueMax = b**FLOOR( (Max - m + 1) * 0.5d0 )
|
|
OV = HUGE(0.0d0)
|
|
*
|
|
* .. Vector X ..
|
|
X(1) = TINY(0.0d0) * b**( DBLE(1-m) )
|
|
X(2) = TINY(0.0d0)
|
|
X(3) = OV
|
|
X(4) = b**( DBLE(Max-1) )
|
|
*
|
|
* .. Then modify X using the step ..
|
|
stepX(1) = 2.0
|
|
stepX(2) = 2.0
|
|
stepX(3) = 0.0
|
|
stepX(4) = 0.5
|
|
*
|
|
* .. Up to the value ..
|
|
limX(1) = X(2)
|
|
limX(2) = 1.0
|
|
limX(3) = 0.0
|
|
limX(4) = 2.0
|
|
*
|
|
* .. Inf entries ..
|
|
aInf = OV * 2
|
|
cInf(1) = DCMPLX( aInf, 0.0d0 )
|
|
cInf(2) = DCMPLX(-aInf, 0.0d0 )
|
|
cInf(3) = DCMPLX( 0.0d0, aInf )
|
|
cInf(4) = DCMPLX( 0.0d0,-aInf )
|
|
cInf(5) = DCMPLX( aInf, aInf )
|
|
*
|
|
* .. NaN entries ..
|
|
aNaN = aInf / aInf
|
|
cNaN(1) = DCMPLX( aNaN, 0.0d0 )
|
|
cNaN(2) = DCMPLX( 0.0d0, aNaN )
|
|
cNaN(3) = DCMPLX( aNaN, aNaN )
|
|
|
|
*
|
|
* .. Tests ..
|
|
*
|
|
if( debug ) then
|
|
print *, '# X :=', X
|
|
print *, '# Blue min constant :=', blueMin
|
|
print *, '# Blue max constant :=', blueMax
|
|
endif
|
|
*
|
|
Xj = X(1)
|
|
if( Xj .eq. 0.0d0 ) then
|
|
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
|
|
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
|
|
print *, "!! fl( subnormal ) may be 0"
|
|
endif
|
|
else
|
|
do 100 i = 1, N
|
|
Xj = X(i)
|
|
if( Xj .eq. 0.0d0 ) then
|
|
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
|
|
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
|
|
print *, "!! fl( subnormal ) may be 0"
|
|
endif
|
|
endif
|
|
100 continue
|
|
endif
|
|
*
|
|
* Test (a) y = x + 0 * I, y/y = 1
|
|
do 10 i = 1, N
|
|
Xj = X(i)
|
|
if( Xj .eq. 0.0d0 ) then
|
|
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
|
|
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
|
|
print *, "!! [a] fl( subnormal ) may be 0"
|
|
endif
|
|
else
|
|
do while( Xj .ne. limX(i) )
|
|
nTests = nTests + 1
|
|
Y = DCMPLX( Xj, 0.0d0 )
|
|
R = Y / Y
|
|
if( R .ne. 1.0D0 ) then
|
|
caseAFails = caseAFails + 1
|
|
if( caseAFails .eq. 1 ) then
|
|
print *, "!! Some (x+0*I)/(x+0*I) differ from 1"
|
|
endif
|
|
WRITE( 0, FMT = 9999 ) 'a',i, Xj,
|
|
$ '(x+0*I)/(x+0*I)', R, 1.0D0
|
|
endif
|
|
Xj = Xj * stepX(i)
|
|
end do
|
|
endif
|
|
10 continue
|
|
*
|
|
* Test (b) y = 0 + x * I, y/y = 1
|
|
do 20 i = 1, N
|
|
Xj = X(i)
|
|
if( Xj .eq. 0.0d0 ) then
|
|
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
|
|
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
|
|
print *, "!! [b] fl( subnormal ) may be 0"
|
|
endif
|
|
else
|
|
do while( Xj .ne. limX(i) )
|
|
nTests = nTests + 1
|
|
Y = DCMPLX( 0.0d0, Xj )
|
|
R = Y / Y
|
|
if( R .ne. 1.0D0 ) then
|
|
caseBFails = caseBFails + 1
|
|
if( caseBFails .eq. 1 ) then
|
|
print *, "!! Some (0+x*I)/(0+x*I) differ from 1"
|
|
endif
|
|
WRITE( 0, FMT = 9999 ) 'b',i, Xj,
|
|
$ '(0+x*I)/(0+x*I)', R, 1.0D0
|
|
endif
|
|
Xj = Xj * stepX(i)
|
|
end do
|
|
endif
|
|
20 continue
|
|
*
|
|
* Test (c) y = x + x * I, y/y = 1
|
|
do 30 i = 1, N
|
|
Xj = X(i)
|
|
if( Xj .eq. 0.0d0 ) then
|
|
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
|
|
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
|
|
print *, "!! [c] fl( subnormal ) may be 0"
|
|
endif
|
|
else
|
|
do while( Xj .ne. limX(i) )
|
|
nTests = nTests + 1
|
|
Y = DCMPLX( Xj, Xj )
|
|
R = Y / Y
|
|
if( R .ne. 1.0D0 ) then
|
|
caseCFails = caseCFails + 1
|
|
if( caseCFails .eq. 1 ) then
|
|
print *, "!! Some (x+x*I)/(x+x*I) differ from 1"
|
|
endif
|
|
WRITE( 0, FMT = 9999 ) 'c',i, Xj,
|
|
$ '(x+x*I)/(x+x*I)', R, 1.0D0
|
|
endif
|
|
Xj = Xj * stepX(i)
|
|
end do
|
|
endif
|
|
30 continue
|
|
*
|
|
* Test (d) y1 = 0 + x * I, y2 = x + 0 * I, y1/y2 = I
|
|
do 40 i = 1, N
|
|
Xj = X(i)
|
|
if( Xj .eq. 0.0d0 ) then
|
|
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
|
|
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
|
|
print *, "!! [d] fl( subnormal ) may be 0"
|
|
endif
|
|
else
|
|
do while( Xj .ne. limX(i) )
|
|
nTests = nTests + 1
|
|
Y = DCMPLX( 0.0d0, Xj )
|
|
Y2 = DCMPLX( Xj, 0.0d0 )
|
|
R = Y / Y2
|
|
if( R .ne. DCMPLX(0.0D0,1.0D0) ) then
|
|
caseDFails = caseDFails + 1
|
|
if( caseDFails .eq. 1 ) then
|
|
print *, "!! Some (0+x*I)/(x+0*I) differ from I"
|
|
endif
|
|
WRITE( 0, FMT = 9999 ) 'd',i, Xj, '(0+x*I)/(x+0*I)',
|
|
$ R, DCMPLX(0.0D0,1.0D0)
|
|
endif
|
|
Xj = Xj * stepX(i)
|
|
end do
|
|
endif
|
|
40 continue
|
|
*
|
|
* Test (e) y1 = 0 + x * I, y2 = x + 0 * I, y2/y1 = -I
|
|
do 50 i = 1, N
|
|
Xj = X(i)
|
|
if( Xj .eq. 0.0d0 ) then
|
|
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
|
|
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
|
|
print *, "!! [e] fl( subnormal ) may be 0"
|
|
endif
|
|
else
|
|
do while( Xj .ne. limX(i) )
|
|
nTests = nTests + 1
|
|
Y = DCMPLX( 0.0d0, Xj )
|
|
Y2 = DCMPLX( Xj, 0.0d0 )
|
|
R = Y2 / Y
|
|
if( R .ne. DCMPLX(0.0D0,-1.0D0) ) then
|
|
caseEFails = caseEFails + 1
|
|
if( caseEFails .eq. 1 ) then
|
|
print *,"!! Some (x+0*I)/(0+x*I) differ from -I"
|
|
endif
|
|
WRITE( 0, FMT = 9999 ) 'e',i, Xj, '(x+0*I)/(0+x*I)',
|
|
$ R, DCMPLX(0.0D0,-1.0D0)
|
|
endif
|
|
Xj = Xj * stepX(i)
|
|
end do
|
|
endif
|
|
50 continue
|
|
*
|
|
* Test (f) y = x + x * I, y/conj(y) = I
|
|
do 60 i = 1, N
|
|
Xj = X(i)
|
|
if( Xj .eq. 0.0d0 ) then
|
|
subnormalTreatedAs0 = subnormalTreatedAs0 + 1
|
|
if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
|
|
print *, "!! [f] fl( subnormal ) may be 0"
|
|
endif
|
|
else
|
|
do while( Xj .ne. limX(i) )
|
|
nTests = nTests + 1
|
|
Y = DCMPLX( Xj, Xj )
|
|
R = Y / DCONJG( Y )
|
|
if( R .ne. DCMPLX(0.0D0,1.0D0) ) then
|
|
caseFFails = caseFFails + 1
|
|
if( caseFFails .eq. 1 ) then
|
|
print *, "!! Some (x+x*I)/(x-x*I) differ from I"
|
|
endif
|
|
WRITE( 0, FMT = 9999 ) 'f',i, Xj, '(x+x*I)/(x-x*I)',
|
|
$ R, DCMPLX(0.0D0,1.0D0)
|
|
endif
|
|
Xj = Xj * stepX(i)
|
|
end do
|
|
endif
|
|
60 continue
|
|
*
|
|
* Test (g) Infs
|
|
do 70 i = 1, nInf
|
|
nTests = nTests + 3
|
|
Y = cInf(i)
|
|
R = czero / Y
|
|
if( (R .ne. czero) .and. (R .eq. R) ) then
|
|
caseInfFails = caseInfFails + 1
|
|
WRITE( *, FMT = 9998 ) 'ia',i, czero, Y, R, 'NaN and 0'
|
|
endif
|
|
R = cone / Y
|
|
if( (R .ne. czero) .and. (R .eq. R) ) then
|
|
caseInfFails = caseInfFails + 1
|
|
WRITE( *, FMT = 9998 ) 'ib',i, cone, Y, R, 'NaN and 0'
|
|
endif
|
|
R = Y / Y
|
|
if( R .eq. R ) then
|
|
caseInfFails = caseInfFails + 1
|
|
WRITE( *, FMT = 9998 ) 'ic',i, Y, Y, R, 'NaN'
|
|
endif
|
|
70 continue
|
|
*
|
|
* Test (h) NaNs
|
|
do 80 i = 1, nNaN
|
|
nTests = nTests + 3
|
|
Y = cNaN(i)
|
|
R = czero / Y
|
|
if( R .eq. R ) then
|
|
caseNaNFails = caseNaNFails + 1
|
|
WRITE( *, FMT = 9998 ) 'na',i, czero, Y, R, 'NaN'
|
|
endif
|
|
R = cone / Y
|
|
if( R .eq. R ) then
|
|
caseNaNFails = caseNaNFails + 1
|
|
WRITE( *, FMT = 9998 ) 'nb',i, cone, Y, R, 'NaN'
|
|
endif
|
|
R = Y / Y
|
|
if( R .eq. R ) then
|
|
caseNaNFails = caseNaNFails + 1
|
|
WRITE( *, FMT = 9998 ) 'nc',i, Y, Y, R, 'NaN'
|
|
endif
|
|
80 continue
|
|
*
|
|
* If any test fails, displays a message
|
|
nFailingTests = caseAFails + caseBFails + caseCFails + caseDFails
|
|
$ + caseEFails + caseFFails + caseInfFails
|
|
$ + caseNaNFails
|
|
if( nFailingTests .gt. 0 ) then
|
|
print *, "# ", nTests-nFailingTests, " tests out of ", nTests,
|
|
$ " pass for complex division,", nFailingTests," fail."
|
|
else
|
|
print *, "# All tests pass for complex division."
|
|
endif
|
|
*
|
|
* If anything was written to stderr, print the message
|
|
if( (caseAFails .gt. 0) .or. (caseBFails .gt. 0) .or.
|
|
$ (caseCFails .gt. 0) .or. (caseDFails .gt. 0) .or.
|
|
$ (caseEFails .gt. 0) .or. (caseFFails .gt. 0) )
|
|
$ print *, "# Please check the failed divisions in [stderr]"
|
|
*
|
|
* .. Formats ..
|
|
9998 FORMAT( '[',A2,I1, '] ', (ES24.16E3,SP,ES24.16E3,"*I"), ' * ',
|
|
$ (ES24.16E3,SP,ES24.16E3,"*I"), ' = ',
|
|
$ (ES24.16E3,SP,ES24.16E3,"*I"), ' differs from ', A10 )
|
|
*
|
|
9999 FORMAT( '[',A2,I1, '] X = ', ES24.16E3, ' : ', A15, ' = ',
|
|
$ (ES24.16E3,SP,ES24.16E3,"*I"), ' differs from ',
|
|
$ (ES24.16E3,SP,ES24.16E3,"*I") )
|
|
*
|
|
* End of zdiv
|
|
*
|
|
END |