Program g02qgfe

!     G02QGF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: g02qgf, g02zkf, g02zlf, g05kff, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: lseed = 1, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: df, rvalue
      Integer                          :: genid, i, ifail, ip, ivalue, j, l,   &
                                          ldbl, lddat, ldres, liopts, lopts,   &
                                          lstate, lwt, m, n, ntau, optype,     &
                                          sorder, subid, tdch
      Character (1)                    :: c1, weight
      Character (30)                   :: cvalue, semeth
      Character (100)                  :: optstr
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: b(:,:), bl(:,:), bu(:,:), ch(:,:,:), &
                                          dat(:,:), opts(:), res(:,:), tau(:), &
                                          wt(:), y(:)
      Integer, Allocatable             :: info(:), iopts(:), isx(:), state(:)
      Integer                          :: seed(lseed)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: count, len_trim, min
!     .. Executable Statements ..
      Write (nout,*) 'G02QGF Example Program Results'
      Write (nout,*)
      Flush (nout)

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

!     Read in the problem size
      Read (nin,*) sorder, c1, weight, n, m, ntau

!     Read in the data
      If (weight=='W' .Or. weight=='w') Then
        lwt = n
      Else
        lwt = 0
      End If
      Allocate (wt(lwt),isx(m),y(n),tau(ntau))

      If (sorder==1) Then
!       DAT(N,M)
        lddat = n
        Allocate (dat(lddat,m))
        If (lwt==0) Then
          Read (nin,*)(dat(i,1:m),y(i),i=1,n)
        Else
          Read (nin,*)(dat(i,1:m),y(i),wt(i),i=1,n)
        End If
      Else
!       DAT(M,N)
        lddat = m
        Allocate (dat(lddat,n))
        If (lwt==0) Then
          Read (nin,*)(dat(1:m,i),y(i),i=1,n)
        Else
          Read (nin,*)(dat(1:m,i),y(i),wt(i),i=1,n)
        End If
      End If

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

!     Calculate IP
      ip = count(isx(1:m)==1)
      If (c1=='Y' .Or. c1=='y') Then
        ip = ip + 1
      End If

!     Read in the quantiles required
      Read (nin,*) tau(1:ntau)

      liopts = 100
      lopts = 100
      Allocate (iopts(liopts),opts(lopts))

!     Initialize the optional argument array
      ifail = 0
      Call g02zkf('INITIALIZE = G02QGF',iopts,liopts,opts,lopts,ifail)

c_lp: Do
!       Read in any optional arguments. Reads in to the end of
!       the input data, or until a blank line is reached
        ifail = 1
        Read (nin,99994,Iostat=ifail) optstr
        If (ifail/=0) Then
          Exit c_lp
        Else If (len_trim(optstr)==0) Then
          Exit c_lp
        End If

!       Set the supplied option
        ifail = 0
        Call g02zkf(optstr,iopts,liopts,opts,lopts,ifail)
      End Do c_lp

!     Assume that no intervals or output matrices are required
!     unless the optional argument state differently
      ldbl = 0
      tdch = 0
      ldres = 0
      lstate = 0

!     Query the optional arguments to see what output is required
      ifail = 0
      Call g02zlf('INTERVAL METHOD',ivalue,rvalue,cvalue,optype,iopts,opts,    &
        ifail)
      semeth = cvalue
      If (semeth/='NONE') Then
!       Require the intervals to be output
        ldbl = ip

        If (semeth=='BOOTSTRAP XY') Then
!         Need to find the length of the state array for the random
!         number generator

!         Read in the generator ID and a seed
          Read (nin,*) genid, subid, seed(1)

!         Query the length of the state array
          Allocate (state(lstate))
          ifail = 0
          Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!         Deallocate STATE so that it can reallocated later
          Deallocate (state)
        End If

        ifail = 0
        Call g02zlf('MATRIX RETURNED',ivalue,rvalue,cvalue,optype,iopts,opts,  &
          ifail)
        If (cvalue=='COVARIANCE') Then
          tdch = ntau
        Else If (cvalue=='H INVERSE') Then
          If (semeth=='BOOTSTRAP XY' .Or. semeth=='IID') Then
!           NB: If we are using bootstrap or IID errors then any request for
!           H INVERSE is ignored
            tdch = 0
          Else
            tdch = ntau + 1
          End If
        End If

        ifail = 0
        Call g02zlf('RETURN RESIDUALS',ivalue,rvalue,cvalue,optype,iopts,opts, &
          ifail)
        If (cvalue=='YES') Then
          ldres = n
        End If
      End If

!     Allocate memory for output arrays
      Allocate (b(ip,ntau),info(ntau),bl(ldbl,ntau),bu(ldbl,ntau),             &
        ch(ldbl,ldbl,tdch),state(lstate),res(ldres,ntau))

      If (lstate>0) Then
!       Doing bootstrap, so initialize the RNG
        ifail = 0
        Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)
      End If

!     Call the model fitting routine
      ifail = -1
      Call g02qgf(sorder,c1,weight,n,m,dat,lddat,isx,ip,y,wt,ntau,tau,df,b,bl, &
        bu,ch,res,iopts,opts,state,info,ifail)
      If (ifail/=0) Then
        If (ifail==231) Then
          Write (nout,*) 'Additional error information (INFO): ', info(1:ntau)
        Else
          Go To 100
        End If
      End If

!     Display the parameter estimates
      Do l = 1, ntau
        Write (nout,99999) 'Quantile: ', tau(l)
        Write (nout,*)
        If (ldbl>0) Then
          Write (nout,*) '        Lower   Parameter   Upper'
          Write (nout,*) '        Limit   Estimate    Limit'
        Else
          Write (nout,*) '      Parameter'
          Write (nout,*) '      Estimate'
        End If
        Do j = 1, ip
          If (ldbl>0) Then
            Write (nout,99998) j, bl(j,l), b(j,l), bu(j,l)
          Else
            Write (nout,99998) j, b(j,l)
          End If
        End Do
        Write (nout,*)
        If (tdch==ntau) Then
          Write (nout,*) 'Covariance matrix'
          Do i = 1, ip
            Write (nout,99997) ch(1:i,i,l)
          End Do
          Write (nout,*)
        Else If (tdch==ntau+1) Then
          Write (nout,*) 'J'
          Do i = 1, ip
            Write (nout,99997) ch(1:i,i,1)
          End Do
          Write (nout,*)
          Write (nout,*) 'H inverse'
          Do i = 1, ip
            Write (nout,99997) ch(1:i,i,l+1)
          End Do
          Write (nout,*)
        End If
        Write (nout,*)
      End Do

      If (ldres>0) Then
        Write (nout,*) 'First 10 Residuals'
        Write (nout,*) '                              Quantile'
        Write (nout,99995) 'Obs.', tau(1:ntau)
        Do i = 1, min(n,10)
          Write (nout,99996) i, res(i,1:ntau)
        End Do
      Else
        Write (nout,*) 'Residuals not returned'
      End If
      Write (nout,*)

100   Continue

99999 Format (1X,A,F6.3)
99998 Format (1X,I3,3(3X,F7.3))
99997 Format (1X,10(E10.3,3X))
99996 Format (2X,I3,10(1X,F10.5))
99995 Format (1X,A,10(3X,F6.3,2X))
99994 Format (A100)
    End Program g02qgfe