Program g02dcfe

!     G02DCF Example Program Text

!     Mark 25 Release. NAG Copyright 2014.

!     .. Use Statements ..
      Use nag_library, Only: g02daf, g02dcf, g02ddf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: rss, tol, wt, y
      Integer                          :: i, idf, ifail, ip, irank, ix, ldq,   &
                                          ldxm, lwt, m, n
      Logical                          :: svd
      Character (1)                    :: mean, update, weight
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: b(:), cov(:), h(:), p(:), q(:,:),    &
                                          res(:), se(:), wk(:), wtm(:), x(:),  &
                                          xm(:,:), ym(:)
      Integer, Allocatable             :: isx(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: count
!     .. Executable Statements ..
      Write (nout,*) 'G02DCF Example Program Results'
      Write (nout,*)

!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n, m, weight, mean

      If (weight=='W' .Or. weight=='w') Then
        lwt = n
      Else
        lwt = 0
      End If
      ldxm = n
      Allocate (xm(ldxm,m),ym(n),wtm(lwt),isx(m),x(m))

!     Read in data
      If (lwt>0) Then
        Read (nin,*)(xm(i,1:m),ym(i),wtm(i),i=1,n)
      Else
        Read (nin,*)(xm(i,1:m),ym(i),i=1,n)
      End If

!     Read in variable inclusion flags
      Read (nin,*) isx(1:m)

!     Calculate IP
      ip = count(isx>0)
      If (mean=='M' .Or. mean=='m') Then
        ip = ip + 1
      End If

      ldq = n
      Allocate (b(ip),cov((ip*ip+ip)/2),h(n),p(ip*(ip+ &
        2)),q(ldq,ip+1),res(n),se(ip),wk(ip*ip+5*(ip-1)))

!     Use suggested value for tolerance
      tol = 0.000001E0_nag_wp

!     Fit initial model using G02DAF
      ifail = 0
      Call g02daf(mean,weight,n,xm,ldxm,m,isx,ip,ym,wtm,rss,idf,b,se,cov,res, &
        h,q,ldq,svd,irank,p,tol,wk,ifail)

!     Display results from initial model fit
      Write (nout,*) 'Results from initial model fit using G02DAF'
      If (svd) Then
        Write (nout,*)
        Write (nout,*) 'Model not of full rank'
      End If
      Write (nout,99999) 'Residual sum of squares = ', rss
      Write (nout,99998) 'Degrees of freedom = ', idf
      Write (nout,*)
      Write (nout,*) 'Variable   Parameter estimate   ', 'Standard error'
      Write (nout,*)
      Write (nout,99997)(i,b(i),se(i),i=1,ip)

!     Updating data is held in X consecutively
      ix = 1

!     Add or delete observations supplied in the data file
u_lp: Do
        Read (nin,*,Iostat=ifail) update
        If (ifail/=0) Then
          Exit u_lp
        End If

        If (lwt>0) Then
          Read (nin,*) x(1:m), y, wt
        Else
          Read (nin,*) x(1:m), y
        End If

!       Update the regression
        ifail = 0
        Call g02dcf(update,mean,weight,m,isx,q,ldq,ip,x,ix,y,wt,rss,wk,ifail)

!       Display titles and update observation count
        Write (nout,*)
        Select Case (update)
        Case ('a','A')
          Write (nout,*) 'Results from adding an observation using G02DCF'
          n = n + 1
        Case ('d','D')
          Write (nout,*) 'Results from dropping an observation using G02DCF'
          n = n - 1
        Case Default
          Write (nout,*) 'Unknown update flag read in from data file'
          Go To 100
        End Select

!       Recalculate the parameter estimates etc
        ifail = 0
        Call g02ddf(n,ip,q,ldq,rss,idf,b,se,cov,svd,irank,p,tol,wk,ifail)

!       Display updated results                        
        Write (nout,99999) 'Residual sum of squares = ', rss
        Write (nout,99998) 'Degrees of freedom = ', idf
        Write (nout,*)
        Write (nout,*) 'Variable   Parameter estimate   ', 'Standard error'
        Write (nout,*)
        Write (nout,99997)(i,b(i),se(i),i=1,ip)
      End Do u_lp

100   Continue

99999 Format (1X,A,E12.4)
99998 Format (1X,A,I4)
99997 Format (1X,I6,2E20.4)
    End Program g02dcfe