C      ALGORITHM 714, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 19, NO. 1, MARCH, 1993, PP. 1-21.
C----------------------------------------------------------------------
C  Program to test CABS
C
C  Accuracy tests are based on the identities
C
C            cabs((3x,4x)) = 5x
C     and
C            cabs((5x,12x)) = 13x
C
C  Data required:
C
C     None
C
C  Subprograms required from this package:
C
C     MACHAR - An environmental inquiry program providing
C              information on the floating point arithmetic
C              system.  Note that the call to MACHAR can
C              be deleted provided the following six
C              parameters are assigned the values indicated:
C
C              IBETA  - the radix of the floating point system;
C              IT     - the number of base-IBETA digits in the
C                       significand of a floating-point number;
C              XMIN   - The smallest non-vanishing floating-point
C                       power of the radix;
C              XMAX   - The largest finite floating-point no.
C
C     REN(K) - a function subprogram returning random real
C              numbers uniformly distributed over (0,1)
C
C
C  Generic Fortran subprograms required:
C
C     ABS, LOG, MAX, REAL (or DBLE), SQRT
C
C  The program as given here must be modified before compiling.
C   If a single (double) precision version is desired, change all
C   occurrences of CS (CD) in columns 1 and 2 to blanks.  For a
C   calibration version, change all occurrences of CS and of CC
C   in columns 1 and 2 to blanks.
C
C  Latest revision - November 7, 1990
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C----------------------------------------------------------------------
      INTEGER I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,
     1        MAXEXP,MINEXP,N,NDUM,NEGEP,NGRD
CS    REAL
      DOUBLE PRECISION
     1   A,AIT,ALL9,ALBETA,B,BETA,CONV,DEL,EIGHT,EPS,EPSNEG,FIVE,FOUR,
     2   ONE,REN,R6,R7,SIXTEN,THREE,TWELVE,TWENTY,W,X,XL,XMAX,XMIN,XN,
     3   X1,Y,Z,ZERO
CS    COMPLEX
      COMPLEX*16
     1   CONVC,CX
CC   2   ,FABS
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
      DATA IOUT/6/
CS    DATA ONE,ZERO,TWENTY,ALL9/1.0E0,0.0E0,2.0E1,-999.0E0/
CS    DATA THREE,FOUR,FIVE,TWELVE/3.0E0,4.0E0,5.0E0,12.0E0/
CS    DATA EIGHT,SIXTEN/8.0E0,16.0E0/
      DATA ONE,ZERO,TWENTY,ALL9/1.0D0,0.0D0,2.0D1,-999.0D0/
      DATA THREE,FOUR,FIVE,TWELVE/3.0D0,4.0D0,5.0D0,12.0D0/
      DATA EIGHT,SIXTEN/8.0D0,16.0D0/
C----------------------------------------------------------------------
C  Statement functions for conversion between integer and float
C----------------------------------------------------------------------
CS    CONV(NDUM) = REAL(NDUM)
CS    CONVC(X,Y) = CMPLX(X,Y)
      CONV(NDUM) = DBLE(NDUM)
      CONVC(X,Y) = DCMPLX(X,Y)
C----------------------------------------------------------------------
      CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1            MAXEXP,EPS,EPSNEG,XMIN,XMAX)
      BETA = CONV(IBETA)
      ALBETA = LOG(BETA)
      AIT = CONV(IT)
      A = ONE
      B = TWENTY
      N = 2000
      XN = CONV(N)
      DEL = (B - A) / XN
      I1 = 0
C-----------------------------------------------------------------
C  Random argument accuracy tests
C-----------------------------------------------------------------
      DO 300 J = 1, 2
         K1 = 0
         K3 = 0
         X1 = ZERO
         R6 = ZERO
         R7 = ZERO
         XL = A
         DO 200 I = 1, N
            X = DEL * REN(I1) + XL
            IF (J .EQ. 1) THEN
                  Y = X * EIGHT
                  Z = X + Y
                  X = Z - Y
                  Y = X * FIVE
                  CX = CONVC(THREE*X,FOUR*X)
               ELSE
                  Y = X * (SIXTEN + SIXTEN)
                  Z = X + Y
                  X = Z - Y
                  Y = X * TWELVE + X
                  CX = CONVC(FIVE*X,TWELVE*X)
            END IF
            Z = ABS(CX)
CC          Z = FABS(CX)
            W = (Z - Y) / Z
C-----------------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------------
            IF (W .GT. ZERO) THEN
                  K1 = K1 + 1
               ELSE IF (W .LT. ZERO) THEN
                  K3 = K3 + 1
            END IF
            W = ABS(W)
            IF (W .GT. R6) THEN
               R6 = W
               X1 = X
            END IF
            R7 = R7 + W * W
            XL = XL + DEL
  200    CONTINUE
C-----------------------------------------------------------------------
C  Process and output statistics
C-----------------------------------------------------------------------
         K2 = N - K1 - K3
         R7 = SQRT(R7/XN)
         IF (J .EQ. 1) THEN
               WRITE (IOUT,1000)
            ELSE
               WRITE (IOUT,1005)
         END IF
         WRITE (IOUT,1010) N,A,B
         WRITE (IOUT,1011) K1,K2,K3
         WRITE (IOUT,1020) IT,IBETA
         IF (R6 .NE. ZERO) THEN
               W = LOG(ABS(R6))/ALBETA
            ELSE
               W = ALL9
         END IF
         WRITE (IOUT,1021) R6,IBETA,W,X1
         W = MAX(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
         IF (R7 .NE. ZERO) THEN
               W = LOG(ABS(R7))/ALBETA
            ELSE
               W = ALL9
         END IF
         WRITE (IOUT,1023) R7,IBETA,W
         W = MAX(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
  300 CONTINUE
C-----------------------------------------------------------------
C  Special tests
C-----------------------------------------------------------------
      WRITE (IOUT,1040)
      X = XMIN
      CX = CONVC(THREE*X,FOUR*X)
      Y = FIVE*X
      WRITE (IOUT,1041)
      Z = ABS(CX)
      WRITE (IOUT,1042) Z,Y
CC    Z = FABS(CX)
C-----------------------------------------------------------------
C  Test of error returns
C-----------------------------------------------------------------
      WRITE (IOUT,1050)
      X = XMAX/SIXTEN
      CX = CONVC(FIVE*X,TWELVE*X)
      WRITE (IOUT,1051) CX
      Z = ABS(CX)
CC    Z = FABS(CX)
      WRITE (IOUT,1055) Z
      X = XMAX/FOUR
      CX = CONVC(THREE*X,XMAX)
      WRITE (IOUT,1052) CX
      Z = ABS(CX)
CC    Z = FABS(CX)
      WRITE (IOUT,1055) Z
      WRITE (IOUT,1100)
      STOP
C-----------------------------------------------------------------
 1000 FORMAT('1Test of CABS((3X,4X)) vs 5X'//)
 1005 FORMAT(//' Test of CABS((5X,12X)) vs 13X'//)
 1010 FORMAT(I7,' Random values of X were drawn from the interval '/
     1 6X,'(',E15.4,',',E15.4,')'//)
 1011 FORMAT(3X,' CABS was larger',I6,' times,'/
     1     12X,' agreed',I6,' times, and'/
     2     8X,'was smaller',I6,' times.'//)
 1020 FORMAT(' There are',I4,' base',I4,
     1    ' significant digits in a floating-point number'//)
 1021 FORMAT(' The maximum relative error of',E15.4,' = ',I4,' **',
     1  F7.2/4X,'occurred for X =',E17.6)
 1022 FORMAT(' The estimated loss of base',I4,
     1  ' significant digits is',F7.2//)
 1023 FORMAT(' The root mean square relative error was',E15.4,
     1    ' = ',I4,' **',F7.2)
 1040 FORMAT(//' Test of special arguments'//)
 1041 FORMAT(' Test of |(3*XMIN,4*XMIN)| vs 5*XMIN.'/
     1       ' This should not trigger an error message'//)
 1042 FORMAT(' CABS((3*XMIN,4*XMIN)) = ',E15.7/
     2       ' 5*XMIN                = ',E15.7//)
 1050 FORMAT(' Test of error returns'//)
 1051 FORMAT(' CABS will be called with the argument',(E15.4,E15.4)/
     1       ' This should not trigger an error message'//)
 1052 FORMAT(' CABS will be called with the argument',(E15.4,E15.4)/
     1       ' This should trigger an error message'//)
 1055 FORMAT(' CABS returned the value ',E15.4///)
 1100 FORMAT(' This concludes the tests')
C---------- Last card of CABS Test Program ----------
      END

CC    COMPLEX FUNCTION FABS(CX)
C----------------------------------------------------------------------
C  Calibration function for cabs tests.
C
C  Latest revision - October 31, 1990
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C----------------------------------------------------------------------
CC    COMPLEX CX
CC    COMPLEX*16 DCX
C----------------------------------------------------------------------
CC    DCX = CX
CC    FABS = ABS(DCX)
CC    RETURN
C---------- Last line of FABS ----------
CC    END
