!   G22YBF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2017.

    Module g22ybfe_mod
!     G22YBF Example Program Module:
!            Parameters and User-defined Routines

!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: construct_labels, fit_lm, read_line
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6

    Contains
      Subroutine read_line(v1)
!       Read in a line from NIN and remove any comments

!       .. Scalar Arguments ..
        Character (*), Intent (Out)    :: v1
!       .. Local Scalars ..
        Integer                        :: pend
!       .. Intrinsic Procedures ..
        Intrinsic                      :: adjustl, index
!       .. Executable Statements ..
        Continue

        Read (nin,'(A200)') v1
        pend = index(v1,'::')
        If (pend/=0) Then
          v1 = v1(1:pend-1)
        End If
        v1 = adjustl(v1)

        Return
      End Subroutine read_line
      Subroutine construct_labels(ip,plab,vnames,vinfo)
!       Construct a set of parameter labels using the information
!       stored in VINFO

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: ip
!       .. Array Arguments ..
        Integer, Intent (In)           :: vinfo(:)
        Character (*), Allocatable, Intent (Out) :: plab(:)
        Character (*), Intent (In)     :: vnames(:)
!       .. Local Scalars ..
        Integer                        :: b, i, j, k, li, vi
        Character (200)                :: line, term
!       .. Intrinsic Procedures ..
        Intrinsic                      :: adjustl, trim
!       .. Executable Statements ..
        Continue

        Allocate (plab(ip))

        k = 1
        Do j = 1, ip
          b = vinfo(k)
          k = k + 1

          If (b==0) Then
            line = 'Intercept'
          Else
            line = ''
            Do i = 1, b
              vi = vinfo(k)
              li = vinfo(k+1)

              If (li/=0) Then
                Write (term,99999) trim(adjustl(vnames(vi))), li
              Else
                Write (term,99998) trim(adjustl(vnames(vi)))
              End If

              If (i==1) Then
                line = trim(line) // trim(term)
              Else
                line = trim(line) // ' . ' // trim(term)
              End If

!             We are ignoring the contrast identifier when constructing
!             these labels
              k = k + 3
            End Do
          End If
          plab(j) = line
        End Do
99999   Format (A,' (Level ',I0,')')
99998   Format (A)

        Return
      End Subroutine construct_labels
      Subroutine fit_lm(hform,intcpt,nobs,mx,x,isx,ip,y,plab)
!       Perform a multiple linear regression using G02DAF

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
        Use nag_library, Only: g02daf, g22znf, nag_wp
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: hform
        Integer, Intent (In)           :: ip, mx, nobs
        Character (*), Intent (In)     :: intcpt
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x(:,:), y(:)
        Integer, Intent (In)           :: isx(:)
        Character (*), Intent (In)     :: plab(:)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: rss, rvalue, tol
        Integer                        :: i, idf, ifail, irank, ivalue, ldq,   &
                                          ldx, lwt, optype
        Logical                        :: svd
        Character (200)                :: cvalue
        Character (1)                  :: weight
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: b(:), cov(:), h(:), p(:), q(:,:),   &
                                          res(:), se(:), wk(:), wt(:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: repeat, size, trim
!       .. Executable Statements ..
        Continue

        ldx = size(x,1)

!       We are assuming un-weighted data
        weight = 'U'
        lwt = 0
        ldq = nobs
        Allocate (b(ip),cov((ip*ip+ip)/2),h(nobs),p(ip*(ip+                    &
          2)),q(ldq,ip+1),res(nobs),se(ip),wk(ip*ip+5*(ip-1)),wt(lwt))

!       Use suggested value for tolerance
        tol = 0.000001E0_nag_wp

!       Fit a regression model
        ifail = 0
        Call g02daf(intcpt,weight,nobs,x,ldx,mx,isx,ip,y,wt,rss,idf,b,se,cov,  &
          res,h,q,ldq,svd,irank,p,tol,wk,ifail)

!       Get the formula for the model being fit
        ifail = 0
        Call g22znf(hform,'Formula',ivalue,rvalue,cvalue,optype,ifail)

!       Display the results
        Write (nout,*) 'Model: ', trim(cvalue)
        Write (nout,*) '                               Parameter   Standard'
        Write (nout,*) 'Coefficients                   Estimate      Error'
        Write (nout,*) repeat('-',51)
        Do i = 1, ip
          Write (nout,99997) plab(i), b(i), se(i)
        End Do
        Write (nout,*) repeat('-',51)
        Write (nout,99998) 'Residual sum of squares = ', rss
        Write (nout,99999) 'Degrees of freedom      = ', idf

        Return
99999   Format (1X,A,I9)
99998   Format (1X,A,F9.4)
99997   Format (1X,A30,1X,F7.3,5X,F7.3)
      End Subroutine fit_lm
    End Module g22ybfe_mod

    Program g22ybfe

!     .. Use Statements ..
      Use g22ybfe_mod, Only: construct_labels, fit_lm, nin, nout, read_line
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: g22yaf, g22ybf, g22ycf, g22ydf, g22zaf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (c_ptr)                     :: hddesc, hform, hxdesc
      Integer                          :: i, ifail, ip, lddat, ldx, lisx,      &
                                          lplab, lvinfo, lvnames, mx, nobs,    &
                                          nvar, sddat, sdx
      Character (200)                  :: formula
      Character (1)                    :: intcpt
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: dat(:,:), x(:,:), y(:)
      Real (Kind=nag_wp)               :: tx(0,0)
      Integer, Allocatable             :: isx(:), levels(:), vinfo(:)
      Integer                          :: tisx(0), tvinfo(3)
      Character (50), Allocatable      :: plab(:), vnames(:)
      Character (1)                    :: tplab(0)
!     .. Executable Statements ..
      Write (nout,*) 'G22YBF Example Program Results'
      Write (nout,*)

      hform = c_null_ptr
      hddesc = c_null_ptr
      hxdesc = c_null_ptr

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

!     Read in size of the data matrix and number of variable labels supplied
      Read (nin,*) nobs, nvar, lvnames

!     Read in number of levels and names for the variables
      Allocate (levels(nvar),vnames(lvnames))
      Read (nin,*) levels(1:nvar)
      If (lvnames>0) Then
        Read (nin,*) vnames(1:lvnames)
      End If

!     Create a description of the data matrix
      ifail = 0
      Call g22ybf(hddesc,nobs,nvar,levels,lvnames,vnames,ifail)

!     Read in the data matrix and response variable
      lddat = nobs
      sddat = nvar
      Allocate (dat(lddat,sddat),y(nobs))
      Read (nin,*)(dat(i,1:nvar),y(i),i=1,nobs)

!     Read in the model formula, remove comments and parse it
      Call read_line(formula)
      ifail = 0
      Call g22yaf(hform,formula,ifail)

!     Start of constructing the design matrix ...
!     Calculate the size of design matrix
      ldx = 0
      sdx = 0
      ifail = 1
      Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
      If (ifail/=91) Then
!       redisplay the error message
        ifail = 0
        Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
      End If

!     Generate the design matrix
      ldx = nobs
      sdx = mx
      Allocate (x(ldx,sdx))
      ifail = 0
      Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,x,ldx,sdx,mx,ifail)
!     ... End of constructing the design matrix

!     Start of getting the ISX vector and information on parameter labels ...
!     Get size of output arrays used by G22YDF
      lvinfo = 3
      lisx = 0
      lplab = 0
      ifail = 1
      Call g22ydf(hform,hxdesc,intcpt,ip,lisx,tisx,lplab,tplab,lvinfo,tvinfo,  &
        ifail)
      If (ifail/=92) Then
!       redisplay the error message
        ifail = 0
        Call g22ydf(hform,hxdesc,intcpt,ip,lisx,tisx,lplab,tplab,lvinfo,       &
          tvinfo,ifail)
      End If

!     Allocate output arrays (we already know that LISX = MX, but G22YDF
!     returns it just in case)
      lisx = tvinfo(1)
      lvinfo = tvinfo(3)
!     We don't need PLAB as we are constructing our own labels from VINFO
      Allocate (isx(lisx),vinfo(lvinfo))

!     Get the ISX flag and parameter labels
      ifail = 0
      Call g22ydf(hform,hxdesc,intcpt,ip,lisx,isx,lplab,tplab,lvinfo,vinfo,    &
        ifail)

!     Construct some verbose labels for the parameters
      Call construct_labels(ip,plab,vnames,vinfo)
!     ... End of getting the ISX vector and information on parameter labels

!     Fit a regression model and print the results
      Call fit_lm(hform,intcpt,nobs,mx,x,isx,ip,y,plab)

!     Clean-up the G22 handles
      ifail = 0
      Call g22zaf(hform,ifail)
      Call g22zaf(hddesc,ifail)
      Call g22zaf(hxdesc,ifail)

      Deallocate (dat,x,y)
      Deallocate (isx,levels,vinfo)
      Deallocate (plab,vnames)
    End Program g22ybfe