test_zcomplexabs.f

Go to the documentation of this file.

*> \brief zabs tests the robustness and precision of the intrinsic ABS for double complex
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \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. Stop on the first success.
*>     Mind that not all platforms might implement subnormal numbers.
*> (2) x = 2**m, where m = MINEXPONENT, ..., 0. Stop on the first success.
*> (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. Stop on the first success.
*>
*> Tests:
*> (a) y = x + 0 * I, |y| = x
*> (b) y = 0 + x * I, |y| = x
*> (c) y = (3/4)*x + x * I, |y| = (5/4)*x whenever (3/4)*x and (5/4)*x can be exactly stored
*> (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2) whenever (1/2)*x can be exactly stored
*>
*> Special cases:
*>
*> (i) Inf propagation
*>    (1) y = Inf + 0 * I, |y| is Inf.
*>    (2) y =-Inf + 0 * I, |y| is Inf.
*>    (3) y = 0 + Inf * I, |y| is Inf.
*>    (4) y = 0 - Inf * I, |y| is Inf.
*>    (5) y = Inf + Inf * I, |y| is Inf.
*>
*> (n) NaN propagation
*>    (1) y = NaN + 0 * I, |y| is NaN.
*>    (2) y = 0 + NaN * I, |y| is NaN.
*>    (3) y = NaN + NaN * I, |y| is NaN.
*>
*> \endverbatim
*
*> \ingroup auxOTHERauxiliary
*
*  =====================================================================
      program zabs
*
*  -- 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, onehalf
      parameter( threefourth = 3.0d0 / 4,
     $                  fivefourth = 5.0d0 / 4,
     $                  onehalf = 1.0d0 / 2 )
*     ..
*     .. Local Variables ..
      integer           i, min, max, m, subnormaltreatedas0,
     $                  caseafails, casebfails, casecfails, casedfails
      double precision  x( n ), r, answerc,
     $                  answerd, ainf, anan, reldiff, b,
     $                  eps, bluemin, bluemax, xj, stepx(n), limx(n)
      double complex    y, cinf( ninf ), cnan( nnan )
*
*     .. Intrinsic Functions ..
      intrinsic         abs, dble, radix, ceiling, tiny, digits, sqrt,
     $                  maxexponent, minexponent, floor, huge, dcmplx,
     $                  epsilon
 
*
*     .. Initialize error counts ..
      subnormaltreatedas0 = 0
      caseafails = 0
      casebfails = 0
      casecfails = 0
      casedfails = 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 )
*
*     .. Vector X ..
      x(1) = tiny(0.0d0) * b**( dble(1-m) )
      x(2) = tiny(0.0d0)
      x(3) = huge(0.0d0)
      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 = x(3) * 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| = x
      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) )
                y = dcmplx( xj, 0.0d0 )
                r = abs( y )
                if( r .ne. xj ) then
                    caseafails = caseafails + 1
                    if( caseafails .eq. 1 ) then
                        print *, "!! Some ABS(x+0*I) differ from ABS(x)"
                    endif
                    WRITE( 0, fmt = 9999 ) 'a',i, xj, '(1+0*I)', r, xj
                endif
                xj = xj * stepx(i)
            end do
        endif
  10  continue
*
*     Test (b) y = 0 + x * I, |y| = x
      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) )
                y = dcmplx( 0.0d0, xj )
                r = abs( y )
                if( r .ne. xj ) then
                    casebfails = casebfails + 1
                    if( casebfails .eq. 1 ) then
                        print *, "!! Some ABS(0+x*I) differ from ABS(x)"
                    endif
                    WRITE( 0, fmt = 9999 ) 'b',i, xj, '(0+1*I)', r, xj
                endif
                xj = xj * stepx(i)
            end do
        endif
  20  continue
*
*     Test (c) y = (3/4)*x + x * I, |y| = (5/4)*x
      do 30 i = 1, n
        if( i .eq. 3 ) go to 30
        if( i .eq. 1 ) then
            xj = 4*x(i)
        else
            xj = x(i)
        endif
        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) )
                answerc = fivefourth * xj
                y = dcmplx( threefourth * xj, xj )
                r = abs( y )
                if( r .ne. answerc ) then
                    casecfails = casecfails + 1
                    if( casecfails .eq. 1 ) then
                        print *, 
     $              "!! Some ABS(x*(3/4+I)) differ from (5/4)*ABS(x)"
                    endif
                    WRITE( 0, fmt = 9999 ) 'c',i, xj, '(3/4+I)', r,
     $                                      answerc
                endif
                xj = xj * stepx(i)
            end do
        endif
  30  continue
*
*     Test (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2)
      do 40 i = 1, n
        if( i .eq. 1 ) then
            xj = 2*x(i)
        else
            xj = x(i)
        endif
        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) )
                answerd = (onehalf * xj) * sqrt(2.0d0)
                if( answerd .eq. 0.0d0 ) then
                    subnormaltreatedas0 = subnormaltreatedas0 + 1
                    if( debug .or. subnormaltreatedas0 .eq. 1 ) then
                        print *, "!! [d] fl( subnormal ) may be 0"
                    endif
                else
                    y = dcmplx( onehalf * xj, onehalf * xj )
                    r = abs( y )
                    reldiff = abs(r-answerd)/answerd
                    if( reldiff .ge. (0.5*eps) ) then
                        casedfails = casedfails + 1
                        if( casedfails .eq. 1 ) then
                            print *, 
     $                "!! Some ABS(x*(1+I)) differ from sqrt(2)*ABS(x)"
                        endif
                        WRITE( 0, fmt = 9999 ) 'd',i, (onehalf*xj),
     $                                  '(1+1*I)', r, answerd
                    endif
                endif
                xj = xj * stepx(i)
            end do
        endif 
  40  continue
*
*     Test (e) Infs
      do 50 i = 1, ninf
        y = cinf(i)
        r = abs( y )
        if( .not.(r .gt. huge(0.0d0)) ) then
            WRITE( *, fmt = 9997 ) 'i',i, y, r
        endif
  50  continue
*
*     Test (f) NaNs
      do 60 i = 1, nnan
        y = cnan(i)
        r = abs( y )
        if( r .eq. r ) then
            WRITE( *, fmt = 9998 ) 'n',i, y, r
        endif
  60  continue
*
*     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) )
     $      print *, "# Please check the failed ABS(a+b*I) in [stderr]"
*
*     .. Formats ..
 9997 FORMAT( '[',a1,i1, '] ABS(', (es8.1,sp,es8.1,"*I"), ' ) = ',
     $        es8.1, ' differs from Inf' )
*
 9998 FORMAT( '[',a1,i1, '] ABS(', (es8.1,sp,es8.1,"*I"), ' ) = ',
     $        es8.1, ' differs from NaN' )
*
 9999 FORMAT( '[',a1,i1, '] ABS(', es24.16e3, ' * ', a7, ' ) = ',
     $         es24.16e3, ' differs from ', es24.16e3 )
*
*     End of zabs
*
      END