*     E04NGA Example Program Text
*     Mark 20 Release. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NMAX, NCMAX
      PARAMETER        (NMAX=10,NCMAX=10)
      INTEGER          LDA
      PARAMETER        (LDA=NCMAX)
      INTEGER          LIWORK, LWORK
      PARAMETER        (LIWORK=1000,LWORK=10000)
      INTEGER          LCWSAV, LIWSAV, LLWSAV, LRWSAV
      PARAMETER        (LCWSAV=1,LIWSAV=610,LLWSAV=120,LRWSAV=475)
*     .. Local Scalars ..
      DOUBLE PRECISION OBJ
      INTEGER          I, IFAIL, INFORM, ITER, J, LDH, N, NCLIN
      CHARACTER        UPLO
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), AX(NCMAX), BL(NMAX+NCMAX),
     +                 BU(NMAX+NCMAX), CLAMDA(NMAX+NCMAX), CVEC(NMAX),
     +                 H(NMAX*(NMAX+1)/2), RUSER(1), RWSAV(LRWSAV),
     +                 WORK(LWORK), X(NMAX)
      INTEGER          ISTATE(NMAX+NCMAX), IUSER(1), IWORK(LIWORK),
     +                 IWSAV(LIWSAV)
      LOGICAL          LWSAV(LLWSAV)
      CHARACTER*80     CWSAV(LCWSAV)
*     .. External Subroutines ..
      EXTERNAL         E04NFA, E04NGA, E04NHA, E04WBF, QPHESS, X04ABF
*     .. Executable Statements ..
      WRITE (NOUT,*) 'E04NGA Example Program Results'
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) N, NCLIN
      IF (N.LE.NMAX .AND. NCLIN.LE.NCMAX) THEN
*
*        Read CVEC, A, BL, BU, X, UPLO and H from data file
*
         READ (NIN,*) (CVEC(I),I=1,N)
         READ (NIN,*) ((A(I,J),J=1,N),I=1,NCLIN)
         READ (NIN,*) (BL(I),I=1,N+NCLIN)
         READ (NIN,*) (BU(I),I=1,N+NCLIN)
         READ (NIN,*) (X(I),I=1,N)
         READ (NIN,*) UPLO
         IF (UPLO.EQ.'U') THEN
*           Read the upper triangle of H
            READ (NIN,*) ((H(J+(2*N-I)*(I-1)/2),J=I,N),I=1,N)
         ELSE IF (UPLO.EQ.'L') THEN
*           Read the lower triangle of H
            READ (NIN,*) ((H(I+(2*N-J)*(J-1)/2),J=1,I),I=1,N)
         END IF
         LDH = N*(N+1)/2
*
*        Set the unit number for advisory messages to NOUT
*
         CALL X04ABF(1,NOUT)
*
*        Initialise using E04WBF and check for error exits
*
         IFAIL = 1
         CALL E04WBF('E04NFA',CWSAV,LCWSAV,LWSAV,LLWSAV,IWSAV,LIWSAV,
     +               RWSAV,LRWSAV,IFAIL)
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99999) 'E04WBF terminated with IFAIL = ', IFAIL
         ELSE
*
*           Set three options using E04NHA
*
            CALL E04NHA(' Check Frequency = 10 ',LWSAV,IWSAV,RWSAV,
     +                  INFORM)
            IF (INFORM.EQ.0) THEN
               CALL E04NHA(' Crash Tolerance = 0.05 ',LWSAV,IWSAV,RWSAV,
     +                     INFORM)
               IF (INFORM.EQ.0) THEN
                  CALL E04NHA(' Infinite Bound Size = 1.0D+25 ',LWSAV,
     +                        IWSAV,RWSAV,INFORM)
               END IF
            END IF
            IF (INFORM.NE.0) THEN
               WRITE (NOUT,99999) 'E04NHA terminated with INFORM = ',
     +           INFORM
            ELSE
*
*              Read the options file for the remaining options
*
               CALL E04NGA(NIN,LWSAV,IWSAV,RWSAV,INFORM)
               IF (INFORM.NE.0) THEN
                  WRITE (NOUT,99999) 'E04NGA terminated with INFORM =',
     +              INFORM
               END IF
            END IF
            IF (INFORM.EQ.0) THEN
*
*              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)
*
*              Check for error exits
*
               WRITE (NOUT,*)
               IF (IFAIL.EQ.6) THEN
                  WRITE (NOUT,99998)
               ELSE
                  WRITE (NOUT,99997) IFAIL
                  WRITE (NOUT,*)
                  WRITE (NOUT,99996)
                  WRITE (NOUT,*)
                  DO 20 I = 1, N
                     WRITE (NOUT,99995) I, ISTATE(I), X(I), CLAMDA(I)
   20             CONTINUE
                  IF (NCLIN.GT.0) THEN
                     WRITE (NOUT,*)
                     WRITE (NOUT,*)
                     WRITE (NOUT,99994)
                     WRITE (NOUT,*)
                     DO 40 I = N + 1, N + NCLIN
                        J = I - N
                        WRITE (NOUT,99993) J, ISTATE(I), AX(J),
     +                    CLAMDA(I)
   40                CONTINUE
                  END IF
                  WRITE (NOUT,*)
                  WRITE (NOUT,*)
                  WRITE (NOUT,99992) OBJ
                  WRITE (NOUT,*)
                  WRITE (NOUT,*)
                  WRITE (NOUT,99991) ITER
               END IF
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (1X,A,I3)
99998 FORMAT (1X,'An input parameter is invalid')
99997 FORMAT (1X,'E04NFA returned with IFAIL = ',I4)
99996 FORMAT (1X,'Varbl',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99995 FORMAT (1X,'V',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99994 FORMAT (1X,'L Con',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99993 FORMAT (1X,'L',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99992 FORMAT (1X,'Final objective value = ',G15.7)
99991 FORMAT (1X,'Exit from problem after',1X,I6,1X,'iterations.')
      END
*
      SUBROUTINE QPHESS(N,JTHCOL,HESS,LDHESS,X,HX,IUSER,RUSER,IWSAV)
*     In this version of QPHESS, the lower triangle of the matrix H is
*     stored in packed form (by columns) in the one-dimensional array
*     HESS. More precisely, the lower triangle of H must be stored with
*     element H(i,j) in HESS(i+(2*N-j)*(j-1)/2) for i .ge. j.
*     Note that storing the lower triangle of H in packed form (by
*     columns) is equivalent to storing the upper triangle of H in
*     packed form (by rows).
*     Note also that LDHESS is used to define the length of HESS, and
*     must therefore be at least N*(N+1)/2.
*     .. Scalar Arguments ..
      INTEGER           JTHCOL, LDHESS, N
*     .. Array Arguments ..
      DOUBLE PRECISION  HESS(LDHESS), HX(N), RUSER(*), X(N)
      INTEGER           IUSER(*), IWSAV(*)
*     .. Local Scalars ..
      DOUBLE PRECISION  S
      INTEGER           I, INC, J, L, LP1
*     .. Executable Statements ..
      IF (JTHCOL.NE.0) THEN
*        Special case -- extract one column of  H.
         L = JTHCOL
         INC = N
         DO 20 I = 1, JTHCOL
            HX(I) = HESS(L)
            INC = INC - 1
            L = L + INC
   20    CONTINUE
         L = L - INC + 1
         IF (JTHCOL.LT.N) THEN
            LP1 = L
            DO 40 I = JTHCOL + 1, N
               HX(I) = HESS(LP1)
               LP1 = LP1 + 1
   40       CONTINUE
         END IF
      ELSE
*        Normal case.
         L = 0
         DO 80 I = 1, N
            S = 0.0D0
            DO 60 J = I, N
               L = L + 1
               S = S + HESS(L)*X(J)
   60       CONTINUE
            HX(I) = S
   80    CONTINUE
         L = 0
         DO 120 J = 1, N - 1
            L = L + 1
            DO 100 I = J + 1, N
               L = L + 1
               HX(I) = HX(I) + HESS(L)*X(J)
  100       CONTINUE
  120    CONTINUE
      END IF
      RETURN
      END
