      CALL TDSQRT(6)
      END
      SUBROUTINE  TDSQRT (IOUT)
*     Need this to allow private dsqrt() testing
      EXTERNAL DSQRT
      DOUBLE PRECISION DSQRT
      EXTERNAL           DRAN
      DOUBLE PRECISION   DRAN
C     PROGRAM TO TEST DSQRT
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      DRANDL(X) - A FUNCTION SUBPROGRAM RETURNING LOGARITHMICALLY
C                 DISTRIBUTED RANDOM DOUBLE PRECISION NUMBERS.
C                 IN PARTICULAR,
C                        A * DRANDL(DLOG(B/A))
C                 IS LOGARITHMICALLY DISTRIBUTED OVER (A,B)
C
C        DRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM DOUBLE
C                 PRECISION NUMBERS UNIFORMLY DISTRIBUTED OVER
C                 (0,1)
C
C
C     STANDARD FORTRAN SUBPROGRAMS REQUIRED
C
C         DABS, DLOG, DMAX1, DFLOAT, DSQRT
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
      DOUBLE PRECISION A,AIT,ALBETA,B,BETA,C,EPS,EPSNEG,ONE,DRANDL,R6,
     1     R7,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 = DBLE(FLOAT(IBETA))
      SQBETA = DSQRT(BETA)
      ALBETA = DLOG(BETA)
      AIT = DBLE(FLOAT(IT))
      ONE = 1.0D+00
      ZERO = 0.0D+00
      A = ONE / SQBETA
      B = ONE
      N = 2000
      XN = DBLE(FLOAT(N))
C-----------------------------------------------------------------
C     RANDOM ARGUMENT ACCURACY TESTS
C-----------------------------------------------------------------
      DO 300 J = 1, 2
         C = DLOG(B/A)
         K1 = 0
         K3 = 0
         X1 = ZERO
         R6 = ZERO
         R7 = ZERO
C
         DO 200 I = 1, N
            X = A * DRANDL(C)
            Y = X * X
            Z = DSQRT(Y)
            W = (Z - X) / X
            IF (W .GT. ZERO) K1 = K1 + 1
            IF (W .LT. ZERO) K3 = K3 + 1
            W = DABS(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 = DSQRT(R7/XN)
         WRITE (IOUT,1000)
         WRITE (IOUT,1010) N,A,B
         WRITE (IOUT,1011) K1,K2,K3
         WRITE (6,1020) IT,IBETA
         W = -999.0D+00
         IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA
         WRITE (IOUT,1021) R6,IBETA,W,X1
         W = DMAX1(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
         W = -999.0D+00
         IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA
         WRITE (IOUT,1023) R7,IBETA,W
         W = DMAX1(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 = DSQRT(X)
      WRITE (IOUT,1041) X,Y
      X = ONE - EPSNEG
      Y = DSQRT(X)
      WRITE (IOUT,1042) EPSNEG,Y
      X = ONE
      Y = DSQRT(X)
      WRITE (IOUT,1043) X,Y
      X = ONE + EPS
      Y = DSQRT(X)
      WRITE (IOUT,1044) EPS,Y
      X = XMAX
      Y = DSQRT(X)
      WRITE (IOUT,1045) X,Y
C-----------------------------------------------------------------
C     TEST OF ERROR RETURNS
C-----------------------------------------------------------------
      WRITE (IOUT,1050)
      X = ZERO
      WRITE (IOUT,1051) X
      Y = DSQRT(X)
      WRITE (IOUT,1055) Y
      X = -ONE
      WRITE (IOUT,1052) X
      Y = DSQRT(X)
      WRITE (IOUT,1055) Y
      WRITE (IOUT,1100)
      RETURN
 1000 FORMAT(24H1TEST OF DSQRT(X*X) - X //)
 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL /
     1 6X,1H(,D15.4,1H,,D15.4,1H)//)
 1011 FORMAT(20H DSQRT(X) WAS LARGER,I6,7H TIMES, /
     1     13X,7H AGREED,I6,11H TIMES, AND /
     2     9X,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,D15.4,3H = ,I4,3H **,
     1  F7.2/4X,16HOCCURRED FOR X =,D17.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,D15.4,
     1    3H = ,I4,3H **,F7.2)
 1040 FORMAT(26H1TEST OF SPECIAL ARGUMENTS//)
 1041 FORMAT(21H DSQRT(XMIN) = DSQRT(,D25.17,4H) = ,D25.17//)
 1042 FORMAT(27H DSQRT(1-EPSNEG) = DSQRT(1-,D15.7,4H) = ,D25.17//)
 1043 FORMAT(20H DSQRT(1.0) = DSQRT(,D25.17,4H) = ,D25.17//)
 1044 FORMAT(24H DSQRT(1+EPS) = DSQRT(1+,D15.7,4H) = ,D25.17//)
 1045 FORMAT(21H DSQRT(XMAX) = DSQRT(,D25.17,4H) = ,D25.17//)
 1050 FORMAT(22H1TEST OF ERROR RETURNS//)
 1051 FORMAT(39H DSQRT WILL BE CALLED WITH THE ARGUMENT,D15.4/
     1       41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//)
 1052 FORMAT(39H0DSQRT WILL BE CALLED WITH THE ARGUMENT,D15.4/
     1       37H THIS SHOULD TRIGGER AN ERROR MESSAGE//)
 1055 FORMAT(25H DSQRT RETURNED THE VALUE,D15.4///)
 1100 FORMAT(25H THIS CONCLUDES THE TESTS )
C     ---------- LAST CARD OF DSQRT TEST PROGRAM ----------
      END

