      CALL TSQRT(6)
      END
      SUBROUTINE  TSQRT (IOUT)
*     Need this to allow private sqrt() testing
      EXTERNAL     SQRT
      REAL         SQRT
      EXTERNAL     RAN
      REAL         RAN
C     PROGRAM TO TEST SQRT
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                 EPS    - THE SMALLEST POSITIVE FLOATING-POINT
C                          NUMBER SUCH THAT 1.0+EPS .NE. 1.0
C                 EPSNEG - THE SMALLEST POSITIVE FLOATING-POINT
C                          NUMBER SUCH THAT 1.0-EPSNEG .NE. 1.0
C                 XMIN   - THE SMALLEST NON-VANISHING FLOATING-POINT
C                          POWER OF THE RADIX
C                 XMAX   - THE LARGEST FINITE FLOATING-POINT NO.
C
C      RANDL(X) - A FUNCTION SUBPROGRAM RETURNING LOGARITHMICALLY
C                 DISTRIBUTED RANDOM REAL NUMBERS.  IN PARTICULAR,
C                        A * RANDL(ALOG(B/A))
C                 IS LOGARITHMICALLY DISTRIBUTED OVER (A,B)
C
C        RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
C                 NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
C
C
C     STANDARD FORTRAN SUBPROGRAMS REQUIRED
C
C         ABS, ALOG, AMAX1, FLOAT, SQRT
C
C
C     LATEST REVISION - AUGUST 2, 1979
C
C     AUTHOR - W. J. CODY
C              ARGONNE NATIONAL LABORATORY
C
C
      INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,K1,K2,K3,MACHEP,
     1        MAXEXP,MINEXP,N,NEGEP,NGRD
      REAL A,AIT,ALBETA,B,BETA,C,EPS,EPSNEG,ONE,RANDL,R6,R7,
     1     SQBETA,W,X,XMAX,XMIN,XN,X1,Y,Z,ZERO
C
C     IOUT = 6
      CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1            MAXEXP,EPS,EPSNEG,XMIN,XMAX)
      BETA = FLOAT(IBETA)
      SQBETA = SQRT(BETA)
      ALBETA = ALOG(BETA)
      AIT = FLOAT(IT)
      ONE = 1.0E0
      ZERO = 0.0E0
      A = ONE / SQBETA
      B = ONE
      N = 2000
      XN = FLOAT(N)
C-----------------------------------------------------------------
C     RANDOM ARGUMENT ACCURACY TESTS
C-----------------------------------------------------------------
      DO 300 J = 1, 2
         C = ALOG(B/A)
         K1 = 0
         K3 = 0
         X1 = ZERO
         R6 = ZERO
         R7 = ZERO
C
         DO 200 I = 1, N
            X = A * RANDL(C)
            Y = X * X
            Z = SQRT(Y)
            W = (Z - X) / X
            IF (W .GT. ZERO) K1 = K1 + 1
            IF (W .LT. ZERO) K3 = K3 + 1
            W = ABS(W)
            IF (W .LE. R6) GO TO 120
            R6 = W
            X1 = X
  120       R7 = R7 + W * W
  200    CONTINUE
C
         K2 = N - K1 - K3
         R7 = SQRT(R7/XN)
         WRITE (IOUT,1000)
         WRITE (IOUT,1010) N,A,B
         WRITE (IOUT,1011) K1,K2,K3
         WRITE (6,1020) IT,IBETA
         W = -999.0E0
         IF (R6 .NE. ZERO) W = ALOG(ABS(R6))/ALBETA
         WRITE (IOUT,1021) R6,IBETA,W,X1
         W = AMAX1(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
         W = -999.0E0
         IF (R7 .NE. ZERO) W = ALOG(ABS(R7))/ALBETA
         WRITE (IOUT,1023) R7,IBETA,W
         W = AMAX1(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
         A = ONE
         B = SQBETA
  300 CONTINUE
C-----------------------------------------------------------------
C     SPECIAL TESTS
C-----------------------------------------------------------------
      WRITE (IOUT,1040)
      X = XMIN
      Y = SQRT(X)
      WRITE (IOUT,1041) X,Y
      X = ONE - EPSNEG
      Y = SQRT(X)
      WRITE (IOUT,1042) EPSNEG,Y
      X = ONE
      Y = SQRT(X)
      WRITE (IOUT,1043) X,Y
      X = ONE + EPS
      Y = SQRT(X)
      WRITE (IOUT,1044) EPS,Y
      X = XMAX
      Y = SQRT(X)
      WRITE (IOUT,1045) X,Y
C-----------------------------------------------------------------
C     TEST OF ERROR RETURNS
C-----------------------------------------------------------------
      WRITE (IOUT,1050)
      X = ZERO
      WRITE (IOUT,1051) X
      Y = SQRT(X)
      WRITE (IOUT,1055) Y
      X = -ONE
      WRITE (IOUT,1052) X
      Y = SQRT(X)
      WRITE (IOUT,1055) Y
      WRITE (IOUT,1100)
      RETURN
 1000 FORMAT(22H1TEST OF SQRT(X*X) - X //)
 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL /
     1 6X,1H(,E15.4,1H,,E15.4,1H)//)
 1011 FORMAT(19H SQRT(X) WAS LARGER,I6,7H TIMES, /
     1     12X,7H AGREED,I6,11H TIMES, AND /
     2     8X,11HWAS SMALLER,I6,7H TIMES.//)
 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4,
     1    46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER  //)
 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,E15.4,3H = ,I4,3H **,
     1  F7.2/4X,16HOCCURRED FOR X =,E17.6)
 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4,
     1  22H SIGNIFICANT DIGITS IS,F7.2//)
 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,E15.4,
     1    3H = ,I4,3H **,F7.2)
 1040 FORMAT(26H1TEST OF SPECIAL ARGUMENTS//)
 1041 FORMAT(19H SQRT(XMIN) = SQRT(,E15.7,4H) = ,E15.7//)
 1042 FORMAT(25H SQRT(1-EPSNEG) = SQRT(1-,E15.7,4H) = ,E15.7//)
 1043 FORMAT(18H SQRT(1.0) = SQRT(,E15.7,4H) = ,E15.7//)
 1044 FORMAT(22H SQRT(1+EPS) = SQRT(1+,E15.7,4H) = ,E15.7//)
 1045 FORMAT(19H SQRT(XMAX) = SQRT(,E15.7,4H) = ,E15.7//)
 1050 FORMAT(22H1TEST OF ERROR RETURNS//)
 1051 FORMAT(38H SQRT WILL BE CALLED WITH THE ARGUMENT,E15.4/
     1       41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//)
 1052 FORMAT(38H0SQRT WILL BE CALLED WITH THE ARGUMENT,E15.4/
     1       37H THIS SHOULD TRIGGER AN ERROR MESSAGE//)
 1055 FORMAT(24H SQRT RETURNED THE VALUE,E15.4///)
 1100 FORMAT(25H THIS CONCLUDES THE TESTS )
C     ---------- LAST CARD OF SQRT TEST PROGRAM ----------
      END
