! E04NGA Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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 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 (jthcol0) 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