* 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) * 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(*) * .. 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