!   G02HDF Example Program Text
!   Mark 25 Release. NAG Copyright 2014.

    Module g02hdfe_mod

!     G02HDF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                               :: betcal, chi, psi
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter        :: dchi = 1.5_nag_wp
      Real (Kind=nag_wp), Parameter        :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: two = 2.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: zero = 0.0_nag_wp
      Integer, Parameter, Public           :: iset = 1, nin = 5, nout = 6
    Contains
      Function psi(t)

!       .. Function Return Value ..
        Real (Kind=nag_wp)                   :: psi
!       .. Parameters ..
        Real (Kind=nag_wp), Parameter        :: c = 1.5_nag_wp
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
!       .. Intrinsic Procedures ..
        Intrinsic                            :: abs
!       .. Executable Statements ..
        If (t<=-c) Then
          psi = -c
        Else If (abs(t)<c) Then
          psi = t
        Else
          psi = c
        End If
        Return
      End Function psi

      Function chi(t)

!       .. Function Return Value ..
        Real (Kind=nag_wp)                   :: chi
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: ps
!       .. Intrinsic Procedures ..
        Intrinsic                            :: abs
!       .. Executable Statements ..
        ps = dchi
        If (abs(t)<dchi) Then
          ps = t
        End If
        chi = ps*ps/two
        Return
      End Function chi

      Subroutine betcal(n,wgt,beta)
!       Calculate BETA for Schweppe type regression

!       .. Use Statements ..
        Use nag_library, Only: s15abf, x01aaf, x02akf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: beta
        Integer, Intent (In)                 :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: wgt(n)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: amaxex, anormc, b, d2, dc, dw, &
                                                dw2, pc, w2
        Integer                              :: i, ifail
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp, log, real, sqrt
!       .. Executable Statements ..
        amaxex = -log(x02akf())
        anormc = sqrt(x01aaf(zero)*two)
        d2 = dchi*dchi
        beta = zero
        Do i = 1, n
          w2 = wgt(i)*wgt(i)
          dw = wgt(i)*dchi
          ifail = 0
          pc = s15abf(dw,ifail)
          dw2 = dw*dw
          dc = zero
          If (dw2<amaxex) Then
            dc = exp(-dw2/two)/anormc
          End If
          b = (-dw*dc+pc-0.5_nag_wp)/w2 + (one-pc)*d2
          beta = b*w2/real(n,kind=nag_wp) + beta
        End Do
        Return
      End Subroutine betcal
    End Module g02hdfe_mod
    Program g02hdfe
!     G02HDF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: g02hdf, nag_wp, x04abf
      Use g02hdfe_mod, Only: betcal, chi, iset, nin, nout, psi
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: beta, eps, psip0, sigma, tol
      Integer                              :: i, ifail, indw, isigma, k, ldx,  &
                                              m, maxit, n, nadv, nit, nitmon
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: rs(:), theta(:), wgt(:), wk(:),  &
                                              x(:,:), y(:)
!     .. Executable Statements ..
      Write (nout,*) 'G02HDF Example Program Results'
      Write (nout,*)

!     Skip heading in data file
      Read (nin,*)

!     Read in the problem size
      Read (nin,*) n, m

      ldx = n
      Allocate (x(ldx,m),y(n),wgt(n),theta(m),rs(n),wk((m+4)*n))

!     Read in data
      Read (nin,*)(x(i,2:m),y(i),wgt(i),i=1,n)

!     Set first column of X to 1 for the constant term
      x(1:n,1) = 1.0E0_nag_wp

!     Set BETA
      Call betcal(n,wgt,beta)

!     Read in value for PSI(0)
      Read (nin,*) psip0

!     Read in control parameters
      Read (nin,*) indw, isigma, nitmon, maxit, tol, eps

!     Set the advisory channel to NOUT for monitoring information
      If (nitmon/=0) Then
        nadv = nout
        Call x04abf(iset,nadv)
      End If

!     Read in initial values
      Read (nin,*) sigma
      Read (nin,*) theta(1:m)

!     Perform bounded influence regression
      ifail = -1
      Call g02hdf(chi,psi,psip0,beta,indw,isigma,n,m,x,ldx,y,wgt,theta,k, &
        sigma,rs,tol,eps,maxit,nitmon,nit,wk,ifail)
      If (ifail==7) Then
        Write (nout,*) 'Some of the following results may be unreliable'
      Else If (ifail/=0) Then
        Go To 100
      End If

!     Display results
      Write (nout,99999) 'G02HDF required ', nit, ' iterations to converge'
      Write (nout,99999) '                   K = ', k
      Write (nout,99998) '               Sigma = ', sigma
      Write (nout,*) '    THETA'
      Write (nout,99997)(theta(i),i=1,m)
      Write (nout,*)
      Write (nout,*) '  Weights  Residuals'
      Write (nout,99996)(wgt(i),rs(i),i=1,n)

100   Continue

99999 Format (1X,A,I0,A)
99998 Format (1X,A,F9.4)
99997 Format (1X,F9.4)
99996 Format (1X,2F9.4)
    End Program g02hdfe