!   E04NGA Example Program Text
!   Mark 27.3 Release. NAG Copyright 2021.

    Module e04ngae_mod

!     E04NGA 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, lcwsav = 1, liwsav = 610,  &
                                          llwsav = 120, lrwsav = 475, nin = 5, &
                                          ninopt = 7, nout = 6
    Contains
      Subroutine qphess(n,jthcol,h,ldh,x,hx,iuser,ruser,iwsav)
!       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)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Integer, Intent (Inout)        :: iuser(*), iwsav(610)
!       .. 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 e04ngae_mod
    Program e04ngae

!     E04NGA Example Main Program

!     .. Use Statements ..
      Use e04ngae_mod, Only: iset, lcwsav, liwsav, llwsav, lrwsav, nin,        &
                             ninopt, nout, qphess
      Use nag_library, Only: e04nfa, e04nga, e04nha, e04wbf, nag_wp, x04abf,   &
                             x04acf, x04baf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Character (*), Parameter         :: fname = 'e04ngae.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(:)
      Real (Kind=nag_wp)               :: ruser(1), rwsav(lrwsav)
      Integer, Allocatable             :: istate(:), iwork(:)
      Integer                          :: iuser(1), iwsav(liwsav)
      Logical                          :: lwsav(llwsav)
      Character (80)                   :: cwsav(lcwsav)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..
      Write (rec,99992) 'E04NGA 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)

!     Initialise E04NFA

      ifail = 0
      Call e04wbf('E04NFA',cwsav,lcwsav,lwsav,llwsav,iwsav,liwsav,rwsav,       &
        lrwsav,ifail)

!     Set three options using E04NHA

      Call e04nha(' Check Frequency = 10 ',lwsav,iwsav,rwsav,inform)

      If (inform==0) Then

        Call e04nha(' Crash Tolerance = 0.05 ',lwsav,iwsav,rwsav,inform)

        If (inform==0) Then

          Call e04nha(' Infinite Bound Size = 1.0D+25 ',lwsav,iwsav,rwsav,     &
            inform)

        End If

      End If

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

!     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 e04nga(ninopt,lwsav,iwsav,rwsav,inform)

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

!     Solve the problem

      ifail = -1
      Call e04nfa(n,nclin,a,lda,bl,bu,cvec,h,ldh,qphess,istate,x,iter,obj,ax,  &
        clamda,iwork,liwork,work,lwork,iuser,ruser,lwsav,iwsav,rwsav,ifail)

      Select Case (ifail)
      Case (0:5,7:)
        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,99998)
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)

        Do i = 1, n
          Write (rec,99997) i, istate(i), x(i), clamda(i)
          Call x04baf(nout,rec)
        End Do

        If (nclin>0) Then
          Write (rec,'()')
          Call x04baf(nout,rec)
          Write (rec,'()')
          Call x04baf(nout,rec)
          Write (rec,99996)
          Call x04baf(nout,rec)
          Write (rec,'()')
          Call x04baf(nout,rec)

          Do i = n + 1, n + nclin
            j = i - n
            Write (rec,99995) j, istate(i), ax(j), clamda(i)
            Call x04baf(nout,rec)
          End Do

        End If

        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,99994) obj
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,99993) iter
        Call x04baf(nout,rec)
      End Select

100   Continue

99999 Format (1X,A,I5)
99998 Format (1X,'Varbl',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99997 Format (1X,'V',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99996 Format (1X,'L Con',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99995 Format (1X,'L',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99994 Format (1X,'Final objective value = ',G15.7)
99993 Format (1X,'Exit from problem after',1X,I6,1X,'iterations.')
99992 Format (1X,A)
    End Program e04ngae