!   E04NGA Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    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
!      .. Parameters ..
       INTEGER, PARAMETER                  :: iset = 1, lcwsav = 1,            &
                                              liwsav = 610, llwsav = 120,      &
                                              lrwsav = 475, nin = 5,           &
                                              ninopt = 7, nout = 6
    CONTAINS
!dec$ attributes offload : mic :: qphess
       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.

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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

!dec$ attributes offload : mic :: qphess
       END SUBROUTINE qphess
    END MODULE e04ngae_mod
    PROGRAM e04ngae

!      E04NGA Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : e04nfa, e04nga, e04nha, e04wbf, nag_wp, x04abf, &
                               x04acf, x04baf
       USE e04ngae_mod, ONLY : iset, lcwsav, liwsav, llwsav, lrwsav, nin,      &
                               ninopt, nout, qphess
!      .. 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 Functions ..
       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 20
       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 20
       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

20     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
