!   E04NGF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    Module e04ngfe_mod

!     E04NGF 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                           :: qphess
!     .. Parameters ..
      Integer, Parameter, Public       :: iset = 1, nin = 5, ninopt = 7,       &
                                          nout = 6
    Contains
      Subroutine qphess(n,jthcol,h,ldh,x,hx)
!       In this version of QPHESS, the lower triangle of matrix H is
!       stored in packed form (by columns) in array H.
!       More precisely, the lower triangle of matrix H must be stored with
!       matrix element H(i,j) in array element H(i+(2*N-j)*(j-1)/2,1),
!       for i .ge. j.
!       Note that storing the lower triangle of matrix H in packed form (by
!       columns) is equivalent to storing the upper triangle of matrix H in
!       packed form (by rows).
!       Note also that LDH is used to define the size of array H, and
!       must therefore be at least N*(N+1)/2.

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: jthcol, ldh, n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: h(ldh,*), x(n)
        Real (Kind=nag_wp), Intent (Out) :: hx(n)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: s
        Integer                        :: i, inc, j, l, lp1
!       .. Executable Statements ..
        If (jthcol/=0) Then

!         Special case -- extract one column of  H.

          l = jthcol
          inc = n

          Do i = 1, jthcol
            hx(i) = h(l,1)
            inc = inc - 1
            l = l + inc
          End Do

          l = l - inc + 1

          If (jthcol<n) Then
            lp1 = l

            Do i = jthcol + 1, n
              hx(i) = h(lp1,1)
              lp1 = lp1 + 1
            End Do

          End If

        Else

!         Normal case.

          l = 0

          Do i = 1, n
            s = 0.0E0_nag_wp

            Do j = i, n
              l = l + 1
              s = s + h(l,1)*x(j)
            End Do

            hx(i) = s
          End Do

          l = 0

          Do j = 1, n - 1
            l = l + 1

            Do i = j + 1, n
              l = l + 1
              hx(i) = hx(i) + h(l,1)*x(j)
            End Do

          End Do

        End If

        Return

      End Subroutine qphess
    End Module e04ngfe_mod
    Program e04ngfe

!     E04NGF Example Main Program

!     .. Use Statements ..
      Use e04ngfe_mod, Only: iset, nin, ninopt, nout, qphess
      Use nag_library, Only: e04nff, e04ngf, e04nhf, nag_wp, x04abf, x04acf,   &
                             x04baf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Character (*), Parameter         :: fname = 'e04ngfe.opt'
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: obj
      Integer                          :: i, ifail, inform, iter, j, lda, ldh, &
                                          liwork, lwork, mode, n, nclin,       &
                                          outchn, sda
      Character (80)                   :: rec
      Character (1)                    :: uplo
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), ax(:), bl(:), bu(:),         &
                                          clamda(:), cvec(:), h(:,:), work(:), &
                                          x(:)
      Integer, Allocatable             :: istate(:), iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..
      Write (rec,99998) 'E04NGF Example Program Results'
      Call x04baf(nout,rec)

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

      Read (nin,*) n, nclin
      liwork = 2*n + 3
      lda = max(1,nclin)

      If (nclin>0) Then
        sda = n
      Else
        sda = 1
      End If

!     This particular example problem is of type QP2 with a nondefault QPHESS,
!     so we allocate CVEC(N) and H(LDH,1), and define LDH and LWORK as below

      ldh = n*(n+1)/2

      If (nclin>0) Then
        lwork = 2*n**2 + 8*n + 5*nclin
      Else
        lwork = n**2 + 8*n
      End If

      Allocate (istate(n+nclin),ax(max(1,nclin)),iwork(liwork),h(ldh,1),bl(n+  &
        nclin),bu(n+nclin),cvec(n),x(n),a(lda,sda),clamda(n+nclin),            &
        work(lwork))

      Read (nin,*) cvec(1:n)
      Read (nin,*)(a(i,1:n),i=1,nclin)
      Read (nin,*) bl(1:(n+nclin))
      Read (nin,*) bu(1:(n+nclin))
      Read (nin,*) x(1:n)
      Read (nin,*) uplo

      If (uplo=='U') Then

!       Read the upper triangle of H

        Read (nin,*)((h(j+(2*n-i)*(i-1)/2,1),j=i,n),i=1,n)
      Else If (uplo=='L') Then

!       Read the lower triangle of H

        Read (nin,*)((h(i+(2*n-j)*(j-1)/2,1),j=1,i),i=1,n)
      End If

      ldh = n*(n+1)/2

!     Set the unit number for advisory messages to OUTCHN

      outchn = nout
      Call x04abf(iset,outchn)

!     Set four options using E04NHF

      Call e04nhf(' Print Level = 1 ')

      Call e04nhf(' Check Frequency = 10 ')

      Call e04nhf(' Crash Tolerance = 0.05 ')

      Call e04nhf(' Infinite Bound Size = 1.0D+25 ')

!     Open the options file for reading

      mode = 0

      ifail = 0
      Call x04acf(ninopt,fname,mode,ifail)

!     Read the options file for the remaining options

      Call e04ngf(ninopt,inform)

      If (inform/=0) Then
        Write (rec,99999) 'E04NGF terminated with INFORM =', inform
        Call x04baf(nout,rec)
        Go To 100
      End If

!     Solve the problem

      ifail = 0
      Call e04nff(n,nclin,a,lda,bl,bu,cvec,h,ldh,qphess,istate,x,iter,obj,ax,  &
        clamda,iwork,liwork,work,lwork,ifail)

100   Continue

99999 Format (1X,A,I5)
99998 Format (1X,A)
    End Program e04ngfe