*BARTLT
      SUBROUTINE BARTLT (VAR,DF,N,CH2,CH2CDF,VARP,DFP,IFLAG)
C
C-----------------------------------------------------------------------
C   BARTLT   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: PERFORMING BARTLETT'S TEST FOR HOMOGENEITY OF VARIANCES ON
C        THREE OR MORE VARIANCES (THE F TEST SHOULD BE USED IN THE CASE
C        OF TWO VARIANCES).  IF THE INPUT PARAMETERS ARE NOT VALID AN 
C        ERROR FLAG IS SET AND NOTHING FURTHER IS COMPUTED, OTHERWISE 
C        THE FOLLOWING ARE COMPUTED: 
C
C           1) THE CHI-SQUARED STATISTIC (CH2),
C           2) THE CUMULATIVE DISTRIBUTION FUNCTION OF THE CHI-SQUARED
C              DISTRIBUTION EVALUATED AT CH2 (CH2CDF), AND
C           3) THE POOLED VARIANCE (VARP) AND ITS CORRESPONDING
C              DEGREES OF FREEDOM (DFP) 
C
C        THE VALUES IN 3) MAY BE USEFUL ONLY IF THE VARIANCES ARE
C        DETERMINED TO BE EQUAL.  THE VALUE OF CH2CDF IS GOOD TO SIX
C        DECIMAL PLACES.
C
C   SUBPROGRAMS CALLED: CDFGAM (GAMMA CUMULATIVE DISTRIBUTION FUNCTION)
C
C   CURRENT VERSION COMPLETED FEBRUARY 3, 1987
C
C   REFERENCES: 
C
C   1) SNEDECOR, GEORGE W. AND COCHRAN, WILLIAM G., 'STATISTICAL
C      METHODS', 6TH EDITION, IOWA STATE UNIVERSITY PRESS, PP. 296-298.
C
C   2) BROWNLEE, K.A., 'STATISTICAL THEORY AND METHODOLOGY IN SCIENCE 
C      AND ENGINEERING', JOHN WILEY & SONS, 1960, PP. 225-227.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C    * VAR = VECTOR (LENGTH N) OF VARIANCES (REAL)
C
C     * DF = VECTOR (LENGTH N) OF DEGREES OF FREEDOM CORRESPONDING
C            TO THE VARIANCES (REAL)
C
C      * N = NUMBER OF VARIANCES [>2] (INTEGER)
C
C      CH2 = THE CHI-SQUARED STATISTIC ASSOCIATED WITH BARTLETT'S TEST
C            (REAL) 
C
C   CH2CDF = THE CUMULATIVE DISTRIBUTION FUNCTION OF THE CHI-SQUARED
C            DISTRIBUTION WITH N-1 DEGREES OF FREEDOM EVALUATED AT CH2
C            (REAL) 
C
C     VARP = THE POOLED VARIANCE DETERMINED FROM THE N VARIANCES (REAL)
C
C      DFP = THE DEGREES OF FREEDOM ASSOCIATED WITH THE POOLED
C            VARIANCE (REAL)
C
C    IFLAG = THE ERROR FLAG ON OUTPUT (INTEGER)   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED
C            1,2 -> ERROR FLAGS FROM SUBROUTINE CDFGAM
C              3 -> N<3
C              4 -> AT LEAST ONE DF(I) IS <= 0.0
C              5 -> AT LEAST ONE VARIANCE(I) IS < 0.0
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION VAR(*),DF(*)
C
C--- CHECK VALIDITY OF INPUT PARAMETERS 
C
      IF (N.LT.3) THEN
         IFLAG = 3
         RETURN
      ENDIF
      DO 10 I = 1, N
         IF (DF(I).LE.0.0) THEN
            IFLAG = 4
            RETURN
         ENDIF
         IF (VAR(I).LT.0.0) THEN
            IFLAG = 5
            RETURN
         ENDIF
   10 CONTINUE
C
C--- COMPUTE NEEDED SUMMATIONS
C
      A = 0.0
      VARP = 0.0
      C = 0.0
      CH2 = 0.0
      DO 20 I = 1, N
         A = A+DF(I)
         VARP = VARP+DF(I)*VAR(I)
         C = C+1.0/DF(I)
         CH2 = CH2+DF(I)*ALOG(VAR(I))
   20 CONTINUE
C
C--- COMPUTE THE POOLED VARIANCE AND ITS DEGREES OF FREEDOM 
C
      VARP = VARP/A 
      DFP = A
C
C--- COMPUTE THE CHI-SQUARED STATISTIC
C
      CH2 = A*ALOG(VARP)-CH2
      A = 1.0+(C-1.0/A)/(3.0*REAL(N-1)) 
      CH2 = CH2/A
C
C--- COMPUTE THE C.D.F. AT CH2
C
      X = 0.5*CH2
      ALPHA = 0.5*REAL(N-1)
      EPS = 0.000001
      CALL CDFGAM (X,ALPHA,EPS,IFLAG,CH2CDF)
      RETURN
      END 
*CDFBET
      SUBROUTINE CDFBET (X,P,Q,EPS,IFLAG,CDFX)
C
C-----------------------------------------------------------------------
C   CDFBET   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE BETA
C        DISTRIBUTION (ALSO KNOWN AS THE INCOMPLETE BETA RATIO) TO A
C        SPECIFIED ACCURACY (TRUNCATION ERROR IN THE INFINITE SERIES).
C        THE ALGORITHM, DESCRIBED IN REFERENCE 2, IS A MODIFICATION OF
C        THE ALGORITHM OF REFERENCE 1.  THREE FEATURES HAVE BEEN ADDED: 
C
C        1) A PRECISE METHOD OF MEETING THE TRUNCATION ACCURACY,
C        2) A CONSTANT W USED IN DETERMINING FOR WHICH X VALUES THE
C           RELATION I(X,P,Q) = 1 - I(1-X,Q,P) IS TO BE USED, AND
C        3) A CONSTANT UFLO >= THE UNDERFLOW LIMIT ON THE COMPUTER.
C
C   SUBPROGRAMS CALLED: DGAMLN (LOG OF GAMMA FUNCTION)
C
C   CURRENT VERSION COMPLETED OCTOBER 24, 1986
C
C   REFERENCES: 
C
C   1) MAJUMDER, K.L. AND BHATTACHARJEE, G.P., 'THE INCOMPLETE BETA
C      INTEGRAL', ALGORITHM AS 63, APPLIED STATISTICS, VOL. 22, NO. 3,
C      1973, PP. 409-411.
C
C   2) REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE BETA C.D.F. 
C      TO A SPECIFIED ACCURACY', STATISTICAL ENGINEERING DIVISION
C      NOTE 86-3, OCTOBER 1986.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C      * X = VALUE AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL)
C
C      * P = FIRST PARAMETER OF THE BETA FUNCTION (>0) (REAL)
C
C      * Q = SECOND PARAMETER OF THE BETA FUNCTION (>0) (REAL)
C
C    * EPS =  THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (>0) (REAL) 
C
C    IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C            0 -> NO ERRORS DETECTED
C            1 -> EITHER P OR Q OR EPS IS <= UFLO 
C            2 -> NUMBER OF TERMS EVALUATED IN THE INFINITE SERIES
C                 EXCEEDS JMAX
C
C     CDFX = THE C.D.F. EVALUATED AT X (REAL)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      LOGICAL LL
      DOUBLE PRECISION DP,DQ,DGAMLN
CCCCC DATA JMAX,W,UFLO / 5000,20.0,1E-100 /
      DATA JMAX,W,UFLO / 5000,20.0,1E-30 /
      CDFX = 0.0
C
C--- CHECK FOR VALIDITY OF ARGUMENTS P, Q, AND EPS
C
      IF (P.LE.UFLO.OR.Q.LE.UFLO.OR.EPS.LE.UFLO) THEN
         IFLAG = 1
         RETURN
      ENDIF
      IFLAG = 0
C
C--- CHECK FOR SPECIAL CASES OF X
C
      IF (X.LE.0.0) RETURN
      IF (X.GE.1.0) THEN
         CDFX = 1.0 
      ELSE
C
C--- SWITCH ARGUMENTS IF NECESSARY
C
         LL = P+W.GE.(P+Q+2.0*W)*X
         IF (LL) THEN
            XY = X
            YX = 1.0-XY
            PQ = P
            QP = Q
         ELSE
            YX = X
            XY = 1.0-YX
            QP = P
            PQ = Q
         ENDIF
C
C--- EVALUATE THE BETA P.D.F. AND CHECK FOR UNDERFLOW
C
         DP = DBLE(PQ-1.0)*DLOG(DBLE(XY))-DGAMLN(PQ)
         DQ = DBLE(QP-1.0)*DLOG(DBLE(YX))-DGAMLN(QP)
         PDFL = SNGL(DGAMLN(PQ+QP)+DP+DQ)
         IF (PDFL.LT.ALOG(UFLO)) THEN
         ELSE
            U = EXP(PDFL)*XY/PQ
            R = XY/YX
   10       IF (QP.LE.1.0) GO TO 20
C
C--- INCREMENT PQ AND DECREMENT QP
C
            IF (U.LE.EPS*(1.0-(PQ+QP)*XY/(PQ+1.0))) GO TO 40
            CDFX = CDFX+U
            PQ = PQ+1.0
            QP = QP-1.0
            U = QP*R*U/PQ
            GO TO 10
   20       V = YX*U
            YXEPS = YX*EPS
C
C--- INCREMENT PQ
C
            DO 30 J = 0, JMAX 
               IF (V.LE.YXEPS) GO TO 40 
               CDFX = CDFX+V
               PQ = PQ+1.0
               V = (PQ+QP-1.0)*XY*V/PQ
   30       CONTINUE
            IFLAG = 2
         ENDIF
   40    IF (.NOT.LL) CDFX = 1.0-CDFX
      ENDIF
      RETURN
      END 
*CDFDNF
      SUBROUTINE CDFDNF (X,DF1,DF2,ALAMB1,ALAMB2,EPS,IFLAG,CDFX)
C
C-----------------------------------------------------------------------
C   CDFDNF   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE DOUBLY 
C        NONCENTRAL F DISTRIBUTION TO A SPECIFIED ACCURACY (TRUNCATION
C        ERROR IN THE INFINITE SERIES REPRESENTATION GIVEN BY EQUATION
C        2.2 IN REFERENCE 1 BELOW).  THE BETA C.D.F. ROUTINE IS CALLED
C        AT MOST TWO TIMES.  FURTHER VALUES OF THE BETA C.D.F. ARE
C        OBTAINED FROM RECURRENCE RELATIONS GIVEN IN REFERENCE 2.
C        REFERENCE 3 GIVES A DETAILED DESCRIPTION OF THE ALGORITHM
C        HEREIN.
C
C        THIS PROGRAM MAY ALSO BE EFFICIENTLY USED TO COMPUTE THE
C        CUMULATIVE DISTRIBUTION FUNCTIONS OF THE SINGLY NONCENTRAL
C        AND CENTRAL F DISTRIBUTIONS BY SETTING THE APPROPRIATE
C        NONCENTRALITY PARAMETERS EQUAL TO ZERO.
C
C        CHECKS ARE MADE TO ASSURE THAT ALL PASSED PARAMETERS ARE
C        WITHIN VALID RANGES AS GIVEN BELOW.  NO UPPER LIMIT IS SET
C        FOR THE NONCENTRALITY PARAMETERS, BUT VALUES UP TO ABOUT
C        10,000 CAN BE HANDLED WITH THE CURRENT DIMENSION LIMITS.  THE
C        COMPUTED VALUE CDFX IS VALID ONLY IF IFLAG=0 ON RETURN.
C
C   NOTE: IN EQUATION 2.2 OF REFERENCE 1 THE AUTHOR HAS MISTAKENLY
C         REVERSED THE ARGUMENTS OF THE INCOMPLETE BETA FUNCTION.
C         THEY SHOULD READ [(M/2)+R,(N/2+S)] WHERE M AND N ARE THE
C         DEGREES OF FREEDOM ASSOCIATED WITH THE NUMERATOR AND
C         DENOMINATOR RESPECTIVELY OF THE F STATISTIC.  TO FURTHER
C         CONFUSE THE ISSUE, THE AUTHOR HAS REVERSED THE USAGE OF
C         M AND N IN SECTION 1 OF THE PAPER.
C
C   NOTE: IN SUBROUTINE EDGEF THE DOUBLE PRECISION CONSTANT DEUFLO IS 
C         THE EXPONENTIAL UNDERFLOW LIMIT WHOSE CURRENT VALUE IS SET
C         AT -69D0.  ON A COMPUTER WHERE DEXP(-69D0) CAUSES UNDERFLOW 
C         THIS LIMIT SHOULD BE CHANGED. 
C
C   SUBPROGRAMS CALLED: CDFBET (BETA C.D.F.)
C                       DGAMLN (DOUBLE PRECISION LOG OF GAMMA FUNCTION)
C                       POISSF, EDGEF (ATTACHED)
C
C   CURRENT VERSION COMPLETED SEPTEMBER 29, 1988
C
C   REFERENCES: 
C
C   1. BULGREN, W.G., 'ON REPRESENTATIONS OF THE DOUBLY NONCENTRAL F
C      DISTRIBUTION', JOURNAL OF THE AMERICAN STATISTICAL ASSOCIATION,
C      MARCH 1971, VOLUME 66, NO. 333, PP. 184-186.
C
C   2. ABRAMOWITZ, MILTON, AND STEGUN, IRENE A., 'HANDBOOK OF
C      MATHEMATICAL FUNCTIONS', NATIONAL BUREAU OF STANDARDS APPLIED
C      MATHEMATICS SERIES 55, NOVEMBER 1970, P. 944.
C
C   3. REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE DOUBLY
C      NONCENTRAL F C.D.F. TO A SPECIFIED ACCURACY', STATISTICAL
C      ENGINEERING DIVISION NOTE 86-4, NOVEMBER 1986.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C        * X = VALUE (>=0) AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL)
C
C      * DF1 = DEGREES OF FREEDOM (>0) IN THE NUMERATOR (REAL)
C
C      * DF2 = DEGREES OF FREEDOM (>0) IN THE DENOMINATOR (REAL)
C
C   * ALAMB1 = THE NONCENTRALITY PARAMETER (>=0) FOR THE NUMERATOR
C              (REAL) [EQUAL TO ZERO FOR THE CENTRAL F DISTRIBUTION]
C
C   * ALAMB2 = THE NONCENTRALITY PARAMETER (>=0) FOR THE DENOMINATOR
C              (REAL) [EQUAL TO ZERO FOR THE SINGLY NONCENTRAL F AND
C              CENTRAL F DISTRIBUTIONS] 
C
C      * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (REAL)
C              [1 >= EPS >= 10**(-10)]
C
C      IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION:  
C                0 -> NO ERRORS DETECTED
C              1,2 -> ERROR FLAGS FROM SUBROUTINE CDFBET
C                3 -> EITHER ALAMB1 OR ALAMB2 IS < 0
C                4 -> EITHER DF1 OR DF2 IS <= 0
C                5 -> EPS IS OUTSIDE THE RANGE [10**(-10),1]
C                6 -> VECTOR DIMENSIONS ARE TOO SMALL - INCREASE NX
C
C       CDFX = THE DOUBLY NONCENTRAL F C.D.F. EVALUATED AT X (REAL)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      PARAMETER (NX=1000)
      DIMENSION BFI(NX),BFJ(NX),POI(NX),POJ(NX)
      CDFX = 0.0
C
C--- CHECK VALIDITY OF ARGUMENTS
C
      IF (ALAMB1.LT.0.0.OR.ALAMB2.LT.0.0) THEN
         IFLAG = 3
         RETURN
      ENDIF
      IF (DF1.LE.0.0.OR.DF2.LE.0.0) THEN
         IFLAG = 4
         RETURN
      ENDIF
      IF (EPS.GT.1.0.OR.EPS.LT.1.0E-10) THEN
         IFLAG = 5
         RETURN
      ENDIF
      IFLAG = 0
C
C--- SET ERROR CRITERION FOR THE BETA C.D.F. (PECULIAR TO CDFBET)
C
      EPS3 = 0.001*EPS
C
      FA = 0.5*ALAMB1
      GA = 0.5*ALAMB2
      FB = 0.5*DF1
      GB = 0.5*DF2
      YY = DF2/(DF2+DF1*X)
      IF (YY.GE.1.0) RETURN
      XX = 1.0-YY
      IF (XX.GE.1.0) THEN
         CDFX = 1.0 
         RETURN
      ENDIF
C
C--- COMPUTE POISSON PROBABILITIES IN VECTORS POI AND POJ
C
      CALL POISSF (FA,EPS,IMIN,NI,POI,NX,IFLAG)
      IF (IFLAG.NE.0) RETURN
      FC = FB+REAL(IMIN)
      CALL POISSF (GA,EPS,JMIN,NJ,POJ,NX,IFLAG)
      IF (IFLAG.NE.0) RETURN
      GC = GB+REAL(JMIN)
C
C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN I=IMIN AND J=JMIN TO JMAX 
C
      CALL EDGEF (NJ,GC,FC,YY,XX,BFJ,CDFX,POJ,POI,EPS3,IFLAG,1)
      IF (NI.LE.1.OR.IFLAG.NE.0) RETURN 
C
C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN J=JMIN AND I=IMIN TO IMAX 
C
      BFI(1) = BFJ(1)
      CALL EDGEF (NI,FC,GC,XX,YY,BFI,CDFX,POI,POJ,EPS3,IFLAG,2)
      IF (NJ.LE.1.OR.IFLAG.NE.0) RETURN 
C
C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN I>IMIN AND J>JMIN
C
      DO 20 I = 2, NI
         BFJ(1) = BFI(I)
         DO 10 J = 2, NJ
            BFJ(J) = XX*BFJ(J)+YY*BFJ(J-1)
            CDFX = CDFX+POI(I)*POJ(J)*BFJ(J)
   10    CONTINUE
   20 CONTINUE
      RETURN
      END 
C
      SUBROUTINE POISSF (ALAMB,EPS,L,NSPAN,V,NV,IFLAG)
C
C--- COMPUTE THE POISSON(ALAMB) PROBABILITIES OVER THE RANGE [L,K]
C--- WHERE THE TOTAL TAIL PROBABILITY IS LESS THAN EPS/2, SUM THE
C--- PROBABILITIES IN DOUBLE PRECISION, AND SHIFT THEM TO THE
C--- BEGINNING OF VECTOR V.
C
      DIMENSION V(*)
      DOUBLE PRECISION DAL,DK,DLIMIT,DSUM,DGAMLN
      DLIMIT = 1.0D0-0.5D0*DBLE(EPS)
      K = INT(ALAMB)
      L = K+1
      IF (ALAMB.EQ.0.0) THEN
         PL = 1.0
      ELSE
         DAL = DBLE(ALAMB)
         DK = DBLE(K)
         PL = SNGL(DEXP(DK*DLOG(DAL)-DAL-DGAMLN(REAL(K+1))))
      ENDIF
      PK = ALAMB*PL/REAL(L)
      NK = NV/2
      NL = NK+1
      DSUM = 0.0
   10 IF (PL.LT.PK) THEN
         NK = NK+1
         IF (NK.GT.NV) THEN
            IFLAG = 6
            RETURN
         ENDIF
         V(NK) = PK 
         DSUM = DSUM+DBLE(PK) 
         K = K+1
         IF (DSUM.GE.DLIMIT) GO TO 20
         PK = ALAMB*PK/REAL(K+1)
      ELSE
         NL = NL-1
         V(NL) = PL 
         DSUM = DSUM+DBLE(PL) 
         L = L-1
         IF (DSUM.GE.DLIMIT) GO TO 20
         PL = REAL(L)*PL/ALAMB
      ENDIF
      GO TO 10
   20 INC = NL-1
      DO 30 I = NL, NK
         V(I-INC) = V(I)
   30 CONTINUE
      NSPAN = NK-INC
      RETURN
      END 
C
      SUBROUTINE EDGEF (NK,FC,GC,XX,YY,BFK,CDFX,POI,POJ,EPS3,IFLAG,L) 
C
C--- COMPUTE THE BETA C.D.F.'S BY A RECURRENCE RELATION ALONG THE EDGES
C--- I = IMIN AND J = JMIN OF A GRID.  THE CORRESPONDING COMPONENTS OF
C--- THE F" C.D.F. ARE INCLUDED IN THE SUMMATION.  TERMS WHICH MIGHT
C--- CAUSE UNDERFLOW ARE SET TO ZERO.
C
      DIMENSION BFK(*),POI(*),POJ(*)
      DOUBLE PRECISION DARG,DEUFLO,DGAMLN
      DATA DEUFLO / -69.0D0 / 
      FD = FC-1.0
      K = MAX0(L,MIN0(NK,INT((GC-1.0)*XX/YY-FD))) 
      FK = FD+REAL(K)
      CALL CDFBET (XX,FK,GC,EPS3,IFLAG,BFK(K))
      IF (IFLAG.NE.0) RETURN
      IF (L.EQ.1) BFK(K) = 1.0-BFK(K)
      IF (NK.EQ.1) GO TO 40
      DARG = DBLE(FK)*DLOG(DBLE(XX))+DBLE(GC)*DLOG(DBLE(YY))-
     *   DLOG(DBLE(FK))+DGAMLN(FK+GC)-DGAMLN(FK)-DGAMLN(GC) 
      IF (DARG.LT.DEUFLO) THEN
         DK = 0.0
      ELSE
         DK = SNGL(DEXP(DARG))*(-1.0)**L
      ENDIF
      IF (K.GE.NK) GO TO 20
      BFK(K+1) = BFK(K)-DK
      DI = DK
      KFLAG = 1
      DO 10 I = K+1, NK-1
         IF (KFLAG.EQ.1) THEN 
            DI = DI*(FD+GC+REAL(I-1))*XX/(FD+REAL(I))
            IF (DK+DI.EQ.DK) THEN
               KFLAG = 0
               DI = 0.0
            ENDIF
         ENDIF
         BFK(I+1) = BFK(I)-DI 
   10 CONTINUE
   20 DI = DK
      KFLAG = 1
      DO 30 I = K-1, L, -1
         IF (KFLAG.EQ.1) THEN 
            DI = DI*(FC+REAL(I))/((FD+GC+REAL(I))*XX)
            IF (DK+DI.EQ.DK) THEN
               KFLAG = 0
               DI = 0.0
            ENDIF
         ENDIF
         BFK(I) = BFK(I+1)+DI 
   30 CONTINUE
   40 DO 50 I = L, NK
         CDFX = CDFX+POI(I)*POJ(1)*BFK(I)
   50 CONTINUE
      RETURN
      END 
*CDFDNT
      SUBROUTINE CDFDNT (X,DF,DELTA,ALAMB,EPS,IFLAG,CDFX)
C
C-----------------------------------------------------------------------
C   CDFDNT   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE DOUBLY 
C        NONCENTRAL T DISTRIBUTION TO A SPECIFIED ACCURACY (TRUNCATION
C        ERROR IN THE INFINITE SERIES REPRESENTATION GIVEN BY EQUATION
C        4 IN REFERENCE 1 BELOW).  WHEN X<0 THE C.D.F. IS COMPUTED
C        FROM CDF(X,DF,DELTA,ALAMB) = 1 - CDF(-X,DF,-DELTA,ALAMB).
C        THE BETA C.D.F. ROUTINE IS CALLED AT MOST FOUR TIMES.  FURTHER
C        VALUES OF THE BETA C.D.F. ARE OBTAINED FROM RECURRENCE
C        RELATIONS GIVEN IN REFERENCE 2.  REFERENCE 3 GIVES A DETAILED
C        DESCRIPTION OF THE ALGORITHM HEREIN.
C
C        THIS PROGRAM MAY ALSO BE EFFICIENTLY USED TO COMPUTE THE
C        CUMULATIVE DISTRIBUTION FUNCTIONS OF THE SINGLY NONCENTRAL
C        AND CENTRAL T DISTRIBUTIONS BY SETTING THE APPROPRIATE
C        NONCENTRALITY PARAMETERS EQUAL TO ZERO.
C
C        CHECKS ARE MADE TO ASSURE THAT ALL PASSED PARAMETERS ARE
C        WITHIN VALID RANGES AS GIVEN BELOW.  NO UPPER LIMIT IS SET
C        FOR THE NONCENTRALITY PARAMETERS, BUT VALUES UP TO ABOUT 100 
C        FOR DELTA AND 10,000 FOR LAMBDA CAN BE HANDLED WITH THE
C        CURRENT DIMENSION LIMITS.  THE COMPUTED VALUE CDFX IS VALID
C        ONLY IF IFLAG=0 ON RETURN.
C
C   NOTE: IN SUBROUTINE EDGET THE DOUBLE PRECISION CONSTANT DEUFLO IS 
C         THE EXPONENTIAL UNDERFLOW LIMIT WHOSE CURRENT VALUE IS SET
C         AT -69D0.  ON A COMPUTER WHERE DEXP(-69D0) CAUSES UNDERFLOW 
C         THIS LIMIT SHOULD BE CHANGED. 
C
C   SUBPROGRAMS CALLED: CDFBET (BETA C.D.F.)
C                       DGAMLN (DOUBLE PRECISION LOG OF GAMMA FUNCTION)
C                       POISST, EDGET, GRID (ATTACHED)
C
C   CURRENT VERSION COMPLETED SEPTEMBER 29, 1988
C
C   REFERENCES: 
C
C   1. KRISHNAN, MARAKATHA, 'SERIES REPRESENTATIONS OF THE DOUBLY
C      NONCENTRAL T DISTRIBUTION', JOURNAL OF THE AMERICAN STATISTICAL
C      ASSOCIATION, SEPTEMBER 1968, VOLUME 63, NO. 323, PP. 1004-1012.
C
C   2. ABRAMOWITZ, MILTON, AND STEGUN, IRENE A., 'HANDBOOK OF
C      MATHEMATICAL FUNCTIONS', NATIONAL BUREAU OF STANDARDS APPLIED
C      MATHEMATICS SERIES 55, NOVEMBER 1970, P. 944.
C
C   3. REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE DOUBLY
C      NONCENTRAL T C.D.F. TO A SPECIFIED ACCURACY', STATISTICAL
C      ENGINEERING DIVISION NOTE 86-5, DECEMBER 1986.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C       * X = VALUE AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL)
C
C      * DF = DEGREES OF FREEDOM (>0) IN THE DENOMINATOR (REAL)
C
C   * DELTA = THE NONCENTRALITY PARAMETER FOR THE NUMERATOR (REAL)
C             [EQUAL TO ZERO FOR THE CENTRAL T DISTRIBUTION]
C
C   * ALAMB = THE NONCENTRALITY PARAMETER (>=0) FOR THE DENOMINATOR
C             (REAL) [EQUAL TO ZERO FOR THE SINGLY NONCENTRAL T AND
C             CENTRAL T DISTRIBUTIONS]
C
C     * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (REAL)
C             [1 >= EPS >= 10**(-10)]
C
C     IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C               0 -> NO ERRORS DETECTED 
C             1,2 -> ERROR FLAGS FROM SUBROUTINE CDFBET
C               3 -> ALAMB IS < 0
C               4 -> DF IS <= 0
C               5 -> EPS IS OUTSIDE THE RANGE [10**(-10),1] 
C               6 -> VECTOR DIMENSIONS ARE TOO SMALL - INCREASE NX
C
C      CDFX = THE DOUBLY NONCENTRAL T C.D.F. EVALUATED AT X (REAL)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      PARAMETER (NX=1000)
      DIMENSION BFI(NX),BFJ(NX),POI(NX),POJ(NX)
      DOUBLE PRECISION DARG,DFA
      LOGICAL LL
      CDFX = 0.0
C
C--- CHECK VALIDITY OF ARGUMENTS
C
      IF (ALAMB.LT.0.0) THEN
         IFLAG = 3
         RETURN
      ENDIF
      IF (DF.LE.0.0) THEN
         IFLAG = 4
         RETURN
      ENDIF
      IF (EPS.GT.1.0.OR.EPS.LT.1.0E-10) THEN
         IFLAG = 5
         RETURN
      ENDIF
      IFLAG = 0
C
C--- SET ERROR CRITERION FOR THE BETA C.D.F. (PECULIAR TO CDFBET)
C
      EPS3 = 0.001*EPS
C
      DELSQ = DELTA**2
      FA = 0.5*DELSQ
      print *,'delta,delsq,fa=',delta,delsq,fa
      GA = 0.5*ALAMB
      GB = 0.5*DF
      YY = DF/(DF+X*X)
      XX = 1.0-YY
C
C--- IF X<0 SET LL=.TRUE., REVERSE SIGN OF DELTA, AND USE THE
C--- IDENTITY DESCRIBED UP FRONT FOR COMPUTING THE C.D.F.
C
      LL = X.LT.0.0 
      IF (XX.GE.1.0) THEN
         CDFX = 1.0 
         GO TO 50
      ENDIF
      SDELTA = DELTA
      IF (LL) SDELTA = -DELTA 
C
C--- COMPUTE POISSON PROBABILITIES IN VECTOR POI
C
      print *,'before poisst, fa=',fa
      CALL POISST (FA,EPS,IMIN,NI,POI,NX,IFLAG)
      print *,'after poisst, fa=',fa
      IF (IFLAG.NE.0) RETURN
      IF (YY.GE.1.0) GO TO 10 
      FC = 0.5+REAL(IMIN)
C
C--- COMPUTE POISSON PROBABILITIES IN VECTOR POJ
C
      CALL POISST (GA,EPS,JMIN,NJ,POJ,NX,IFLAG)
      IF (IFLAG.NE.0) RETURN
      GC = GB+REAL(JMIN)
C
C--- SUM THE TERMS CORRESPONDING TO 'EVEN' VALUES OF INDEX I
C
      CALL GRID (NI,NJ,FC,GC,BFI,BFJ,POI,POJ,XX,YY,EPS3,CDFX,IFLAG)
      IF (IFLAG.NE.0) RETURN
   10 IF (DELTA.EQ.0.0) THEN
         NI = 0
         SUM = 0.0
         IF (YY.GE.1.0) GO TO 40
      ELSE
C
C--- COMPUTE 'POISSON-LIKE' PROBABILITIES IN VECTOR POI
C
         print *,'before k, fa=',fa
         K = INT(FA)
         IF (IMIN.GT.0) THEN
            IMIN = IMIN-1
            NI = NI+1
         ENDIF
         DFA = DBLE(FA)
         DARG = (DBLE(K)+0.5D0)*DLOG(DFA)-DFA-DGAMLN(REAL(K)+1.5)
         print *,'k,dfa,darg=',k,dfa,darg
         L = K-IMIN+1
         POI(L) = SIGN(SNGL(DEXP(DARG)),SDELTA)
         SUM = POI(L)
         DO 20 I = K-1, IMIN, -1
            L = L-1 
            POI(L) = POI(L+1)*(REAL(I)+1.5)/FA
            SUM = SUM+POI(L)
   20    CONTINUE
         L = K-IMIN+1
         DO 30 I = K+1, IMIN+NI-1
            L = L+1 
            POI(L) = POI(L-1)*FA/(REAL(I)+0.5)
            SUM = SUM+POI(L)
   30    CONTINUE
         IF (YY.GE.1.0) GO TO 40
         FC = 1.0+REAL(IMIN)
C
C--- SUM THE TERMS CORRESPONDING TO 'ODD' VALUES OF INDEX I 
C
         CALL GRID (NI,NJ,FC,GC,BFI,BFJ,POI,POJ,XX,YY,EPS3,CDFX,IFLAG)
         print *,'after call grid, cdfx,iflag=',cdfx,iflag
         IF (IFLAG.NE.0) RETURN
      ENDIF
C
C--- COMPUTE THE NORMAL C.D.F. AT -SDELTA
C
   40 PHI = 0.5*(1.0-SUM)
C
C--- COMPUTE THE DOUBLY NONCENTRAL T C.D.F. AT X, USING AN IDENTITY
C--- IF X<0
C
      CDFX = 0.5*CDFX+PHI
   50 IF (LL) CDFX = 1.0-CDFX 
      RETURN
      END 
C
      SUBROUTINE POISST (ALAMB,EPS,L,NSPAN,V,NV,IFLAG)
C
C--- COMPUTE THE POISSON(ALAMB) PROBABILITIES OVER THE RANGE [L,K]
C--- WHERE THE TOTAL TAIL PROBABILITY IS LESS THAN EPS/3, SUM THE
C--- PROBABILITIES IN DOUBLE PRECISION, AND SHIFT THEM TO THE
C--- BEGINNING OF VECTOR V.
C
      DIMENSION V(*)
      DOUBLE PRECISION DAL,DK,DLIMIT,DSUM,DGAMLN
      DLIMIT = 1.0D0-2.0D0*DBLE(EPS)/3.0D0
      K = INT(ALAMB)
      L = K+1
      IF (ALAMB.EQ.0.0) THEN
         PL = 1.0
      ELSE
         DAL = DBLE(ALAMB)
         DK = DBLE(K)
         PL = SNGL(DEXP(DK*DLOG(DAL)-DAL-DGAMLN(REAL(K+1))))
      ENDIF
      PK = ALAMB*PL/REAL(L)
      NK = NV/2
      NL = NK+1
      DSUM = 0.0
   10 IF (PL.LT.PK) THEN
         NK = NK+1
         IF (NK.GT.NV) THEN
            IFLAG = 6
            RETURN
         ENDIF
         V(NK) = PK 
         DSUM = DSUM+DBLE(PK) 
         K = K+1
         IF (DSUM.GE.DLIMIT) GO TO 20
         PK = ALAMB*PK/REAL(K+1)
      ELSE
         NL = NL-1
         V(NL) = PL 
         DSUM = DSUM+DBLE(PL) 
         L = L-1
         IF (DSUM.GE.DLIMIT) GO TO 20
         PL = REAL(L)*PL/ALAMB
      ENDIF
      GO TO 10
   20 INC = NL-1
      DO 30 I = NL, NK
         V(I-INC) = V(I)
   30 CONTINUE
      NSPAN = NK-INC
      RETURN
      END 
C
      SUBROUTINE EDGET (NK,FC,GC,XX,YY,BFK,CDFX,POI,POJ,EPS3,IFLAG,L) 
C
C--- COMPUTE THE BETA C.D.F.'S BY A RECURRENCE RELATION ALONG THE EDGES
C--- I = IMIN AND J = JMIN OF A GRID.  THE CORRESPONDING COMPONENTS OF
C--- THE T" C.D.F. ARE INCLUDED IN THE SUMMATION.  TERMS WHICH MIGHT
C--- CAUSE UNDERFLOW ARE SET TO ZERO.
C
      DIMENSION BFK(*),POI(*),POJ(*)
      DOUBLE PRECISION DARG,DEUFLO,DGAMLN
      DATA DEUFLO / -69.0D0 / 
      FD = FC-1.0
      K = MAX0(L,MIN0(NK,INT((GC-1.0)*XX/YY-FD))) 
      FK = FD+REAL(K)
      CALL CDFBET (XX,FK,GC,EPS3,IFLAG,BFK(K))
      IF (IFLAG.NE.0) RETURN
      IF (L.EQ.1) BFK(K) = 1.0-BFK(K)
      IF (NK.EQ.1) GO TO 40
      DARG = DBLE(FK)*DLOG(DBLE(XX))+DBLE(GC)*DLOG(DBLE(YY))-
     *   DLOG(DBLE(FK))+DGAMLN(FK+GC)-DGAMLN(FK)-DGAMLN(GC) 
      IF (DARG.LT.DEUFLO) THEN
         DK = 0.0
      ELSE
         DK = SNGL(DEXP(DARG))*(-1.0)**L
      ENDIF
      IF (K.GE.NK) GO TO 20
      BFK(K+1) = BFK(K)-DK
      DI = DK
      KFLAG = 1
      DO 10 I = K+1, NK-1
         IF (KFLAG.EQ.1) THEN 
            DI = DI*(FD+GC+REAL(I-1))*XX/(FD+REAL(I))
            IF (DK+DI.EQ.DK) THEN
               KFLAG = 0
               DI = 0.0
            ENDIF
         ENDIF
         BFK(I+1) = BFK(I)-DI 
   10 CONTINUE
   20 DI = DK
      KFLAG = 1
      DO 30 I = K-1, L, -1
         IF (KFLAG.EQ.1) THEN 
            DI = DI*(FC+REAL(I))/((FD+GC+REAL(I))*XX)
            IF (DK+DI.EQ.DK) THEN
               KFLAG = 0
               DI = 0.0
            ENDIF
         ENDIF
         BFK(I) = BFK(I+1)+DI 
   30 CONTINUE
   40 DO 50 I = L, NK
         CDFX = CDFX+POI(I)*POJ(1)*BFK(I)
   50 CONTINUE
      RETURN
      END 
C
      SUBROUTINE GRID (NI,NJ,FC,GC,BFI,BFJ,POI,POJ,XX,YY,EPS3,CDFX,IFLAG
     *   )
C
C--- COMPUTE DOUBLE SUMMATION OF COMPONENTS OF THE T" C.D.F. OVER THE 
C--- GRID I=IMIN TO IMAX AND J=JMIN TO JMAX
C
      DIMENSION BFI(*),BFJ(*),POI(*),POJ(*)
C
C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN I=IMIN, J=JMIN TO JMAX
C
      CALL EDGET (NJ,GC,FC,YY,XX,BFJ,CDFX,POJ,POI,EPS3,IFLAG,1)
      IF (NI.LE.1.OR.IFLAG.NE.0) RETURN 
C
C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN J=JMIN, I=IMIN TO IMAX
C
      BFI(1) = BFJ(1)
      CALL EDGET (NI,FC,GC,XX,YY,BFI,CDFX,POI,POJ,EPS3,IFLAG,2)
      IF (NJ.LE.1.OR.IFLAG.NE.0) RETURN 
C
C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN I>IMIN, J>JMIN
C
      DO 20 I = 2, NI
         BFJ(1) = BFI(I)
         DO 10 J = 2, NJ
            BFJ(J) = XX*BFJ(J)+YY*BFJ(J-1)
            CDFX = CDFX+POI(I)*POJ(J)*BFJ(J)
   10    CONTINUE
   20 CONTINUE
      RETURN
      END 
*CDFF
      SUBROUTINE CDFF (X,DF1,DF2,EPS,IFLAG,CDFX)
C
C-----------------------------------------------------------------------
C   CDFF   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C          DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C          MARYLAND  20899
C
C   FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE F
C        DISTRIBUTION FROM THE CUMULATIVE DISTRIBUTION FUNCTION OF
C        THE BETA DISTRIBUTION.  THE RELATIONSHIP BETWEEN THE TWO
C        DISTRIBUTIONS IS GIVEN IN THE REFERENCE BELOW.
C
C   SUBPROGRAMS CALLED: CDFBET (BETA C.D.F.)
C
C   CURRENT VERSION COMPLETED AUGUST 15, 1987
C
C   REFERENCE: REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE BETA
C              C.D.F. TO A SPECIFIED ACCURACY', STATISTICAL ENGINEERING
C              DIVISION NOTE 86-3, OCTOBER 1986.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C      * X = VALUE AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL)
C
C    * DF1 = FIRST 'DEGREES OF FREEDOM' PARAMETER (>0) (REAL)
C
C    * DF2 = SECOND 'DEGREES OF FREEDOM' PARAMETER (>0) (REAL)
C
C    * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (>0) (REAL)
C
C    IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED
C            1,2 -> ERROR INDICATORS FROM THE BETA C.D.F. ROUTINE
C              3 -> DF1 AND/OR DF2 IS NON-POSITIVE
C
C     CDFX = THE C.D.F. EVALUATED AT X (REAL)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      CDFX = 0.0
C
C--- CHECK FOR NON-POSITIVE DEGREES OF FREEDOM
C
      IF (AMIN1(DF1,DF2).LE.0.0) THEN
         IFLAG = 3
         RETURN
C
      ENDIF
      H = DF1*X
      Y = H/(H+DF2) 
      P = 0.5*DF1
      Q = 0.5*DF2
      CALL CDFBET (Y,P,Q,EPS,IFLAG,CDFX)
      RETURN
C
      END 
*CDFGAM
      SUBROUTINE CDFGAM (X,ALPHA,EPS,IFLAG,CDFX)
C
C-----------------------------------------------------------------------
C   CDFGAM   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MD  20899
C
C   FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE GAMMA
C        DISTRIBUTION (ALSO KNOWN AS THE INCOMPLETE GAMMA RATIO) TO A 
C        SPECIFIED ACCURACY (TRUNCATION ERROR IN THE INFINITE SERIES).
C        THE ALGORITHM, DESCRIBED IN REFERENCE 2, IS A MODIFICATION OF
C        THE ALGORITHM OF REFERENCE 1.  THREE FEATURES HAVE BEEN ADDED: 
C
C        1) A PRECISE METHOD OF MEETING THE TRUNCATION ACCURACY,
C        2) COMPUTATION OF THE UPPER TAIL AREA BY DECREMENTING ALPHA
C           WHEN THAT METHOD IS MORE EFFICIENT, AND
C        3) A CONSTANT UFLO >= THE UNDERFLOW LIMIT ON THE COMPUTER.
C
C   SUBPROGRAMS CALLED: DGAMLN (LOG OF GAMMA FUNCTION)
C
C   CURRENT VERSION COMPLETED OCTOBER 29, 1986
C
C   REFERENCES: 
C
C   1) LAU, CHI-LEUNG, 'A SIMPLE SERIES FOR THE INCOMPLETE GAMMA
C      INTEGRAL', ALGORITHM AS 147, APPLIED STATISTICS, VOL. 29,
C      NO. 1, 1980, PP. 113-114.
C
C   2) REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE GAMMA C.D.F.
C      TO A SPECIFIED ACCURACY', STATISTICAL ENGINEERING DIVISION
C      NOTE 86-2, OCTOBER 1986.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C       * X = VALUE AT WHICH THE C.D.F IS TO BE COMPUTED (REAL)
C
C   * ALPHA = PARAMETER OF THE GAMMA FUNCTION (>0) (REAL)
C
C     * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F (>0) (REAL)
C
C     IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C             0 -> NO ERRORS DETECTED
C             1 -> EITHER ALPHA OR EPS IS <= UFLO 
C             2 -> NUMBER OF TERMS EVALUATED IN THE INFINITE SERIES
C                  EXCEEDS IMAX.
C
C      CDFX = THE C.D.F. EVALUATED AT X (REAL)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      LOGICAL LL
      DOUBLE PRECISION DX,DGAMLN
      DATA IMAX,UFLO / 5000,1.0E-100 /
      CDFX = 0.0
C
C--- CHECK FOR VALIDITY OF ARGUMENTS ALPHA AND EPS
C
      IF (ALPHA.LE.UFLO.OR.EPS.LE.UFLO) THEN
         IFLAG = 1
         RETURN
      ENDIF
      IFLAG = 0
C
C--- CHECK FOR SPECIAL CASE OF X
C
      IF (X.LE.0.0) RETURN
C
C--- EVALUATE THE GAMMA P.D.F. AND CHECK FOR UNDERFLOW
C
      DX = DBLE(X)
      PDFL = SNGL(DBLE(ALPHA-1.0)*DLOG(DX)-DX-DGAMLN(ALPHA))
      IF (PDFL.LT.ALOG(UFLO)) THEN
         IF (X.GE.ALPHA) CDFX = 1.0
      ELSE
         P = ALPHA
         U = EXP(PDFL)
C
C--- DETERMINE WHETHER TO INCREMENT OR DECREMENT ALPHA (A.K.A. P)
C
         LL = .TRUE.
         IF (X.GE.P) THEN
            K = INT(P)
            IF (P.LE.REAL(K)) K = K-1
            ETA = P-REAL(K)
            BL = SNGL(DBLE(ETA-1.0)*DLOG(DX)-DX-DGAMLN(ETA))
            LL = BL.GT.ALOG(EPS)
         ENDIF
         EPSX=EPS/X 
         IF (LL) THEN
C
C--- INCREMENT P
C
            DO 10 I = 0, IMAX 
               IF (U.LE.EPSX*(P-X)) RETURN
               U = X*U/P
               CDFX = CDFX+U
               P = P+1.0
   10       CONTINUE
            IFLAG = 2
         ELSE
C
C--- DECREMENT P
C
            DO 20 J = 1, K
               P = P-1.0
               IF (U.LE.EPSX*(X-P)) GO TO 30
               CDFX = CDFX+U
               U = P*U/X
   20       CONTINUE
   30       CDFX = 1.0-CDFX
         ENDIF
      ENDIF
      RETURN
      END 
*CDFNOR
      SUBROUTINE CDFNOR (Z,EPS,IFLAG,CDFZ)
C
C-----------------------------------------------------------------------
C   CDFNOR   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE STANDARD
C        NORMAL DISTRIBUTION TO A SPECIFIED ACCURACY.  THE C.D.F. OF
C        THE GAMMA DISTRIBUTION IS FIRST COMPUTED, THEN TRANSFORMED TO
C        THE C.D.F. OF THE NORMAL DISTRIBUTION.
C
C   NOTE: TIMING TESTS ON THE CYBER 180/855 COMPUTER AT N.I.S.T.
C         INDICATE THAT THE AVERAGE CPU TIME FOR COMPUTING ONE C.D.F. 
C         IS ABOUT 0.0004 SECOND.
C
C   SUBPROGRAMS CALLED: CDFGAM (C.D.F. OF THE GAMMA DISTRIBUTION)
C
C   CURRENT VERSION COMPLETED APRIL 7, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C     * Z = THE VALUE FOR WHICH THE NORMAL C.D.F. IS TO BE COMPUTED
C           [REAL]
C
C   * EPS = THE ABSOLUTE ACCURACY REQUIREMENT FOR THE C.D.F. [REAL]
C
C   IFLAG = ERROR INDICATOR ON OUTPUT [INTEGER]   INTERPRETATION: 
C           0 -> NO ERRORS DETECTED.
C           1 -> ERROR FLAG FROM CDFGAM 
C           2 -> ERROR FLAG FROM CDFGAM 
C
C    CDFZ = THE C.D.F. OF THE STANDARD NORMAL DISTRIBUTION EVALUATED
C           AT X [REAL]
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DEL = 2.0*EPS 
      IF (Z.EQ.0.0) THEN
         CDFZ = 0.0 
      ELSE
         X = 0.5*Z*Z
         CALL CDFGAM (X,0.5,DEL,IFLAG,CDFX)
         IF (IFLAG.NE.0) RETURN
         IF (Z.GT.0.0) THEN
            CDFZ = 0.5+0.5*CDFX
         ELSE
            CDFZ = 0.5-0.5*CDFX
         ENDIF
      ENDIF
      RETURN
C
      END 
*CDFT
      SUBROUTINE CDFT (X,DF,EPS,IFLAG,CDFX)
C
C-----------------------------------------------------------------------
C   CDFT   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C          DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C          MARYLAND  20899
C
C   FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE T
C        DISTRIBUTION FROM THE CUMULATIVE DISTRIBUTION FUNCTION OF
C        THE BETA DISTRIBUTION.  THE RELATIONSHIP BETWEEN THE TWO
C        DISTRIBUTIONS IS GIVEN IN THE REFERENCE BELOW.
C
C   SUBPROGRAMS CALLED: CDFBET (BETA C.D.F.)
C
C   CURRENT VERSION COMPLETED AUGUST 15, 1987
C
C   REFERENCE: REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE BETA
C              C.D.F. TO A SPECIFIED ACCURACY', STATISTICAL ENGINEERING
C              DIVISION NOTE 86-3, OCTOBER 1986.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C      * X = VALUE AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL)
C
C     * DF = 'DEGREES OF FREEDOM' PARAMETER (>0) (REAL)
C
C    * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (>0) (REAL)
C
C    IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED
C            1,2 -> ERROR INDICATORS FROM THE BETA C.D.F. ROUTINE
C              3 -> DF IS NON-POSITIVE
C
C     CDFX = THE C.D.F. EVALUATED AT X (REAL)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      CDFX = 0.0
C
C--- CHECK FOR NON-POSITIVE DEGREES OF FREEDOM
C
      IF (DF.LE.0.0) THEN
         IFLAG = 3
         RETURN
C
      ENDIF
      H = X*X
      Y = H/(H+DF)
      P = 0.5
      Q = 0.5*DF
      CALL CDFBET (Y,P,Q,EPS,IFLAG,CDFY)
      CDFX = 0.5*(1.0+SIGN(CDFY,X))
      RETURN
C
      END 
*CENSCL
      SUBROUTINE CENSCL (X,W,N1,N2,IOPT,ITER,P1TAIL,CENTER,SCALE)
C
C    *  *  *  *  *  *  *  *  *  * CPR242 *  *  *  *  *  *  *  *  *  *  *
C    CENSCL   WRITTEN BY CHARLES P. REEVE
C    FOR: COMPUTING STATISTICAL ESTIMATES OF THE CENTER AND SCALE
C         OF ENTRIES N1 THROUGH N2 (INCLUSIVE) IN A SET OF DATA.
C         CLASSICAL OR ROBUST ESTIMATES MAY BE SPECIFIED.
C    EXTERNAL SUBPROGRAMS CALLED: SORT1 
C    INTERNAL SUBPROGRAMS CALLED: MAD, MEANSD, MEDIAN, VCOPY
C    CURRENT VERSION COMPLETED MARCH 2, 1983
C    *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
C
C    DEFINITION OF PASSED PARAMETERS: 
C
C    * X(*) =  INPUT: VECTOR (DIMENSION >= N2) CONTAINING DATA
C             OUTPUT: NO CHANGE
C
C      W(*) = WORKSPACE VECTOR (DIMENSION >= N2)
C
C      * N1 = FIRST ENTRY OF INTEREST IN DATA VECTOR X
C
C      * N2 = LAST ENTRY OF INTEREST IN DATA VECTOR X
C
C    * IOPT = 1 TO COMPUTE MEAN AND STANDARD DEVIATION
C             2 TO COMPUTE MEDIAN AND MEDIAN ABSOLUTE DEVIATION
C             3 TO COMPUTE TRIMMED MEAN AND STANDARD DEVIATION
C             4 TO COMPUTE WINSORIZED MEAN AND STANDARD DEVIATION
C             5 TO COMPUTE M(T) AND S(T)  (REEVE) 
C             6 TO COMPUTE M(R) AND S(R)  (REEVE) 
C             7 TO COMPUTE BIWEIGHT MEAN AND STD. DEV. (GROSS)
C
C    * ITER = MAXIMUM NUMBER OF ITERATIONS TO BE USED IN COMPUTATION
C             OF BIWEIGHT ESTIMATE OF LOCATION.  A POSITIVE VALUE WILL
C             CAUSE A DIAGNOSTIC MESSAGE TO BE PRINTED WHEN 'ITER'
C             ITERATIONS HAVE OCCURRED.  A NEGATIVE VALUE WILL CAUSE
C             NO DIAGNOSTIC MESSAGE TO BE PRINTED.  IN EITHER CASE, THE
C             LATEST COMPUTED VALUE OF THE BIWEIGHT ESTIMATE IS RETURNED
C             WHEN 'ABS(ITER)' ITERATIONS HAVE OCCURRED.  CONVERGENCE TO
C             'X(K)' WILL OCCUR EARLIER ON ITERATION 'K' IF 
C
C               ABS( X(K)-X(K-1) )  <  0.0001 * ABS( X(K) ) + .000001 
C
C
C  * P1TAIL = PROPORTION OF ORDERED DATA AT EACH EXTREME TO BE
C             CONSIDERED THE TAIL REGION (USED IF IOPT=3 OR 4)
C
C    CENTER = MEAN                                    IF IOPT=1
C             MEDIAN                                  IF IOPT=2
C             TRIMMED MEAN                            IF IOPT=3
C             WINSORIZED MEAN                         IF IOPT=4
C             M(T)  (REEVE)                           IF IOPT=5
C             M(R)  (REEVE)                           IF IOPT=6
C             BIWEIGHT MEAN                           IF IOPT=7
C
C     SCALE = STANDARD DEVIATION                      IF IOPT=1
C             MEDIAN ABSOLUTE DEVIATION               IF IOPT=2
C             TRIMMED STANDARD DEVIATION              IF IOPT=3
C             WINSORIZED STANDARD DEVIATION           IF IOPT=4
C             S(T)  (REEVE)                           IF IOPT=5
C             S(R)  (REEVE)                           IF IOPT=6
C             BIWEIGHT STD. DEV.  (GROSS)             IF IOPT=7
C
C    * DENOTES VARIABLES REQUIRING INPUT VALUES
C
C    *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
C
      DIMENSION X(*),W(*)
C
C--- COMPUTE NUMBER OF DATA POINTS USED 
C
      N = N2-N1+1
C
C--- CHECK FOR AT LEAST FOUR (4) DATA POINTS -- OTHERWISE, STOP!
C
      IF (N.GE.4) GO TO 10
      PRINT *,' *** NUMBER OF DATA POINTS IS < 4' 
      PRINT *,' *** EXECUTION STOPPED IN SUBROUTINE CENSCL' 
      STOP
C
C--- SET ACCIDENTAL NEGATIVE TAIL PROPORTION TO ZERO
C
   10 P1TAIL = AMAX1(P1TAIL,0.)
C
C--- BRANCH TO CLASSICAL OR ROBUST ESTIMATES
C
      GO TO (20,30,40,50,70,70,130), IOPT
C
C****************************************************************
C*** COMPUTE CLASSICAL ESTIMATES: MEAN AND STANDARD DEVIATION ***
C****************************************************************
C
   20 CALL MEANSD (X,N1,N2,CENTER,SCALE)
      RETURN
C
C**********************************************************************
C*** COMPUTE ROBUST ESTIMATES: MEDIAN AND MEDIAN ABSOLUTE DEVIATION ***
C**********************************************************************
C
   30 CALL VCOPY (X,W,N1,N2)
      CALL MEDIAN (W,N1,N2,CENTER)
      CALL MAD (W,N1,N2,CENTER,SCALE)
      RETURN
C
C*********************************************************************
C*** COMPUTE ROBUST ESTIMATES: TRIMMED MEAN AND STANDARD DEVIATION ***
C*********************************************************************
C
   40 CALL VCOPY (X,W,N1,N2)
      CALL SORT1 (W,N1,N2)
      NTRIM = MIN0((N-1)/2,NINT(P1TAIL*REAL(N)))
      NLO = N1+NTRIM
      NHI = N2-NTRIM
      CALL MEANSD (W,NLO,NHI,CENTER,SCALE)
      RETURN
C
C***************************************************************
C*** COMPUTE ROBUST ESTIMATES: WINSORIZED MEAN AND STD. DEV. ***
C***************************************************************
C
   50 CALL VCOPY (X,W,N1,N2)
      CALL SORT1 (W,N1,N2)
      NTRIM = MIN0((N-1)/2,NINT(P1TAIL*REAL(N)))
      DO 60 I = 1, NTRIM
         W(N1+I-1) = W(N1+NTRIM)
         W(N2-I+1) = W(N2-NTRIM)
   60 CONTINUE
      CALL MEANSD (W,N1,N2,CENTER,SCALE)
      RETURN
C
C****************************************************************
C*** COMPUTE ROBUST ESTIMATES: M(T) AND S(T) OR M(R) AND S(R) ***
C****************************************************************
C
   70 CALL VCOPY (X,W,N1,N2)
      CALL SORT1 (W,N1,N2)
C
C--- SET CONSTANT WHICH RELATES SHORTEST INTERVAL TO STANDARD
C--- DEVIATION OF THE NORMAL DISTRIBUTION
C
      CNORM = 1.348979
C
C--- SET REJECTION LIMITS AS A MULTIPLE OF THE STANDARD DEVIATION
C
      CRHO = 3.0
      RMIN = 1E30
      L = N/2
      NR = (N-1)/2
      B = REAL(L-NR+1)
      A = 4.-B
C
C--- FIND THE SHORTEST INTERVAL COVERING HALF THE POINTS
C
      DO 80 I = 1, NR
         J = N1+I-1 
         K = J+L
         R = (A*(W(K+1)-W(J))+B*(W(K)-W(J+1)))/4. 
         IF (R.GE.RMIN) GO TO 80
         RMIN = R
         K1 = J
   80 CONTINUE
      K2 = K1+1
      K3 = K1+L
      K4 = K3+1
      CENTER = (A*(W(K1)+W(K4))+B*(W(K2)+W(K3)))/8.
      SCALE = RMIN/CNORM
C
C--- RETURN M(R) AND S(R) IF IOPT=6
C
      IF (IOPT.EQ.6) RETURN
      BOUND = CRHO*SCALE
C
C--- TRIM OUTLYING OBSERVATIONS AT LOWER END
C
      NLO = N1
   90 IF (CENTER-W(NLO).LE.BOUND) GO TO 100
      NLO = NLO+1
      GO TO 90
C
  100 NHI = N2
C
C--- TRIM OUTLYING OBSERVATIONS AT UPPER END
C
  110 IF (W(NHI)-CENTER.LE.BOUND) GO TO 120
      NHI = NHI-1
      GO TO 110
C
  120 CALL MEANSD (W,NLO,NHI,CENTER,SCALE)
C
C--- RETURN M(T) AND S(T) FOR IOPT=5
C
      RETURN
C
C**********************************************************************
C*** COMPUTE ROBUST ESTIMATES: BIWEIGHT MEAN AND STANDARD DEVIATION ***
C**********************************************************************
C--- DEFINE CONVERGENCE CONSTANTS
C
  130 C1 = 1E-4
      C2 = 1E-6
      ITERA = IABS(ITER)
C
C--- COUNT NUMBER OF ITERATIONS FOR BIWEIGHT MEAN COMPUTATION
C
      KOUNT = 0
      CALL VCOPY (X,W,N1,N2)
      CALL MEDIAN (W,N1,N2,CENTER)
  140 CALL MAD (W,N1,N2,CENTER,SCALE)
      S1 = 0.
      S2 = 0.
      DO 150 I = N1, N2
         U = (W(I)-CENTER)/(6.*SCALE)
         IF (ABS(U).GE.1.0) GO TO 150
         U = (1.-U*U)**2
         S1 = S1+U*W(I)
         S2 = S2+U
  150 CONTINUE
      COLD = CENTER 
      CENTER = S1/S2
      IF (ABS(CENTER-COLD).LT.C1*ABS(CENTER)+C2) GO TO 160
      KOUNT = KOUNT+1
      IF (KOUNT.LE.ITERA) GO TO 140
      IF (ITER.LT.0) GO TO 160
      N = N2-N1+1
      WRITE (6,1010) ITER,COLD,CENTER,N 
  160 S1 = 0.
      S2 = 0.
      DO 170 I = N1, N2
         U = (W(I)-CENTER)/(9.*SCALE)
         IF (ABS(U).GE.1.0) GO TO 170
         V = 1.-U*U 
         S1 = S1+(W(I)-CENTER)**2*V**4
         S2 = S2+V*(1.-5.*U*U)
  170 CONTINUE
      SCALE = SQRT(REAL(N)*S1/(S2*(S2-1.)))
      RETURN
 1010 FORMAT (1X,'*** ',I3,
     *   ' ITERATION LIMIT ON BIWEIGHT ***    OLD CENTER =',G12.5,3X, 
     *   'NEW CENTER =',G12.5,3X,'SAMPLE SIZE =',I6)
C
      END 
C
C
      SUBROUTINE MEANSD (X,NLO,NHI,AMEAN,SD)
C
C--- COMPUTE CLASSICAL MEAN AND STANDARD DEVIATION OF DATA BETWEEN
C--- POSITIONS X(NLO) AND X(NHI) INCLUSIVE
C
      DIMENSION X(*)
      AMEAN = 0.
      DO 10 I = NLO, NHI
         AMEAN = AMEAN+X(I)
   10 CONTINUE
      AMEAN = AMEAN/REAL(NHI-NLO+1)
      SD = 0.
      DO 20 I = NLO, NHI
         SD = SD+(X(I)-AMEAN)**2
   20 CONTINUE
      SD = SQRT(SD/REAL(NHI-NLO))
      RETURN
C
      END 
C
C
      SUBROUTINE VCOPY (A,B,NLO,NHI)
C
C--- COPY SEGMENT OF VECTOR A INTO SAME SEGMENT OF VECTOR B 
C
      DIMENSION A(*),B(*)
      DO 10 I = NLO, NHI
         B(I) = A(I)
   10 CONTINUE
      RETURN
C
      END 
C
C
      SUBROUTINE MEDIAN (X,NLO,NHI,AMED)
C
C--- COMPUTE MEDIAN OF DATA BETWEEN POSITIONS X(NLO) AND X(NHI)
C--- INCLUSIVE.  DATA IN THIS REGION OF VECTOR X ARE ALTERED.
C
      DIMENSION X(*)
C
C--- SORT REGION OF INTEREST IN VECTOR X
C
      CALL SORT1 (X,NLO,NHI)
C
C--- COMPUTE MEDIAN 
C
      I = (NLO+NHI)/2
      J = (NLO+NHI+1)/2
      AMED = (X(I)+X(J))/2.
      RETURN
C
      END 
C
C
      SUBROUTINE MAD (X,NLO,NHI,C,AMAD) 
C
C--- COMPUTE MEDIAN ABSOLUTE DEVIATION (MAD) OF X FROM C, AN ESTIMATE 
C--- OF THE CENTER OF DATA, BETWEEN X(NLO) AND X(NHI) INCLUSIVE.  VECTOR
C--- X IS EXPECTED TO BE SORTED AND IS UNCHANGED BY THIS SUBROUTINE.
C--- NOTE: IF THE NUMBER OF ENTRIES OF INTEREST, N, IS EVEN THE MAD IS
C--- THE N/2 LARGEST DEVIATION, A SLIGHT OVERESTIMATE.
C
      DIMENSION X(*)
      MLO = NLO
      MHI = NHI
      K = (MHI-MLO)/2
      DO 20 I = 1, K
         IF (X(MHI)+X(MLO).GT.2.0*C) GO TO 10
         MLO = MLO+1
         GO TO 20
C
   10    MHI = MHI-1
   20 CONTINUE
      AMAD = MAX(ABS(X(MHI)-C),ABS(X(MLO)-C))
      RETURN
C
      END 
*CIELIP
      SUBROUTINE CIELIP (X,Y,N,P,YP,M,A,B,XBAR,S,XP,XL,XU,IOPT)
C
C-----------------------------------------------------------------------
C   CIELIP   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, MD 
C
C   FOR: INVERSE PREDICTION FOR A STRAIGHT LINE FIT TO DATA.  THE
C        DEPENDENT VARIABLE IS Y (WITH ERROR) AND THE INDEPENDENT
C        VARIABLE IS X (WITHOUT ERROR).  GIVEN N PAIRS (X,Y) A
C        LEAST SQUARES FIT OF THE MODEL 
C
C              Y(I) = ALPHA + BETA * (X - MEAN(X)) + ERROR
C
C        IS MADE.  FOR A USER SPECIFIED VALUE Y', (THE MEAN OF M
C        MEASUREMENTS WITH THE SAME VARIANCE AS THE Y-VALUES)
C        A PREDICTED X-VALUE, X', IS COMPUTED ALONG WITH ITS
C        100(1-P)% CONFIDENCE INTERVAL.  THE METHOD IS DUE TO
C        CHURCHILL EISENHART [SEE REFERENCE 1 BELOW, ESPECIALLY
C        EQUATIONS (25) AND (26)].  THE CONFIDENCE INTERVAL IS
C        ASYMMETRIC ABOUT X'. 
C
C   NOTE: THE CONFIDENCE INTERVAL DEFINED BY XL AND XU SHOULD BE
C         INTERPRETED AS FOLLOWS: 
C
C         1) IF XL < XU THEN IT IS THE REGION BETWEEN XL AND XU.
C
C         2) IF XL > XU THEN IT IS THE UNION OF THE TWO
C            SEMI-INFINITE RAYS (-INF,XU) AND (XL,+INF).
C
C         3) IF XL = XU = 0 THEN NO CONFIDENCE INTERVAL EXISTS.
C            THIS IS CAUSED BY HIGHLY SCATTERED DATA WHICH IS OF
C            DUBIOUS VALUE FOR INVERSE PREDICTION.
C
C   SUBPROGRAMS CALLED: MDSTI (IMSL) - T PERCENT POINT FUNCTION
C
C   CURRENT VERSION COMPLETED JULY 13, 1984
C
C   REFERENCES: 
C
C   1. EISENHART, CHURCHILL, 'THE INTERPRETATION OF CERTAIN REGRESSION
C      METHODS AND THEIR USE IN BIOLOGICAL AND INDUSTRIAL RESEARCH',
C      THE ANNALS OF MATHEMATICAL STATISTICS, VOLUME X, NO. 2, JUNE
C      1939, PP. 162-186.
C
C   2. REEVE, CHARLES P., 'A COMPARISON OF SOME INVERSE PREDICTION
C      METHODS FOR STRAIGHT LINE DATA', STATISTICAL ENGINEERING
C      DIVISION NOTE 84-1, JUNE 1984.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C      * X(*) = VECTOR (LENGTH N) OF INDEPENDENT VARIABLE
C
C      * Y(*) = VECTOR (LENGTH N) OF DEPENDENT VARIABLE
C
C         * N = NUMBER OF DATA VALUES
C
C         * P = SIGNIFICANCE LEVEL (0.05 FOR 95% C.I., ETC.)
C
C        * YP = Y-VALUE FOR WHICH PREDICTED X-VALUE IS TO BE COMPUTED 
C
C         * M = NUMBER OF MEASUREMENTS INVOLVED IN OBTAINING YP.  IF
C               YP IS EXACT, SET M=0
C
C           A = ESTIMATE OF ALPHA (CONSTANT TERM) 
C
C           B = ESTIMATE OF BETA (SLOPE)
C
C        XBAR = MEAN OF THE N X-VALUES
C
C           S = RESIDUAL STANDARD DEVIATION
C
C          XP = PREDICTED X-VALUE CORRESPONDING TO YP
C
C          XL = LOWER BOUND OF CONFIDENCE INTERVAL FOR XP
C
C          XU = UPPER BOUND OF CONFIDENCE INTERVAL FOR XP
C
C      * IOPT = 0 IF THE QUANTITIES A, B, XBAR, AND S ARE TO BE COMPUTED
C                 FROM A LEAST SQUARES FIT USING INPUT VALUES OF X(*),
C                 Y(*), AND N.
C             = 1 IF INPUT VALUES OF THE QUANTITIES A, B, XBAR, S, AND
C                 N ARE TO BE USED
C
C     * DENOTES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION X(*),Y(*)
      RN = REAL(N)
C
C--- BRANCH ACCORDING TO IOPT 
C
      IF (IOPT.GT.0) GO TO 40 
C
C--- COMPUTE LEAST SQUARES ESTIMATES OF COEFFICIENTS
C
      SX = 0.0
      SY = 0.0
      DO 10 I = 1, N
         SX = SX+X(I)
         SY = SY+Y(I)
   10 CONTINUE
      A = SY/RN
      XBAR = SX/RN
      SXX = 0.0
      SXY = 0.0
      DO 20 I = 1, N
         SXX = SXX+(X(I)-XBAR)**2
         SXY = SXY+(X(I)-XBAR)*Y(I)
   20 CONTINUE
      R = 1.0/SXX
      B = R*SXY
C
C--- COMPUTE RESIDUAL STANDARD DEVIATION
C
      S = 0.0
      DO 30 I = 1, N
         S = S+(Y(I)-A-B*(X(I)-XBAR))**2
   30 CONTINUE
      DF = REAL(N-2)
      S = SQRT(S/DF)
C
C--- COMPUTE T CRITICAL VALUE 
C
CALAN CALL MDSTI (P,DF,T,IER) 
C
C--- COMPUTE INVERSE OF M (UNLESS M=0)
C
   40 IF (M.EQ.0) THEN
         AIM = 0.0
      ELSE
         AIM = 1.0/ABS(REAL(M))
      ENDIF
C
C--- COMPUTE PREDICTED X-VALUE FOR YP (SET XP = 0 IF SLOPE = 0)
C
      IF (B.EQ.0.0) THEN
         XP = 0.0
      ELSE
         XP = XBAR+(YP-A)/B
      ENDIF
C
C--- COMPUTE CONFIDENCE INTERVAL FOR PREDICTED X-VALUE
C
      Q = B**2-R*(T*S)**2
      RADICL = R*(YP-A)**2+(AIM+1.0/RN)*Q
      IF (RADICL.LT.0.0) THEN 
         XU = 0.0
         XL = 0.0
      ELSE
         TERM1 = XBAR+B*(YP-A)/Q
         TERM2 = T*S*SQRT(RADICL)/Q
         XL = TERM1-TERM2
         XU = TERM1+TERM2
      ENDIF
      RETURN
C
      END 
*DGAMLN
      DOUBLE PRECISION FUNCTION DGAMLN (X)
C
C-----------------------------------------------------------------------
C   DGAMLN   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING THE DOUBLE PRECISION LOG OF THE GAMMA FUNCTION WITH
C        SINGLE PRECISION PARAMETER X>0.  THE MAXIMUM TRUNCATION ERROR
C        IN THE INFINITE SERIES (SEE REFERENCE 1) IS DETERMINED BY THE
C        CONSTANT XMIN.  WHEN X<XMIN A RECURRENCE RELATION IS USED IN 
C        ORDER TO ACHIEVE THE REQUIRED ABSOLUTE ACCURACY.  THE TABLE
C        BELOW GIVES THE MINIMUM VALUE OF X WHICH YIELDS THE CORRES-
C        PONDING ABSOLUTE ACCURACY IN DGAMLN(X) ASSUMING THE MACHINE
C        CARRIES ENOUGH DIGITS WHEN THOSE TO THE LEFT OF THE DECIMAL
C        ARE CONSIDERED (SEE REFERENCE 2 FOR FURTHER DISCUSSION).  IF 
C        THE LATTER CONDITION IS NOT MET, AN ERROR MESSAGE IS PRINTED.
C
C        THE CYBER 180/855 AT NBS CARRIES ABOUT 15 DIGITS IN SINGLE
C        PRECISION, THEREFORE THE PRE-SET VALUE OF ABSACC IS 10**(-15)
C        AND THE CORRESPONDING VALUE OF XMIN IS 6.894.  ON A DIFFERENT
C        MACHINE THESE CONSTANTS SHOULD BE CHANGED ACCORDINGLY.
C
C           XMIN    ACCURACY    XMIN    ACCURACY    XMIN    ACCURACY
C          ------   --------   ------   --------   ------   --------
C           1.357     1E-3      4.592     1E-12    15.539     1E-21
C           2.037     1E-6      6.894     1E-15    23.330     1E-24
C           3.059     1E-9     10.351     1E-18    35.025     1E-27
C
C   NOTE: THIS IS EXACTLY THE SAME SOFTWARE AS SUBROUTINE DGAMLN
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED MAY 1, 1989
C
C   REFERENCES: 
C
C   1) ABRAMOWITZ, MILTON AND STEGUN, IRENE, 'HANDBOOK OF MATHEMATICAL
C      FUNCTIONS', NBS APPLIED MATHEMATICS SERIES 55, NOV. 1970,
C      EQ. 6.1.40, P 257.
C
C   2) REEVE, CHARLES P., 'ACCURATE COMPUTATION OF THE LOG OF THE GAMMA
C      FUNCTION', STATISTICAL ENGINEERING DIVISION NOTE 86-1, OCTOBER 
C      1986.
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION ABSACC,B1,B2,B3,B4,B5,B6,B7,B8,C,DX,Q,R,XMIN,XN
      DATA XMIN,ABSACC / 6.894D0,1.0D-15 /
      DATA C / 0.918938533204672741780329736D0 /
      DATA B1 / 0.833333333333333333333333333D-1 /
      DATA B2 / -0.277777777777777777777777778D-2 /
      DATA B3 / 0.793650793650793650793650794D-3 /
      DATA B4 / -0.595238095238095238095238095D-3 /
      DATA B5 / 0.841750841750841750841750842D-3 /
      DATA B6 / -0.191752691752691752691752692D-2 /
      DATA B7 / 0.641025641025641025641025641D-2 /
      DATA B8 / -0.295506535947712418300653595D-1 /
C
C--- TERMINATE EXECUTION IF X<=0.0
C
      IF (X.LE.0.0) STOP  '*** X<=0.0 IN FUNCTION DGAMLN ***'
      DX = DBLE(X)
      N = MAX0(0,INT(XMIN-DX+1.0D0))
      XN = DX+DBLE(N)
      R = 1.0D0/XN
      Q = R*R
      DGAMLN = R*(B1+Q*(B2+Q*(B3+Q*(B4+Q*(B5+Q*(B6+Q*(B7+Q*B8)))))))+C+
     *   (XN-0.5D0)*DLOG(XN)-XN
C
C--- USE RECURRENCE RELATION WHEN N>0 (X<XMIN)
C
      IF (N.GT.0) THEN
         Q = 1.0D0
         DO 20 I  =  0, N-1
            Q = Q*(DX+DBLE(I))
   20    CONTINUE
         DGAMLN=DGAMLN-DLOG(Q)
      ENDIF
C
C--- PRINT WARNING IF ABSOLUTE ACCURACY HAS NOT BEEN ATTAINED
C
      IF (DGAMLN+ABSACC.EQ.DGAMLN) THEN 
CCCCC    PRINT *,' ********* WARNING FROM FUNCTION DGAMLN *********'
CCCCC    PRINT *,' REQUIRED ABSOLUTE ACCURACY NOT ATTAINED FOR X = ',X
      ENDIF
      RETURN
      END 
*DNCMLN
      DOUBLE PRECISION FUNCTION DNCMLN (N,M)
C
C-----------------------------------------------------------------------
C   DNCMLN   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: COMPUTING THE DOUBLE PRECISION NATURAL LOGARITHM OF
C        'N CHOOSE M' = N!/[M!(N-M)!] WHERE N>=0, M>=0, AND M<=N.  THAT
C        VALUE IS RETURNED IN DNCMLN.
C
C   SUBPROGRAMS CALLED: DGAMLN (NATURAL LOGARITHM OF GAMMA FUNCTION)
C
C   CURRENT VERSION COMPLETED MAY 1, 1989
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION D1,D2,D3,DGAMLN
      IF (N.LT.0) STOP '*** STOP: N<0 IN SUBROUTINE DNCMLN ***'
      IF (M.LT.0) STOP '*** STOP: M<0 IN SUBROUTINE DNCMLN ***'
      IF (M.GT.N) STOP '*** STOP: M>N IN SUBROUTINE DNCMLN ***'
      R1 = REAL(N+1)
      R2 = REAL(M+1)
      R3 = REAL(N-M+1)
      D1 = DGAMLN(R1)
      D2 = DGAMLN(R2)
      D3 = DGAMLN(R3)
      DNCMLN = D1+D2-D3
      RETURN
C
      END 
*DPR1LN
      DOUBLE PRECISION FUNCTION DPR1LN (N1,M1,N2,M2)
C
C-----------------------------------------------------------------------
C   DPR1LN   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: COMPUTING THE DOUBLE PRECISION NATURAL LOGARITHM OF
C
C                  N1!/[M1!(N1-M1)!] N2!/[M2!(N2-M2)!]
C                  -----------------------------------
C                   (N1+N2)!/[(M1+M2)!(N1+N2-M1-M2)!]
C
C        WHERE N1>=0, M1>=0, M1<=N1, N2>=0, M2>=0, AND M2<=N2.  THAT
C        VALUE IS RETURNED IN DPR1LN.
C
C   SUBPROGRAMS CALLED: DNCMLN (NATURAL LOGARITHM OF 'N CHOOSE M')
C
C   CURRENT VERSION COMPLETED MAY 1, 1989
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION D1,D2,D3,DNCMLN
      IF (N1.LT.0) STOP '*** STOP: N1<0 IN SUBROUTINE DPR1LN ***'
      IF (M1.LT.0) STOP '*** STOP: M1<0 IN SUBROUTINE DPR1LN ***'
      IF (M1.GT.N1) STOP '*** STOP: M1>N1 IN SUBROUTINE DPR1LN ***'
      IF (N2.LT.0) STOP '*** STOP: N2<0 IN SUBROUTINE DPR1LN ***'
      IF (M2.LT.0) STOP '*** STOP: M2<0 IN SUBROUTINE DPR1LN ***'
      IF (M2.GT.N2) STOP '*** STOP: M2>N2 IN SUBROUTINE DPR1LN ***'
      D1 = DNCMLN(N1,M1)
      D2 = DNCMLN(N2,M2)
      D3 = DNCMLN(N1+N2,M1+M2)
      DPR1LN = D1+D2-D3
      RETURN
C
      END 
*DWESD
      SUBROUTINE DWESD (X,LDX,N,M,U,LDU,V,LDV,E,SD)
C
C-----------------------------------------------------------------------
C   DWESD   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C           DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C           MARYLAND  20899
C
C   FOR: COMPUTING THE EXPECTED VALUE AND VARIANCE OF THE DURBIN-WATSON
C        STATISTIC D AS APPLIED TO THE RESIDUALS FROM THE LEAST SQUARES
C        FIT OF THE MODEL  Y = XB + ERROR.  FOR A GIVEN MATRIX X, EXACT
C        FORMULAS FOR E(D) AND VAR(D) ARE GIVEN IN THE REFERENCES BELOW.
C
C   SUBPROGRAMS CALLED: MATXXI (STSPAC) 
C                       XTAX, XTA2X, TRACES (ATTACHED)
C
C   CURRENT VERSION COMPLETED NOVEMBER 6, 1985
C
C   REFERENCES: 
C
C   1. DURBIN, J. AND WATSON, G.S., 'TESTING FOR SERIAL CORRELATION IN
C      LEAST SQUARES REGRESSION. I', BIOMETRIKA 37, 1950, PP. 409-428 
C
C   2. DURBIN, J. AND WATSON, G.S., 'TESTING FOR SERIAL CORRELATION IN
C      LEAST SQUARES REGRESSION. II', BIOMETRIKA 38, 1951, PP. 159-178
C
C   3. DURBIN, J. AND WATSON, G.S., 'TESTING FOR SERIAL CORRELATION IN
C      LEAST SQUARES REGRESSION. III', BIOMETRIKA 58, 1971, PP. 1-19
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C     * X = MATRIX (SIZE NXM) IN THE MODEL Y = XB + ERROR
C
C   * LDX = LEADING DIMENSION OF X (LDX >=N+2)
C
C     * N = NUMBER OF EQUATIONS AND NUMBER OF ROWS OF X
C
C     * M = NUMBER OF UNKNOWNS AND NUMBER OF COLUMNS OF X
C
C     * U = MATRIX (SIZE MXM) FOR INTERMEDIATE COMPUTATIONS 
C
C   * LDU = LEADING DIMENSION OF U (LDU >= M)
C
C     * V = MATRIX (SIZE MXM) FOR INTERMEDIATE COMPUTATIONS 
C
C   * LDV = LEADING DIMENSION OF V (LDV >= M)
C
C       E = EXPECTED VALUE OF DURBIN-WATSON STATISTIC FOR X 
C
C      SD = STANDARD DEVIATION OF DURBIN-WATSON STATISTIC FOR X
C
C   * DENOTES VARIABLES REQUIRING INPUT VALUES
C-----------------------------------------------------------------------
C
      DIMENSION X(LDX,*),U(LDU,*),V(LDV,*)
C
C--- CHECK FOR PARAMETERS OUT OF RANGE
C
      IF (M.GE.N.OR.N.LT.5) THEN
         WRITE (6,1000)
         STOP
C
      ENDIF
C
C--- COMPUTE U = X'AX
C
      CALL XTAX (X,LDX,U,LDU,N,M)
C
C--- COMPUTE V = X'AAX
C
      CALL XTA2X (X,LDX,V,LDV,N,M)
C
C--- COMPUTE W = INV(X'X)  (STORE IN X) 
C
      CALL MATXXI (X,LDX,N,M) 
C
C--- COMPUTE TRACE(UW) AND TRACE(UWUW)
C
      CALL TRACES (U,LDU,X,LDX,M,M,TR1,TR2,2)
C
C--- COMPUTE TRACE(VW)
C
      CALL TRACES (V,LDV,X,LDX,M,M,TR3,DUM,1)
C
C--- COMPUTE CONSTANTS P AND Q
C
      P = REAL(2*N-2)-TR1
      Q = REAL(6*N-8)-2.0*TR3+TR2
C
C--- COMPUTE EXPECTED VALUE AND STANDARD DEVIATION OF D
C
      E = P/REAL(N-M)
      SD = SQRT(2.0*(Q-P*E)/REAL((N-M)*(N-M+2)))
      RETURN
C
C
 1000 FORMAT (/1X,'*** FATAL ERROR -- MUST HAVE N>4 AND N>M IN DWESD'/)
C
      END 
C
      SUBROUTINE XTAX (X,LDX,Z,LDZ,N,M) 
C
C-----------------------------------------------------------------------
C   XTAX   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C          DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C           MARYLAND  20899
C
C   FOR: COMPUTING X'AX WHERE X IS NXM AND A IS AN NXN TRIDIAGONAL
C        MATRIX ASSOCIATED WITH THE DURBIN-WATSON STATISTIC.  MATRIX
C        A HAS THE FORM
C
C               1 -1  0  0  ...  0  0
C              -1  2 -1  0  ...  0  0
C               0 -1  2 -1  ...  0  0
C                    .      . 
C                    .       .
C                    .        .
C               0  0  0  0  ...  2 -1
C               0  0  0  0  ... -1  1 . 
C
C        THE COMPUTATION TAKES THE FORM Z = (L'X)'(L'X) WHERE L IS
C        THE CHOLESKY FACTORIZATION OF A (A=LL'). 
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED NOVEMBER 6, 1985
C-----------------------------------------------------------------------
C
      DIMENSION X(LDX,*),Z(LDZ,*)
      DO 30 I = 1, M
         DO 20 J = 1, I
            S = 0.0 
            DO 10 K = 2, N
               S = S+(X(K-1,I)-X(K,I))*(X(K-1,J)-X(K,J))
   10       CONTINUE
            Z(I,J) = S
            Z(J,I) = S
   20    CONTINUE
   30 CONTINUE
      RETURN
C
      END 
C
      SUBROUTINE XTA2X (X,LDX,Z,LDZ,N,M)
C
C-----------------------------------------------------------------------
C   XTA2X   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C           DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C           MARYLAND  20899
C
C   FOR: COMPUTING X'AAX WHERE X IS NXM AND AA IS AN NXN PENTADIAGONAL
C        MATRIX ASSOCIATED WITH THE DURBIN-WATSON STATISTIC.  MATRIX
C        AA HAS THE FORM
C
C               2 -3  1  0  0  ...  0  0  0
C              -3  6 -4  1  0  ...  0  0  0
C               1 -4  6 -4  1  ...  0  0  0
C               0  1 -4  6 -4  ...  0  0  0
C                     .        .
C                     .         .
C                     .          .
C               0  0  0  0  0  ...  6 -4  1
C               0  0  0  0  0  ... -4  6 -3
C               0  0  0  0  0  ...  1 -3  2 .
C
C        THE COMPUTATION TAKES THE FORM Z = (AX)'(AX) SINCE A=A'.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED NOVEMBER 6, 1985
C-----------------------------------------------------------------------
C
      DIMENSION X(LDX,*),Z(LDZ,*)
      DO 30 I = 1, M
         DO 20 J = 1, I
            S = (X(1,I)-X(2,I))*(X(1,J)-X(2,J))+(X(N,I)-X(N-1,I))*(X(N,J
     *         )-X(N-1,J))
            DO 10 K = 2, N-1
               S = S+(-X(K-1,I)+2.0*X(K,I)-X(K+1,I))*(-X(K-1,J)+2.0*X(K,
     *            J)-X(K+1,J))
   10       CONTINUE
            Z(I,J) = S
            Z(J,I) = S
   20    CONTINUE
   30 CONTINUE
      RETURN
C
      END 
C
      SUBROUTINE TRACES (U,LDU,V,LDV,N,M,TR1,TR2,IND)
C
C-----------------------------------------------------------------------
C   TRACES   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING THE TRACES OF THE NXN MATRICES UV AND UVUV WHERE
C        U IS NXM AND V IS MXN.  NOTE THAT NEITHER UV NOR UVUV
C        ACTUALLY NEEDS TO BE COMPUTED. 
C
C   NOTE: IF IND=1 THEN ONLY TRACE(UV) IS COMPUTED
C         IF IND=2 THEN BOTH TRACE(UV) AND TRACE(UVUV) ARE COMPUTED
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED NOVEMBER 5, 1985
C-----------------------------------------------------------------------
C
      DIMENSION U(LDU,*),V(LDV,*)
      TR1 = 0.0
      DO 20 I = 1, N
         S = 0.0
         DO 10 K = 1, M
            S = S+U(I,K)*V(K,I)
   10    CONTINUE
         TR1 = TR1+S
   20 CONTINUE
      TR2 = 0.0
      IF (IND.NE.2) RETURN
      DO 50 I = 1, N
         DO 40 J = 1, N
            S1 = 0.0
            S2 = 0.0
            DO 30 K = 1, M
               S1 = S1+U(I,K)*V(K,J)
               S2 = S2+U(J,K)*V(K,I)
   30       CONTINUE
            TR2 = TR2+S1*S2
   40    CONTINUE
   50 CONTINUE
      RETURN
C
      END 
C
*ELLPTS
      SUBROUTINE ELLPTS (E,N,X,Y,IFLAG) 
C
C-----------------------------------------------------------------------
C   ELLPTS   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899. 
C
C   FOR: COMPUTING N REGULARLY-SPACED POINTS ON THE PERIMETER OF THE
C        ELLIPSE DEFINED BY
C
C                               [E(1)  E(3)] [X-E(4)]
C              [X-E(4)  Y-E(5)] [          ] [      ] = E(6) .
C                               [E(3)  E(2)] [Y-E(5)]
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED JANUARY 5, 1987
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C     * E = VECTOR (LENGTH 6) OF CONSTANTS DEFINING ELLIPSE (REAL)
C
C     * N = NUMBERS OF POINTS TO BE GENERATED ON PERIMETER OF ELLIPSE 
C           (INTEGER)
C
C       X = VECTOR (LENGTH N) OF X COORDINATES OF POINTS ON PERIMETER 
C           OF ELLIPSE (REAL) 
C
C       Y = VECTOR (LENGTH N) OF Y COORDINATES OF POINTS ON PERIMETER 
C           OF ELLIPSE (REAL) 
C
C   IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C           0 -> NO ERRORS DETECTED
C           1 -> N<1
C           2 -> E(1) OR E(2) OR E(6) IS NON-POSITIVE
C           3 -> E(1)*E(2)<=E(3)**2
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION E(*),X(*),Y(*)
      IF (N.LE.0) THEN
         IFLAG = 1
      ELSEIF (AMIN1(E(1),E(2),E(6)).LE.0.0) THEN
         IFLAG = 2
      ELSEIF (E(1)*E(2).LE.E(3)**2) THEN
         IFLAG = 3
      ELSE
         U = 2.0*E(3)
         V = E(1)-E(2)
         IF (U.EQ.0.0.AND.V.EQ.0.0) THEN
            THETA = 0.0
         ELSE
            THETA = 0.5*ATAN2(U,V)
         ENDIF
         CT = COS(THETA)
         ST = SIN(THETA)
         Q = SQRT(V**2+U**2)
         U = 2.0*E(6)
         V = E(1)+E(2)
         A = SQRT(U/(V+Q))
         B = SQRT(U/(V-Q))
         C = 8.0*ATAN(1.0)/REAL(N)
         DO 10 I = 1, N
            PHI = C*REAL(I)
            CP = COS(PHI)
            SP = SIN(PHI)
            X(I) = A*CT*CP-B*ST*SP+E(4) 
            Y(I) = A*ST*CP+B*CT*SP+E(5) 
   10    CONTINUE
         IFLAG = 0
      ENDIF
      RETURN
      END 
*FACTOR
      SUBROUTINE FACTOR (N,NFAC,L,IFLAG)
C
C-----------------------------------------------------------------------
C   FACTOR   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: FINDING THE PRIME FACTORS OF THE ABSOLUTE VALUE OF N, AN
C        INTEGER OF 10 DIGITS OR LESS.  THE FACTORS ARE RETURNED IN
C        THE VECTOR NFAC WITH L INDICATING THE NUMBER OF FACTORS.
C
C        **********************************************************
C        *** IN THE CALLING PROGRAM THE VECTOR NFAC SHOULD BE   ***
C        *** DIMENSIONED AT LEAST 33 TO ASSURE THAT IT WILL NOT ***
C        *** BE OVER-INDEXED IN THIS ROUTINE (2**34 > 10**10).  ***
C        **********************************************************
C
C        UPON THE FIRST CALL TO THIS SUBROUTINE, ALL PRIME INTEGERS IN
C        THE RANGE [2,100003], IN ASCENDING ORDER, ARE READ FROM A FILE
C        NAMED 'PRIMES' INTO THE VECTOR IPRM.  N IS THEN 'FACTORED BY 
C        DIVISION' USING THESE PRIMES.  THE ALGORITHM USED IS A SLIGHT
C        MODIFICATION OF ALGORITHM 'A' IN THE REFERENCE.  THE RETURNED
C        VALUES FROM THIS ROUTINE ARE VALID ONLY IF THE ERROR FLAG
C        (IFLAG) IS ZERO ON RETURN.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   DATA FILES CALLED: PRIMES 
C
C   CURRENT VERSION COMPLETED FEBRUARY 4, 1987
C
C   REFERENCE: KNUTH, DONALD E., 'THE ART OF COMPUTER PROGRAMMING',
C              VOL. 2/SEMINUMERICAL ALGORITHMS, 2ND EDITION, ADDISON- 
C              WESLEY PUBLISHING CO., 1981, PP. 364-398.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C     * N = AN INTEGER OF 10 DIGITS OR LESS.  THE ABSOLUTE VALUE OF N 
C           WILL BE FACTORED INTO PRIMES.
C
C    NFAC = A VECTOR (LENGTH L) OF PRIME FACTORS OF ABS(N) (INTEGER)
C
C       L = THE NUMBER OF PRIME FACTORS OF N (INTEGER)
C
C   IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C           0 -> NO ERRORS DETECTED
C           1 -> N HAS MORE THAN 10 DIGITS
C           2 -> AN END-OF-FILE ENCOUNTERED IN READING FILE 'PRIMES'
C           3 -> A FORMAT ERROR ENCOUNTERED IN READING FILE 'PRIMES'
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      PARAMETER (JX=9593)
      DIMENSION IPRM(JX),NFAC(*)
      DATA NCALL / 0 /
      M = IABS(N)
      IFLAG = 0
      L = 0
      IF (M.EQ.0) RETURN
      IF (M.GT.9999999999) THEN
         IFLAG = 1
      ELSEIF (M.EQ.1) THEN
         L = 1
         NFAC(1) = 1
      ELSE
         IF (NCALL.EQ.0) THEN 
C
C--- READ ALL PRIMES FROM 2 TO 100,003
C
            OPEN (UNIT=4,FILE='PRIMES') 
            REWIND 4
            READ (4,*,ERR=50,END=40) (IPRM(J),J=1,JX)
         ENDIF
C
C--- COUNT NUMBER OF CALLS TO THIS SUBROUTINE
C
         NCALL = NCALL+1
C
C--- IMPLEMENT ALGORITHM A (FACTORING BY PRIMES) ON PP. 364-5 OF THE
C--- REFERENCE
C
         K = 1
   10    IF (M.EQ.1) RETURN
   20    IQUO = M/IPRM(K)
         IREM = MOD(M,IPRM(K))
         IF (IREM.NE.0) GO TO 30
         L = L+1
         NFAC(L) = IPRM(K)
         M = IQUO
         GO TO 10
   30    IF (IQUO.GT.IPRM(K)) THEN
            K = K+1 
            GO TO 20
         ENDIF
         L = L+1
         NFAC(L) = M
      ENDIF
      RETURN
C
C--- SET ERROR FLAG FOR 'END-OF-FILE' ON UNIT 4
C
   40 IFLAG = 2
      RETURN
C
C--- SET ERROR FLAG FOR 'READ ERROR' ON UNIT 4
C
   50 IFLAG = 3
      RETURN
      END 
*FIBMIN
      SUBROUTINE FIBMIN (FUNC,XLO,XHI,TOL,XMIN,FXMIN)
C-----------------------------------------------------------------------
C   FIBMIN   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND 20899
C
C   FOR: COMPUTING THE X VALUE (XMIN) AT WHICH THE FUNCTION FUNC(X) IS
C        MINIMAL IN THE INTERVAL [XLO,XHI].  IN THAT INTERVAL FUNC IS 
C        ASSUMED TO BE UNIMODAL.  THE TOLERANCE ON XMIN (TOL) MUST BE 
C        SPECIFIED.  THE FUNCTION VALUE AT XMIN (FXMIN) IS RETURNED.
C        THE MINIMUM IS FOUND USING A FIBONACCI SEARCH ALGORITHM.
C
C   NOTE: IF TOL < [XHI-XLO]/[THE JX(TH) FIBONACCI NUMBER], THEN JX
C         SHOULD BE INCREASED.  THE 100(TH) FIBONACCI NUMBER IS ABOUT 
C         3.5E+20.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED AUGUST 12, 1988
C-----------------------------------------------------------------------
      PARAMETER (JX=100)
      DIMENSION IFIB(JX)
      DATA IFIB(1),IFIB(2) / 1,1 /
      XRANGE = XHI-XLO
      BOUND = XRANGE/TOL
      J = 2
   10 J = J+1
      IF (J.GT.JX) STOP ' J>JX IN SUBROUTINE FIBMIN - STOP' 
      IFIB(J) = IFIB(J-2)+IFIB(J-1)
      IF (REAL(IFIB(J)).LT.BOUND) GO TO 10
      DELTA = XRANGE/REAL(IFIB(J))
      C1 = XLO
      C2 = DELTA
      K = 0
      M = IFIB(J-2) 
      N = IFIB(J-1) 
      L = IFIB(J)
      F1 = FUNC(XLO)
      F2 = FUNC(XHI)
      I = J-2
   20 I = I-1
      IF (I.LT.1) GO TO 30
      F1 = FUNC(C1+C2*REAL(M))
      F2 = FUNC(C1+C2*REAL(N))
      IF (F1.EQ.F2) THEN
         K = M
         L = N
         I = I-2
         IF (I.LT.1) GO TO 30 
         M = K+IFIB(I)
         N = K+L-M
         F3 = FUNC(C1+C2*REAL(M))
         IF (F3.GE.F1) GO TO 30
         F4 = FUNC(C1+C2*REAL(N))
         IF (F4.GE.F1) GO TO 30
      ELSEIF (F1.GT.F2) THEN
         K = M
         M = N
         N = K+L-M
      ELSE
         L = N
         N = M
         M = K+L-N
      ENDIF
      GO TO 20
   30 FBMIN = 0.5*REAL(K+L)
      XMIN = C1+C2*FBMIN
      FXMIN = FUNC(XMIN)
      RETURN
C
      END 
*IGCD
      INTEGER FUNCTION IGCD (M,N)
C
C-----------------------------------------------------------------------
C   IGCD   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C          DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C          MARYLAND  20899
C
C   FOR: COMPUTING THE GREATEST COMMON DIVISOR OF THE INTEGERS M AND N.
C        IF M=N=0 THEN THE G.C.D. IS DEFINED TO BE ZERO.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED FEBRUARY 3, 1987
C-----------------------------------------------------------------------
C
      I = IABS(N)
      IGCD = IABS(M)
      IF (IGCD.NE.0) GO TO 10 
      IGCD = I
      RETURN
   10 K = MOD(I,IGCD)
      IF (K.EQ.0) RETURN
      I = IGCD
      IGCD = K
      GO TO 10
      END 
*IGCDM
      INTEGER FUNCTION IGCDM (M,LDM,NR,NC)
C
C-----------------------------------------------------------------------
C   IGCDM   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C           DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C           MARYLAND  20899
C
C   FOR: COMPUTING THE GREATEST COMMON DIVISOR OF THE NR*NC INTEGERS IN
C        THE MATRIX M.  IF M(1,1)=M(1,2)=...=M(NR,NC)=0 THEN THE G.C.D.
C        IS DEFINED TO BE ZERO (NR IS THE NUMBER OF ROWS OF M, NC IS
C        THE NUMBER OF COLUMNS OF M, AND LDM IS THE LEADING DIMENSION 
C        OF M AS DIMENSIONED IN THE CALLING PROGRAM).  NO CHECK IS MADE
C        FOR THE ERROR CONDITION NR>LDM.
C
C   SUBPROGRAMS CALLED: IGCD
C
C   CURRENT VERSION COMPLETED FEBRUARY 3, 1987
C-----------------------------------------------------------------------
C
      DIMENSION M(LDM,*)
      IGCDM = M(1,1)
      DO 20 J = 1, NC
         DO 10 I = 1, NR
            L = IABS(M(I,J))
            IGCDM = IGCD(IGCDM,L)
            IF (IGCDM.EQ.1) RETURN
   10    CONTINUE
   20 CONTINUE
      RETURN
      END 
*IGCDV
      INTEGER FUNCTION IGCDV (M,N)
C
C-----------------------------------------------------------------------
C   IGCDV   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C           DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C           MARYLAND  20899
C
C   FOR: COMPUTING THE GREATEST COMMON DIVISOR OF THE N INTEGERS IN THE
C        VECTOR M.  IF M(1)=M(2)=...=M(N)=0 THEN THE G.C.D. IS DEFINED
C        TO BE ZERO.
C
C   SUBPROGRAMS CALLED: IGCD
C
C   CURRENT VERSION COMPLETED FEBRUARY 3, 1987
C-----------------------------------------------------------------------
C
      DIMENSION M(*)
      IGCDV = M(1)
      DO 10 I = 1, N
         L = IABS(M(I))
         IGCDV = IGCD(IGCDV,L)
         IF (IGCDV.EQ.1) RETURN
   10 CONTINUE
      RETURN
      END 
*LINSYS
      SUBROUTINE LINSYS (A,LDA,N,B,IPVT,IFLAG)
C
C-----------------------------------------------------------------------
C   LINSYS   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: SOLVING THE NXN LINEAR SYSTEM OF EQUATIONS AX=B FOR THE VECTOR
C        X.  GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING IS USED TO
C        FACTOR A INTO LU.  THE TRIANGULAR SYSTEM LY=B IS THEN SOLVED 
C        FOR Y, AND SIMILARLY UX=Y IS SOLVED FOR X.  THIS ALGORITHM IS
C        DESCRIBED ON PAGES 1.11 AND 1.14 OF THE REFERENCE BELOW.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED JANUARY 6, 1987
C
C   REFERENCE: DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, 
C              G.W., "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C      * A = MATRIX (SIZE NXN) OF COEFFICIENTS ON INPUT AND LU
C            FACTORIZATION ON OUTPUT (REAL)
C
C    * LDA = THE LEADING DIMENSION OF MATRIX A (>=N) (INTEGER)
C
C      * N = NUMBER OF EQUATIONS AND UNKNOWNS (INTEGER)
C
C      * B = VECTOR (LENGTH N) CONTAINING THE RIGHT HAND SIDE OF THE
C            SYSTEM ON INPUT AND THE SOLUTION VECTOR X ON OUTPUT (REAL)
C
C     IPVT = VECTOR (LENGTH N) FOR PIVOTING INFORMATION (INTEGER)
C
C    IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C            0 -> NO ERRORS DETECTED
C            K -> THE PIVOT ELEMENT AT THE KTH STEP OF THE ELIMINATION
C                 (1<=K<=N) IS ZERO.  THE LU FACTORIZATION OF A IS
C                 VALID, BUT THE MATRIX U IS SINGULAR AND CANNOT BE USED
C                 FOR SOLVING THE LINEAR SYSTEM.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION A(LDA,*),B(*),IPVT(*)
      IFLAG = 0
C
C--- COMPUTE LU FACTORIZATION OF A BY GAUSSIAN ELIMINATION WITH
C--- PARTIAL PIVOTING.
C
      DO 50 K = 1, N-1
         AMAX = 0.0 
         DO 10 I = K, N
            IF (ABS(A(I,K)).GE.AMAX) THEN
               AMAX = ABS(A(I,K))
               L = I
            ENDIF
   10    CONTINUE
         IPVT(K) = L
         IF (A(L,K).EQ.0.0) THEN
            IFLAG = K
            GO TO 50
         ENDIF
         IF (L.NE.K) THEN
            Q = A(K,K)
            A(K,K) = A(L,K)
            A(L,K) = Q
         ENDIF
         DO 20 I = K+1, N
            A(I,K) = -A(I,K)/A(K,K)
   20    CONTINUE
         DO 40 J = K+1, N
            IF (L.NE.K) THEN
               Q = A(K,J)
               A(K,J) = A(L,J)
               A(L,J) = Q
            ENDIF
            DO 30 I = K+1, N
               A(I,J) = A(I,J)+A(I,K)*A(K,J)
   30       CONTINUE
   40    CONTINUE
   50 CONTINUE
      IF (A(N,N).EQ.0.0) IFLAG = N
      IF (IFLAG.NE.0) RETURN
C
C--- SOLVE LY=B FOR Y BY BACK SUBSTITUTION
C
      DO 70 K = 1, N-1
         L = IPVT(K)
         IF (L.NE.K) THEN
            Q = B(K)
            B(K) = B(L)
            B(L) = Q
         ENDIF
         DO 60 I = K+1, N
            B(I) = B(I)+A(I,K)*B(K)
   60    CONTINUE
   70 CONTINUE
C
C--- SOLVE UX=Y FOR X BY BACK SUBSTITUTION
C
      DO 90 K = N, 1, -1
         B(K) = B(K)/A(K,K)
         DO 80 I = 1, K-1
            B(I) = B(I)-A(I,K)*B(K)
   80    CONTINUE
   90 CONTINUE
      RETURN
      END 
*LISYPD
      SUBROUTINE LISYPD (A,LDA,N,B,IFLAG)
C
C-----------------------------------------------------------------------
C   LISYPD   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: SOLVING THE NXN LINEAR SYSTEM OF EQUATIONS AX=B FOR THE VECTOR
C        X WHERE A IS SYMMETRIC AND POSITIVE DEFINITE.  THE CHOLESKY
C        FACTORIZATION OF A YIELDS A LOWER TRIANGULAR MATRIX L SUCH
C        THAT A=LL'.  THE TRIANGULAR SYSTEM LY=B IS THEN SOLVED FOR
C        Y, AND SIMILARLY L'X=Y IS SOLVED FOR X.  THE SOLUTION VECTOR 
C        X IS RETURNED IN B.  THIS ALGORITHM IS DESCRIBED ON PAGES
C        3.10 AND 3.11 OF THE REFERENCE BELOW.  NOTE THAT ONLY THE
C        LOWER TRIANGLE OF A IS USED.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED NOVEMBER 29, 1988
C
C   REFERENCE: DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, 
C              G.W., "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C      * A = MATRIX (SIZE NXN) OF COEFFICIENTS ON INPUT AND LU
C            FACTORIZATION ON OUTPUT (REAL)
C
C    * LDA = THE LEADING DIMENSION OF MATRIX A (>=N) (INTEGER)
C
C      * N = NUMBER OF EQUATIONS AND UNKNOWNS (INTEGER)
C
C      * B = VECTOR (LENGTH N) CONTAINING THE RIGHT HAND SIDE OF THE
C            SYSTEM ON INPUT AND THE SOLUTION VECTOR X ON OUTPUT (REAL)
C
C    IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C            -3 -> N<1
C             0 -> NO ERRORS DETECTED
C             K -> THE KTH DIAGONAL ELEMENT OF THE CHOLESKY FACTOR-
C                  IZATION IS <0, THUS THE FACTORIZATION CANNOT BE
C                  COMPLETED.  INDICATES A IS NOT POSITIVE DEFINITE.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION A(LDA,*),B(*) 
      IFLAG = 0
      IF (N.LT.1) THEN
         IFLAG = -3 
         RETURN
C
      ENDIF
C
C--- COMPUTE THE CHOLESKY FACTORIZATION A=LL'
C
      DO 40 K = 1, N
         DO 20 I = 1, K-1
            S = A(K,I)
            DO 10 J = 1, I-1
               S = S-A(I,J)*A(K,J)
   10       CONTINUE
            A(K,I) = S/A(I,I) 
   20    CONTINUE
         S = A(K,K) 
         DO 30 J = 1, K-1
            S = S-A(K,J)**2
   30    CONTINUE
         IF (S.LE.0.0) THEN
            IFLAG = K
            RETURN
C
         ENDIF
         A(K,K) = SQRT(S)
   40 CONTINUE
C
C--- SOLVE LY=B (STORE Y IN B)
C
      DO 60 I = 1, N
         S = B(I)
         DO 50 J = 1, I-1
            S = S-A(I,J)*B(J) 
   50    CONTINUE
         B(I) = S/A(I,I)
   60 CONTINUE
C
C--- SOLVE L'X=Y (STORE X IN B)
C
      DO 80 I = N, 1, -1
         S = B(I)
         DO 70 J = I+1, N
            S = S-A(J,I)*B(J) 
   70    CONTINUE
         B(I) = S/A(I,I)
   80 CONTINUE
      RETURN
C
      END 
*LSQSVD
      SUBROUTINE LSQSVD (X,Y,B,WORK,S,E,V,N,M,LDX,MX,K,IFLAG)
C
C-----------------------------------------------------------------------
C   LSQSVD   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: SOLVING THE LINEAR SYSTEM Y = XB FOR B GIVEN Y AND X WHERE
C        Y IS OF LENGTH N AND X IS AN N BY M MATRIX (N>=M).  IF X IS
C        OF FULL RANK (M) THEN B IS THE (UNWEIGHTED) LEAST SQUARES
C        SOLUTION.  IF X IS NOT OF FULL RANK THEN B IS THE MINIMUM
C        NORM LEAST SQUARES SOLUTION.  A LINPACK ROUTINE IS USED TO
C        PERFORM A SINGULAR VALUE FACTORIZATION OF X.
C
C   NOTE: RNDERR IS A MACHINE DEPENDENT CONSTANT WHICH IS THE MACHINE 
C         ROUNDING ERROR (OR A LITTLE LARGER).  IT IS USED TO DETERMINE
C         WHEN A SINGULAR VALUE IS ZERO.  SEE THE DISCUSSION OF THIS
C         POINT ON PAGE 11.2 OF REFERENCE 1.
C
C   SUBPROGRAMS CALLED: SSVDC (LINPACK) 
C
C   CURRENT VERSION COMPLETED SEPTEMBER 28, 1987
C
C   REFERENCES: 
C
C   1) DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, G.W.,
C      "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979, CH. 11.
C
C   2) LAWSON, CHARLES L. AND HANSON, RICHARD J., "SOLVING LEAST
C      SQUARES PROBLEMS", PRENTICE-HALL, INC., CH. 7.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   * X(LDX,*) = MATRIX (SIZE N BY M).  THE SECOND DIMENSION OF X,
C                DEFINED IN THE CALLING PROGRAM, MUST BE >=M [REAL]
C
C       * Y(*) = VECTOR (LENGTH N) [REAL]
C
C         B(*) = VECTOR (LENGTH M) CONTAINING THE LEAST SQUARES
C                SOLUTION IF RANK(X)=M, OR THE MINIMUM NORM LEAST
C                SQUARES SOLUTION IF RANK(X)<M, ON OUTPUT [REAL]
C
C      WORK(*) = VECTOR (LENGTH N) USED AS WORKSPACE [REAL] 
C
C         S(*) = VECTOR (LENGTH M) OF SINGULAR VALUES IN DESCENDING
C                ORDER ON RETURN, PROVIDED INFO=0 ON RETURN [REAL]
C
C         E(*) = VECTOR (LENGTH M) OF ZEROS ON RETURN, PROVIDED INFO=0
C                ON RETURN [REAL]
C
C      V(MX,*) = MATRIX (SIZE M BY M) USED FOR INTERMEDIATE 
C                COMPUTATIONS [REAL]
C
C          * N = NUMBER OF ROWS IN MATRIX X [INTEGER]
C
C          * M = NUMBER OF COLUMNS IN MATRIX X (M<=N) [INTEGER]
C
C        * LDX = LEADING DIMENSION OF MATRIX X (LDX>=N) [INTEGER]
C
C         * MX = LEADING DIMENSION OF MATRIX V (MX>=M) [INTEGER]
C
C            K = NUMBER OF NONZERO SINGULAR VALUES ON RETURN [INTEGER]
C
C        IFLAG = ERROR INDICATOR ON OUTPUT [INTEGER]  INTERPRETATION: 
C                1 -> N>LDX OR M>MX.
C                2 -> N<M.
C                3 -> INFO<>0 RETURNED FROM SSVDC (SINGULAR VALUES WERE
C                     NOT COMPUTED CORRECTLY, THUS MOORE-PENROSE
C                     PSEUDO-INVERSE NOT COMPUTED).
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION X(LDX,*),Y(*),B(*),WORK(*),S(*),E(*),V(MX,*)
      DATA RNDERR / 1.0E-14 / 
      IFLAG = 0
      IF (N.GT.LDX.OR.M.GT.MX) THEN
         IFLAG = 1
         RETURN
C
      ENDIF
      IF (N.LT.M) THEN
         IFLAG = 2
         RETURN
C
      ENDIF
C
C--- COMPUTE LARGEST ELEMENT OF X (IN ABSOLUTE VALUE)
C
      XMAX = 0.0
      DO 20 I = 1, N
         DO 10 J = 1, M
            XMAX = AMAX1(XMAX,ABS(X(I,J)))
   10    CONTINUE
   20 CONTINUE
C
C--- COMPUTE CUTOFF POINT FOR A SINGULAR VALUE BEING ZERO
C
      CUTOFF = 10.0*RNDERR*XMAX
C
C--- PERFORM SINGULAR VALUE FACTORIZATION USING LINPACK
C
      CALL SSVDC (X,LDX,N,M,S,E,X,LDX,V,MX,WORK,21,INFO)
C
C--- CHECK WHETHER SINGULAR VALUES HAVE BEEN COMPUTED CORRECTLY
C
      IF (INFO.NE.0) THEN
         IFLAG = 3
         RETURN
C
      ENDIF
C
C--- DETERMINE NUMBER OF NONZERO SINGULAR VALUES
C
      K = 0
      DO 30 J = 1, M
         IF (ABS(S(J)).GT.CUTOFF) THEN
            K = K+1 
         ELSE
            GO TO 80
C
         ENDIF
   30 CONTINUE
C
C--- COMPUTE WORK = INV(S)U'Y 
C
   80 DO 50 I = 1, K
         T = 0.0
         DO 40 J = 1, N
            T = T+X(J,I)*Y(J) 
   40    CONTINUE
         WORK(I) = T/S(I)
   50 CONTINUE
C
C--- COMPUTE B = VW 
C
      DO 70 I = 1, M
         T = 0.0
         DO 60 J = 1, K
            T = T+V(I,J)*WORK(J)
   60    CONTINUE
         B(I) = T
   70 CONTINUE
      END 
*LSTSQR
      SUBROUTINE LSTSQR (X,LDX,N,M,Y,RES,SRES,HAT,UNK,SDUNK,LU,IND,ID)
C
C-----------------------------------------------------------------------
C   LSTSQR   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: ANALYZING THE STATISTICAL MODEL  Y = X*B + ERROR  BY UNWEIGHTED
C        LINEAR LEAST SQUARES.  THREE SUBROUTINES IN LINPACK ARE CALLED.
C        COMPUTED QUANTITIES ARE: 
C
C           1) ESTIMATES OF UNKNOWNS (B)
C           2) STANDARD ERRORS OF ESTIMATES OF UNKNOWNS
C           3) T-VALUES FOR UNKNOWNS
C           4) RESIDUAL STANDARD DEVIATION
C           5) PREDICTED VALUES OF OBSERVATIONS
C           6) RESIDUALS
C           7) STANDARDIZED RESIDUALS
C           8) DIAGONAL ELEMENTS OF 'HAT MATRIX'
C           9) DURBIN-WATSON STATISTIC FOR SEQUENTIAL CORRELATION
C              OF RESIDUALS
C
C        THESE QUANTITIES ARE AUTOMATICALLY PRINTED AT USER'S OPTION. 
C        ALL BUT 3 AND 5 ARE RETURNED IN THE PASSED ARGUMENTS.  AN
C        OPTION EXISTS TO SUPPRESS THE COMPUTATION OF 7 AND 8 SINCE
C        THE TIME REQUIREMENTS FOR THAT ARE ON THE ORDER OF N**2.  TO 
C        ILLUSTRATE THE TIME REQUIREMENTS, A FIFTH-DEGREE POLYNOMIAL
C        WAS FIT TO N POINTS WITH AND WITHOUT COMPUTING 7 AND 8.  THE 
C        CPU TIMES WERE AS FOLLOWS ON THE CYBER 180/855 AT N.I.S.T.:  
C
C                             WITHOUT 7,8  WITH 7,8
C                        N      (IND=0)    (IND=10)
C                      ----   -----------  --------
C                        20      0.005       0.017
C                       100      0.013       0.153
C                       500      0.050       3.045
C                      2500      0.277      83.763
C
C        THE RESIDUAL STANDARD DEVIATION IS RETURNED IN X(1,1) AND THE
C        DURBIN-WATSON STATISTIC IS RETURNED IN X(2,1).  THE VECTOR OF
C        OBSERVATIONS, Y, IS NOT CHANGED BY THIS SUBROUTINE.
C
C   NOTE: THE VECTORS WHICH APPEAR IN THE ARGUMENT LIST ARE NECESSARY 
C         FOR WORK SPACE IN THIS SUBROUTINE.  THE USER SHOULD NOT
C         ATTEMPT TO REMOVE ANY OF THEM.
C
C   NOTE: IN THE CALLING PROGRAM VECTORS Y, RES, SRES, AND HAT MUST
C         BE DIMENSIONED AT LEAST N.  VECTORS UNK AND SDUNK MUST BE
C         DIMENSIONED AT LEAST M.  A SAMPLE MAIN PROGRAM MIGHT BE: 
C
C         PARAMETER (NX=200,MX=10)
C         DIMENSION X(NX,MX)
C         DIMENSION Y(NX),RES(NX),SRES(NX),HAT(NX)
C         DIMENSION UNK(MX),SDUNK(MX)
C         CHARACTER*60 ID
C         ID = '<UP TO 60 LINES IDENTIFYING DATA SET>'
C         LU = 6
C         N = <IN THE RANGE [2,NX]>
C         M = <IN THE RANGE [1,MX]>
C               .
C               .
C               .
C
C         CALL LSTSQR(X,NX,N,M,Y,RES,SRES,HAT,UNK,SDUNK,LU,IND,ID)
C               .
C               .
C               .
C         END
C
C   SUBPROGRAMS CALLED: SQRDC, SQRSL, STRSL (LINPACK)
C
C   CURRENT VERSION COMPLETED APRIL 21, 1989
C
C   REFERENCE: DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, 
C              G.W., "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979,
C              CH. 9.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C    * X(LDX,*) = MATRIX (SIZE N BY M) IN THE MODEL Y = XB + ERROR.
C                 THE SECOND DIMENSION OF X, DEFINED IN THE CALLING
C                 PROGRAM, MUST BE >=M.  THIS MATRIX IS CHANGED
C                 CONSIDERABLY ON OUTPUT. [REAL]
C
C         * LDX = LEADING DIMENSION OF MATRIX X (LDX>=N) [INTEGER]
C
C           * N = NUMBER OF ROWS OF MATRIX X (N<=LDX) [INTEGER]
C
C           * M = NUMBER OF COLUMNS OF MATRIX X (M<=N) [INTEGER]
C
C        * Y(*) = VECTOR (LENGTH N) OF OBSERVATIONS.  NOT CHANGED ON
C                 OUTPUT. [REAL]
C
C        RES(*) = VECTOR (LENGTH N) USED FOR WORKSPACE.  ON OUTPUT
C                 CONTAINS RESIDUALS [REAL]
C
C       SRES(*) = VECTOR (LENGTH N) USED FOR WORKSPACE.  ON OUTPUT
C                 CONTAINS STANDARDIZED RESIDUALS [REAL]
C
C        HAT(*) = VECTOR (LENGTH N) USED FOR WORKSPACE.  ON OUTPUT
C                 CONTAINS DIAGONAL ELEMENTS OF THE HAT MATRIX [REAL] 
C
C        UNK(*) = VECTOR (LENGTH M) CONTAINING ESTIMATES OF UNKNOWNS
C                 ON OUTPUT [REAL]
C
C      SDUNK(*) = VECTOR (LENGTH M) USED FOR WORKSPACE.  ON OUTPUT
C                 CONTAINS STANDARD ERRORS OF ESTIMATES OF UNKNOWNS
C                 [REAL]
C
C          * LU = LOGICAL UNIT OF PRINT FILE WHERE AUTOMATIC PRINTOUT 
C                 OCCURS, IF REQUESTED [INTEGER]
C
C         * IND = INTEGER INDICATOR VARIABLE OF THE FORM AB WHERE
C                   A=0 TO SUPPRESS COMPUTING OF HAT MATRIX DIAGONALS 
C                   A=1 TO COMPUTE HAT MATRIX DIAGONALS
C                   B=0 FOR NO PRINTOUT TO LOGICAL UNIT LU
C                   B=1 FOR MINIMAL PRINTOUT TO LOGICAL UNIT LU (NO
C                       RESIDUALS, STANDARDIZED RESIDUALS, OR HAT
C                       MATRIX DIAGONAL ELEMENTS PRINTED)
C                   B=2 FOR COMPLETE PRINTOUT TO LOGICAL UNIT LU
C
C          * ID = CHARACTER VARIABLE OF LENGTH 60 FOR PASSING THE
C                 IDENTIFICATION OF THE DATA SET BEING FITTED.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION X(LDX,*),Y(*),RES(*),SRES(*),HAT(*),UNK(*),SDUNK(*)
      CHARACTER*1 STAR
      CHARACTER*60 ID
C
C--- DETERMINE INDICATORS FOR PRINTOUT AND COMPUTATION OF HAT
C--- MATRIX DIAGONAL
C
      IHAT = IND/10 
      IPRT = IND-10*IHAT
C
C--- CHECK FOR ILLEGAL INPUT ARGUMENTS
C
      IF (M.LT.1) THEN
         IFLAG = 1
      ELSEIF (N.LE.M) THEN
         IFLAG = 2
      ELSEIF (LDX.LT.N) THEN
         IFLAG = 3
      ELSEIF (IHAT.LT.0.OR.IHAT.GT.1.OR.IPRT.LT.0.OR.IPRT.GT.2) THEN
         IFLAG = 4
      ELSE
         IF (IPRT.GE.1) WRITE (LU,1000) ID
C
C--- COMPUTE QR FACTORIZATION OF X
C
         CALL SQRDC (X,LDX,N,M,SDUNK,JDUM,DUM,0)
C
C--- COMPUTE DIAGONAL OF HAT MATRIX, X*INV(X'X)*X'.  IF COMPUTATION OF
C--- THIS QUANTITY IS SUPPRESSED, SET EACH ELEMENT TO M/N.
C
         RMN = REAL(M)/REAL(N)
         HATLIM = 2.0*RMN
         IF (IHAT.EQ.1) THEN
            DO 10 I = 1, N
               SRES(I) = 0.0
   10       CONTINUE
            DO 30 I = 1, N
               SRES(I) = 1.0
               IF (I.GT.1) SRES(I-1) = 0.0
               CALL SQRSL (X,LDX,N,M,SDUNK,SRES,DUM,RES,DUM,DUM,DUM,
     *            01000,INFO) 
               S = 0.0
               DO 20 J = 1, M 
                  S = S+RES(J)**2
   20          CONTINUE
               HAT(I) = S
   30       CONTINUE
         ELSE
            DO 40 I = 1, N
               HAT(I) = RMN
   40       CONTINUE
         ENDIF
C
C--- COMPUTE ESTIMATES OF UNKNOWNS AND RESIDUALS
C
         CALL SQRSL (X,LDX,N,M,SDUNK,Y,DUM,SRES,UNK,RES,DUM,00110,INFO)
C
C--- COMPUTE NUMERATOR OF DURBIN-WATSON STATISTIC 
C
         DW = 0.0
         DO 50 I = 2, N
            DW = DW+(RES(I)-RES(I-1))**2
   50    CONTINUE
C
C--- COMPUTE RESIDUAL SUM-OF-SQUARES
C
         S = 0.0
         DO 60 I = 1, N
            S = S+RES(I)**2
   60    CONTINUE
C
C--- COMPUTE DURBIN-WATSON STATISTIC, DEGREES OF FREEDOM, AND
C--- RESIDUAL STANDARD DEVIATION
C
         DW = DW/S
         IDF = N-M
         RSD = SQRT(S/REAL(IDF))
         IF (IPRT.GE.1) WRITE (LU,1030) N,M,IDF,RSD
         IF (IPRT.GE.1) WRITE (LU,1010) 
C
C--- COMPUTE STANDARD ERRORS AND T-VALUES FOR THE UNKNOWNS
C
         DO 90 J = 1, M
            DO 70 K = 1, M
               SRES(K) = 0.0
   70       CONTINUE
            SRES(J) = 1.0
            CALL STRSL (X,LDX,M,SRES,11,INFO)
            S = 0.0 
            DO 80 K = 1, M
               S = S+SRES(K)**2
   80       CONTINUE
            SDUNK(J) = RSD*SQRT(S)
            T = UNK(J)/SDUNK(J)
            IF (IPRT.GE.1) WRITE (LU,1050) J,UNK(J),SDUNK(J),T
   90    CONTINUE
         IF (IPRT.GE.1) WRITE (6,1080) DW
         IF (IPRT.GE.2) WRITE (6,1020)
C
C--- COMPUTE STANDARDIZED RESIDUALS AND INDICATE LARGE HAT MATRIX
C--- DIAGONAL ELEMENTS, IF ANY
C
         DO 100 I = 1, N
            IF (HAT(I).GT.HATLIM) THEN
               STAR = '*'
            ELSE
               STAR = ' '
            ENDIF
            PRED = Y(I)-RES(I)
            IF (HAT(I).GE.1.0) THEN
               SRES(I) = 0.0
            ELSE
               SRES(I) = RES(I)/(RSD*SQRT(1.0-HAT(I)))
            ENDIF
            IF (IPRT.GE.2) WRITE (LU,1040) I,Y(I),PRED,RES(I),SRES(I),
     *         HAT(I),STAR
  100    CONTINUE
C
C--- WRITE MESSAGE ABOUT HAT MATIX DIAGONAL ELEMENTS
C
         IF (IPRT.GE.2) THEN
            IF (IHAT.EQ.1) THEN
               WRITE (LU,1060) HATLIM
            ELSE
               WRITE (LU,1070)
            ENDIF
         ENDIF
C
C--- RETURN RESIDUAL STANDARD DEVIATION AND DURBIN-WATSON STATISTIC
C--- IN MATRIX X (WHICH HAS ALREADY BEEN CHANGED) 
C
         X(1,1) = RSD
         X(2,1) = DW
      ENDIF
      RETURN
C
C
 1000 FORMAT ('1',18X,'UNWEIGHTED LINEAR LEAST SQUARES (LSTSQR)'///7X,
     *   'ID: ',A60)
 1010 FORMAT (//15X,'UNKNOWN',8X,'STD. DEV.',8X,'T-VALUE'/) 
 1020 FORMAT (//50X,'STANDARDIZED',2X,'HAT MATRIX'/11X,'OBSERVED',6X, 
     *   'PREDICTED',5X,'RESIDUAL',5X,'RESIDUAL',5X,'DIAGONAL'/)
 1030 FORMAT (//10X,'NUMBER OF OBSERVATIONS = ',I5/14X,'NUMBER OF ',
     *   'UNKNOWNS = ',I5/8X,'RESIDUAL DEG. OF FREEDOM = ',I5,3X,
     *   'RESIDUAL STD. DEV. = ',G12.5) 
 1040 FORMAT (1X,I5,2X,2G15.8,G12.5,F9.3,3X,F9.4,1X,A1)
 1050 FORMAT (5X,I5,1X,G16.8,1X,G13.5,2X,G14.3)
 1060 FORMAT (//10X,'DIAGONAL ELEMENTS OF HAT MATRIX WHICH EXCEED ',
     *   'RULE-OF-THUMB'/10X,'LIMIT 2*(# UNKNOWNS)/(# EQUATIONS) = ', 
     *   F5.3,' ARE INDICATED BY *'//)
 1070 FORMAT (//10X,'DIAGONAL ELEMENTS OF HAT MATRIX WERE NOT COMPU', 
     *   'TED.  THEY WERE'/10X,'SET TO THEIR AVERAGE VALUE, (# UNKNOWN',
     *   'S)/(# EQUATIONS).  AS A'/10X,'RESULT, THE STANDARDIZED RESI',
     *   'DUALS ARE APPROXIMATE.'//)
 1080 FORMAT (//10X,'DURBIN-WATSON STATISTIC = ',F6.3)
C
      END 
*MATCNO
      SUBROUTINE MATCNO (A,LDA,N,CONDNO,IFLAG)
C
C-----------------------------------------------------------------------
C   MATCNO   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: COMPUTING THE CONDITION NUMBER OF THE N BY N MATRIX A.  THE
C        CONDITION NUMBER IS KAPPA = NORM(A)*NORM[INV(A)].  HERE THE
C        INFINITY-NORM IS USED.  THE CONDITION NUMBER MEASURES HOW
C        "ILL-CONDITIONED" THE MATRIX A IS, AND THUS HOW ACCURATELY THE
C        SYSTEM A*X = B MAY BE SOLVED.  IF A IS SINGULAR THEN KAPPA IS
C        INFINITE AND AN ERROR FLAG IS RETURNED.  SEE THE REFERENCE FOR
C        A GOOD DISCUSSION OF THE PROBLEM.
C
C   SUBPROGRAMS CALLED: MATINV (STSPAC) 
C                       ANORMI (ATTACHED)
C
C   CURRENT VERSION COMPLETED JUNE 28, 1989
C
C   REFERENCE: GOLUB, G.H. AND VAN LOAN, C.F., "MATRIX COMPUTATIONS", 
C              JOHNS HOPKINS UNIVERSITY PRESS, BALTIMORE, 1985, P. 72.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   * A(LDA,*) = MATRIX (SIZE N BY N) WHOSE CONDITION NUMBER IS TO BE 
C                DETERMINED.  THE SECOND DIMENSION OF A, DEFINED IN THE
C                CALLING PROGRAM, MUST BE >=N [REAL]
C
C        * LDA = THE LEADING DIMENSION OF A (>=N) [INTEGER] 
C
C          * N = SIZE OF MATRIX A (N<=LDA) [INTEGER]
C
C       CONDNO = THE CONDITION NUMBER OF A (KAPPA) ON RETURN WHEN
C                IFLAG=0 [REAL]
C
C        IFLAG = ERROR INDICATOR ON OUTPUT.  THE OUTPUT VALUE OF IFLAG
C                FROM SUBROUTINE MATINV IS PASSED DIRECTLY THROUGH THIS
C                SUBROUTINE.  IF IFLAG=0, THEN ALL IS WELL.  OTHERWISE,
C                SEE DOCUMENTATION OF MATINV FOR INTERPRETATION OF THE
C                THE ERROR FLAG IFLAG [INTEGER]
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION A(LDA,*)
      FAC1 = ANORMI(A,LDA,N)
      CALL MATINV (A,LDA,N,IFLAG)
      IF (IFLAG.NE.0) THEN
         CONDNO = 0.999999E29 
      ELSE
         FAC2 = ANORMI(A,LDA,N)
         CONDNO = FAC1*FAC2
      ENDIF
      RETURN
C
      END 
C
C=======================================================================
C
      FUNCTION ANORMI (A,LDA,N)
C
C--- COMPUTE THE INFINITY-NORM OF THE N BY N MATRIX A
C
      DIMENSION A(LDA,*)
      ANORMI = 0.0
      DO 20 I = 1, N
         T = 0.0
         DO 10 J = 1, N
            T = T+ABS(A(I,J)) 
   10    CONTINUE
         ANORMI = AMAX1(ANORMI,T)
   20 CONTINUE
      RETURN
C
      END 
*MATHDI
      SUBROUTINE MATHDI (X,LDX,N,M,Y,QTY,HAT,QRAUX)
C
C-----------------------------------------------------------------------
C   MATHDI   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: COMPUTING THE DIAGONAL OF THE 'HAT MATRIX' CORRESPONDING TO
C        THE 'DESIGN MATRIX' X IN THE (UNWEIGHTED) LINEAR LEAST SQUARES
C        MODEL  Y = XB + ERROR.  THE HAT MATRIX IS DEFINED TO BE (SEE 
C        REFERENCES 2 AND 3 BELOW)  H = X*INV(X'X)*X'.  WHEN X HAS N
C        ROWS AND M COLUMNS (M<=N) THEN H IS N BY N.  THE DIAGONAL
C        ENTRIES OF H MEASURE THE 'LEVERAGE' OR 'INFLUENCE' OF THE
C        CORRESPONDING OBSERVATIONS.  WHEN X IS OF FULL RANK, THE SUM OF
C        THE DIAGONAL ENTRIES OF H IS M, THEREFORE, THE AVERAGE VALUE 
C        OF THE ENTRIES IS M/N.  A RULE-OF-THUMB GIVEN IN THE REFERENCES
C        IS THAT, FOR A GOOD DESIGN MATRIX, NO DIAGONAL ELEMENT OF H
C        EXCEEDS 2*M/N.
C
C        THE COMPUTATIONAL PROCEDURE BELOW INCORPORATES THE STANDARD
C        QR FACTORIZATION OF X AS IMPLEMENTED IN LINPACK.  WHEN X IS
C        NOT OF FULL RANK, THE RANK DEFICIENCY IS EXPRESSED IN THE
C        FACTOR R.  SINCE ONLY THE FACTOR Q IS USED IN COMPUTING THE
C        HAT MATRIX, THE PROGRAM WILL RETURN RESULTS WHEN X IS NOT OF 
C        FULL RANK.  THOSE RESULTS WILL BE ERRONEOUS, HOWEVER.  FOR THE
C        SAKE OF SIMPLICITY, NO CHECK IS MADE ON THE RANK OF X.  IT IS
C        THEREFORE UP TO THE USER TO KNOW THAT X IS OF FULL RANK BEFORE
C        ACCEPTING THE RESULTS OF THIS ROUTINE.
C
C   SUBPROGRAMS CALLED: SQRDC, SQRSL (LINPACK)
C
C   CURRENT VERSION COMPLETED JANUARY 25, 1989
C
C   REFERENCES: 
C
C   1) DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, G.W.,
C      "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979, CH. 9.
C
C   2) HOAGLIN, D.C. AND WELSCH, R.E., "THE HAT MATRIX IN REGRESSION
C      AND ANOVA", THE AMERICAN STATISTICIAN, VOL. 32, NO. 1, FEBRUARY
C      1978, PP. 17-22.
C
C   3) VELLEMAN, P.F. AND WELSCH, R.E., "EFFICIENT COMPUTING OF
C      REGRESSION DIAGNOSTICS", THE AMERICAN STATISTICIAN, VOL. 35,
C      NO. 4, NOVEMBER 1981, PP. 234-242.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C    * X(LDX,*) = MATRIX (SIZE N BY M) IN THE MODEL Y = XB + ERROR.
C                 THE SECOND DIMENSION OF X, DEFINED IN THE CALLING
C                 PROGRAM, MUST BE >=M [REAL]
C
C         * LDX = LEADING DIMENSION OF MATRIX X (LDX>=N) [INTEGER]
C
C           * N = NUMBER OF ROWS OF MATRIX X (N<=LDX) [INTEGER]
C
C           * M = NUMBER OF COLUMNS OF MATRIX X (M<=N) [INTEGER]
C
C          Y(*) = VECTOR (LENGTH N) USED FOR WORKSPACE [REAL]
C
C        QTY(*) = VECTOR (LENGTH N) USED FOR WORKSPACE [REAL]
C
C        HAT(*) = VECTOR (LENGTH N) CONTAINING DIAGONAL ENTRIES OF
C                 THE 'HAT MATRIX' ON OUTPUT [REAL]
C
C      QRAUX(*) = VECTOR (LENGTH M) USED FOR WORKSPACE [REAL]
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION X(LDX,*),Y(*),QTY(*),HAT(*),QRAUX(*)
C
C--- INITIALIZE UNIT DIRECTIONAL VECTOR 
C
      DO 10 I = 1, N
         Y(I) = 0.0 
   10 CONTINUE
C
C--- COMPUTE QR FACTORIZATION OF X (USING LINPACK)
C
      CALL SQRDC (X,LDX,N,M,QRAUX,JDUM,DUM,0)
C
C--- COMPUTE HAT MATRIX DIAGONAL FOR DESIGN MATRIX X
C
      DO 30 I = 1, N
         Y(I) = 1.0 
         IF (I.GT.1) Y(I-1) = 0.0
         CALL SQRSL (X,LDX,N,M,QRAUX,Y,DUM,QTY,DUM,DUM,DUM,01000,INFO)
         S = 0.0
         DO 20 J = 1, M
            S = S+QTY(J)**2
   20    CONTINUE
         HAT(I) = S 
   30 CONTINUE
      RETURN
C
      END 
*MATINV
      SUBROUTINE MATINV (A,LDA,N,IFLAG) 
C
C-----------------------------------------------------------------------
C   MATINV   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING THE INVERSE OF A GENERAL N BY N MATRIX IN PLACE,
C        I.E., THE INVERSE OVERWRITES THE ORIGINAL MATRIX.  THE STEPS 
C        OF THE ALGORITHM ARE DESCRIBED BELOW AS THEY OCCUR.  ROW
C        INTERCHANGES ARE DONE AS NEEDED IN ORDER TO INCREASE THE
C        ACCURACY OF THE INVERSE MATRIX.  WITHOUT INTERCHANGES THIS
C        ALGORITHM WILL FAIL WHEN ANY OF THE LEADING PRINCIPAL
C        SUBMATRICES ARE SINGULAR OR WHEN THE MATRIX ITSELF IS
C        SINGULAR.  WITH INTERCHANGES THIS ALGORITHM WILL FAIL ONLY
C        WHEN THE MATRIX ITSELF IS SINGULAR.  THE LEADING PRINCIPAL
C
C                                   [A B C]
C        SUBMATRICES OF THE MATRIX  [D E F]  ARE  [A]  AND  [A B] .
C                                   [G H I]                 [D E]
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED JANUARY 15, 1987
C
C   REFERENCE: STEWART, G.W., 'INTRODUCTION TO MATRIX COMPUTATIONS',
C              ACADEMIC PRESS, INC., 1973
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS
C
C     * A = MATRIX (SIZE NXN) TO BE INVERTED (REAL)
C
C   * LDA = LEADING DIMENSION OF MATRIX A [LDA>=N] (INTEGER)
C
C     * N = NUMBER OF ROWS AND COLUMNS OF MATRIX A (INTEGER)
C
C   IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C           -2 -> TOO MANY ROW INTERCHANGES NEEDED - INCREASE MX
C           -1 -> N>LDA
C            0 -> NO ERRORS DETECTED
C            K -> MATRIX A FOUND TO BE SINGULAR AT THE KTH STEP OF
C                 THE CROUT REDUCTION (1<=K<=N)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      PARAMETER (MX=100)
      DIMENSION A(LDA,*),IEX(MX,2)
      IFLAG = 0
C
C--- CHECK CONSISTENCY OF PASSED PARAMETERS
C
      IF (N.GT.LDA) THEN
         IFLAG = -1 
         RETURN
      ENDIF
C
C--- COMPUTE A = LU BY THE CROUT REDUCTION WHERE L IS LOWER TRIANGULAR
C--- AND U IS UNIT UPPER TRIANGULAR (ALGORITHM 3.4, P. 138 OF THE
C--- REFERENCE)
C
      NEX = 0
      DO 70 K = 1, N
         DO 20 I = K, N
            S = A(I,K)
            DO 10 L = 1, K-1
               S = S-A(I,L)*A(L,K)
   10       CONTINUE
            A(I,K) = S
   20    CONTINUE
C
C--- INTERCHANGE ROWS IF NECESSARY
C
         Q = 0.0
         L = 0
         DO 30 I = K, N
            R = ABS(A(I,K))
            IF (R.GT.Q) THEN
               Q = R
               L = I
            ENDIF
   30    CONTINUE
         IF (L.EQ.0) THEN
            IFLAG = K
            RETURN
         ENDIF
         IF (L.NE.K) THEN
            NEX = NEX+1
            IF (NEX.GT.MX) THEN
               IFLAG = -2
               RETURN
            ENDIF
            IEX(NEX,1) = K
            IEX(NEX,2) = L
            DO 40 J = 1, N
               Q = A(K,J)
               A(K,J) = A(L,J)
               A(L,J) = Q
   40       CONTINUE
         ENDIF
C
C--- END ROW INTERCHANGE SECTION
C
         DO 60 J = K+1, N
            S = A(K,J)
            DO 50 L = 1, K-1
               S = S-A(K,L)*A(L,J)
   50       CONTINUE
            A(K,J) = S/A(K,K) 
   60    CONTINUE
   70 CONTINUE
C
C--- INVERT THE LOWER TRIANGLE L IN PLACE (SIMILAR TO ALGORITHM 1.5,
C--- P. 110 OF THE REFERENCE) 
C
      DO 100 K = N, 1, -1
         A(K,K) = 1.0/A(K,K)
         DO 90 I = K-1, 1, -1 
            S = 0.0 
            DO 80 J = I+1, K
               S = S+A(J,I)*A(K,J)
   80       CONTINUE
            A(K,I) = -S/A(I,I)
   90    CONTINUE
  100 CONTINUE
C
C--- INVERT THE UPPER TRIANGLE U IN PLACE (ALGORITHM 1.5, P. 110 OF
C--- THE REFERENCE) 
C
      DO 130 K = N, 1, -1
         DO 120 I = K-1, 1, -1
            S = A(I,K)
            DO 110 J = I+1, K-1
               S = S+A(I,J)*A(J,K)
  110       CONTINUE
            A(I,K) = -S
  120    CONTINUE
  130 CONTINUE
C
C--- COMPUTE INV(A) = INV(U)*INV(L)
C
      DO 160 I = 1, N
         DO 150 J = 1, N
            IF (J.GT.I) THEN
               S = 0.0
               L = J
            ELSE
               S = A(I,J)
               L = I+1
            ENDIF
            DO 140 K = L, N
               S = S+A(I,K)*A(K,J)
  140       CONTINUE
            A(I,J) = S
  150    CONTINUE
  160 CONTINUE
C
C--- INTERCHANGE COLUMNS OF INV(A) TO REVERSE EFFECT OF ROW 
C--- INTERCHANGES OF A
C
      DO 180 I = NEX, 1, -1
         K = IEX(I,1)
         L = IEX(I,2)
         DO 170 J = 1, N
            Q = A(J,K)
            A(J,K) = A(J,L)
            A(J,L) = Q
  170    CONTINUE
  180 CONTINUE
      RETURN
      END 
*MATIPD
      SUBROUTINE MATIPD (A,LDA,N,IFLAG) 
C
C-----------------------------------------------------------------------
C   MATIPD   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING THE INVERSE OF A SYMMETRIC POSITIVE DEFINITE N BY N
C        MATRIX IN PLACE, I.E., THE INVERSE OVERWRITES THE ORIGINAL
C        MATRIX.  THE STEPS OF THE ALGORITHM ARE DESCRIBED BELOW AS
C        THEY OCCUR.  ONLY THE LOWER TRIANGLE (INCLUDING THE DIAGONAL)
C        OF THE INPUT MATRIX IS USED IN THE INVERSION.  THE COMPLETE
C        INVERSE MATRIX IS RETURNED ON OUTPUT.  IF THE INPUT MATRIX IS
C        NOT POSITIVE DEFINITE THE INVERSE CANNOT BE COMPUTED BY THE
C        CHOLESKY FACTORIZATION USED HEREIN, THUS AN ERROR FLAG IS SET.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED JANUARY 15, 1987
C
C   REFERENCE: STEWART, G.W., 'INTRODUCTION TO MATRIX COMPUTATIONS',
C              ACADEMIC PRESS, INC., 1973
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS
C
C     * A = SYMMETRIC POSITIVE DEFINITE MATRIX (SIZE NXN) TO BE INVERTED
C           (REAL)
C
C   * LDA = LEADING DIMENSION OF MATRIX A [LDA>=N] (INTEGER)
C
C     * N = NUMBER OF ROWS AND COLUMNS OF MATRIX A (INTEGER)
C
C   IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C           -1 -> N>LDA
C            0 -> NO ERRORS DETECTED
C            K -> MATRIX A FOUND NOT TO BE POSITIVE DEFINITE AT THE
C                 KTH STEP OF THE CHOLESKY FACTORIZATION (1<=K<=N)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION A(LDA,*)
      IFLAG = 0
C
C--- CHECK CONSISTENCY OF PASSED PARAMETERS
C
      IF (N.GT.LDA) THEN
         IFLAG = -1 
         RETURN
      ENDIF
C
C--- COMPUTE THE LOWER TRIANGULAR MATRIX L OF THE CHOLESKY
C--- FACTORIZATION A=LL' (ALGORITHM 3.9, P. 142 OF THE REFERENCE)
C
      DO 40 K = 1, N
         DO 20 I = 1, K-1
            S = A(K,I)
            DO 10 J = 1, I-1
               S = S-A(I,J)*A(K,J)
   10       CONTINUE
            A(K,I) = S/A(I,I) 
   20    CONTINUE
         S = A(K,K) 
         DO 30 J = 1, K-1
            S = S-A(K,J)**2
   30    CONTINUE
         IF (S.LE.0.0) THEN
            IFLAG = K
            RETURN
         ENDIF
         A(K,K) = SQRT(S)
   40 CONTINUE
C
C--- INVERT THE LOWER TRIANGLE L IN PLACE (SIMILAR TO ALGORITHM 1.5,
C--- P. 110 OF THE REFERENCE) 
C
      DO 70 K = N, 1, -1
         A(K,K) = 1.0/A(K,K)
         DO 60 I = K-1, 1, -1 
            S = 0.0 
            DO 50 J = I+1, K
               S = S+A(J,I)*A(K,J)
   50       CONTINUE
            A(K,I) = -S/A(I,I)
   60    CONTINUE
   70 CONTINUE
C
C--- COMPUTE INV(A)=INV(L)'*INV(L)
C
      DO 100 I = 1, N
         DO 90 J = I, N
            S = 0.0 
            DO 80 K = J, N
               S = S+A(K,I)*A(K,J)
   80       CONTINUE
            A(I,J) = S
            A(J,I) = S
   90    CONTINUE
  100 CONTINUE
      RETURN
      END 
*MATMPI
      SUBROUTINE MATMPI (X,WORK,S,E,V,N,M,NX,MX,K,IFLAG)
C
C-----------------------------------------------------------------------
C   MATMPI   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING THE MOORE-PENROSE PSEUDO-INVERSE OF AN NXM MATRIX
C        X WHERE N >= M.  THE TRANSPOSE OF THE INVERSE IS RETURNED IN 
C        THE ORIGINAL MATRIX.  A LINPACK ROUTINE IS USED TO PERFORM
C        A SINGULAR VALUE DECOMPOSITION OF X FROM WHICH THE INVERSE
C        X+ IS COMPUTED.  IF X IS OF FULL RANK THEN X+ = INV(X'X)*X'. 
C
C   NOTE: RNDERR IS A MACHINE DEPENDENT CONSTANT WHICH IS THE MACHINE 
C         ROUNDING ERROR (OR A LITTLE LARGER).  IT IS USED TO DETERMINE
C         WHEN A SINGULAR VALUE IS ZERO (IN CASES WHERE THE USER HAS
C         REQUESTED AUTOMATIC DETERMINATION OF RANK BY INPUTTING K=0).
C         SEE DISCUSSION OF THIS POINT ON PAGE 11.2 OF REFERENCE 1.
C
C   SUBPROGRAMS CALLED: SSVDC (LINPACK) 
C
C   CURRENT VERSION COMPLETED DECEMBER 14, 1989
C
C   REFERENCES: 
C
C   1) DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, G.W.,
C      "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979, CH. 11.
C
C   2) LAWSON, CHARLES L. AND HANSON, RICHARD J., "SOLVING LEAST
C      SQUARES PROBLEMS", PRENTICE-HALL, INC., CH. 7.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C    * X(NX,*) = MATRIX (SIZE N BY M) WHOSE PSEUDO-INVERSE IS TO BE
C                COMPUTED.  THE SECOND DIMENSION OF X, DEFINED IN THE 
C                CALLING PROGRAM, MUST BE >=M. [REAL]
C
C      WORK(*) = VECTOR (LENGTH N) USED AS WORKSPACE [REAL] 
C
C         S(*) = VECTOR (LENGTH M) OF SINGULAR VALUES IN DESCENDING
C                ORDER ON RETURN, PROVIDED INFO=0 ON RETURN [REAL]
C
C         E(*) = VECTOR (LENGTH M) OF ZEROS ON RETURN, PROVIDED INFO=0
C                ON RETURN [REAL]
C
C      V(MX,*) = MATRIX (SIZE M BY M) USED FOR INTERMEDIATE 
C                COMPUTATIONS [REAL]
C
C          * N = NUMBER OF ROWS IN MATRIX X [INTEGER]
C
C          * M = NUMBER OF COLUMNS IN MATRIX X (M<=N) [INTEGER]
C
C         * NX = LEADING DIMENSION OF MATRIX X (NX>=N) [INTEGER]
C
C         * MX = LEADING DIMENSION OF MATRIX V (MX>=M) [INTEGER]
C
C          * K = ON INPUT: K>0 INDICATES KNOWN RANK OF MATRIX X.
C                          K=0 INDICATES RANK SHOULD BE AUTOMATICALLY 
C                              DETERMINED BY PROGRAM.
C                ON OUTPUT: = UNCHANGED IF K>0 ON INPUT.
C                           = COMPUTED RANK OF X IF K=0 ON INPUT.
C
C        IFLAG = ERROR INDICATOR ON OUTPUT [INTEGER]  INTERPRETATION: 
C                1 -> N>NX OR M>MX.
C                2 -> N<M.
C                3 -> INFO<>0 RETURNED FROM SSVDC (SINGULAR VALUES WERE
C                     NOT COMPUTED CORRECTLY, THUS MOORE-PENROSE
C                     PSEUDO-INVERSE NOT COMPUTED).
C                4 -> K<0 ON INPUT.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION X(NX,*),WORK(*),S(*),E(*),V(MX,*) 
      DATA RNDERR / 1.0E-14 / 
      IFLAG = 0
      IF (N.GT.NX.OR.M.GT.MX) THEN
         IFLAG = 1
         RETURN
C
      ENDIF
      IF (N.LT.M) THEN
         IFLAG = 2
         RETURN
C
      ENDIF
      IF (K.LT.0) THEN
         IFLAG = 4
         RETURN
C
      ENDIF
      IF (K.EQ.0) THEN
C
C--- COMPUTE LARGEST ELEMENT OF X (IN ABSOLUTE VALUE)
C
         XMAX = 0.0 
         DO 20 I = 1, N
            DO 10 J = 1, M
               XMAX = AMAX1(XMAX,ABS(X(I,J)))
   10       CONTINUE
   20    CONTINUE
C
C--- COMPUTE CUTOFF POINT FOR A SINGULAR VALUE BEING ZERO
C
         CUTOFF = 10.0*RNDERR*XMAX
      ENDIF
C
C--- PERFORM SINGULAR VALUE FACTORIZATION USING LINPACK
C
      CALL SSVDC (X,NX,N,M,S,E,X,NX,V,MX,WORK,21,INFO)
C
C--- CHECK WHETHER SINGULAR VALUES HAVE BEEN COMPUTED CORRECTLY
C
      IF (INFO.NE.0) THEN
         IFLAG = 3
         RETURN
C
      ENDIF
      IF (K.EQ.0) THEN
C
C--- DETERMINE NUMBER OF NONZERO SINGULAR VALUES
C
         K = 0
         DO 30 J = 1, M
            IF (ABS(S(J)).GT.CUTOFF) THEN
               K = K+1
            ELSE
               GO TO 40
C
            ENDIF
   30    CONTINUE
      ENDIF
C
C--- COMPUTE THE MOORE-PENROSE PSEUDO-INVERSE OF X (TRANSPOSED)
C
   40 DO 60 J = 1, M
         DO 50 L = 1, K
            V(J,L) = V(J,L)/S(L)
   50    CONTINUE
   60 CONTINUE
      DO 100 I = 1, N
         DO 80 J = 1, M
            T = 0.0 
            DO 70 L = 1, K
               T = T+V(J,L)*X(I,L)
   70       CONTINUE
            E(J) = T
   80    CONTINUE
         DO 90 J = 1, M
            X(I,J) = E(J)
   90    CONTINUE
  100 CONTINUE
      RETURN
C
      END 
*MATXXI
      SUBROUTINE MATXXI (X,LDX,N,M,DET,IFLAG)
C
C-----------------------------------------------------------------------
C   MATXXI   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING THE INVERSE OF X'X AND THE DETERMINANT OF X'X WHERE
C        X IS AN N BY M MATRIX WITH N >= M.  THE INVERSE IS RETURNED
C        IN THE FIRST M ROWS AND M COLUMNS OF THE ORIGINAL MATRIX X.
C        ENTRIES BELOW ROW M ARE LIKELY TO BE CHANGED.  NOTE THAT ROWS
C        N+1 AND N+2 OF MATRIX X ARE USED FOR SCRATCH AREA, THEREFORE 
C        LDX >= N+2.  THE STEPS OF THE ALGORITHM ARE DESCRIBED BELOW
C        AS THEY OCCUR.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED APRIL 8, 1988
C
C   REFERENCE: STEWART, G.W., 'INTRODUCTION TO MATRIX COMPUTATIONS',
C              ACADEMIC PRESS, INC., 1973
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C     * X = MATRIX (SIZE NXM) FOR WHICH INV(X'X) AND DET(X'X) ARE TO
C           BE COMPUTED (REAL)
C
C   * LDX = LEADING DIMENSION OF X [LDX>=N+2] (INTEGER)
C
C     * N = NUMBER OF ROWS IN MATRIX X (INTEGER)
C
C     * M = NUMBER OF COLUMNS IN MATRIX X [M<=N] (INTEGER)
C
C     DET = DETERMINANT OF (X'X) (REAL) 
C
C   IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C           -2 -> M>N
C           -1 -> LDX<N+2
C            0 -> NO ERRORS DETECTED
C            K -> THE KTH DIAGONAL ELEMENT OF THE TRIANGULAR MATRIX R,
C                 FROM THE QR FACTORIZATION, IS ZERO
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION X(LDX,*)
      IFLAG = 0
C
C--- CHECK CONSISTENCY OF PASSED PARAMETERS
C
      IF (LDX.LT.N+2) THEN
         IFLAG = -1 
         RETURN
      ENDIF
      IF (N.LT.M) THEN
         IFLAG = -2 
         RETURN
      ENDIF
C
C--- DEFINE SOME CONSTANTS
C
      N1 = N+1
      N2 = N+2
C
C--- COMPUTE QR FACTORIZATION OF MATRIX X (ALGORITHM 3.8, P. 236 OF
C--- THE REFERENCE) 
C
      L = MIN0(N-1,M)
      DO 70 K = 1, L
         E = 0.0
         DO 10 I = K, N
            Q = ABS(X(I,K))
            IF (Q.GT.E) E = Q 
   10    CONTINUE
         IF (E.EQ.0.0) THEN
            X(N+1,K) = 0.0
         ELSE
            DO 20 I = K, N
               X(I,K) = X(I,K)/E
   20       CONTINUE
            S = 0.0 
            DO 30 I = K, N
               S = S+X(I,K)**2
   30       CONTINUE
            S = SIGN(SQRT(S),X(K,K))
            X(K,K) = X(K,K)+S 
            X(N1,K) = S*X(K,K)
            X(N2,K) = -E*S
            DO 60 J = K+1, M
               T = 0.0
               DO 40 I = K, N 
                  T = T+X(I,K)*X(I,J)
   40          CONTINUE
               T = T/X(N1,K)
               DO 50 I = K, N 
                  X(I,J) = X(I,J)-T*X(I,K)
   50          CONTINUE
   60       CONTINUE
         ENDIF
   70 CONTINUE
      IF (N.EQ.M) X(N2,N) = X(N,N)
C
C--- MOVE DIAGONAL ELEMENTS BACK TO DIAGONAL AND COMPUTE DETERMINANT
C
      DET = 1.0
      DO 80 J = 1, M
         X(J,J) = X(N2,J)
         DET = DET*X(J,J)**2
         IF (X(J,J).EQ.0.0) THEN
            IFLAG = J
            RETURN
         ENDIF
   80 CONTINUE
C
C--- INVERT UPPER TRIANGULAR MATRIX R IN PLACE (ALGORITHM 1.5, P. 110 
C--- OF THE REFERENCE)
C
      DO 110 K = M, 1, -1
         X(K,K) = 1.0/X(K,K)
         DO 100 I = K-1, 1, -1
            S = 0.0 
            DO 90 J = I+1, K
               S = S-X(I,J)*X(J,K)
   90       CONTINUE
            X(I,K) = S/X(I,I) 
  100    CONTINUE
  110 CONTINUE
C
C--- COMPUTE INV(X'X) = INV(R)*INV(R)'
C
      DO 140 I = 1, M
         DO 130 J = I, M
            S = 0.0 
            DO 120 K = J, M
               S = S+X(I,K)*X(J,K)
  120       CONTINUE
            X(I,J) = S
            X(J,I) = S
  130    CONTINUE
  140 CONTINUE
      RETURN
      END 
*NEXPER
      SUBROUTINE NEXPER (RED,PINK,BROWN)
C-----------------------------------------------------------------------
C   NEXPER   COPIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND 20899 FROM THE REFERENCE BELOW (COLOR ADDED)
C
C   FOR: COMPUTING THE NEXT PERMUTATION OF THE INTEGERS 1, 2, ..., N. 
C        THE CALLING SEQUENCE IS
C
C                        CALL NEXPER (IPERM,N,LL) 
C
C        WHERE IPERM IS AN INTEGER VECTOR DIMENSIONED AT LEAST N AND
C        LL IS A LOGICAL VARIABLE.  NONE OF THESE PASSED PARAMETERS
C        NEEDS TO BE DEFINED ON INPUT.  ON OUTPUT IPERM CONTAINS THE
C        CURRENT PERMUTATION OF THE FIRST N INTEGERS AND LL IS .TRUE. 
C        UNLESS THIS PERMUTATION IS THE LAST PERMUTATION OF THE CYCLE 
C        (IN WHICH CASE LL IS .FALSE.). 
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED JANUARY 20, 1987
C
C   REFERENCE: NIJENHUIS, ALBERT AND WILF, HERBERT S., 'COMBINATORIAL 
C              ALGORITHMS', ACADEMIC PRESS, 1975, PP. 49-59.
C-----------------------------------------------------------------------
      IMPLICIT INTEGER (A-Z)
      LOGICAL BROWN 
      DIMENSION RED(*)
      DATA MAROON / 0 /
      IF (PINK.EQ.MAROON) GO TO 40
   10 MAROON = PINK 
      ORANGE = 1
      SILVER = 1
      PURPLE = 1
      DO 20 BLUE = 1, PINK
         PURPLE = PURPLE*BLUE 
         RED(BLUE) = BLUE
   20 CONTINUE
   30 BROWN = ORANGE.NE.PURPLE
      RETURN
   40 IF (.NOT.BROWN) GO TO 10
      GO TO (50,60), SILVER
   50 GOLD = RED(2) 
      RED(2) = RED(1)
      RED(1) = GOLD 
      SILVER = 2
      ORANGE = ORANGE+1
      GO TO 30
   60 YELLOW = 3
      BLACK = ORANGE/2
   70 VIOLET = MOD(BLACK,YELLOW)
      IF (VIOLET.NE.0) GO TO 80
      BLACK = BLACK/YELLOW
      YELLOW = YELLOW+1
      GO TO 70
   80 BLACK = PINK
      GREEN = YELLOW-1
      DO 90 BLUE = 1, GREEN
         WHITE = RED(BLUE)-RED(YELLOW)
         IF (WHITE.LT.0) WHITE = WHITE+PINK
         IF (WHITE.GE.BLACK) GO TO 90
         BLACK = WHITE
         INDIGO = BLUE
   90 CONTINUE
      GOLD = RED(YELLOW)
      RED(YELLOW) = RED(INDIGO)
      RED(INDIGO) = GOLD
      SILVER = 1
      ORANGE = ORANGE+1
      RETURN
      END 
*PERMAN
      SUBROUTINE PERMAN (A,LDA,N,IWORK,WORK,APERM)
C
C-----------------------------------------------------------------------
C   PERMAN   COPIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND 20899 FROM THE REFERENCE BELOW
C
C   FOR: COMPUTING THE PERMANENT OF THE N BY N MATRIX A.  THE PERMANENT
C        OF A MATRIX IS SIMILAR TO THE DETERMINANT EXCEPT THAT THE
C        ALTERNATING SIGN CHANGES FOR THE TERMS ARE EXCLUDED.  FOR
C        EXAMPLE, GIVEN THE MATRIX
C
C                               [A B C] 
C                               [D E F] 
C                               [G H I] 
C
C        THE DETERMINANT IS  AEI+BFG+CDH-CEG-BDI-AFH  WHEREAS THE
C        PERMANENT IS  AEI+BFG+CDH+CEG+BDI+AFH .
C
C   NOTE: THE COMPUTING TIME IS PROPORTIONAL TO 2**N, AND ON THE
C         CYBER 180/855 AT NBS A VALUE OF N=20 REQUIRES ABOUT 30
C         SECONDS OF CPU TIME.
C
C   SUBPROGRAMS CALLED: NEXSUB (ATTACHED)
C
C   CURRENT VERSION COMPLETED JANUARY 20, 1987
C
C   REFERENCE: NIJENHUIS, ALBERT AND WILF, HERBERT S., 'COMBINATORIAL 
C              ALGORITHMS', ACADEMIC PRESS, 1975, CH. 1,19. 
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C     * A = MATRIX (SIZE NXN) WHOSE PERMANENT IS TO BE COMPUTED (REAL)
C
C   * LDA = THE LEADING DIMENSION OF MATRIX A [LDA>=N] (INTEGER)
C
C     * N = THE NUMBER OF ROWS AND COLUMNS IN MATRIX A (INTEGER)
C
C   IWORK = VECTOR (LENGTH N) USED AS SCRATCH AREA (INTEGER)
C
C    WORK = VECTOR (LENGTH N) USED AS SCRATCH AREA (REAL)
C
C   APERM = THE PERMANENT OF MATRIX A (REAL)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
      LOGICAL MTC
      DIMENSION A(LDA,*),IWORK(*),WORK(*)
      P = 0.0
      N1 = N-1
      DO 20 I = 1, N
         SUM = 0.0
         DO 10 J = 1, N
            SUM = SUM+A(I,J)
   10    CONTINUE
         WORK(I) = A(I,N)-SUM/2.0
   20 CONTINUE
      SGN = -1.0
   30 SGN = -SGN
      PROD = SGN
      CALL NEXSUB (N1,IWORK,MTC,NCARD,J)
      IF (NCARD.EQ.0) GO TO 50
      Z = 2.0*REAL(IWORK(J))-1.0
      DO 40 I = 1, N
         WORK(I) = WORK(I)+Z*A(I,J)
   40 CONTINUE
   50 DO 60 I = 1, N
         PROD = PROD*WORK(I)
   60 CONTINUE
      P = P+PROD
      IF (MTC) GO TO 30
      APERM = 2.0*REAL(2*MOD(N,2)-1)*P
      RETURN
      END 
C
C=======================================================================
C
      SUBROUTINE NEXSUB (N,IWORK,MTC,NCARD,J)
      LOGICAL MTC
      DIMENSION IWORK(*)
      DATA NLAST / 0 /
      IF (N.EQ.NLAST) GO TO 30
   10 M = 0
      MTC = .TRUE.
      DO 20 I = 1, N
         IWORK(I) = 0
   20 CONTINUE
      NCARD = 0
      NLAST = N
      RETURN
   30 IF (.NOT.MTC) GO TO 10
      M = M+1
      M1 = M
      J = 0
   40 J = J+1
      IF (MOD(M1,2).EQ.1) GO TO 50
      M1 = M1/2
      GO TO 40
   50 L = IWORK(J)
      IWORK(J) = 1-L
      NCARD = NCARD+1-2*L
      MTC = NCARD.NE.1.OR.IWORK(N).EQ.0 
      RETURN
      END 
*PLOTCR
      SUBROUTINE PLOTCR (X,Y,CH,N,IPR)
C
C-----------------------------------------------------------------------
C   PLOTCR   MODIFIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND 20899 FROM THE PROGRAM PLOTC WRITTEN BY JAMES
C            J. FILLIBEN FOR THE SOFTWARE PACKAGE DATAPAC.
C
C   FOR: CREATING A SIMPLE ONE-PAGE PRINTER PLOT OF Y(I) VS. X(I) WITH
C        PLOT CHARACTERS CONTROLLED BY CH(I).  THE PLOTTING GRID IS
C        45 LINES BY 109 COLUMNS WITH THE ENTIRE PLOT TAKING 125
C        COLUMNS.  THE PLOT CHARACTER USED AT THE Ith PLOTTING POSITION
C        WILL BE: 
C
C              1 IF CH(I) IS BETWEEN  0.5 AND  1.5
C              2 IF CH(I) IS BETWEEN  1.5 AND  2.5
C                               : 
C                               : 
C              9 IF CH(I) IS BETWEEN  8.5 AND  9.5
C              0 IF CH(I) IS BETWEEN  9.5 AND 10.5
C              A IF CH(I) IS BETWEEN 10.5 AND 11.5
C              B IF CH(I) IS BETWEEN 11.5 AND 12.5
C                               : 
C                               : 
C              Y IF CH(I) IS BETWEEN 34.5 AND 35.5
C              Z IF CH(I) IS BETWEEN 35.5 AND 36.5
C              * IF CH(I) IS OUTSIDE  0.5  TO 36.5
C              + TO INDICATE MULTIPLE POINTS.
C
C        THE NUMBER OF CHARACTERS PLOTTED IS N.  IF N <= 0 ON INPUT
C        THEN NO PLOT IS DONE AND NO WARNING IS PRINTED.  NO LEADING
C        OR TRAILING BLANK LINES ARE PRINTED.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED NOVEMBER 25, 1987
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C    X(*) = VECTOR (LENGTH N) OF X VALUES PLOTTED HORIZONTALLY (REAL) 
C
C    Y(*) = VECTOR (LENGTH N) OF Y VALUES PLOTTED VERTICALLY (REAL)
C
C   CH(*) = VECTOR (LENGTH N) OF PLOT CHARACTER CONTROL VALUES (REAL) 
C
C       N = NUMBER OF PAIRS (X(I),Y(I)) TO PLOT (INTEGER)
C
C     IPR = LOGICAL UNIT NUMBER OF OUTPUT PRINT FILE (INTEGER)
C-----------------------------------------------------------------------
C
      DIMENSION X(*),Y(*),CH(*),YLABEL(11)
      CHARACTER*1 GR(45,109),PL(37)
      DATA PL / '1','2','3','4','5','6','7','8','9','0','A','B','C','D',
     *   'E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T'
     *   ,'U','V','W','X','Y','Z','*'/
C
C--- RETURN IF N<=0 
C
      IF (N.LE.0) RETURN
C
C--- DETERMINE THE VALUES TO BE LISTED ON THE LEFT VERTICAL AXIS
C
      YMIN = Y(1)
      YMAX = YMIN
      DO 10 I = 1, N
         YMIN = AMIN1(YMIN,Y(I))
         YMAX = AMAX1(YMAX,Y(I))
   10 CONTINUE
      DO 20 I = 1, 9
         YLABEL(I) = YMAX-(REAL(I-1)/8.0)*(YMAX-YMIN)
   20 CONTINUE
C
C--- DETERMINE THE VALUES TO BE LISTED ON THE BOTTOM HORIZONTAL AXIS
C--- DETERMINE XMIN, XMAX, XMID, X25 (=THE 25% POINT), AND
C--- X75 (=THE 75% POINT)
C
      XMIN = X(1)
      XMAX = XMIN
      DO 30 I = 1, N
         XMIN = AMIN1(XMIN,X(I))
         XMAX = AMAX1(XMAX,X(I))
   30 CONTINUE
      XMID = 0.5*(XMIN+XMAX)
      X25 = 0.75*XMIN+0.25*XMAX
      X75 = 0.25*XMIN+0.75*XMAX
C
C--- BLANK OUT THE GRAPH
C
      DO 50 I = 1, 45
         DO 40 J = 1, 109
            GR(I,J) = ' '
   40    CONTINUE
   50 CONTINUE
C
C--- PRODUCE THE VERTICAL AXES
C
      DO 60 I = 3, 43
         GR(I,5) = 'I'
         GR(I,109) = 'I'
   60 CONTINUE
      DO 70 I = 3, 43, 5
         GR(I,5) = '-'
         GR(I,109) = '-'
   70 CONTINUE
      GR(3,1) = '=' 
      GR(3,2) = 'M' 
      GR(3,3) = 'A' 
      GR(3,4) = 'X' 
      GR(23,1) = '='
      GR(23,2) = 'M'
      GR(23,3) = 'I'
      GR(23,4) = 'D'
      GR(43,1) = '='
      GR(43,2) = 'M'
      GR(43,3) = 'I'
      GR(43,4) = 'N'
C
C--- PRODUCE THE HORIZONTAL AXES
C
      DO 80 J = 7, 107
         GR(1,J) = '-'
         GR(45,J) = '-'
   80 CONTINUE
      DO 90 J = 7, 107, 25
         GR(1,J) = 'I'
         GR(45,J) = 'I'
   90 CONTINUE
      DO 100 J = 20, 107, 25
         GR(1,J) = 'I'
         GR(45,J) = 'I'
  100 CONTINUE
C
C--- DETERMINE THE (X,Y) PLOT POSITIONS 
C
      RATIOY = 40.0/(YMAX-YMIN)
      RATIOX = 100.0/(XMAX-XMIN)
      DO 110 I = 1, N
         MX = 7+NINT(RATIOX*(X(I)-XMIN))
         MY = 43-NINT(RATIOY*(Y(I)-YMIN))
         IARG = NINT(CH(I))
         IF (IARG.LT.1.OR.IARG.GT.36) IARG = 37
         IF (GR(MY,MX).EQ.' ') THEN
            GR(MY,MX) = PL(IARG)
         ELSE
            GR(MY,MX) = '+'
         ENDIF
  110 CONTINUE
C
C--- WRITE OUT THE GRAPH
C
      DO 120 I = 1, 45
         IP2 = I+2
         IFLAG = IP2-(IP2/5)*5
         K = IP2/5
         IF (IFLAG.NE.0) THEN 
            WRITE (IPR,1000) (GR(I,J),J=1,109)
         ELSE
            WRITE (IPR,1010) YLABEL(K),(GR(I,J),J=1,109)
         ENDIF
  120 CONTINUE
      WRITE (IPR,1020) XMIN,X25,XMID,X75,XMAX
      RETURN
C
 1000 FORMAT (1X,15X,109A1)
 1010 FORMAT (1X,G15.4,109A1) 
 1020 FORMAT (1X,7X,G20.4,5X,G20.4,5X,G20.4,5X,G20.4,1X,G20.4)
C
      END 
*PPFBET
      SUBROUTINE PPFBET (PR,P,Q,TOL,IFLAG,X)
C
C-----------------------------------------------------------------------
C   PPFBET   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND 20899
C
C   FOR: EVALUATING THE INVERSE CUMULATIVE DISTRIBUTION FUNCTION
C        (PERCENT POINT FUNCTION) OF THE BETA(P,Q) DISTRIBUTION.
C        FOR A GIVEN PROBABILITY PR.  THE PERCENT POINT X IS COMPUTED 
C        TO A SPECIFIED ABSOLUTE ACCURACY TOL WHEN POSSIBLE. THE
C        METHOD OF BRENT, AS DESCRIBED IN THE REFERENCES BELOW, IS
C        USED TO COMPUTE THE APPROXIMATE ZERO OF I(X,P,Q)-PR WHERE
C        I(X,P,Q) IS THE CUMULATIVE DISTRIBUTION FUNCTION OF THE
C        BETA(P,Q) DISTRIBUTION EVALUATED AT X.  THIS METHOD DOES
C        NOT REQUIRE DERIVATIVES.
C
C   NOTE: THE CONSTANT EPS FOR MACHINE FLOATING POINT PRECISION IS
C         MACHINE DEPENDENT.
C
C   SUBPROGRAMS CALLED: CDFBET (BETA CUMULATIVE DISTRIBUTION FUNCTION)
C
C   CURRENT VERSION COMPLETED OCTOBER 13, 1987
C
C   REFERENCES: 
C
C   1) PRESS, WILLIAM H., FLANNERY, BRIAN P., TEUKOLSKY, SAUL A.,
C      AND VETTERLING, WILLIAM T., 'NUMERICAL RECIPES - THE ART OF
C      SCIENTIFIC COMPUTING', CAMBRIDGE UNIVERSITY PRESS, 1986,
C      PP. 251-254. 
C
C   2) BRENT, RICHARD P., 'ALGORITHMS FOR MINIMIZATION WITHOUT
C      DERIVATIVES', PRENTICE-HALL, 1973, CH. 3-4.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C    * PR = A PROBABILITY VALUE IN THE INTERVAL [0,1] (REAL)
C
C     * P = THE FIRST PARAMETER (>0) OF THE BETA(P,Q) DISTRIBUTION
C           (REAL)
C
C     * Q = THE SECOND PARAMETER (>0) OF THE BETA(P,Q) DISTRIBUTION
C           (REAL)
C
C   * TOL = THE REQUIRED ACCURACY (>=1.0E-8) OF THE PERCENT POINT X
C           (REAL)
C
C   IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C             0 -> NO ERRORS DETECTED
C           1,2 -> ERROR FLAGS FROM SUBROUTINE CDFBET
C             3 -> PR<0 OR PR>1
C             4 -> P<=0 OR Q<=0
C             5 -> TOL<1.0E-8 
C             6 -> THE CDF'S AT THE ENDPOINTS HAVE THE SAME SIGN - NO 
C                  VALUE OF X IS DEFINED (THIS SHOULD NEVER OCCUR)
C             7 -> MAXIMUM ITERATIONS EXCEEDED - CURRENT VALUE OF X
C                  RETURNED
C
C       X = THE COMPUTED PERCENT POINT (REAL)
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
C--- DEFINE MAXIMUM ITERATIONS AND MACHINE FLOATING POINT PRECISION
C
      DATA ITMAX,EPS / 50,1.0E-12 /
C
C--- CHECK VALIDITY OF INPUT ARGUMENTS
C
      IFLAG = 0
      IF (PR.LT.0.0.OR.PR.GT.1.0) THEN
         IFLAG = 3
         RETURN
C
      ENDIF
      IF (AMIN1(P,Q).LE.0.0) THEN
         IFLAG = 4
         RETURN
C
      ENDIF
      IF (TOL.LT.1.0E-8) THEN 
         IFLAG = 5
         RETURN
C
      ENDIF
      A = 0.0
      B = 1.0
      FA = -PR
      FB = 1.0-PR
C
C--- CHECK FOR ZERO BEING BRACKETED
C
      IF (FB*FA.GT.0.0) THEN
         IFLAG = 6
         RETURN
C
      ENDIF
      FC = FB
      DO 10 ITER = 1, ITMAX
         IF (FB*FC.GT.0.0) THEN
C
C--- RENAME A,B,C AND ADJUST BOUNDING INTERVAL D
C
            C = A
            FC = FA 
            D = B-A 
            E = D
         ENDIF
         IF (ABS(FC).LT.ABS(FB)) THEN
            A = B
            B = C
            C = A
            FA = FB 
            FB = FC 
            FC = FA 
         ENDIF
C
C--- CONVERGENCE CHECK
C
         TOL1 = 2.0*EPS*ABS(B)+0.5*TOL
         XM = 0.5*(C-B)
         IF (ABS(XM).LE.TOL1.OR.FB.EQ.0.0) THEN
            X = B
            RETURN
C
         ENDIF
         IF (ABS(E).GE.TOL1.AND.ABS(FA).GT.ABS(FB)) THEN
C
C--- ATTEMPT INVERSE QUADRATIC INTERPOLATION
C
            S = FB/FA
            IF (A.EQ.C) THEN
               U = 2.0*XM*S
               V = 1.0-S
            ELSE
               V = FA/FC
               R = FB/FC
               U = S*(2.0*XM*V*(V-R)-(B-A)*(R-1.0))
               V = (V-1.0)*(R-1.0)*(S-1.0)
            ENDIF
C
C--- CHECK WHETHER IN BOUNDS
C
            IF (U.GT.0.0) V = -V
            U = ABS(U)
            IF (2.0*U.LT.AMIN1(3.0*XM*V-ABS(TOL1*V),ABS(E*V))) THEN
C
C--- ACCEPT INTERPLATION
C
               E = D
               D = U/V
            ELSE
C
C--- INTERPOLATION FAILED, USE BISECTION
C
               D = XM
               E = D
            ENDIF
         ELSE
C
C--- BOUNDS DECREASING TOO SLOWLY, USE BISECTION
C
            D = XM
            E = D
         ENDIF
C
C--- MOVE LAST BEST GUESS TO A
C
         A = B
         FA = FB
C
C--- EVALUATE NEW TRIAL ZERO
C
         IF (ABS(D).GT.TOL1) THEN
            B = B+D 
         ELSE
            B = B+SIGN(TOL1,XM)
         ENDIF
         CALL CDFBET (B,P,Q,EPS,IFLAG,CDF)
         IF (IFLAG.NE.0) RETURN
         FB = CDF-PR
   10 CONTINUE
      IFLAG = 7
      X = B
      RETURN
C
      END 
*QDASIM
      SUBROUTINE QDASIM (F,A,B,MXPOW2,ABSERR,AREA,IFLAG)
C
C-----------------------------------------------------------------------
C   QDASIM   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899. 
C
C   FOR: NUMERICAL INTEGRATION USING AN ADAPTIVE SIMPSON'S RULE.  THE 
C        NUMBER OF PARTITIONS IS DOUBLED UNTIL SUCCESSIVE INTEGRAL
C        VALUES AGREE TO ABSERR OR BETTER.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED JUNE 9, 1986
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C        F = USER-DEFINED REAL FUNCTION WITH ONE ARGUMENT, I.E., F(X) 
C
C        A = LEFT ENDPOINT OF REGION OF INTEGRATION (REAL)
C
C        B = RIGHT ENDPOINT OF REGION OF INTEGRATION (REAL) 
C
C   MXPOW2 = INTEGER INDICATING THAT THE MAXIMUM NUMBER OF PARTITIONS 
C            TO BE USED IS 2**MXPOW2
C
C   ABSERR = ABSOLUTE ERROR TOLERANCE ON THE COMPUTED VALUE OF THE
C            INTEGRAL (REAL)
C
C     AREA = COMPUTED VALUE OF THE INTEGRAL (REAL)
C
C    IFLAG = INTEGER ERROR FLAG ON RETURN.     INTERPRETATION: 
C            0 -> INTEGRAL COMPUTED
C            1 -> INTEGRAL NOT COMPUTED SINCE A>B 
C            2 -> MAXIMUM PARTITION SIZE REACHED WITHOUT CONVERGENCE
C-----------------------------------------------------------------------
C
      IFLAG = 0
      AREA = 0.0
      IF (A.GT.B) THEN
         IFLAG = 1
      ELSE
         H = 0.5*(B-A)
         S = 2.0*F(A+H)
         R = F(A)
         R = R+F(B) 
         C = H/3.0
         AOLD = C*(R+S+S)
         N = 2
         DO 20 I = 1, MXPOW2
            N = N+N 
            H = 0.5*H
            R = R+S 
            S = 0.0 
            DO 10 J = 1, N, 2 
               X = A+H*REAL(J)
               S = S+F(X)
   10       CONTINUE
            S = S+S 
            C = 0.5*C
            AREA = C*(R+S+S)
            IF (ABS(AREA-AOLD).LE.ABSERR) RETURN
            AOLD = AREA
   20    CONTINUE
         IFLAG = 2
      ENDIF
      RETURN
C
      END 
*QDTANH
      SUBROUTINE QDTANH (F,A,B,N,AREA,IFLAG)
C
C-----------------------------------------------------------------------
C   QDTANH   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND 20899.
C
C   FOR: NUMERICAL INTEGRATION USING THE 'TANH RULE' AS DESCRIBED IN
C        THE REFERENCE BELOW.  THIS RULE IS ESPECIALLY GOOD FOR SMOOTH
C        FUNCTIONS WHICH ARE INTEGRABLE BUT HAVE SINGULARITIES OR
C        UNBOUNDED FIRST DERIVATIVES AT THE ENDPOINTS.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED JUNE 9, 1986
C
C   REFERENCE: HABER, SEYMOUR, 'THE TANH RULE FOR NUMERICAL 
C              INTEGRATION', SIAM JOURNAL OF NUMERICAL ANALYSIS,
C              VOL. 14, NO. 4, SEPTEMBER 1977.
C-----------------------------------------------------------------------
C    DEFINITION OF PASSED PARAMETERS: 
C
C        F = USER-DEFINED REAL FUNCTION WITH ONE ARGUMENT, I.E., F(X) 
C
C        A = LEFT ENDPOINT OF REGION OF INTEGRATION (REAL)
C
C        B = RIGHT ENDPOINT OF REGION OF INTEGRATION (REAL) 
C
C        N = AN INTEGER GREATER THAN ZERO WHICH DETERMINES THE NUMBER 
C            OF ABSCISSAE AT WHICH THE FUNCTION F IS EVALUATED (2*N+1)
C
C     AREA = COMPUTED VALUE OF THE INTEGRAL
C
C    IFLAG = INTEGER ERROR FLAG ON RETURN.     INTERPRETATION: 
C            0 -> INTEGRAL COMPUTED
C            1 -> INTEGRAL NOT COMPUTED SINCE N<1 
C            2 -> INTEGRAL NOT COMPUTED SINCE A>B 
C-----------------------------------------------------------------------
C
      IFLAG = 0
      AREA = 0.0
      IF (N.LT.1) THEN
         IFLAG = 1
         RETURN
C
      ENDIF
      IF (A.GT.B) THEN
         IFLAG = 2
      ELSE
         CON = 0.5/REAL(N)
         H = 4.0*ATAN(1.0)*SQRT(CON)-CON
         HWIDTH = 0.5*(B-A)
         X = 0.5*(A+B)
         AREA = F(X)
         DO 10 I = 1, N
            HI = H*REAL(I)
            Y = HWIDTH*TANH(HI)
            FSUM = F(X+Y)
            AREA = AREA+(FSUM+F(X-Y))/COSH(HI)**2 
   10    CONTINUE
         AREA = HWIDTH*H*AREA 
      ENDIF
      RETURN
C
      END 
*QDTRAP
      FUNCTION QDTRAP (F,A,B,N,FDA,FDB,IFLAG)
C
C--- CORRECTED TRAPEZOIDAL RULE FOR INTEGRATION
C--- CPR  JUNE 5, 1986   REF: CONTE & DEBOOR
C
      H = (B-A)/REAL(N)
      FA = F(A)
      QDTRAP = 0.5*(FA+F(B))
      DO 10 I = 1, N-1
         X = A+H*REAL(I)
         QDTRAP = QDTRAP+F(X) 
   10 CONTINUE
      QDTRAP = H*QDTRAP
      IF (IFLAG.EQ.1) QDTRAP = QDTRAP+(FDA-FDB)*H*H/12.0
      RETURN
C
      END 
*QSMNMX
      SUBROUTINE QSMNMX (A,B,C,D,E,F,XA,XB,YA,YB,XMIN,YMIN,ZMIN,XMAX, 
     *   YMAX,ZMAX) 
C
C-----------------------------------------------------------------------
C   QSMNMX   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: FINDING THE MINIMUM (ZMIN) AND MAXIMUM (ZMAX) VALUES OF THE
C        QUADRATIC SURFACE DEFINED BY THE FUNCTION
C
C             G(X,Y) = A*X**2 + B*X*Y + C*Y**2 + D*X + E*Y + F
C
C        OVER THE RECTANGLE XA<=X<=XB AND YA<=Y<=YB.  THE X AND Y
C        VALUES CORRESPONDING TO ZMIN (XMIN,YMIN) AND ZMAX (XMAX,YMAX)
C        ARE RETURNED.  IF XA>XB OR YA>YB ON INPUT THEN NOTHING IS
C        DONE.  ALL PASSED PARAMETERS ARE OF TYPE REAL.
C
C   SUBPROGRAMS CALLED: MINMAX (ATTACHED)
C
C   CURRENT VERSION COMPLETED JANUARY 15, 1987
C-----------------------------------------------------------------------
C
      G(X,Y) = X*(A*X+B*Y+D)+Y*(C*Y+E)+F
C
C--- DO NOTHING IF XA>XB OR YA>YB
C
      IF (XA.GT.XB.OR.YA.GT.YB) RETURN
C
C--- INITIALIZE MINIMUM AND MAXIMUM VALUES OF FUNCTION G
C
      ZMIN = 1.0E30 
      ZMAX = -ZMIN
C
C--- EVALUATE FUNCTION AT CORNERS (XA,YA), (XB,YA), (XB,YA), (XB,YB)
C
      Z = G(XA,YA)
      CALL MINMAX (XA,YA,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
      Z = G(XA,YB)
      CALL MINMAX (XA,YB,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
      Z = G(XB,YA)
      CALL MINMAX (XB,YA,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
      Z = G(XB,YB)
      CALL MINMAX (XB,YB,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
C
C--- CHECK BOUNDARIES Y=YA AND Y=YB FOR CANDIDATE POINTS
C
      IF (A.NE.0.0) THEN
         X = -0.5*(B*YA+D)/A
         IF (X.GT.XA.AND.X.LT.XB) THEN
            Z = G(X,YA)
            CALL MINMAX (X,YA,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
         ENDIF
         X = -0.5*(B*YB+D)/A
         IF (X.GT.XA.AND.X.LT.XB) THEN
            Z = G(X,YB)
            CALL MINMAX (X,YB,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
         ENDIF
      ENDIF
C
C--- CHECK BOUNDARIES X=XA AND X=XB FOR CANDIDATE POINTS
C
      IF (C.NE.0.0) THEN
         Y = -0.5*(B*XA+E)/C
         IF (Y.GT.YA.AND.Y.LT.YB) THEN
            Z = G(XA,Y)
            CALL MINMAX (XA,Y,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
         ENDIF
         Y = -0.5*(B*XB+E)/C
         IF (Y.GT.YA.AND.Y.LT.YB) THEN
            Z = G(XB,Y)
            CALL MINMAX (XB,Y,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
         ENDIF
      ENDIF
C
C--- CHECK INTERIOR XA<X<XB AND YA<Y<YB FOR CANDIDATE POINT 
C
      DEN = B**2-4.0*A*C
      IF (DEN.NE.0.0) THEN
         X = (2.0*C*D-B*E)/DEN
         Y = (2.0*A*E-B*D)/DEN
         IF (X.GT.XA.AND.X.LT.XB.AND.Y.GT.YA.AND.Y.LT.YB) THEN
            Z = G(X,Y)
            CALL MINMAX (X,Y,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
         ENDIF
      ENDIF
      RETURN
      END 
C=======================================================================
      SUBROUTINE MINMAX (X,Y,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX)
C
C--- UPDATE THE MINIMUM OR MAXIMUM VALUE OF Z IF APPROPRIATE
C
      IF (Z.LT.ZMIN) THEN
         XMIN = X
         YMIN = Y
         ZMIN = Z
      ENDIF
      IF (Z.GT.ZMAX) THEN
         XMAX = X
         YMAX = Y
         ZMAX = Z
      ENDIF
      RETURN
      END 
*RDBETA
      FUNCTION RDBETA (P,Q)
C
C-----------------------------------------------------------------------
C   RDBETA   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C
C   FOR: GENERATING A RANDOM DEVIATE FROM THE BETA(P,Q) DISTRIBUTION
C        WHERE P,Q > 0.  IF THIS CONDITION IS NOT MET THEN PROGRAM
C        EXECUTION IS TERMINATED.
C
C        IF MIN(P,Q) > 1 THEN THE CHENG ALGORITHM BB IS USED.
C        IF 0 < MIN(P,Q) <= 1 THEN THE CHENG ALGORITHM BC IS USED.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - UNIFORM(0,1) GENERATOR
C
C   CURRENT VERSION COMPLETED JULY 30, 1984
C
C   REFERENCE: KENNEDY, WILLIAM J. AND GENTLE, JAMES E., "STATISTICAL 
C              COMPUTING", MARCEL DEKKER, INC., 1980, PP. 216-219
C
C   REFERENCE: CHENG, R.C.H., "GENERATING BETA VARIATES WITH
C              NONINTEGRAL SHAPE PARAMETERS", COMMUNICATIONS OF
C              THE ACM, VOL. 21, NO. 4, APRIL 1978, PP. 317-322.
C-----------------------------------------------------------------------
C
      H = AMIN1(P,Q)
      IF (H.GT.1.0) THEN
C
C
C     CHENG ALGORITHM (BB)
C
         A = H
         B = AMAX1(P,Q)
         ALPHA = A+B
         BETA = SQRT((ALPHA-2.)/(2.*A*B-ALPHA))
         GAMMA = A+1./BETA
C
C--- STEP 1
C
   10    U1 = RDUNI(0)
         U2 = RDUNI(0)
         V = BETA*ALOG(U1/(1.-U1))
         W = A*EXP(V)
         Z = U2*U1**2
         R = GAMMA*V-1.3862944
         S = A+R-W
C
C--- STEP 2
C
         IF (S+2.609438.GE.5.*Z) GO TO 20
C
C--- STEP 3
C
         T = ALOG(Z)
         IF (S.GE.T) GO TO 20 
C
C--- STEP 4
C
         IF (R+ALPHA*ALOG(ALPHA/(B+W)).LT.T) GO TO 10
C
C--- STEP 5
C
   20    IF (A.EQ.P) THEN
            RDBETA = W/(B+W)
         ELSE
            RDBETA = B/(B+W)
         ENDIF
      ELSEIF (H.GT.0.0) THEN
C
C
C     CHENG ALGORITHM (BC)
C
         A = AMAX1(P,Q)
         B = H
         ALPHA = A+B
         BETA = 1./B
         DELTA = 1.+A-B
         KAPPA1 = DELTA*(0.013888889+0.041666667*B)/(A*BETA-0.77777778)
         KAPPA2 = 0.25+(0.5+0.25/DELTA)*B
C
C--- STEP 1
C
   30    U1 = RDUNI(0)
         U2 = RDUNI(0)
         IF (U1.GE.0.5) GO TO 40
C
C--- STEP 2
C
         Y = U1*U2
         Z = U1*Y
         IF (0.25*U2+Z-Y.GE.KAPPA1) THEN
            GO TO 30
C
         ELSE
            GO TO 50
C
         ENDIF
C
C--- STEP 3
C
   40    Z = U2*U1**2
         IF (Z.LE.0.25) THEN
            V = BETA*ALOG(U1/(1.-U1))
            W = A*EXP(V)
            GO TO 60
C
         ENDIF
C
C--- STEP 4
C
         IF (Z.GE.KAPPA2) GO TO 30
C
C--- STEP 5
C
   50    V = BETA*ALOG(U1/(1.-U1))
         W = A*EXP(V)
         IF (ALPHA*(ALOG(ALPHA/(B+W))+V)-1.3862944.LT.ALOG(Z)) GO TO 30
C
C--- STEP 6
C
   60    IF (A.EQ.P) THEN
            RDBETA = W/(B+W)
         ELSE
            RDBETA = B/(B+W)
         ENDIF
      ELSE
C
C
C     PARAMETERS OUT-OF-RANGE 
C
C
         PRINT *,' *** EITHER P OR Q IS <= 0'
         PRINT *,' *** EXECUTION STOPPED IN FUNCTION RDBETA'
         STOP
C
      ENDIF
      RETURN
      END 
*RDCHI2
      FUNCTION RDCHI2 (DF)
C
C-----------------------------------------------------------------------
C   RDCHI2   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C
C   FOR: GENERATING A RANDOM DEVIATE FROM THE CHI-SQUARED(DF)
C        DISTRIBUTION WHERE DF > 0 AND DOES NOT HAVE TO BE AN INTEGER.
C        A CHI-SQUARED(DF) DEVIATE IS EQUAL TO TWICE A GAMMA(DF/2)
C        DEVIATE.  SEE REFERENCE BELOW. 
C
C   SUBPROGRAMS CALLED: RDGAMM (STSPAC) 
C                         RDNOR (STSPAC)
C
C   CURRENT VERSION COMPLETED FEBRUARY 28, 1986
C
C   REFERENCE: KENNEDY, WILLIAM J. AND GENTLE, JAMES E., 'STATISTICAL 
C              COMPUTING', MARCEL DEKKER, INC., 1980, PP. 209-216
C-----------------------------------------------------------------------
      RDCHI2 = 2.*RDGAMM(DF/2.)
      RETURN
C
      END 
*RDF
      FUNCTION RDF (DF1,DF2)
C
C-----------------------------------------------------------------------
C   RDF   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C         DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C
C   FOR: GENERATING A RANDOM DEVIATE FROM THE F(DF1,DF2) DISTRIBUTION 
C        WHERE DF1,DF2 > 0 CAN BE NON-INTEGERS.  IF X IS A
C        BETA(DF1/2,DF2/2) DEVIATE THEN DF2*X/(DF1*(1-X)) IS AN
C        F(DF1,DF2) DEVIATE.  SEE REFERENCE BELOW.
C
C   SUBPROGRAMS CALLED: RDBETA (STSPAC) 
C
C   CURRENT VERSION COMPLETED FEBRUARY 28, 1986
C
C   REFERENCE: KENNEDY, WILLIAM J. AND GENTLE, JAMES E., 'STATISTICAL 
C              COMPUTING', MARCEL DEKKER, INC., 1980, PP. 219-221
C-----------------------------------------------------------------------
C
      P = DF1/2.
      Q = DF2/2.
   10 X = RDBETA(P,Q)
      IF (X.GE.1.0) THEN
         PRINT *,' *** A BETA DEVIATE OF 1.0 CANNOT BE TRANSFORMED TO'
         PRINT *,' *** AN F DEVIATE --- GENERATE ANOTHER BETA DEVIATE'
         GO TO 10
C
      ELSE
         RDF = DF2*X/(DF1*(1.-X))
      ENDIF
      RETURN
      END 
*RDGAMM
      FUNCTION RDGAMM (ALPHA) 
C
C-----------------------------------------------------------------------
C   RDGAMM   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, MD 
C
C   FOR: GENERATING A RANDOM DEVIATE FROM THE GAMMA(ALPHA) DISTRIBUTION.
C        ONE OF THREE METHODS IS USED DEPENDING ON THE VALUE OF THE
C        PARAMETER ALPHA (WHICH DOES NOT HAVE TO BE AN INTEGER): 
C
C             VALUE OF A              METHOD USED 
C            -------------         ------------------
C            0 < ALPHA < 1      AHRENS (GS)   [FIRST REFERENCE]
C              ALPHA = 1        LOG TRANSFORMATION
C              ALPHA > 1        DIETER-AHRENS (GD)   [SECOND REFERENCE]
C
C        IF A <= 0 AN ERROR MESSAGE IS PRINTED AND EXECUTION IS
C        TERMINATED.
C
C        DESCRIPTIONS OF EACH OF THESE ALGORITHMS CAN BE FOUND IN
C        THE REFERENCE GIVEN BELOW.
C
C   SUBPROGRAMS CALLED:  RDUNI (STSPAC) - UNIFORM(0,1) GENERATOR
C                       RDNOR (STSPAC) - NORMAL(0,1) GENERATOR
C
C   CURRENT VERSION COMPLETED FEBRUARY 28, 1986
C
C   REFERENCE: KENNEDY, WILLIAM J. AND GENTLE, JAMES E., 'STATISTICAL 
C              COMPUTING', MARCEL DEKKER, INC., 1980, PP. 209-216
C
C   REFERENCE: AHRENS, J.H. AND DIETER, U., 'GENERATING GAMMA VARIATES
C              BY A MODIFIED REJECTION TECHNIQUE', COMMUNICATIONS OF
C              THE ACM, VOL. 25, NO. 1, PP. 47-54, JANUARY 1982.
C-----------------------------------------------------------------------
C
      DIMENSION Q(8),A(8),E(6)
      DATA AP,APP / 0.,0. /
      DATA Q / +0.041666661,+0.020834040,+0.007970958,+0.001686911,-
     *   0.000775882,+0.001274709,-0.000506403,+0.000213943 /
      DATA A / +0.33333332,-0.24999995,+0.20000622,-0.16667748,+
     *   0.14236572,-0.12438558,+0.12337954,-0.11275089 /
      DATA E / 1.,0.50000027,0.1666605,0.04171864,0.00813673,0.00172501
     *   /
      IF (ALPHA.GT.1.0) THEN
C
C
C     DIETER-AHRENS ALGORITHM (GD)
C
C--- STEP 1
C
         IF (ALPHA.NE.AP) THEN
            AP = ALPHA
            S2 = ALPHA-0.5
            S = SQRT(S2)
            D = SQRT(32.0)-12.*S
         ENDIF
C
C--- STEP 2
C
         T = RDNOR(0)
         X = S+0.5*T
         IF (T.GE.0.0) THEN
            RDGAMM = X**2
            RETURN
C
         ENDIF
C
C--- STEP 3
C
         U = RDUNI(0)
         IF (D*U.LE.T**3) THEN
            RDGAMM = X**2
            RETURN
C
         ENDIF
C
C--- STEP 4
C
         IF (ALPHA.NE.APP) THEN
            APP = ALPHA
            Q0 = 0. 
            DO 10 K = 1, 8
               Q0 = Q0+Q(K)*ALPHA**(-K) 
   10       CONTINUE
            IF (ALPHA.LE.3.686) THEN
               B = 0.463+S+0.178*S2
               SIGMA = 1.235
               C = 0.195/S-0.079+0.16*S 
            ELSEIF (ALPHA.LE.13.022) THEN
               B = 1.654+0.0076*S2
               SIGMA = 1.68/S+0.275
               C = 0.062/S+0.024
            ELSE
               B = 1.77
               SIGMA = 0.75
               C = 0.1515/S
            ENDIF
         ENDIF
C
C--- STEP 5
C
         IF (X.LE.0.0) GO TO 30
C
C--- STEP 6
C
         V = T/(S+S)
         IF (ABS(V).GT.0.25) THEN
            QQ = Q0-S*T+0.25*T**2+(S2+S2)*ALOG(1.+V)
         ELSE
            SUM = 0.
            DO 20 K = 1, 8
               SUM = SUM+A(K)*V**K
   20       CONTINUE
            QQ = Q0+(0.5*T**2)*SUM
         ENDIF
C
C--- STEP 7
C
         IF (ALOG(1.-U).LE.QQ) THEN
            RDGAMM = X**2
            RETURN
C
         ENDIF
C
C--- STEP 8
C
   30    EE = -ALOG(RDUNI(0)) 
         U = 2.*RDUNI(0)-1.
         T = B+SIGN(EE*SIGMA,U)
C
C--- STEP 9
C
         IF (T.LE.-0.71874484) GO TO 30 
C
C--- STEP 10
C
         V = T/(S+S)
         IF (ABS(V).GT.0.25) THEN
            QQ = Q0-S*T+0.25*T**2+(S2+S2)*ALOG(1.+V)
         ELSE
            SUM = 0.
            DO 40 K = 1, 8
               SUM = SUM+A(K)*V**K
   40       CONTINUE
            QQ = Q0+(0.5*T**2)*SUM
         ENDIF
C
C--- STEP 11
C
         IF (QQ.LE.0.0) GO TO 30
         ET = EE-0.5*T**2
         IF (QQ.LE.0.5) THEN
            FACTOR = 0.
            DO 50 K = 1, 6
               FACTOR = FACTOR+E(K)*QQ**K
   50       CONTINUE
            EXPROD = FACTOR*EXP(ET)
         ELSE
            EXPROD = EXP(QQ+ET)-EXP(ET) 
         ENDIF
         IF (C*ABS(U).GT.EXPROD) GO TO 30
C
C--- STEP 12
C
         RDGAMM = (S+0.5*T)**2
      ELSEIF (ALPHA.EQ.1.0) THEN
C
C
C     LOG TRANSFORMATION
C
C
         RDGAMM = -ALOG(RDUNI(0))
      ELSEIF (ALPHA.GT.0.0) THEN
C
C
C     AHRENS ALGORITHM (GS)
C
C--- STEP 1
C
         EE = 2.7182818
   60    U = RDUNI(0)
         B = (EE+ALPHA)/EE
         P = B*U
         IF (P.GT.1.0) GO TO 70
C
C--- STEP 2
C
         RDGAMM = P**(1.0/ALPHA)
         V = RDUNI(0)
         IF (V.GT.EXP(-RDGAMM)) GO TO 60
         RETURN
C
C--- STEP 3
C
   70    RDGAMM = -ALOG((B-P)/ALPHA)
         V = RDUNI(0)
         IF (V.GT.RDGAMM**(ALPHA-1.0)) GO TO 60
         RETURN
C
      ELSE
         PRINT *,' *** THE GAMMA PARAMETER MUST BE > 0'
         PRINT *,' *** EXECUTION STOPPED IN FUNCTION RDGAMM'
         STOP
      ENDIF
      RETURN
C
      END 
*RDMNOR
      SUBROUTINE RDMNOR (AMU,SIG,LDSIG,N,LTF,ZM,IFLAG)
C
C-----------------------------------------------------------------------
C   RDMNOR   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: COMPUTING A VECTOR OF PSEUDO-RANDOM MULTIVARIATE NORMAL
C        DEVIATES WITH MEAN AMU AND VARIANCE SIG.  BEFORE THE FIRST
C        CALL WITH A GIVEN SIG, THE LOGICAL VARIABLE LTF MUST BE
C        ASSIGNED THE VALUE .TRUE. SO THAT A CHOLESKY FACTORIZATION
C        OF SIG IS PERFORMED AS IN THE REFERENCE.  THE RESULT IS A
C        LOWER TRIANGULAR MATRIX L, STORED IN THE LOWER TRIANGLE OF
C        SIG, SUCH THAT SIG=LL'.  FURTHER CALL TO RDMNOR USE ONLY THE 
C        LOWER TRIANGLE OF SIG (L) IN COMPUTING THE DEVIATES UNTIL THE
C        VALUE OF LTF IS RESET TO .TRUE., EVEN IS SIG IS REDEFINED.
C
C   NOTE: BEFORE THE FIRST CALL TO THIS ROUTINE THE CALL
C
C                           Z = RDNOR(ISEED)
C
C         SHOULD BE MADE IN ORDER TO INITIALIZE THE NORMAL RANDOM
C         NUMBER GENERATOR WHERE ISEED IS A POSITIVE INTEGER.  THIS
C         ALLOWS THE USER TO ESTABLISH A REPEATABLE SEQUENCE OF
C         DEVIATES. 
C
C   SUBPROGRAMS CALLED: RDNOR (PSEUDO-RANDOM NORMAL GENERATOR)
C
C   CURRENT VERSION COMPLETED MAY 15, 1987
C
C   REFERENCE: STEWART, G.W., 'INTRODUCTION TO MATRIX COMPUTATIONS',
C              ACADEMIC PRESS, ALGORITHM 3.9, P 142.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C     * AMU = MEAN VECTOR (LENGTH N) OF THE MULTIVARIATE NORMAL
C             DEVIATES ZM (REAL)
C
C     * SIG = COVARIANCE MATRIX (SIZE NXN) OF THE MULTIVARIATE NORMAL 
C             DEVIATES ZM (REAL)
C
C   * LDSIG = THE LEADING DIMENSION OF MATRIX SIG (>=N) (INTEGER)
C
C       * N = THE LENGTH OF THE VECTOR OF DEVIATES ZM (INTEGER)
C
C     * LTF = AN INDICATOR VARIABLE FOR PERFORMING A CHOLESKY
C             FACTORIZATION OF A NEW COVARIANCE MATRIX SIG (LOGICAL)
C
C        ZM = A PSEUDO-RANDOM MULTIVARIATE NORMAL VECTOR (LENGTH N)
C             WITH MEAN AMU AND VARIANCE SIG (REAL)
C
C     IFLAG = AN ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C             0 -> NO ERRORS DETECTED
C             1 -> THE MATRIX SIG IS NOT POSITIVE SEMIDEFINITE, THUS
C                  CANNOT BE A COVARIANCE MATRIX - NO DEVIATE GENERATED
C
C   * INDICATES VARIABLES REQUIRING INPUT VALUES
C-----------------------------------------------------------------------
      DIMENSION AMU(*),SIG(LDSIG,*),ZM(*)
      LOGICAL LTF
C
C--- IF NEW MATRIX SIG, PERFORM CHOLESKY FACTORIZATION.  SET ERROR
C--- FLAG IF SIG IS NOT POSITIVE DEFINITE
C
      IF (LTF) THEN 
         DO 40 K = 1, N
            DO 20 I = 1, K-1
               S = 0.0
               DO 10 J = 1, I-1
                  S = S+SIG(I,J)*SIG(K,J)
   10          CONTINUE
               SIG(K,I) = (SIG(K,I)-S)/SIG(I,I)
   20       CONTINUE
            S = 0.0 
            DO 30 J = 1, K-1
               S = S+SIG(K,J)**2
   30       CONTINUE
            Q = SIG(K,K)-S
            IF (Q.LT.0.0) THEN
               IFLAG = 1
               RETURN
            ELSE
               SIG(K,K) = SQRT(Q)
            ENDIF
   40    CONTINUE
         LTF = .FALSE.
      ENDIF
      IFLAG = 0
C
C--- COMPUTE N INDEPENDENT N(0,1) PSEUDO-RANDOM DEVIATES IN ZM
C
      DO 50 I = 1, N
         ZM(I) = RDNOR(0)
   50 CONTINUE
C
C--- COMPUTE THE PSEUDO-RANDOM MULTIVARIATE NORMAL DEVIATES IN ZM
C
      DO 70 I = N, 1, -1
         S = 0.0
         DO 60 J = 1, I
            S = S+SIG(I,J)*ZM(J)
   60    CONTINUE
         ZM(I) = AMU(I)+S
   70 CONTINUE
      RETURN
      END 
*RDNF
      FUNCTION RDNF (DF1,DF2,ALAMB1,ALAMB2)
C
C-----------------------------------------------------------------------
C   RDNF   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C          DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C
C   FOR: GENERATING A DOUBLY NONCENTRAL F DEVIATE AS THE RATIO OF TWO 
C        NONCENTRAL CHI-SQUARED DEVIATES WITH DEGREES OF FREEDOM DF1
C        AND DF2 (BOTH >= 1) AND NONCENTRALITY PARAMETERS ALAMB1 AND
C        ALAMB2 (BOTH >= 0).
C
C   SUBPROGRAMS CALLED: RDCHI2 (STSPAC) - NONCENTRAL CHI-SQUARED DEVIATE
C                       RDNOR (STSPAC) - NORMAL(0,1) DEVIATE
C
C   CURRENT VERSION COMPLETED FEBRUARY 28, 1986
C-----------------------------------------------------------------------
      IF (DF1.GE.1.0.AND.DF2.GE.1.0) GO TO 10
         PRINT *,' *** BOTH DEGREES OF FREEDOM MUST BE >= 1.0'
         PRINT *,' *** EXECUTION STOPPED IN FUNCTION RDNF'
         STOP
10       X1=(RDNOR(0)+SQRT(ALAMB1))**2
         IF(DF1.GT.1.0) X1=X1+RDCHI2(DF1-1.0)
         X2=(RDNOR(0)+SQRT(ALAMB2))**2
         IF(DF2.GT.1.0) X2=X2+RDCHI2(DF2-1.0)
      RDNF = DF2*X1/(DF1*X2)
      RETURN
C
      END 
*RDNOR
      FUNCTION RDNOR (JD)
C
C-----------------------------------------------------------------------
C   RDNOR   OBTAINED BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C           DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C           MD 20899 FROM THE SOFTWARE LIBRARY CMLIB AT NBS 
C
C   FOR: GENERATING NORMAL PSEUDO-RANDOM DEVIATES WITH MEAN ZERO AND
C        VARIANCE ONE.  THE INITIAL CALL TO RDNOR SHOULD BE OF THE
C        FORM Z=RDNOR(K) WHERE K IS A POSITIVE INTEGER.  FURTHER
C        CALLS SHOULD BE OF THE FORM Z=RDNOR(0).  MDIG IS A LOWER
C        BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE FOR 
C        REPRESENTING INTEGERS, INCLUDING THE SIGN BIT.  THIS VALUE
C        MUST BE AT LEAST 16.  TO ALLOW THE LONGEST SEQUENCE OF
C        RANDOM DEVIATES WITHOUT CYCLING, SET MDIG TO ITS MAXIMUM
C        VALUE.
C
C        THE NBS CENTRAL COMPUTERS WERE FOUND TO GENERATE DEVIATES
C        AT THE FOLLOWING RATES USING THE INDICATED MAXIMUM VALUES
C        OF MDIG: 
C
C             CYBER 180/855 ->   54,893 DEV/SEC   (MDIG=49) 
C             CYBER 205/622 ->   89,098 DEV/SEC   (MDIG=48) 
C
C   SUBPROGRAMS CALLED: RDUNI 
C
C   CURRENT VERSION COMPLETED FEBRUARY 9, 1987
C
C   REFERENCE: MARSAGLIA & TSANG, 'FAST, EASILY IMPLEMENTED METHOD
C              FOR SAMPLING FROM DECREASING OR SYMMETRIC UNIMODAL
C              DENSITY FUNCTIONS', SIAM J SISC 1983.
C-----------------------------------------------------------------------
C
      REAL V(65),W(65)
      INTEGER M(17) 
      SAVE I1,J1,M,M1,M2,RMAX 
      DATA AA,B,C,RMAX / 12.37586,.4878992,12.67706,3.0518509E-5 /
      DATA C1,C2,PC,XN / 0.9689279,1.301198,.1958303E-1,2.776994 /
      DATA V / 0.3409450,0.4573146,0.5397793,0.6062427,0.6631691,
     *   0.7136975,0.7596125,0.8020356,0.8417227,0.8792102,0.9148948, 
     *   0.9490791,0.9820005,1.0138492,1.0447810,1.0749254,1.1043917, 
     *   1.1332738,1.1616530,1.1896010,1.2171815,1.2444516,1.2714635, 
     *   1.2982650,1.3249008,1.3514125,1.3778399,1.4042211,1.4305929, 
     *   1.4569915,1.4834526,1.5100121,1.5367061,1.5635712,1.5906454, 
     *   1.6179680,1.6455802,1.6735255,1.7018503,1.7306045,1.7598422, 
     *   1.7896223,1.8200099,1.8510770,1.8829044,1.9155830,1.9492166, 
     *   1.9839239,2.0198430,2.0571356,2.0959930,2.1366450,2.1793713, 
     *   2.2245175,2.2725185,2.3239338,2.3795007,2.4402218,2.5075117, 
     *   2.5834658,2.6713916,2.7769943,2.7769943,2.7769943,2.7769943 /
      DATA W / 0.10405134E-04,0.13956560E-04,0.16473259E-04,0.18501623E-
     *   04,0.20238931E-04,0.21780983E-04,0.23182241E-04,0.24476931E-04,
     *   0.25688121E-04,0.26832186E-04,0.27921226E-04,0.28964480E-04, 
     *   0.29969191E-04,0.30941168E-04,0.31885160E-04,0.32805121E-04, 
     *   0.33704388E-04,0.34585827E-04,0.35451919E-04,0.36304851E-04, 
     *   0.37146564E-04,0.37978808E-04,0.38803170E-04,0.39621114E-04, 
     *   0.40433997E-04,0.41243096E-04,0.42049621E-04,0.42854734E-04, 
     *   0.43659562E-04,0.44465208E-04,0.45272764E-04,0.46083321E-04, 
     *   0.46897980E-04,0.47717864E-04,0.48544128E-04,0.49377973E-04, 
     *   0.50220656E-04,0.51073504E-04,0.51937936E-04,0.52815471E-04, 
     *   0.53707761E-04,0.54616606E-04,0.55543990E-04,0.56492112E-04, 
     *   0.57463436E-04,0.58460740E-04,0.59487185E-04,0.60546402E-04, 
     *   0.61642600E-04,0.62780711E-04,0.63966581E-04,0.65207221E-04, 
     *   0.66511165E-04,0.67888959E-04,0.69353880E-04,0.70922996E-04, 
     *   0.72618816E-04,0.74471933E-04,0.76525519E-04,0.78843526E-04, 
     *   0.81526890E-04,0.84749727E-04,0.84749727E-04,0.84749727E-04, 
     *   0.84749727E-04 /
      DATA M / 30788,23052,2053,19346,10646,19427,23975,19049,10949,
     *   19693,29746,26748,2796,23890,29168,31924,16499 /
      DATA M1,M2,I1,J1 / 32767,256,5,17 /
      DATA MDIG / 49 /
      IF (JD.NE.0) GO TO 40
C
C--- EXECUTE THIS SEGMENT EACH CALL
C
   10 I = M(I1)-M(J1)
      IF (I.LT.0) I = I+M1
      M(J1) = I
      I1 = I1-1
      IF (I1.EQ.0) I1 = 17
      J1 = J1-1
      IF (J1.EQ.0) J1 = 17
      J = MOD(I,64)+1
      RDNOR = I*W(J+1)
      IF (((I/M2)/2)*2.EQ.(I/M2)) RDNOR = -RDNOR
      IF (ABS(RDNOR).LE.V(J)) RETURN
C
C--- SLOW PART - AA IS A*F(0) 
C
      X = (ABS(RDNOR)-V(J))/(V(J+1)-V(J))
      Y = RDUNI(0)
      S = X+Y
      IF (S.GT.C2) GO TO 30
      IF (S.LE.C1) RETURN
      IF (Y.GT.C-AA*EXP(-0.5*(B-B*X)**2)) GO TO 30
      IF (EXP(-0.5*V(J+1)**2)+Y*PC/V(J+1).LE.EXP(-0.5*RDNOR**2)) RETURN
C
C--- TAIL PART - 3.855849 IS 0.5*XN**2
C
   20 S = XN-ALOG(RDUNI(0))/XN
      IF (3.855849+ALOG(RDUNI(0))-XN*S.GT.-0.5*S**2) GO TO 20
      RDNOR = SIGN(S,RDNOR)
      RETURN
   30 RDNOR = SIGN(B-B*X,RDNOR)
      RETURN
C
C--- EXECUTE THIS SEGMENT ONLY IF JD<>0 
C
   40 M1 = 2**(MDIG-2)+(2**(MDIG-2)-1)
      M2 = 2**(MDIG/2)
      JSEED = MIN0(IABS(JD),M1)
      IF (MOD(JSEED,2).EQ.0) JSEED = JSEED-1
      K0 = MOD(9069,M2)
      K1 = 9069/M2
      J0 = MOD(JSEED,M2)
      J1 = JSEED/M2 
      DO 50 I = 1, 17
         JSEED = J0*K0
         J1 = MOD(JSEED/M2+J0*K1+J1*K0,M2/2)
         J0 = MOD(JSEED,M2)
   50    M(I) = J0+M2*J1
      J1 = 17
      I1 = 5
      RMAX = 1.0/REAL(M1)
C
C--- SEED UNIFORM (0,1) GENERATOR (JUST A DUMMY CALL)
C
      RDNOR = RDUNI(JD)
      DO 60 I = 1, 65
   60    W(I) = RMAX*V(I)
      GO TO 10
      END 
*RDNOR3
      FUNCTION RDNOR3 (IX,IY,IZ)
C
C-----------------------------------------------------------------------
C   RDNOR3   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: GENERATING STANDARD NORMAL PSEUDO-RANDOM DEVIATES BY THE
C        POLAR METHOD AS DESCRIBED IN THE REFERENCE BELOW.  UNIFORM
C        PSEUDO-RANDOM DEVIATES IN THE INTERVAL (0,1) ARE SUPPLIED
C        BY A WICHMAN-HILL GENERATOR WHICH REQUIRES THREE INTEGER
C        ARGUMENTS IX, IY, AND IZ.  THE NBS CENTRAL COMPUTERS WERE
C        FOUND TO GENERATE NORMAL DEVIATES AT THE FOLLOWING RATES: 
C
C                  CYBER 180/855 ->  22,546 DEVIATES/SECOND 
C                  CYBER 205/622 ->  42,997 DEVIATES/SECOND  .
C
C        BEFORE THE FIRST CALL TO THIS ROUTINE THE INTEGERS IX, IY, AND
C        IZ SHOULD BE SET TO VALUES WITHIN THE RANGE [1,30000].  AFTER
C        THAT THEIR VALUES SHOULD NOT BE CHANGED EXCEPT BY THIS ROUTINE.
C
C   SUBPROGRAMS CALLED: RDUNWH (WICHMAN-HILL UNIFORM GENERATOR)
C
C   CURRENT VERSION COMPLETED MAY 15, 1987
C
C   REFERENCE: KNUTH, DONALD E., 'THE ART OF COMPUTER PROGRAMMING',
C              VOL. 2/SEMINUMERICAL ALGORITHMS, 2ND EDITION, ADDISON- 
C              WESLEY PUBLISHING CO., 1981, PP. 117-118.
C-----------------------------------------------------------------------
C
      LOGICAL LL
      SAVE LL
      DATA LL / .FALSE. /
      IF (LL) THEN
         RDNOR3 = V*Q
      ELSE
10       U = 2.0*RDUNWH(IX,IY,IZ)-1.0
         V = 2.0*RDUNWH(IX,IY,IZ)-1.0
         S = U*U+V*V
         IF (S.GE.1.0) GO TO 10
         Q = SQRT(-2.0*ALOG(S)/S)
         RDNOR3 = U*Q
      ENDIF
      LL = .NOT.LL
      RETURN
      END 
*RDT
      FUNCTION RDT (DF)
C
C-----------------------------------------------------------------------
C   RDT   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C         DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C
C   FOR: GENERATING A RANDOM DEVIATE FROM THE T(DF) DISTRIBUTION.
C        ONE OF THREE METHODS IS USED DEPENDING ON THE VALUE OF THE
C        PARAMETER DF (WHICH DOES NOT HAVE TO BE AN INTEGER): 
C
C             VALUE OF DF             METHOD USED 
C            -------------         ------------------
C              0 < DF < 1          NORMAL/SQRT(CHI-SQUARED/DF)
C                DF = 1            TANGENT TRANSFORMATION (CST)
C                DF > 1            KINDERMAN-MONAHAN-RAMAGE (TIR)
C
C        IF DF <= 0 AN ERROR MESSAGE IS PRINTED AND EXECUTION IS
C        TERMINATED.
C
C        DESCRIPTIONS OF EACH OF THESE ALGORITHMS CAN BE FOUND IN
C        THE REFERENCE GIVEN BELOW.
C
C   SUBPROGRAMS CALLED:  RDUNI (STSPAC) - UNIFORM(0,1) GENERATOR
C                        RDNOR (STSPAC) - NORMAL(0,1) GENERATOR
C                       RDCHI2 (STSPAC) - CHI-SQUARED GENERATOR
C
C   CURRENT VERSION COMPLETED FEBRUARY 28, 1986
C
C   REFERENCE: KINDERMAN, A.J., MONAHAN, J.F., AND RAMAGE, J.G.,
C              "COMPUTER METHODS FOR SAMPLING FROM STUDENT'S T
C              DISTRIBUTION", MATHEMATICS OF COMPUTATION, VOLUME 31,
C              NUMBER 140, OCTOBER 1977, PP. 1009-1018
C-----------------------------------------------------------------------
C
      F(X,A) = (1.0+X*X/A)**(-(A+1.0)/2.0)
      IF (DF.GT.1.0) THEN
C
C
C   KINDERMAN-MONAHAN-RAMAGE ALGORITHM (TIR)
C
C--- STEP 1
C
   10    U = RDUNI(0)
         IF (U.GE.0.23079283) GO TO 20
         RDT = 4.0*U-0.46158566
C
C--- STEP 2
C
         V = RDUNI(0)
         IF (V.LE.1.0-0.5*ABS(RDT)) RETURN
         IF (V.LE.F(RDT,DF)) RETURN
         GO TO 10
C
C--- STEP 3
C
   20    IF (U.GE.0.5) GO TO 40
         S = 4.0*U-1.46158566 
         RDT = SIGN(ABS(S)+0.46158566,S)
         V = RDUNI(0)
C
C--- STEP 4
C
   30    IF (V.LE.1.0-0.5*ABS(RDT)) RETURN
         IF (V.GE.1.2130613/(1.0+RDT*RDT)) GO TO 10
         IF (V.LE.F(RDT,DF)) RETURN
         GO TO 10
C
C--- STEP 5
C
   40    IF (U.GE.0.75) GO TO 50
         S = 8.0*U-5.0
         RDT = 2.0/SIGN(ABS(S)+1.0,S)
         V = RDUNI(0)/(RDT*RDT)
         GO TO 30
C
C--- STEP 6
C
   50    RDT = 2.0/(8.0*U-7.0)
         V = RDUNI(0)
         IF (V.LT.RDT*RDT*F(RDT,DF)) RETURN
         GO TO 10
C
      ELSEIF (DF.EQ.1.0) THEN 
C
C
C   SYNTHETIC TANGENT ALGORITHM (CST)
C
C--- STEP 1
C
   60    U = RDUNI(0)
         V = 2.0*RDUNI(0)-1.0 
C
C--- STEP 2
C
         IF (U*U+V*V.GT.1.0) GO TO 60
         RDT = V/U
         RETURN
C
      ELSEIF (DF.GT.0.0) THEN 
C
C
C   RATIO OF STANDARD NORMAL AND SQUARE ROOT OF
C   CHI-SQUARED DIVIDED BY ITS DEGREES OF FREEDOM 
C
C
         D = SQRT(RDCHI2(DF)/DF)
         RDT = RDNOR(0)/D
      ELSE
         PRINT *,' *** DEGREES OF FREEDOM MUST BE > 0'
         PRINT *,' *** EXECUTION STOPPED IN FUNCTION RDT'
         STOP
C
      ENDIF
      RETURN
      END 
*RDUNI
      FUNCTION RDUNI (JD)
C
C-----------------------------------------------------------------------
C   RDUNI   OBTAINED BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C           DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C           MD 20899 FROM THE SOFTWARE LIBRARY CMLIB AT NBS 
C
C   FOR: GENERATING UNIFORM PSEUDO-RANDOM DEVIATES IN THE INTERVAL
C        [0,1) USING A LAGGED FIBONACCI GENERATOR.  THE INITIAL CALL
C        TO RDUNI SHOULD BE OF THE FORM Z=RDUNI(K) WHERE K IS A
C        POSITIVE INTEGER.  FURTHER CALLS SHOULD BE OF THE FORM
C        Z=RDUNI(0).  MDIG IS A LOWER BOUND ON THE NUMBER OF BINARY
C        DIGITS AVAILABLE FOR REPRESENTING INTEGERS, INCLUDING THE
C        SIGN BIT.  THIS VALUE MUST BE AT LEAST 16.  TO ALLOW THE
C        LONGEST SEQUENCE OF RANDOM NUMBERS WITHOUT CYCLING, SET
C        MDIG TO ITS MAXIMUM VALUE.
C
C        THE NBS CENTRAL COMPUTERS WERE FOUND TO GENERATE DEVIATES
C        AT THE FOLLOWING RATES USING THE INDICATED MAXIMUM VALUES
C        OF MDIG: 
C
C             CYBER 180/855 ->   100,050 DEV/SEC   (MDIG=49)
C             CYBER 205/622 ->   105,471 DEV/SEC   (MDIG=48)
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED FEBRUARY 9, 1987
C
C   REFERENCE: MARSAGLIA, GEORGE, 'COMMENTS ON THE PERFECT UNIFORM
C              RANDOM NUMBER GENERATOR', UNPUBLISHED NOTES, 
C              WASHINGTON STATE UNIVERSITY.
C-----------------------------------------------------------------------
C
      INTEGER M(17) 
      SAVE I,J,M,M1,M2
      DATA M / 30788,23052,2053,19346,10646,19427,23975,19049,10949,
     *   19693,29746,26748,2796,23890,29168,31924,16499 /
      DATA M1,M2,I,J / 32767,256,5,17 / 
      DATA MDIG / 49 /
      IF (JD.EQ.0) GO TO 20
C
C--- EXECUTE THIS SEGMENT ONLY IF JD<>0 
C
      M1 = 2**(MDIG-2)+(2**(MDIG-2)-1)
      M2 = 2**(MDIG/2)
      JSEED = MIN0(IABS(JD),M1)
      IF (MOD(JSEED,2).EQ.0) JSEED = JSEED-1
      K0 = MOD(9069,M2)
      K1 = 9069/M2
      J0 = MOD(JSEED,M2)
      J1 = JSEED/M2 
      DO 10 I = 1, 17
         JSEED = J0*K0
         J1 = MOD(JSEED/M2+J0*K1+J1*K0,M2/2)
         J0 = MOD(JSEED,M2)
         M(I) = J0+M2*J1
   10 CONTINUE
      I = 5
      J = 17
C
C--- EXECUTE THIS SEGMENT EACH CALL
C
   20 K = M(I)-M(J) 
      IF (K.LT.0) K = K+M1
      M(J) = K
      I = I-1
      IF (I.EQ.0) I = 17
      J = J-1
      IF (J.EQ.0) J = 17
      RDUNI = REAL(K)/REAL(M1)
      RETURN
      END 
*RDUNLL
      FUNCTION RDUNLL (ISEED) 
C
C-----------------------------------------------------------------------
C   RDUNLL   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: GENERATING UNIFORM PSEUDO-RANDOM DEVIATES IN THE INTERVAL
C        (0,1) USING A CONGRUENTIAL GENERATOR.  THIS GENERATOR WAS
C        DEVELOPED AT LAWRENCE LIVERMORE LABORATORIES UNDER THE NAME
C        'LLRANDOMI' AND IS CURRENTLY USED IN THE IMSL LIBRARY.  IT
C        HAS A CYCLE LENGTH OF 2**31-2 = 2,147,483,646.  BEFORE THE
C        FIRST CALL TO THIS ROUTINE ISEED SHOULD BE ASSIGNED A VALUE
C        IN THE RANGE [1,2147483646].  THEREAFTER THE VALUE OF ISEED
C        SHOULD NOT BE CHANGED.
C
C        THE NBS CENTRAL COMPUTERS WERE FOUND TO GENERATE DEVIATES
C        AT THE FOLLOWING RATES USING THIS ROUTINE: 
C
C                 CYBER 180/855 ->   125,707 DEV/SEC
C                 CYBER 205/622 ->   260,373 DEV/SEC
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED FEBRUARY 10, 1987
C-----------------------------------------------------------------------
C
      DATA M,RM / 2147483647,2147483647.0 /
      ISEED = MOD(16807*ISEED,M)
      RDUNLL = REAL(ISEED)/RM 
      RETURN
      END 
*RDUNWH
      FUNCTION RDUNWH (IX,IY,IZ)
C
C-----------------------------------------------------------------------
C   RDUNWH   COPIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND 20899 FROM THE REFERENCE BELOW WITH SLIGHT
C            MODIFICATIONS
C
C   FOR: GENERATING UNIFORM PSEUDO-RANDOM DEVIATES IN THE INTERVAL
C        [0,1) BY THE WICHMAN-HILL METHOD.  AN EXACT VALUE OF ZERO
C        IS POSSIBLE BUT NOT LIKELY.  THE CYCLE LENGTH EXCEEDS 6.95E12.
C        TWO METHODS, WHICH GENERATE IDENTICAL DEVIATES, ARE AVAILABLE
C        FOR USE DEPENDING ON THE LIMITS OF THE INTEGER ARITHMETIC.
C        THE FIRST METHOD MAY BE USED ON COMPUTERS WITH 16-BIT WORDS. 
C        THE METHOD NOT BEING USED CAN BE 'COMMENTED' OUT.  USING THE 
C        SECOND METHOD, THE NBS CENTRAL COMPUTERS WERE FOUND TO
C        GENERATE DEVIATES AT THE FOLLOWING RATES: 
C
C                 CYBER 180/855 ->   48,995 DEVIATES/SECOND 
C                 CYBER 205/622 ->  126,564 DEVIATES/SECOND .
C
C        BEFORE THE FIRST CALL TO THIS ROUTINE THE INTEGERS IX, IY, AND
C        IZ SHOULD BE SET TO VALUES WITHIN THE RANGE [1,30000].  AFTER
C        THAT THEIR VALUES SHOULD NOT BE CHANGED EXCEPT BY THIS ROUTINE.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED MAY 15, 1987
C
C   REFERENCE: GRIFFITHS, P. AND HILL, I.D. (EDS.), 'APPLIED STATISTICS
C              ALGORITHMS', THE ROYAL STATISTICAL SOCIETY/ELLIS HARWOOD
C              LIMITED, 1985, ALGORITHM AS 183, PP. 238-242.
C-----------------------------------------------------------------------
C
C--- THE NEXT SIX LINES REQUIRE INTEGER ARITHMETIC UP TO 30,323
C
C     IX = 171*MOD(IX,177)-2*(IX/177)
C     IY = 172*MOD(IY,176)-35*(IY/176)
C     IZ = 170*MOD(IZ,178)-63*(IZ/178)
C     IF (IX.LT.0) IX = IX+30269
C     IF (IY.LT.0) IY = IY+30307
C     IF (IZ.LT.0) IZ = IZ+30323
C
C--- THE NEXT THREE LINES REQUIRE INTEGER ARITHMETIC UP TO 5,212,632
C
      IX = MOD(171*IX,30269)
      IY = MOD(172*IY,30307)
      IZ = MOD(170*IZ,30323)
      SUM = REAL(IX)/30269.0+REAL(IY)/30307.0+REAL(IZ)/30323.0
      RDUNWH = AMOD(SUM,1.0)
      RETURN
      END 
*RDCONS
      SUBROUTINE RDCONS (P,LDP,N,H,PSI,IFLAG)
C
C-----------------------------------------------------------------------
C   RDCONS   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE
C        SURFACE OF A CONE WITH HEIGHT H AND HALF-ANGLE PSI.  NO
C        POINTS WILL BE GENERATED UNLESS H>0 AND 0<PSI<PI/2.  THE CONE
C        TIP IS AT THE ORIGIN, AND ITS AXIS COINCIDES WITH THE Z AXIS.
C        THE BASE LIES IN THE PLANE Z=H.  ON OUTPUT, EACH OF THE N
C        COLUMNS OF MATRIX P CONTAINS THE X, Y, AND Z COORDINATES OF
C        THE GENERATED POINTS.  THE LEADING DIMENSION OF P (LDP) MUST 
C        BE 3 OR MORE IN THE CALLING PROGRAM.
C
C        IF IFLAG=0 ON INPUT, THEN POINTS WILL BE GENERATED ONLY ON THE
C        CURVED SURFACE.  IF IFLAG<>0 ON INPUT, THEN POINTS WILL BE
C        GENERATED ON BOTH THE CURVED SURFACE AND THE BASE. 
C
C        POINTS ARE GENERATED ON THE CURVED SURFACE BY TRANSFORMATIONS,
C        AND ON THE BASE (IF REQUESTED) BY A REJECTION PROCEDURE
C        WHICH REJECTS 1 - PI/4 = 0.2146 OF THE POINTS.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED MAY 26, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED ON THE CONE SURFACE.  THE LEADING
C              DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION,
C              DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N.
C              [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C        * H = HEIGHT OF THE CONE (>0) [REAL]
C
C      * PSI = HALF-ANGLE OF THE CONE (>0) [REAL] 
C
C    * IFLAG = INDICATOR FOR WHERE TO GENERATE POINTS ON INPUT
C              (0 -> CURVED SURFACE ONLY; OTHERWISE -> CURVED SURFACE 
C              AND BASE) [INTEGER]
C
C              ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> H<=0 OR PSI<=0 OR PSI>=PI/2
C              2 -> N<1.
C              3 -> LDP<3.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
      TWOPI = 8.0*ATAN(1.0)
      HALFPI = 0.25*TWOPI
      IF (H.LE.0.0.OR.PSI.LE.0.0.OR.PSI.GE.HALFPI) THEN
         IFLAG = 1
      ELSEIF (N.LT.1) THEN
         IFLAG = 2
      ELSEIF (LDP.LT.3) THEN
         IFLAG = 3
      ELSE
         TANPSI = TAN(PSI)
         UC = 1.0/(1.0+SIN(PSI))
         RR = H*TANPSI
         IF (IFLAG.EQ.0) THEN 
            DO 10 I = 1, N
               Z = H*SQRT(RDUNI(0))
               R = Z*TANPSI
               PHI = TWOPI*RDUNI(0)
               P(1,I) = R*COS(PHI)
               P(2,I) = R*SIN(PHI)
               P(3,I) = Z
   10       CONTINUE
         ELSE
            DO 30 I = 1, N
               U1 = RDUNI(0)
               IF (U1.LE.UC) THEN
                  Z = H*SQRT(U1/UC)
                  R = Z*TANPSI
                  PHI = TWOPI*RDUNI(0)
                  P(1,I) = R*COS(PHI)
                  P(2,I) = R*SIN(PHI)
                  P(3,I) = Z
               ELSE 
   20             U1 = 2.0*RDUNI(0)-1.0 
                  U2 = 2.0*RDUNI(0)-1.0 
                  Q = U1**2+U2**2
                  IF (Q.GT.1.0) GO TO 20
                  P(1,I) = RR*U1
                  P(2,I) = RR*U2
                  P(3,I) = H
               ENDIF
   30       CONTINUE
         ENDIF
         IFLAG = 0
      ENDIF
      RETURN
C
      END 
*RDCONV
      SUBROUTINE RDCONV (P,LDP,N,H,PSI,IFLAG)
C
C-----------------------------------------------------------------------
C   RDCONV   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS WITHIN THE
C        CONE WITH HEIGHT H AND HALF-ANGLE PSI.  NO POINTS WILL BE
C        GENERATED UNLESS H>0 AND 0<PSI<PI/2.  THE CONE TIP IS AT THE 
C        ORIGIN, AND ITS AXIS COINCIDES WITH THE Z AXIS.  THE BASE
C        LIES IN THE PLANE Z=H.  ON OUTPUT, EACH OF THE N COLUMNS OF
C        MATRIX P CONTAINS THE X, Y, AND Z COORDINATES OF THE GENERATED
C        POINTS.  THE LEADING DIMENSION OF P (LDP) MUST BE 3 OR MORE
C        IN THE CALLING PROGRAM.
C
C        POINTS ARE GENERATED BY A REJECTION PROCEDURE WHEREBY 0.2146 
C        = 1-PI/4 OF THE POINTS ARE REJECTED.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED MAY 30, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED INSIDE THE CONE.  THE LEADING
C              DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION,
C              DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N.
C              [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C        * H = HEIGHT OF THE CONE (>0) [REAL]
C
C      * PSI = HALF-ANGLE OF THE CONE (>0) [REAL] 
C
C      IFLAG = ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> H<=0 OR PSI<=0 OR H<PSI.
C              2 -> N<1.
C              3 -> LDP<3.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
      TWOPI = 8.0*ATAN(1.0)
      HALFPI = 0.25*TWOPI
      IF (H.LE.0.0.OR.PSI.LE.0.0.OR.PSI.GT.HALFPI) THEN
         IFLAG = 1
      ELSEIF (N.LT.1) THEN
         IFLAG = 2
      ELSEIF (LDP.LT.3) THEN
         IFLAG = 3
      ELSE
         IFLAG = 0
         CUBERT = 1.0/3.0
         TANPSI = TAN(PSI)
         DO 10 I = 1, N
            Z = H*RDUNI(0)**CUBERT
            R = Z*TANPSI*SQRT(RDUNI(0)) 
            PHI = TWOPI*RDUNI(0)
            P(1,I) = R*COS(PHI)
            P(2,I) = R*SIN(PHI)
            P(3,I) = Z
   10    CONTINUE
      ENDIF
      RETURN
C
      END 
*RDCYLS
      SUBROUTINE RDCYLS (P,LDP,N,H,R,IFLAG)
C
C-----------------------------------------------------------------------
C   RDCYLS   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE
C        SURFACE OF THE CYLINDER WITH HEIGHT H AND RADIUS R.  NO
C        POINTS WILL BE GENERATED UNLESS H>=0 AND R>=0.  THE AXIS
C        OF THE CYLINDER COINCIDES WITH THE POSITIVE Z AXIS, AND THE
C        ENDS OF THE CYLINDER ARE FORMED BY THE PLANES Z=0 AND Z=H.
C        ON OUTPUT, EACH OF THE N COLUMNS OF MATRIX P CONTAINS THE X, 
C        Y, AND Z COORDINATES OF THE GENERATED POINTS.  THE LEADING
C        DIMENSION OF P (LDP) MUST BE 3 OR MORE IN THE CALLING PROGRAM.
C
C        AT THE USERS OPTION, POINTS MAY BE GENERATED ON THE CURVED
C        SURFACE ONLY OR ON BOTH THE CURVED AND FLAT SURFACES.  POINTS
C        ARE GENERATED ON THE CURVED SURFACE BY A TRANSFORMATION, AND 
C        POINTS ARE GENERATED ON THE FLAT SURFACES BY A REJECTION
C        PROCEDURE WHEREBY THE PROPORTION OF POINTS REJECTED IS 0.2146
C        = 1-PI/4.
C
C        NOTE THAT POINTS CAN BE GENERATED ON A STRAIGHT LINE BY
C        SETTING H>0 AND R=0.  POINTS CAN BE GENERATED ON A CIRCLE BY 
C        SETTING H=0 AND R>0. 
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED MAY 30, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED ON THE CYLINDER SURFACE.  THE LEADING 
C              DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION,
C              DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N.
C              [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C        * H = HEIGHT OF THE CYLINDER (>=0) [REAL]
C
C        * R = RADIUS OF THE CYLINDER (>=0) [REAL]
C
C    * IFLAG = INDICATOR FOR WHERE TO GENERATE POINTS ON INPUT
C              (0 -> CURVED SURFACE ONLY; OTHERWISE -> CURVED SURFACE 
C              AND BASE) [INTEGER]
C
C              ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> H<0.
C              2 -> R<0.
C              3 -> LDP<3.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
      TWOPI = 8.0*ATAN(1.0)
      IF (H.LT.0.0) THEN
         IFLAG = 1
      ELSEIF (R.LT.0.0) THEN
         IFLAG = 2
      ELSE
         IF (IFLAG.EQ.0) THEN 
C
C--- GENERATE POINTS ON THE CURVED SURFACE ONLY
C
            DO 10 I = 1, N
               P(3,I) = H*RDUNI(0)
               U = TWOPI*RDUNI(0)
               P(1,I) = R*COS(U)
               P(2,I) = R*SIN(U)
   10       CONTINUE
         ELSE
C
C--- GENERATE POINTS ON BOTH THE CURVED AND FLAT SURFACES
C
            UC1 = H/(H+R)
            UC2 = 0.5*(1.0+UC1)
            DO 30 I = 1, N
               U = RDUNI(0)
               IF (U.LE.UC1) THEN
                  P(3,I) = H*U/UC1
                  U = TWOPI*RDUNI(0)
                  P(1,I) = R*COS(U)
                  P(2,I) = R*SIN(U)
               ELSE 
   20             U1 = 2.0*RDUNI(0)-1.0 
                  U2 = 2.0*RDUNI(0)-1.0 
                  Q = U1*U1+U2*U2
                  IF (Q.GT.1.0) GO TO 20
                  P(1,I) = R*U1
                  P(2,I) = R*U2
                  IF (U.LE.UC2) THEN
                     P(3,I) = 0.0
                  ELSE
                     P(3,I) = H
                  ENDIF
               ENDIF
   30       CONTINUE
         ENDIF
         IFLAG = 0
      ENDIF
      RETURN
C
      END 
*RDCYLV
      SUBROUTINE RDCYLV (P,LDP,N,H,R,IFLAG)
C
C-----------------------------------------------------------------------
C   RDCYLV   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS WITHIN THE
C        CYLINDER WITH HEIGHT H AND RADIUS R.  NO POINTS WILL BE
C        GENERATED UNLESS H>=0 AND R>=0.  THE AXIS OF THE CYLINDER
C        COINCIDES WITH THE POSITIVE Z AXIS, AND THE ENDS OF THE
C        CYLINDER ARE FORMED BY THE PLANES Z=0 AND Z=H.  ON OUTPUT,
C        EACH OF THE N COLUMNS OF MATRIX P CONTAINS THE X, Y, AND Z
C        COORDINATES OF THE GENERATED POINTS.  THE LEADING DIMENSION
C        OF P (LDP) MUST BE 3 OR MORE IN THE CALLING PROGRAM.
C
C        POINTS ARE GENERATED BY A REJECTION PROCEDURE WHEREBY THE
C        PROPORTION OF POINTS REJECTED IS 0.2146 = 1-PI/4.
C
C        NOTE THAT POINTS CAN BE GENERATED (RATHER INEFFICIENTLY) ON
C        A STRAIGHT LINE BY SETTING H>0 AND R=0.  POINTS CAN BE
C        GENERATED ON A (TWO-DIMENSIONAL) DISK BY SETTING H=0 AND R>0.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED MAY 30, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED INSIDE THE CYLINDER.  THE LEADING
C              DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION,
C              DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N.
C              [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C        * H = HEIGHT OF THE CYLINDER (>=0) [REAL]
C
C        * R = RADIUS OF THE CYLINDER (>=0) [REAL]
C
C      IFLAG = ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> H<0.
C              2 -> R<0.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
      IF (H.LT.0.0) THEN
         IFLAG = 1
      ELSEIF (R.LT.0.0) THEN
         IFLAG = 2
      ELSE
         IFLAG = 0
         DO 20 I = 1, N
            P(3,I) = H*RDUNI(0)
   10       U = 2.0*RDUNI(0)-1.0
            V = 2.0*RDUNI(0)-1.0
            Q = U*U+V*V
            IF (Q.GT.1.0) GO TO 10
            P(1,I) = R*U
            P(2,I) = R*V
   20    CONTINUE
      ENDIF
      RETURN
C
      END 
*RDELLS
      SUBROUTINE RDELLS (P,LDP,N,M,A,IFLAG)
C
C-----------------------------------------------------------------------
C   RDELLS   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE
C        SURFACE OF THE M-DIMENSIONAL ELLIPSOID DEFINED BY
C
C                            M
C                           SUM (X /A )**2 = 1
C                           J=1   J  J
C
C        WHERE EACH A(J)>0.  THE CENTER OF THE ELLIPSOID IS AT THE
C        ORIGIN.  NOTE THAT M STANDARD NORMAL DEVIATES AND 1 STANDARD 
C        UNIFORM DEVIATE ARE REQUIRED TO GENERATE EACH POINT.  IF THE 
C        A(J) ARE NOT ALL EQUAL, THERE IS A POSITIVE PROBABILITY THAT 
C        THE POINT WILL BE REJECTED.
C
C   SUBPROGRAMS CALLED: RDNOR (STSPAC) - N(0,1) RANDOM DEVIATES
C                       RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED JANUARY 8, 1990
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE M BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED ON THE SURFACE OF THE ELLIPSOID.  THE 
C              LEADING DIMENSION OF P IS LDP (>=M) AND THE SECOND
C              DIMENSION, DEFINED IN THE CALLING PROGRAM, MUST BE AT
C              LEAST N. [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=M) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C        * M = DIMENSION OF THE ELLIPSOID (>0) [INTEGER]
C
C     * A(*) = VECTOR (LENGTH M) OF ELLIPSOID SEMIAXES (>0) [REAL]
C
C      IFLAG = ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> N<1 OR M<1.
C              2 -> M>LDP.
C              3 -> AT LEAST ONE A(J)<=0.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION A(*),P(LDP,*) 
      IF (MIN0(N,M).LE.0) THEN
         IFLAG = 1
      ELSEIF (M.GT.LDP) THEN
         IFLAG = 2
      ELSE
         AMIN = A(1)
         DO 10 J = 2, M
            AMIN = AMIN1(AMIN,A(J))
   10    CONTINUE
         IF (AMIN.LE.0.0) THEN
            IFLAG = 3
         ELSE
            IFLAG = 0
            DO 50 I = 1, N
   20          Q = 0.0
               R = 0.0
               DO 30 J = 1, M 
                  P(J,I) = RDNOR(0)
                  Q = Q+P(J,I)**2
                  R = R+(P(J,I)/A(J))**2
   30          CONTINUE
               Q = SQRT(Q)
               R = SQRT(R)
               U1 = RDUNI(0)
               IF (AMIN*R.LT.Q*U1) GO TO 20
               DO 40 J = 1, M 
                  P(J,I) = A(J)*P(J,I)/Q
   40          CONTINUE
   50       CONTINUE
         ENDIF
      ENDIF
      RETURN
C
      END 
*RDELLV
      SUBROUTINE RDELLV (P,LDP,N,M,A,IFLAG)
C
C-----------------------------------------------------------------------
C   RDELLV   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS INSIDE THE
C        M-DIMENSIONAL ELLIPSOID DEFINED BY
C
C                            M
C                           SUM (X /A )**2 = 1
C                           J=1   J  J
C
C        WHERE EACH A(J)>0.  THE CENTER OF THE ELLIPSOID IS AT THE
C        ORIGIN.  NOTE THAT M STANDARD NORMAL DEVIATES AND ONE
C        UNIFORM DEVIATE ARE REQUIRED TO GENERATE EACH POINT.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C                       RDNOR (STSPAC) - N(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED MARCH 16, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE M BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED INSIDE THE ELLIPSOID.  THE
C              LEADING DIMENSION OF P IS LDP (>=M) AND THE SECOND
C              DIMENSION, DEFINED IN THE CALLING PROGRAM, MUST BE AT
C              LEAST N. [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=M) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C        * M = DIMENSION OF THE ELLIPSOID (>0) [INTEGER]
C
C     * A(*) = VECTOR (LENGTH M) OF ELLIPSOID SEMIAXES (>0) [REAL]
C
C      IFLAG = ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> N<1 OR M<1.
C              2 -> M>LDP.
C              3 -> AT LEAST ONE A(J)<=0.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION A(*),P(LDP,*) 
      IF (MIN0(N,M).LE.0) THEN
         IFLAG = 1
      ELSEIF (M.GT.LDP) THEN
         IFLAG = 2
      ELSE
         AMIN = A(1)
         DO 10 J = 2, M
            AMIN = AMIN1(AMIN,A(J))
   10    CONTINUE
         IF (AMIN.LE.0.0) THEN
            IFLAG = 3
         ELSE
            IFLAG = 0
            POW = 1.0/REAL(M) 
            DO 40 I = 1, N
               Q = 0.0
               DO 20 J = 1, M 
                  P(J,I) = RDNOR(0)
                  Q = Q+P(J,I)**2
   20          CONTINUE
               R = RDUNI(0)**POW
               FAC = R/SQRT(Q)
               DO 30 J = 1, M 
                  P(J,I) = FAC*A(J)*P(J,I)
   30          CONTINUE
   40       CONTINUE
         ENDIF
      ENDIF
      RETURN
C
      END 
*RDHLX
      SUBROUTINE RDHLX (P,LDP,N,R,D,H,PHI,KAPPA,IFLAG)
C
C-----------------------------------------------------------------------
C   RDHLX   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C           DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, 
C           GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON A HELIX
C        IN THREE-DIMENSIONAL SPACE.  THE HELIX MAY BE 'RIGHT-HANDED' 
C        OR 'LEFT-HANDED' AT THE USERS OPTION.  THE AXIS OF THE HELIX 
C        COINCIDES WITH THE Z-AXIS BETWEEN Z=0 AND Z=H.  THE HELIX
C        RADIUS IS R, AND THE SEPARATION BETWEEN ADJACENT POINTS WITH 
C        THE SAME (X,Y) COORDINATES IS D.  THE INITIAL HELIX ANGLE
C        (AT Z=0) IS PHI.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED JANUARY 9, 1990
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED ON THE HELIX.  THE LEADING
C              DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION,
C              DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N.
C              [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED (N>0) [INTEGER]
C
C        * R = THE RADIUS OF THE HELIX (R>=0) [REAL]
C
C        * D = THE Z-DISTANCE BETWEEN ADJACENT POINTS WHICH HAVE THE
C              SAME (X,Y) COORDINATES (D>0) [REAL]
C
C        * H = THE HEIGHT OF THE HELIX (H>0) [REAL]
C
C      * PHI = INITIAL HELIX ANGLE WHEN Z=0 [REAL]
C
C    * KAPPA =  1 FOR A RIGHT-HANDED HELIX,
C              -1 FOR A LEFT-HANDED HELIX [INTEGER]
C
C      IFLAG = ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> N<1 OR D<=0 OR H<=0.
C              2 -> R<0.
C              3 -> ABS(KAPPA)<>1.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
      IF (N.LT.1.OR.D.LE.0.0.OR.H.LE.0.0) THEN
         IFLAG = 1
      ELSEIF (R.LT.0.0) THEN
         IFLAG = 2
      ELSEIF (IABS(KAPPA).NE.1) THEN
         IFLAG = 3
      ELSE
         IFLAG = 0
         CON = 8.0*ATAN(1.0)*REAL(KAPPA)/D
         DO 10 I=1,N
            U1 = RDUNI(0)
            Z = H*U1
            THETA = PHI+CON*Z 
            P(1,I) = R*COS(THETA)
            P(2,I) = R*SIN(THETA)
            P(3,I) = Z
10       CONTINUE
      ENDIF
      RETURN
C
      END 
*RDRECS
      SUBROUTINE RDRECS (P,LDP,N,A,B,C,IFLAG)
C
C-----------------------------------------------------------------------
C   RDRECS   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE
C        SURFACE OF A 3-DIMENSIONAL RECTANGULAR SOLID WITH SIDES OF
C        LENGTH A, B, AND C.  NO POINTS WILL BE GENERATED UNLESS
C        MIN(A,B,C)>=0.  ONE CORNER OF THE RECTANGULAR SOLID HAS
C        COORDINATES (0,0,0) AND THE OPPOSITE CORNER HAS COORDINATES
C        (A,B,C).  ON OUTPUT, EACH OF THE N COLUMNS OF MATRIX P
C        CONTAINS THE X, Y, AND Z COORDINATES OF THE GENERATED POINTS.
C        THE LEADING DIMENSION OF P (LDP) MUST BE 3 OR MORE IN THE
C        CALLING PROGRAM.
C
C        NOTE THAT POINTS MAY BE GENERATED IN A PLANE BY SETTING ONE
C        OF (A,B,C) TO ZERO, AND POINTS MAY BE GENERATED ON A LINE BY 
C        SETTING TWO OF (A,B,C,) TO ZERO.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED MAY 31, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED ON THE RECTANGULAR SOLID SURFACE.  THE
C              LEADING DIMENSION OF P IS LDP (>=3) AND THE SECOND
C              DIMENSION, DEFINED IN THE CALLING PROGRAM, MUST BE AT
C              LEAST N. [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C        * A = LENGTH OF RECTANGULAR SOLID IN X DIRECTION (>=0) [REAL]
C
C        * B = LENGTH OF RECTANGULAR SOLID IN Y DIRECTION (>=0) [REAL]
C
C        * C = LENGTH OF RECTANGULAR SOLID IN Z DIRECTION (>=0) [REAL]
C
C      IFLAG = ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> MIN(A,B,C)<0.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
      IF (AMIN1(A,B,C).LT.0.0) THEN
         IFLAG = 1
      ELSE
         SA = A*(B+C)+B*C
         C1 = A*B/SA
         C2 = A*(B+C)/SA
C
C--- GENERATE POINTS ON THE SURFACE OF THE RECTANGULAR SOLID
C
         DO 10 I = 1, N
            U = RDUNI(0)
            V = RDUNI(0)
            IF (U.LT.C1) THEN 
               P(1,I) = A*U/C1
               IF (V.LE.0.5) THEN
                  P(2,I) = 2.0*B*V
                  P(3,I) = 0.0
               ELSE 
                  P(2,I) = 2.0*B*(V-0.5)
                  P(3,I) = C
               ENDIF
            ELSEIF (U.LT.C2) THEN
               P(1,I) = A*(U-C1)/(C2-C1)
               IF (V.LE.0.5) THEN
                  P(3,I) = 2.0*C*V
                  P(2,I) = 0.0
               ELSE 
                  P(3,I) = 2.0*C*(V-0.5)
                  P(2,I) = B
               ENDIF
            ELSE
               P(2,I) = B*(U-C2)/(1.0-C2)
               IF (V.LE.0.5) THEN
                  P(3,I) = 2.0*C*V
                  P(1,I) = 0.0
               ELSE 
                  P(3,I) = 2.0*C*(V-0.5)
                  P(1,I) = A
               ENDIF
            ENDIF
   10    CONTINUE
         IFLAG = 0
      ENDIF
      RETURN
C
      END 
*RDRECV
      SUBROUTINE RDRECV (P,LDP,N,A,B,C,IFLAG)
C
C-----------------------------------------------------------------------
C   RDRECV   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS INSIDE A
C        3-DIMENSIONAL RECTANGULAR SOLID WITH SIDES OF LENGTH A, B,
C        AND C.  NO POINTS WILL BE GENERATED UNLESS MIN(A,B,C)>=0.  ONE
C        CORNER OF THE RECTANGULAR SOLID HAS COORDINATES (0,0,0) AND
C        THE OPPOSITE CORNER HAS COORDINATES (A,B,C).  ON OUTPUT,
C        EACH OF THE N COLUMNS OF MATRIX P CONTAINS THE X, Y, AND Z
C        COORDINATES OF THE GENERATED POINTS.  THE LEADING DIMENSION
C        OF P (LDP) MUST BE 3 OR MORE IN THE CALLING PROGRAM.
C
C        NOTE THAT POINTS MAY BE GENERATED IN A PLANE BY SETTING ONE
C        OF (A,B,C) TO ZERO, AND POINTS MAY BE GENERATED ON A LINE BY 
C        SETTING TWO OF (A,B,C,) TO ZERO.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED MAY 31, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED INSIDE THE RECTANGULAR SOLID.  THE
C              LEADING DIMENSION OF P IS LDP (>=3) AND THE SECOND
C              DIMENSION, DEFINED IN THE CALLING PROGRAM, MUST BE AT
C              LEAST N. [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C        * A = LENGTH OF RECTANGULAR SOLID IN X DIRECTION (>=0) [REAL]
C
C        * B = LENGTH OF RECTANGULAR SOLID IN Y DIRECTION (>=0) [REAL]
C
C        * C = LENGTH OF RECTANGULAR SOLID IN Z DIRECTION (>=0) [REAL]
C
C      IFLAG = ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> MIN(A,B,C)<0.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
      IF (AMIN1(A,B,C).LT.0.0) THEN
         IFLAG = 1
      ELSE
C
C--- GENERATE POINTS INSIDE A RECTANGULAR SOLID
C
         DO 10 I = 1, N
            P(1,I) = A*RDUNI(0)
            P(2,I) = B*RDUNI(0)
            P(3,I) = C*RDUNI(0)
   10    CONTINUE
         IFLAG = 0
      ENDIF
      RETURN
C
      END 
*RDSPSH
      SUBROUTINE RDSPSH (P,LDP,L,N,R1,R2,IFLAG)
C
C-----------------------------------------------------------------------
C   RDSPSH   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING N RANDOM POINTS 'UNIFORMLY' SPACED IN A REGION
C        BETWEEN THE SURFACES OF TWO SPHERES OF RADIUS R1 AND R2.
C        THE ONLY REQUIREMENTS ARE THAT N>0, L>0, R1>=0, AND R2>=0.
C        NOTE THAT SETTING R1=0 GENERATES POINTS WITHIN A SPHERE OF
C        RADIUS R2. 
C
C        THE COMPUTATIONAL PROCEDURE INVOLVES FIRST GENERATING POINTS 
C        UNIFORMLY ON THE SURFACE OF AN L-DIMENSIONAL UNIT SPHERE BY
C        THE METHOD DESCRIBED IN THE REFERENCE BELOW.  THEN THE POINTS
C        ARE SCALED BY THE RANDOM VARIABLE R SO THAT THEIR LENGTHS ARE
C        IN THE INTERVAL [R1,R2].  WHEN R1=0 AND R2=1, R IS DISTRIBUTED
C        AS THE LTH ROOT OF A UNIFORM(0,1) RANDOM VARIABLE. 
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C                       RDNOR (STSPAC) - N(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED FEBRUARY 9, 1989
C
C   REFERENCE: MULLER, M.E., 'A NOTE ON A METHOD FOR GENERATING POINTS
C              UNIFORMLY ON N-DIMENSIONAL SPHERES', COMMUNICATIONS OF 
C              THE ACM, 2, PP. 19-20, 1959.
C-----------------------------------------------------------------------
C   DEFINITIONS OF PASSED PARAMETERS: 
C
C  P(LDP,*) = MATRIX (SIZE L BY N) OF POINT COORDINATES.  ON OUTPUT
C             THE JTH COLUMN CONTAINS THE L COORDINATES OF THE JTH
C             POINT (J=1,2,...,N) [REAL]
C
C     * LDP = LEADING DIMENSION OF MATRIX P AS DEFINED IN THE CALLING 
C             PROGRAM (LDP>=L) [INTEGER]
C
C       * L = DIMENSION OF THE POINTS AND SPHERES (L<=LDP) [INTEGER]
C
C       * N = NUMBER OF RANDOM POINTS TO BE GENERATED.  THE SECOND
C             DIMENSION OF P MUST BE >=N IN CALLING PROGRAM [INTEGER] 
C
C      * R1 = THE RADIUS OF THE FIRST SPHERE (R1>=0) [REAL] 
C
C      * R2 = THE RADIUS OF THE SECOND SPHERE (R2>=0) [REAL]
C
C     IFLAG = ERROR INDICATOR [INTEGER]   ON OUTPUT: 
C             0 -> NO ERRORS DETECTED, N POINTS GENERATED.
C             1 -> EITHER L OR N IS NONPOSITIVE.
C             2 -> L>LDP.
C             3 -> EITHER R1<0 OR R2<0. 
C
C  * INDICATES PARAMETERS REQUIRING INPUT VALUES
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
C
C--- CHECK FOR ILLEGAL VALUES OF INPUT PARAMETERS AND SET ERROR FLAG
C
      IF (MIN0(L,N).LT.1) THEN
         IFLAG = 1
      ELSEIF (L.GT.LDP) THEN
         IFLAG = 2
      ELSEIF (AMIN1(R1,R2).LT.0.0) THEN 
         IFLAG = 3
      ELSE
         IFLAG = 0
C
C--- COMPUTE THREE CONSTANTS
C
         CON1 = R1**L
         CON2 = R2**L
         POWER = 1.0/REAL(L)
C
C--- BEGIN LOOP FOR POINTS
C
         DO 30 J = 1, N
C
C--- COMPUTE ONE UNIFORM DEVIATE (U) FOR JTH POINT, THEN COMPUTE A
C--- DEVIATE R IN THE INTERVAL [R1,R2]
C
            U = RDUNI(0)
            R = (U*CON2+(1.0-U)*CON1)**POWER
C
C--- BEGIN LOOP FOR COORDINATES OF JTH POINT
C
            S = 0.0 
            DO 10 I = 1, L
C
C--- COMPUTE ONE NORMAL DEVIATE FOR ITH COORDINATE OF JTH POINT.  S IS
C--- THE SQUARED DISTANCE OF THE JTH POINT FROM THE ORIGIN. 
C
               Z = RDNOR(0)
               S = S+Z*Z
               P(I,J) = Z
   10       CONTINUE
C
C--- COMPUTE 'COMBINED SCALE FACTOR', THEN SCALE EACH COORDINATE OF THE
C--- JTH POINT SO THAT THE POINT LIES IN THE REGION BETWEEN SPHERES OF
C--- RADIUS R1 AND R2.
C
            SCALE = R/SQRT(S) 
            DO 20 I = 1, L
               P(I,J) = SCALE*P(I,J)
   20       CONTINUE
   30    CONTINUE
      ENDIF
      RETURN
C
      END 
*RDTORS
      SUBROUTINE RDTORS (P,LDP,N,R1,R2,IFLAG)
C
C-----------------------------------------------------------------------
C   RDTORS   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE
C        SURFACE OF THE TORUS WHOSE MAJOR RADIUS IS R1 AND MINOR RADIUS
C        IS R2.  NO POINTS WILL BE GENERATED UNLESS R1>0, R2>0, R1>=R2,
C        AND N>0.  THE TORUS LIES IN THE (X,Y) PLANE CENTERED AT THE
C        ORIGIN.  ON OUTPUT, EACH OF THE N COLUMNS OF MATRIX P CONTAINS
C        THE X, Y, AND Z COORDINATES OF THE GENERATED POINTS.  THE
C        LEADING DIMENSION OF P (LDP) MUST BE 3 OR MORE IN THE CALLING
C        PROGRAM.
C
C        THE ANGLE THETA (ASSOCIATED WITH R2) IS GENERATED BY A
C        REJECTION METHOD REQUIRING AN AVERAGE OF 2+2*(R2/R1) UNIFORM 
C        RANDOM DEVIATES.  THE ANGLE PHI (ASSOCIATED WITH R1) REQUIRES
C        A SINGLE UNIFORM RANDOM DEVIATE.  THEREFORE, EACH GENERATED
C        POINT REQUIRES AN AVERAGE OF 3+2*(R2/R1) UNIFORM DEVIATES.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED MARCH 23, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED ON THE TORUS SURFACE.  THE LEADING
C              DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION,
C              DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N.
C              [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C       * R1 = MAJOR RADIUS OF THE TORUS (>0) [REAL]
C
C       * R2 = MINOR RADIUS OF THE TORUS (>0) [REAL]
C
C      IFLAG = ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> R1<=0 OR R2<=0 OR R1<R2.
C              2 -> N<1.
C              3 -> LDP<3.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
      IF (AMIN1(R1,R2).LE.0.0.OR.R1.LT.R2) THEN
         IFLAG = 1
      ELSEIF (N.LT.1) THEN
         IFLAG = 2
      ELSEIF (LDP.LT.3) THEN
         IFLAG = 3
      ELSE
         IFLAG = 0
         TWOPI = 8.0*ATAN(1.0)
         RSUM = R1+R2
         DO 20 I = 1, N
   10       THETA = TWOPI*RDUNI(0)
            U = RSUM*RDUNI(0) 
            R = R1+R2*COS(THETA)
            IF (U.GT.R) GO TO 10
            PHI = TWOPI*RDUNI(0)
            P(1,I) = R*COS(PHI)
            P(2,I) = R*SIN(PHI)
            P(3,I) = R2*SIN(THETA)
   20    CONTINUE
      ENDIF
      RETURN
C
      END 
*RDTORV
      SUBROUTINE RDTORV (P,LDP,N,R1,R2,IFLAG)
C
C-----------------------------------------------------------------------
C   RDTORV   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MARYLAND  20899
C
C   FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS WITHIN THE
C        TORUS WHOSE MAJOR RADIUS IS R1 AND MINOR RADIUS IS R2.  NO
C        POINTS WILL BE GENERATED UNLESS R1>0, R2>=0, R1>=R2 AND N>0. 
C        THE TORUS LIES IN THE (X,Y) PLANE CENTERED AT THE ORIGIN.
C        ON OUTPUT, EACH OF THE N COLUMNS OF MATRIX P CONTAINS THE
C        X, Y, AND Z COORDINATES OF THE GENERATED POINTS.  THE LEADING
C        DIMENSION OF P (LDP) MUST BE 3 OR MORE IN THE CALLING PROGRAM.
C
C        A REJECTION PROCEDURE IS USED WHEREBY POINTS ARE FIRST
C        GENERATED UNIFORMLY INSIDE A TORUS-LIKE FIGURE WHOSE
C        CROSS-SECTIONAL AREA IS SQUARE (EACH SIDE LENGTH 2*R2).
C        POINTS WHICH LIE OUTSIDE THE INSCRIBED TORUS ARE REJECTED.
C        REGARDLESS OF R1 AND R2, THE PROPORTION OF POINTS REJECTED
C        IS 0.2146 = 1 - PI/4.  THEREFORE, THE EXPECTED NUMER OF
C        UNIFORM DEVIATES REQUIRED TO GENERATE EACH POINT IS 3.5465.
C
C   SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES
C
C   CURRENT VERSION COMPLETED MARCH 8, 1989
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C   P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM
C              POINTS GENERATED INSIDE THE TORUS.  THE LEADING
C              DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION,
C              DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N.
C              [REAL]
C
C      * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER]
C
C        * N = NUMBER OF POINTS TO BE GENERATED [INTEGER]
C
C       * R1 = MAJOR RADIUS OF THE TORUS (>0) [REAL]
C
C       * R2 = MINOR RADIUS OF THE TORUS (>0) [REAL]
C
C      IFLAG = ERROR FLAG ON OUTPUT.   INTERPRETATION: 
C              0 -> NO ERRORS DETECTED. 
C              1 -> R1<=0 OR R2<0 OR R1<R2.
C              2 -> N<1.
C              3 -> LDP<3.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION P(LDP,*)
      IF (R1.LE.0.0.OR.R2.LT.0.0.OR.R1.LT.R2) THEN
         IFLAG = 1
      ELSEIF (N.LT.1) THEN
         IFLAG = 2
      ELSEIF (LDP.LT.3) THEN
         IFLAG = 3
      ELSE
         IFLAG = 0
         SUM = R1+R2
         DIFF = R1-R2
         R2SQ = R2**2
         TWOPI = 8.0*ATAN(1.0)
         DO 20 I = 1, N
   10       U = RDUNI(0)
            R = SQRT(U*SUM**2+(1.0-U)*DIFF**2)
            Z = R2*(2.0*RDUNI(0)-1.0)
            IF ((R-R1)**2+Z**2.GT.R2SQ) GO TO 10
            PHI = TWOPI*RDUNI(0)
            P(1,I) = R*COS(PHI)
            P(2,I) = R*SIN(PHI)
            P(3,I) = Z
   20    CONTINUE
      ENDIF
      RETURN
C
      END 
*REJ1
      SUBROUTINE REJ1 (Y,N,NPASS,NREJ,SMEAN,SSD,IFLAG)
C
C-----------------------------------------------------------------------
C   REJ1   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C          DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY,
C          GAITHERSBURG, MARYLAND  20899
C
C   FOR: COMPUTING THE MEAN AND STANDARD DEVIATION OF A SAMPLE OF
C        'NORMAL' DATA IN WHICH OUTLIERS MAY BE PRESENT.  OUTLIERS ARE
C        FIRST REJECTED BY A PROCEDURE BASED ON THE SHORTEST INTERVAL 
C        COVERING HALF THE POINTS.  THE PROCEDURE IS ITERATED UNTIL
C
C           1) A USER-SPECIFIED NUMBER OF PASSES OCCURS, OR 
C           2) THE PROPORTION OF VALUES REJECTED IN A GIVEN PASS IS
C              0.01 OR LESS.
C
C        SIMULATION STUDIES ON NORMAL DATA WERE USED TO DETERMINE
C        THE APPROPRIATE VALUES OF CONSTANTS IN THIS PROGRAM.  THEY
C        WERE CHOSEN SO THAT, ON THE FIRST PASS, THE EXPECTED PROPOR- 
C        TION OF 'GOOD' VALUES REJECTED WAS 0.01 REGARDLESS OF SAMPLE 
C        SIZE.  WHEN THE NUMBER OF PASSES ARE NOT LIMITED, THE ACTUAL 
C        PROPORTION OF VALUES REJECTED WAS FOUND TO BE 0.010 TO 0.012 
C        FOR ALL SAMPLE SIZES.
C
C        THE PROCEDURE WAS ORIGINALLY DESIGNED FOR USE ON LARGE SETS
C        OF 'NORMAL' DATA ASYMMETRICALLY CONTAMINATED WITH UP TO 50%
C        OUTLIERS.  ITS BEHAVIOR WAS EXAMINED FOR SAMPLE SIZES OF 15
C        TO 10,000 AND IT APPEARS TO WORK WELL.  WHEN THE SAMPLE SIZE 
C        IS 25 OR LESS, HOWEVER, THE USER MAY WANT TO CONSIDER THE
C        WELL-ESTABLISHED DIXON TEST AS AN ALTERNATIVE.  THAT TEST IS 
C        DISCUSSED IN MANY STATISTICS BOOKS.
C
C   NOTE: ON INPUT, THE SAMPLE SIZE MUST BE AT LEAST 15 AND THE
C         USER-SPECIFIED MAXIMUM NUMBER OF PASSES MUST BE AT LEAST ONE.
C         ON RETURN, THE NUMBER OF 'GOOD' VALUES CAN BE COMPUTED AS:  
C
C              N - NREJ(1) - NREJ(2) - ... - NREJ(NPASS).
C
C   SUBPROGRAMS CALLED: SORT1 (STSPAC)
C
C   CURRENT VERSION COMPLETED OCTOBER 2, 1989
C
C   REFERENCE: REEVE, CHARLES P., 'AN OUTLIER REJECTION PROCEDURE,
C              BASED ON RANGES, FOR ASYMMETRICALLY-CONTAMINATED NORMAL
C              DATA', STATISTICAL ENGINEERING DIVISION NOTE 89-2,
C              DECEMBER 1989. 
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C     * Y(*) = VECTOR (LENGTH N) OF SAMPLE VALUES.  MUST BE 
C              DIMENSIONED AT LEAST N IN CALLING PROGRAM.  [REAL]
C
C        * N = SAMPLE SIZE (N>=15) [INTEGER]
C
C    * NPASS = ON INPUT, MAXIMUM NUMBER OF PASSES (ITERATIONS) OF
C              THE OUTLIER REJECTION PROCEDURE.  ON OUTPUT, THE
C              ACTUAL NUMBER OF PASSES.  [INTEGER]
C
C    NREJ(*) = VECTOR (LENGTH NPASS) CONTAINING, ON OUTPUT, THE
C              NUMBER OF VALUES REJECTED AT EACH PASS.  IT MUST
C              BE DIMENSIONED AT LEAST NPASS IN THE CALLING 
C              PROGRAM.  [INTEGER]
C
C      SMEAN = THE COMPUTED SAMPLE MEAN OF THE NON-REJECTED VALUES
C              [REAL]
C
C        SSD = THE COMPUTED SAMPLE STANDARD DEVIATION OF THE
C              NON-REJECTED VALUES [REAL]
C
C      IFLAG = ERROR INDICATOR ON OUTPUT [INTEGER]   INTERPRETATION:  
C              0 -> NO ERRORS DETECTED. 
C              1 -> N<15 ON INPUT.
C              2 -> NPASS<1 ON INPUT.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION NREJ(*),Y(*)
      LOGICAL LG
      DATA CNORM / 1.349 /
      DATA C1,C2,C3,C4 / 2.576,9.573,-3.013,-0.6989 /
      DATA C5,C6,C7,C8 / 2.576,7.889,1.687,-0.6729 /
C
C--- CHECK FOR VALIDITY OF INPUT VALUES 
C
      IF (N.LT.15) THEN
         IFLAG = 1
      ELSEIF (NPASS.LT.1) THEN
         IFLAG = 2
      ELSE
         IFLAG = 0
C
C--- SORT Y-VALUES
C
         CALL SORT1 (Y,1,N)
C
C--- DETERMINE OUTLIERS BY FINDING THE SHORTEST INTERVAL COVERING
C--- HALF THE POINTS
C
         NADJ = N
         NIT = 0
   10    NIT = NIT+1
         RNA = REAL(NADJ)
         IF (MOD(NADJ,2).EQ.0) THEN
            SIGMLT = C1+C2*(RNA+C3)**C4 
         ELSE
            SIGMLT = C5+C6*(RNA+C7)**C8 
         ENDIF
         L = (NADJ+1)/2
         RMIN = Y(N)-Y(1)
         DO 20 I = 1, N-L
            R = Y(I+L)-Y(I)
            IF (R.LE.RMIN) THEN
               RMIN = R
               K = I
            ENDIF
   20    CONTINUE
         SMEAN = 0.5*(Y(K)+Y(K+L))
         BOUND = SIGMLT*RMIN/CNORM
C
C--- TRIM OUTLIERS AT LOWER END
C
         NLO = 1
   30    IF (SMEAN-Y(NLO).GT.BOUND) THEN
            NLO = NLO+1
            GO TO 30
C
         ENDIF
C
C--- TRIM OUTLIERS AT UPPER END
C
         NHI = N
   40    IF (Y(NHI)-SMEAN.GT.BOUND) THEN
            NHI = NHI-1
            GO TO 40
C
         ENDIF
         NGOOD = NHI-NLO+1
         NREJ(NIT) = NADJ-NGOOD
         LG = NIT.EQ.NPASS.OR.NGOOD.LT.15
         IF (NREJ(NIT).LE.1+INT(0.01*RNA).OR.LG) THEN
            NPASS = NIT
C
C--- COMPUTE MEAN AND STANDARD DEVIATION OF NON-REJECTED VALUES
C
            SMEAN = 0.0
            DO 50 I = NLO, NHI
               SMEAN = SMEAN+Y(I)
   50       CONTINUE
            SMEAN = SMEAN/REAL(NGOOD)
            SSD = 0.0
            DO 60 I = NLO, NHI
               SSD = SSD+(Y(I)-SMEAN)**2
   60       CONTINUE
            SSD = SQRT(SSD/REAL(NGOOD-1))
         ELSE
            NADJ = NGOOD
            GO TO 10
C
         ENDIF
      ENDIF
      RETURN
C
      END 
*SKEKUR
      SUBROUTINE SKEKUR (Y,NLO,NHI,YSKEW,YKURT,IOPT)
C
C-----------------------------------------------------------------------
C   SKEKUR   WRITTEN BY CHARLES P. REEVE
C
C   FOR: COMPUTING SKEWNESS AND KURTOSIS FOR ENTRIES NLO THROUGH NHI
C        IN VECTOR Y.  THE VALUES MAY BE CENTERED ABOUT EITHER THE
C        MEAN (IOPT <> 0) OR ABOUT ZERO (IOPT = 0).  THE TRADITIONAL
C        DIVISIOR OF N (NOT N-1) IS USED WHEN THE MEAN IS ESTIMATED.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED FEBRUARY 28, 1986
C-----------------------------------------------------------------------
C
      DIMENSION Y(*)
      RN = REAL(NHI-NLO+1)
      IF (IOPT.EQ.0) THEN
         S = 0.
      ELSE
         S = 0.
         DO 10 I = NLO, NHI
            S = S+Y(I)
   10    CONTINUE
         S = S/RN
      ENDIF
      T2 = 0.
      T3 = 0.
      T4 = 0.
      DO 20 I = NLO, NHI
         D = Y(I)-S 
         T2 = T2+D**2
         T3 = T3+D**3
         T4 = T4+D**4
   20 CONTINUE
      YSKEW=SQRT(RN)*T3/T2**1.5
      YKURT=RN*T4/T2**2
      RETURN
C
      END 
*SORT1
      SUBROUTINE SORT1 (X,M,N)
C
C-----------------------------------------------------------------------
C   SORT1   OBTAINED BY CHARLES P. REEVE FROM DR. D.A. ZAHN AT THE
C           FLORIDA STATE UNIVERSITY, TALLAHASSEE, FLORIDA UNDER THE
C           NAME 'FTASORT'.  SEVERAL MINOR CORRECTIONS WERE MADE TO
C           THE ORIGINAL VERSION AFTER LINE NUMBER 90.
C
C   FOR: SORTING THE SEGMENT OF A REAL ARRAY BETWEEN ENTRIES M AND N
C        FROM SMALLEST TO LARGEST.  THE ARRAYS LP(K) AND LQ(K) PERMIT 
C        SORTING UP TO 2**(K+1)-1 ELEMENTS, THUS K=25 PERMITS SORTING 
C        UP TO 67,108,863 ELEMENTS.
C
C    SUBROUTINES CALLED: -NONE-
C
C    CURRENT VERSION COMPLETED JANUARY 30, 1978
C-----------------------------------------------------------------------
C
      DIMENSION LP(25),LQ(25),X(*)
      IF (M.GE.N) RETURN
      I = 1
      KM = M
      KN = N
      IF (KM.GE.KN) GO TO 70
   10 J = KM
      K = (KN+KM)/2 
      A = X(K)
      IF (X(KM).LE.A) GO TO 20
      X(K) = X(KM)
      X(KM) = A
      A = X(K)
   20 L = KN
      IF (X(KN).GE.A) GO TO 40
      X(K) = X(KN)
      X(KN) = A
      A = X(K)
      IF (X(KM).LE.A) GO TO 40
      X(K) = X(KM)
      X(KM) = A
      A = X(K)
      GO TO 40
   30 X(L) = X(J)
      X(J) = B
   40 L = L-1
      IF (X(L).GT.A) GO TO 40 
      B = X(L)
   50 J = J+1
      IF (X(J).LT.A) GO TO 50 
      IF (J.LE.L) GO TO 30
      IF (L-KM.LE.KN-J) GO TO 60
      LP(I) = KM
      LQ(I) = L
      KM = J
      I = I+1
      GO TO 80
   60 LP(I) = J
      LQ(I) = KN
      KN = L
      I = I+1
      GO TO 80
   70 I = I-1
      IF (I.EQ.0) RETURN
      KM = LP(I)
      KN = LQ(I)
   80 IF (KN-KM.GE.11) GO TO 10
      KM = KM-1
   90 KM = KM+1
      IF (KM.EQ.KN) GO TO 70
      A = X(KM+1)
      IF (X(KM).LE.A) GO TO 90
      J = KM
  100 X(J+1) = X(J) 
      J = J-1
      IF (J.EQ.M-1) GO TO 110 
      IF (A.LT.X(J)) GO TO 100
  110 X(J+1) = A
      GO TO 90
      END 
*SORT2
      SUBROUTINE SORT2 (X,Y,M,N)
C
C-----------------------------------------------------------------------
C    SORT2   OBTAINED BY CHARLES P. REEVE FROM DR. D. A. ZAHN AT
C            THE FLORIDA STATE UNIVERSITY, TALLAHASSEE, FLORIDA
C            UNDER THE NAME *FTASORT*.  SEVERAL MINOR CORRECTIONS
C            WERE MADE TO THE ORIGINAL VERSION AFTER LINE NUMBER 90.
C
C   FOR: SORTING THE SEGMENT OF REAL ARRAY X BETWEEN ENTRIES M AND
C        N FROM SMALLEST TO LARGEST.  A SECOND ARRAY Y IS CARRIED
C        ALONG AND SORTED IDENTICALLY AS THE FIRST.  THIS FEATURE
C        IS A MODIFICATION TO THE SUBROUTINE SORT1.
C
C    SUBROUTINES CALLED: -NONE-
C
C    CURRENT VERSION COMPLETED OCTOBER 30, 1978
C
C    NOTE: ARRAYS LP(K) AND LQ(K) PERMIT SORTING UP TO 2**(K+1)-1
C          ELEMENTS, I.E., FOR K=25 YOU MAY SORT 67,108,863 ELEMENTS. 
C-----------------------------------------------------------------------
C
      DIMENSION LP(25),LQ(25),X(1),Y(1) 
      I = 1
      KM = M
      KN = N
      IF (KM.GE.KN) GO TO 70
   10 J = KM
      K = (KN+KM)/2 
      A = X(K)
      B = Y(K)
      IF (X(KM).LE.A) GO TO 20
      X(K) = X(KM)
      Y(K) = Y(KM)
      X(KM) = A
      Y(KM) = B
      A = X(K)
      B = Y(K)
   20 L = KN
      IF (X(KN).GE.A) GO TO 40
      X(K) = X(KN)
      Y(K) = Y(KN)
      X(KN) = A
      Y(KN) = B
      A = X(K)
      B = Y(K)
      IF (X(KM).LE.A) GO TO 40
      X(K) = X(KM)
      Y(K) = Y(KM)
      X(KM) = A
      Y(KM) = B
      A = X(K)
      B = Y(K)
      GO TO 40
C
   30 X(L) = X(J)
      Y(L) = Y(J)
      X(J) = G
      Y(J) = H
   40 L = L-1
      IF (X(L).GT.A) GO TO 40 
      G = X(L)
      H = Y(L)
   50 J = J+1
      IF (X(J).LT.A) GO TO 50 
      IF (J.LE.L) GO TO 30
      IF (L-KM.LE.KN-J) GO TO 60
      LP(I) = KM
      LQ(I) = L
      KM = J
      I = I+1
      GO TO 80
C
   60 LP(I) = J
      LQ(I) = KN
      KN = L
      I = I+1
      GO TO 80
C
   70 I = I-1
      IF (I.EQ.0) RETURN
      KM = LP(I)
      KN = LQ(I)
   80 IF (KN-KM.GE.11) GO TO 10
      KM = KM-1
   90 KM = KM+1
      IF (KM.EQ.KN) GO TO 70
      A = X(KM+1)
      B = Y(KM+1)
      IF (X(KM).LE.A) GO TO 90
      J = KM
  100 X(J+1) = X(J) 
      Y(J+1) = Y(J) 
      J = J-1
      IF (J.EQ.M-1) GO TO 110 
      IF (A.LT.X(J)) GO TO 100
  110 X(J+1) = A
      Y(J+1) = B
      GO TO 90
C
      END 
C
*TRISYS
      SUBROUTINE TRISYS (SUB,DIAG,SUP,B,N,IFLAG)
C
C-----------------------------------------------------------------------
C   TRISYS   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL INSTITUTE FOR STANDARDS AND TECHNOLOGY,
C            GAITHERSBURG, MD  20899
C
C   FOR: SOLVING THE TRIDIAGONAL SYSTEM OF LINEAR EQUATIONS 
C
C             SUB X    + DIAG X  + SUP X    = B
C                I I-1       I I      I I+1    I
C
C        WHERE I = 1,2,...,N.  SUB(1) AND SUP(N) ARE ZERO BY DEFINITION,
C        THEREFORE THEIR PASSED VALUES ARE IGNORED.  THE CONTENTS OF
C        DIAG ARE ALTERED AND THE SOLUTION VECTOR (X) IS RETURNED IN B.
C        SEE ALGORITHM 3.3 IN THE REFERENCE.
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED SEPTEMBER 23, 1988
C
C   REFERENCE: CONTE, S.D. AND DEBOOR, CARL, 'ELEMENTARY NUMERICAL
C              ANALYSIS', MCGRAW-HILL BOOK COMPANY, 1972, PP. 121-122.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C    * SUB(*) = VECTOR (LENGTH N) OF SUB-DIAGONAL ELEMENTS.  THE VALUE
C               IN SUB(1) IS NOT USED. [REAL]
C
C   * DIAG(*) = VECTOR (LENGTH N) OF DIAGONAL ELEMENTS.  ALL VALUES
C               ARE CHANGED ON OUTPUT. [REAL]
C
C    * SUP(*) = VECTOR (LENGTH N) OF SUPER-DIAGONAL ELEMENTS.  THE
C               VALUE IN SUP(N) IS NOT USED. [REAL]
C
C      * B(*) = VECTOR (LENGTH N) CONTAINING THE RIGHT HAND SIDE OF
C               THE SYSTEM ON INPUT, AND THE SOLUTION VECTOR ON OUTPUT.
C               [REAL]
C
C         * N = THE NUMBER OF EQUATIONS AND UNKNOWNS [INTEGER]
C
C       IFLAG = ERROR INDICATOR ON OUTPUT.   INTERPRETATION: 
C               -1 -> N<1, THEREFORE NO SOLUTION RETURNED.
C                0 -> NO ERRORS DETECTED.
C                K -> AT SOME POINT DIAG(K)=0.0, THEREFORE FAILURE.
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
      DIMENSION SUB(*),DIAG(*),SUP(*),B(*)
      IFLAG = 0
C
C--- CHECK FOR N<1
C
      IF (N.LT.1) THEN
         IFLAG = -1 
         RETURN
C
      ENDIF
C
C--- SOLVE TRIDIAGONAL SYSTEM 
C
      DO 10 K = 2, N
         IF (DIAG(K-1).EQ.0.0) THEN
            IFLAG = K-1
            RETURN
C
         ELSE
            RATIO = -SUB(K)/DIAG(K-1)
            DIAG(K) = DIAG(K)+RATIO*SUP(K-1)
            B(K) = B(K)+RATIO*B(K-1)
         ENDIF
   10 CONTINUE
      IF (DIAG(N).EQ.0.0) THEN
         IFLAG = N
         RETURN
C
      ELSE
         B(N) = B(N)/DIAG(N)
      ENDIF
      DO 20 K = N-1, 1, -1
         IF (DIAG(K).EQ.0.0) THEN
            IFLAG = K
            RETURN
C
         ELSE
            B(K) = (B(K)-SUP(K)*B(K+1))/DIAG(K)
         ENDIF
   20 CONTINUE
      RETURN
C
      END 
*UPDATE
      SUBROUTINE UPDATE (OLDMN,OLDSD,NEWMN,NEWSD,YNEW,N)
C
C-----------------------------------------------------------------------
C   UPDATE   WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND  20899
C
C   FOR: TAKING A PREVIOUSLY COMPUTED MEAN (OLDMN) AND STANDARD
C        DEVIATION (OLDSD) BASED ON N-1 VALUES AND COMPUTING AN
C        UPDATED MEAN (NEWMN) AND STANDARD DEVIATION (OLDSD) BASED
C        ON THE N-1 PREVIOUS VALUES PLUS THE NTH VALUE (YNEW).
C
C   NOTE: ALL PASSED PARAMETERS OTHER THAN N ARE OF TYPE REAL.
C
C   SPECIAL CASES:  
C
C        IF N < 1 THEN A MESSAGE IS PRINTED AND EXECUTION STOPPED.
C        IF N = 1 THEN INPUT VALUES OF OLDMN AND OLDSD ARE IGNORED
C                 WHILE COMPUTING NEWMN = YNEW AND NEWSD = 0.0 .
C
C   SUBPROGRAMS CALLED: -NONE-
C
C   CURRENT VERSION COMPLETED MARCH 2, 1986
C-----------------------------------------------------------------------
C
      REAL NEWMN,NEWSD
      IF(N.GE.2) THEN
      RN = REAL(N)
      RN1 = REAL(N-1)
      NEWMN = (RN1*OLDMN+YNEW)/RN
      NEWSD = SQRT(REAL(N-2)*OLDSD**2/RN1+(YNEW-OLDMN)**2/RN)
      ELSE IF(N.EQ.1) THEN
      NEWMN=YNEW
      NEWSD=0.0
      ELSE
      PRINT *,' *** N MUST BE A POSITIVE INTEGER' 
      PRINT *,' *** EXECUTION STOPPED IN SUBROUTINE UPDATE' 
      STOP
      END IF
      RETURN
      END 
*ZEROBR
      SUBROUTINE ZEROBR (FUNC,X1,X2,TOL,X0,IFLAG) 
C
C-----------------------------------------------------------------------
C   ZEROBR   COPIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING
C            DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
C            MARYLAND 20899 FROM REFERENCE 1 BELOW WITH SLIGHT
C            MODIFICATIONS.
C
C   FOR: FINDING THE REAL ZERO X0 OF THE REAL FUNCTION FUNC KNOWN TO
C        LIE BETWEEN X1 AND X2 USING BRENT'S METHOD AS DESCRIBED IN
C        THE REFERENCES.  THE ZERO IS REFINED UNTIL ITS ABSOLUTE
C        ACCURACY IS TOL.
C
C   NOTE: THE CONSTANT EPS FOR MACHINE FLOATING POINT PRECISION IS
C         MACHINE DEPENDENT.
C
C   SUBPROGRAMS CALLED: FUNC (USER DEFINED)
C
C   CURRENT VERSION COMPLETED OCTOBER 13, 1987
C
C   REFERENCES: 
C
C   1) PRESS, WILLIAM H., FLANNERY, BRIAN P., TEUKOLSKY, SAUL A.,
C      AND VETTERLING, WILLIAM T., 'NUMERICAL RECIPES - THE ART OF
C      SCIENTIFIC COMPUTING', CAMBRIDGE UNIVERSITY PRESS, 1986,
C      PP. 251-254. 
C
C   2) BRENT, RICHARD P., 'ALGORITHMS FOR MINIMIZATION WITHOUT
C      DERIVATIVES', PRENTICE-HALL, 1973, CH. 3-4.
C-----------------------------------------------------------------------
C   DEFINITION OF PASSED PARAMETERS: 
C
C    FUNC = AN EXTERNAL FUNCTION WITH ONE ARGUMENT (REAL)
C
C    * X1 = THE LEFT ENDPOINT OF AN INTERVAL CONTAINING A ZERO (REAL) 
C
C    * X2 = THE RIGHT ENDPOINT OF AN INTERVAL CONTAINING A ZERO (REAL)
C
C   * TOL = THE ABSOLUTE TOLERANCE ON THE COMPUTED ZERO (REAL)
C
C      X0 = THE COMPUTED ZERO OF FUNC IN THE INTERVAL (X1,X2) (REAL)
C
C   IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER)   INTERPRETATION: 
C           0 -> NO ERRORS DETECTED
C           1 -> FUNC(X1) AND FUNC(X2) HAVE THE SAME SIGN - THE ZERO
C                IS NOT DEFINED
C           2 -> MAXIMUM ITERATIONS EXCEEDED - CURRENT VALUE OF ZERO
C                RETURNED
C
C   * INDICATES PARAMETERS REQUIRING INPUT VALUES 
C-----------------------------------------------------------------------
C
C--- DEFINE MAXIMUM ITERATIONS AND MACHINE FLOATING POINT PRECISION
C
      DATA ITMAX,EPS / 50,1.0E-14 /
      A = X1
      B = X2
      FA = FUNC(A)
      FB = FUNC(B)
C
C--- CHECK FOR ZERO BEING BRACKETED
C
      IF (FB*FA.GT.0.0) THEN
         IFLAG = 1
         RETURN
C
      ENDIF
      FC = FB
      DO 10 ITER = 1, ITMAX
         IF (FB*FC.GT.0.0) THEN
C
C--- RENAME A,B,C AND ADJUST BOUNDING INTERVAL D
C
            C = A
            FC = FA 
            D = B-A 
            E = D
         ENDIF
         IF (ABS(FC).LT.ABS(FB)) THEN
            A = B
            B = C
            C = A
            FA = FB 
            FB = FC 
            FC = FA 
         ENDIF
C
C--- CONVERGENCE CHECK
C
         TOL1 = 2.0*EPS*ABS(B)+0.5*TOL
         XM = 0.5*(C-B)
         IF (ABS(XM).LE.TOL1.OR.FB.EQ.0.0) THEN
            X0 = B
            RETURN
C
         ENDIF
         IF (ABS(E).GE.TOL1.AND.ABS(FA).GT.ABS(FB)) THEN
C
C--- ATTEMPT INVERSE QUADRATIC INTERPOLATION
C
            S = FB/FA
            IF (A.EQ.C) THEN
               P = 2.0*XM*S
               Q = 1.0-S
            ELSE
               Q = FA/FC
               R = FB/FC
               P = S*(2.0*XM*Q*(Q-R)-(B-A)*(R-1.0))
               Q = (Q-1.0)*(R-1.0)*(S-1.0)
            ENDIF
C
C--- CHECK WHETHER IN BOUNDS
C
            IF (P.GT.0.0) Q = -Q
            P = ABS(P)
            IF (2.0*P.LT.AMIN1(3.0*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
C
C--- ACCEPT INTERPLATION
C
               E = D
               D = P/Q
            ELSE
C
C--- INTERPOLATION FAILED, USE BISECTION
C
               D = XM
               E = D
            ENDIF
         ELSE
C
C--- BOUNDS DECREASING TOO SLOWLY, USE BISECTION
C
            D = XM
            E = D
         ENDIF
C
C--- MOVE LAST BEST GUESS TO A
C
         A = B
         FA = FB
C
C--- EVALUATE NEW TRIAL ZERO
C
         IF (ABS(D).GT.TOL1) THEN
            B = B+D 
         ELSE
            B = B+SIGN(TOL1,XM)
         ENDIF
         FB = FUNC(B)
   10 CONTINUE
      IFLAG = 2
      X0 = B
      RETURN
C
      END