*     G02HDF Example Program Text
*     Mark 14 Revised. NAG Copyright 1989.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NMAX, MMAX
      PARAMETER        (NMAX=9,MMAX=3)
*     .. Local Scalars ..
      DOUBLE PRECISION BETA, EPS, PSIP0, SIGMA, TOL
      INTEGER          I, IFAIL, INDW, ISIGMA, IX, J, K, M, MAXIT, N,
     +                 NIT, NITMON
*     .. Local Arrays ..
      DOUBLE PRECISION RS(NMAX), THETA(MMAX), WGT(NMAX),
     +                 WK(NMAX*(MMAX+4)), X(NMAX,MMAX), Y(NMAX)
*     .. External Functions ..
      DOUBLE PRECISION CHI, PSI
      EXTERNAL         CHI, PSI
*     .. External Subroutines ..
      EXTERNAL         BETCAL, G02HDF, X04ABF
*     .. Executable Statements ..
      WRITE (NOUT,*) 'G02HDF Example Program Results'
*     Skip heading in data file
      READ (NIN,*)
      CALL X04ABF(1,NOUT)
*     Read in the dimensions of X
      READ (NIN,*) N, M
      IF ((N.LE.NMAX) .AND. (M.LE.MMAX)) THEN
*        Read in the X matrix, the Y values and set X(i,1) to 1 for the
*        constant term
         DO 20 I = 1, N
            READ (NIN,*) (X(I,J),J=2,M), Y(I)
            X(I,1) = 1.0D0
   20    CONTINUE
*        Read in weights
         DO 40 I = 1, N
            READ (NIN,*) WGT(I)
   40    CONTINUE
         CALL BETCAL(N,WGT,BETA)
*        Set other parameter values
         IX = NMAX
         MAXIT = 50
         TOL = 0.5D-4
         EPS = 0.5D-5
         PSIP0 = 1.0D0
*        Set value of ISIGMA and initial value of SIGMA
         ISIGMA = 1
         SIGMA = 1.0D0
*        Set initial value of THETA
         DO 60 J = 1, M
            THETA(J) = 0.0D0
   60    CONTINUE
*        * Change NITMON to a positive value if monitoring information
*          is required *
         NITMON = 0
*        Schweppe type regression
         INDW = 1
         IFAIL = -1
*
         CALL G02HDF(CHI,PSI,PSIP0,BETA,INDW,ISIGMA,N,M,X,IX,Y,WGT,
     +               THETA,K,SIGMA,RS,TOL,EPS,MAXIT,NITMON,NIT,WK,IFAIL)
*
         WRITE (NOUT,*)
         IF (IFAIL.NE.0 .AND. IFAIL.NE.7) THEN
            WRITE (NOUT,99999) 'G02HDF fails, IFAIL = ', IFAIL
         ELSE
            IF (IFAIL.EQ.7) THEN
               WRITE (NOUT,99999) 'G02HDF returned IFAIL = ', IFAIL
               WRITE (NOUT,*)
     +           'Some of the following results may be unreliable'
            END IF
            WRITE (NOUT,99998) 'G02HDF required ', NIT,
     +        ' iterations to converge'
            WRITE (NOUT,99998) '                   K = ', K
            WRITE (NOUT,99997) '               Sigma = ', SIGMA
            WRITE (NOUT,*) '    THETA'
            DO 80 J = 1, M
               WRITE (NOUT,99996) THETA(J)
   80       CONTINUE
            WRITE (NOUT,*)
            WRITE (NOUT,*) '  Weights  Residuals'
            DO 100 I = 1, N
               WRITE (NOUT,99995) WGT(I), RS(I)
  100       CONTINUE
         END IF
      END IF
      STOP
*
99999 FORMAT (1X,A,I2)
99998 FORMAT (1X,A,I4,A)
99997 FORMAT (1X,A,F9.4)
99996 FORMAT (1X,F9.4)
99995 FORMAT (1X,2F9.4)
      END
*
      DOUBLE PRECISION FUNCTION PSI(T)
*     .. Parameters ..
      DOUBLE PRECISION              C
      PARAMETER                     (C=1.5D0)
*     .. Scalar Arguments ..
      DOUBLE PRECISION              T
*     .. Intrinsic Functions ..
      INTRINSIC                     ABS
*     .. Executable Statements ..
      IF (T.LE.-C) THEN
         PSI = -C
      ELSE IF (ABS(T).LT.C) THEN
         PSI = T
      ELSE
         PSI = C
      END IF
      RETURN
      END
*
      DOUBLE PRECISION FUNCTION CHI(T)
*     .. Parameters ..
      DOUBLE PRECISION              DCHI
      PARAMETER                     (DCHI=1.5D0)
*     .. Scalar Arguments ..
      DOUBLE PRECISION              T
*     .. Local Scalars ..
      DOUBLE PRECISION              PS
*     .. Intrinsic Functions ..
      INTRINSIC                     ABS
*     .. Executable Statements ..
      PS = DCHI
      IF (ABS(T).LT.DCHI) PS = T
      CHI = PS*PS/2.0D0
      RETURN
      END
*
      SUBROUTINE BETCAL(N,WGT,BETA)
*     Calculate BETA for Schweppe type regression
*     .. Parameters ..
      DOUBLE PRECISION  DCHI
      PARAMETER         (DCHI=1.5D0)
*     .. Scalar Arguments ..
      DOUBLE PRECISION  BETA
      INTEGER           N
*     .. Array Arguments ..
      DOUBLE PRECISION  WGT(N)
*     .. Local Scalars ..
      DOUBLE PRECISION  AMAXEX, ANORMC, B, D2, DC, DW, DW2, PC, W2
      INTEGER           I, IFAIL
*     .. External Functions ..
      DOUBLE PRECISION  S15ABF, X01AAF, X02AKF
      EXTERNAL          S15ABF, X01AAF, X02AKF
*     .. Intrinsic Functions ..
      INTRINSIC         DBLE, EXP, LOG, SQRT
*     .. Executable Statements ..
      AMAXEX = -LOG(X02AKF())
      ANORMC = SQRT(X01AAF(0.0D0)*2.0D0)
      D2 = DCHI*DCHI
      BETA = 0.0D0
      DO 20 I = 1, N
         W2 = WGT(I)*WGT(I)
         DW = WGT(I)*DCHI
         IFAIL = 0
         PC = S15ABF(DW,IFAIL)
         DW2 = DW*DW
         DC = 0.0D0
         IF (DW2.LT.AMAXEX) DC = EXP(-DW2/2.0D0)/ANORMC
         B = (-DW*DC+PC-0.5D0)/W2 + (1.0D0-PC)*D2
         BETA = B*W2/DBLE(N) + BETA
   20 CONTINUE
      RETURN
      END