#include "../include/dprec.fh" !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine ql0001 here] subroutine ql0001(m,me,mmax,n,nmax,mnn,c,d,a,b,xl,xu, & x,u,iout,ifail,iprint,war,lwar,iwar,liwar,eps) implicit none !************************************************************************** ! SOLUTION OF QUADRATIC PROGRAMMING PROBLEMS ! QL0001 SOLVES THE QUADRATIC PROGRAMMING PROBLEM ! MINIMIZE .5*X'*C*X + D'*X ! SUBJECT TO A(J)*X + B(J) = 0 , J=1,...,ME ! A(J)*X + B(J) >= 0 , J=ME+1,...,M ! XL <= X <= XU ! HERE C MUST BE AN N BY N SYMMETRIC AND POSITIVE MATRIX, D AN N-DIMENSIONAL ! VECTOR, A AN M BY N MATRIX AND B AN M-DIMENSIONAL VECTOR. THE ABOVE ! SITUATION IS INDICATED BY IWAR(1)=1. ALTERNATIVELY, I.E. IF IWAR(1)=0, ! THE OBJECTIVE FUNCTION MATRIX CAN ALSO BE PROVIDED IN FACTORIZED FORM. ! IN THIS CASE, C IS AN UPPER TRIANGULAR MATRIX. ! THE SUBROUTINE REORGANIZES SOME DATA SO THAT THE PROBLEM CAN BE SOLVED ! BY A MODIFICATION OF AN ALGORITHM PROPOSED BY POWELL (1983). ! USAGE: ! QL0001(M,ME,MMAX,N,NMAX,MNN,C,D,A,B,XL,XU,X,U,IOUT,IFAIL,IPRINT, ! WAR,LWAR,IWAR,LIWAR) ! DEFINITION OF THE PARAMETERS: ! M : TOTAL NUMBER OF CONSTRAINTS. ! ME : NUMBER OF EQUALITY CONSTRAINTS. ! MMAX : ROW DIMENSION OF A. MMAX MUST BE AT LEAST ONE AND GREATER ! THAN M. ! N : NUMBER OF VARIABLES. ! NMAX : ROW DIMENSION OF C. NMAX MUST BE GREATER OR EQUAL TO N. ! MNN : MUST BE EQUAL TO M + N + N. ! C(NMAX,NMAX): OBJECTIVE FUNCTION MATRIX WHICH SHOULD BE SYMMETRIC AND ! POSITIVE DEFINITE. IF IWAR(1) = 0, C IS SUPPOSED TO BE THE ! CHOLESKEY-FACTOR OF ANOTHER MATRIX, I.E. C IS UPPER ! TRIANGULAR. ! D(NMAX) : CONTAINS THE CONSTANT VECTOR OF THE OBJECTIVE FUNCTION. ! A(MMAX,NMAX): CONTAINS THE DATA MATRIX OF THE LINEAR CONSTRAINTS. ! B(MMAX) : CONTAINS THE CONSTANT DATA OF THE LINEAR CONSTRAINTS. ! XL(N),XU(N): CONTAIN THE LOWER AND UPPER BOUNDS FOR THE VARIABLES. ! X(N) : ON RETURN, X CONTAINS THE OPTIMAL SOLUTION VECTOR. ! U(MNN) : ON RETURN, U CONTAINS THE LAGRANGE MULTIPLIERS. THE FIRST ! M POSITIONS ARE RESERVED FOR THE MULTIPLIERS OF THE M ! LINEAR CONSTRAINTS AND THE SUBSEQUENT ONES FOR THE ! MULTIPLIERS OF THE LOWER AND UPPER BOUNDS. ON SUCCESSFUL ! TERMINATION, ALL VALUES OF U WITH RESPECT TO INEQUALITIES ! AND BOUNDS SHOULD BE GREATER OR EQUAL TO ZERO. ! IOUT : INTEGER INDICATING THE DESIRED OUTPUT UNIT NUMBER, I.E. ! ALL WRITE-STATEMENTS START WITH 'WRITE(IOUT,... '. ! IFAIL : SHOWS THE TERMINATION REASON. ! IFAIL = 0 : SUCCESSFUL RETURN. ! IFAIL = 1 : TOO MANY ITERATIONS (MORE THAN 40*(N+M)). ! IFAIL = 2 : ACCURACY INSUFFICIENT TO SATISFY CONVERGENCE ! CRITERION. ! IFAIL = 5 : LENGTH OF A WORKING ARRAY IS TOO SHORT. ! IFAIL > 10 : THE CONSTRAINTS ARE INCONSISTENT. ! IPRINT : OUTPUT CONTROL. ! IPRINT = 0 : NO OUTPUT OF QL0001. ! IPRINT > 0 : BRIEF OUTPUT IN ERROR CASES. ! WAR(LWAR) : _REAL_ WORKING ARRAY. THE LENGTH LWAR SHOULD BE GRATER THAN ! 3*NMAX*NMAX/2 + 10*NMAX + 2*MMAX. ! IWAR(LIWAR): INTEGER WORKING ARRAY. THE LENGTH LIWAR SHOULD BE AT ! LEAST N. ! IF IWAR(1)=1 INITIALLY, THEN THE CHOLESKY DECOMPOSITION ! WHICH IS REQUIRED BY THE DUAL ALGORITHM TO GET THE FIRST ! UNCONSTRAINED MINIMUM OF THE OBJECTIVE FUNCTION, IS ! PERFORMED INTERNALLY. OTHERWISE, I.E. IF IWAR(1)=0, THEN ! IT IS ASSUMED THAT THE USER PROVIDES THE INITIAL FAC- ! TORIZATION BY HIMSELF AND STORES IT IN THE UPPER TRIAN- ! GULAR PART OF THE ARRAY C. ! A NAMED COMMON-BLOCK /CMACHE/EPS MUST BE PROVIDED BY THE USER, ! WHERE EPS DEFINES A GUESS FOR THE UNDERLYING MACHINE PRECISION. ! AUTHOR (C): K. SCHITTKOWSKI, ! MATHEMATISCHES INSTITUT, ! UNIVERSITAET BAYREUTH, ! 95440 BAYREUTH, ! GERMANY, F.R. ! VERSION: 1.5 (JUNE, 1991) !********************************************************************* !Passed variables lw ?? integer m,me,nmax,mmax,n,mnn,lwar,liwar,iout,ifail,iprint,iwar(liwar) _REAL_ c(nmax,n),d(n),a(mmax,n),b(mmax), & xl(n),xu(n),x(n),u(mnn),war(lwar) _REAL_ diag,zero,eps,qpeps,ten integer inw1,inw2,j,lw,mn,i,idiag,info,nact,maxit,in logical lql ! INTRINSIC FUNCTIONS: sqrt !common /cmache/eps ! CONSTANT DATA lql=.false. if (iwar(1) == 1) lql=.true. zero=0.0d+0 ten=1.d+1 maxit=40*(m+n) qpeps=eps inw1=1 inw2=inw1+mmax ! PREPARE PROBLEM DATA FOR EXECUTION if (m <= 0) goto 20 in=inw1 do j=1,m war(in)=-b(j) in=in+1 end do 20 lw=3*nmax*nmax/2 + 10*nmax + m if ((inw2+lw) > lwar) goto 80 if (liwar < n) goto 81 if (mnn < m+n+n) goto 82 mn=m+n ! CALL OF QL0002 call ql0002(n,m,me,mmax,mn,mnn,nmax,lql,a,war(inw1), & d,c,xl,xu,x,nact,iwar,maxit,qpeps,info,diag, & war(inw2),lw) ! TEST OF MATRIX CORRECTIONS ifail=0 if (info == 1) goto 40 if (info == 2) goto 90 idiag=0 if ((diag > zero).and.(diag < 1000.0)) idiag=diag if ((iprint > 0).and.(idiag > 0)) & write(iout,1000) idiag if (info < 0) goto 70 ! REORDER MULTIPLIER do j=1,mnn u(j)=zero end do in=inw2-1 if (nact == 0) goto 30 do i=1,nact j=iwar(i) u(j)=war(in+i) end do 30 continue return ! ERROR MESSAGES 70 ifail=-info+10 if ((iprint > 0).and.(nact > 0)) & write(iout,1100) -info,(iwar(i),i=1,nact) return 80 ifail=5 if (iprint > 0) write(iout,1200) return 81 ifail=5 if (iprint > 0) write(iout,1210) return 82 ifail=5 if (iprint > 0) write(iout,1220) return 40 ifail=1 if (iprint > 0) write(iout,1300) maxit return 90 ifail=2 if (iprint > 0) write(iout,1400) return ! FORMAT-INSTRUCTIONS 1000 format(/8x,28h***ql: matrix g was enlarged,i3, & 20h-times by unitmatrix) 1100 format(/8x,18h***ql: constraint ,i5, & 19h not consistent to ,/,(10x,10i5)) 1200 format(/8x,21h***ql: lwar too small) 1210 format(/8x,22h***ql: liwar too small) 1220 format(/8x,20h***ql: mnn too small) 1300 format(/8x,37h***ql: too many iterations (more than,i6,1h)) 1400 format(/8x,50h***ql: accuracy insufficient to attain convergence) end subroutine ql0001 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine ql0002 here] subroutine ql0002(n,m,meq,mmax,mn,mnn,nmax,lql,a,b,grad,g, & xl,xu,x,nact,iact,maxit,vsmall,info,diag,w,lw) implicit none !************************************************************************** ! THIS SUBROUTINE SOLVES THE QUADRATIC PROGRAMMING PROBLEM ! MINIMIZE GRAD'*X + 0.5 * X*G*X ! SUBJECT TO A(K)*X = B(K) K=1,2,...,MEQ, ! A(K)*X >= B(K) K=MEQ+1,...,M, ! XL <= X <= XU ! THE QUADRATIC PROGRAMMING METHOD PROCEEDS FROM AN INITIAL CHOLESKY- ! DECOMPOSITION OF THE OBJECTIVE FUNCTION MATRIX, TO CALCULATE THE ! UNIQUELY DETERMINED MINIMIZER OF THE UNCONSTRAINED PROBLEM. ! SUCCESSIVELY ALL VIOLATED CONSTRAINTS ARE ADDED TO A WORKING SET ! AND A MINIMIZER OF THE OBJECTIVE FUNCTION SUBJECT TO ALL CONSTRAINTS ! IN THIS WORKING SET IS COMPUTED. IT IS POSSIBLE THAT CONSTRAINTS ! HAVE TO LEAVE THE WORKING SET. ! DESCRIPTION OF PARAMETERS: ! N : IS THE NUMBER OF VARIABLES. ! M : TOTAL NUMBER OF CONSTRAINTS. ! MEQ : NUMBER OF EQUALITY CONTRAINTS. ! MMAX : ROW DIMENSION OF A, DIMENSION OF B. MMAX MUST BE AT ! LEAST ONE AND GREATER OR EQUAL TO M. ! MN : MUST BE EQUAL M + N. ! MNN : MUST BE EQUAL M + N + N. ! NMAX : ROW DIEMSION OF G. MUST BE AT LEAST N. ! LQL : DETERMINES INITIAL DECOMPOSITION. ! LQL = .FALSE. : THE UPPER TRIANGULAR PART OF THE MATRIX G ! CONTAINS INITIALLY THE CHOLESKY-FACTOR OF A SUITABLE ! DECOMPOSITION. ! LQL = .TRUE. : THE INITIAL CHOLESKY-FACTORISATION OF G IS TO BE ! PERFORMED BY THE ALGORITHM. ! A(MMAX,NMAX) : A IS A MATRIX WHOSE COLUMNS ARE THE CONSTRAINTS NORMALS. ! B(MMAX) : CONTAINS THE RIGHT HAND SIDES OF THE CONSTRAINTS. ! GRAD(N) : CONTAINS THE OBJECTIVE FUNCTION VECTOR GRAD. ! G(NMAX,N): CONTAINS THE SYMMETRIC OBJECTIVE FUNCTION MATRIX. ! XL(N), XU(N): CONTAIN THE LOWER AND UPPER BOUNDS FOR X. ! X(N) : VECTOR OF VARIABLES. ! NACT : FINAL NUMBER OF ACTIVE CONSTRAINTS. ! IACT(K) (K=1,2,...,NACT): INDICES OF THE FINAL ACTIVE CONSTRAINTS. ! INFO : REASON FOR THE RETURN FROM THE SUBROUTINE. ! INFO = 0 : CALCULATION WAS TERMINATED SUCCESSFULLY. ! INFO = 1 : MAXIMUM NUMBER OF ITERATIONS ATTAINED. ! INFO = 2 : ACCURACY IS INSUFFICIENT TO MAINTAIN INCREASING ! FUNCTION VALUES. ! INFO < 0 : THE CONSTRAINT WITH INDEX ABS(INFO) AND THE CON- ! STRAINTS WHOSE INDICES ARE IACT(K), K=1,2,...,NACT, ! ARE INCONSISTENT. ! MAXIT : MAXIMUM NUMBER OF ITERATIONS. ! VSMALL : REQUIRED ACCURACY TO BE ACHIEVED (E.G. IN THE ORDER OF THE ! MACHINE PRECISION FOR SMALL AND WELL-CONDITIONED PROBLEMS). ! DIAG : ON RETURN DIAG IS EQUAL TO THE MULTIPLE OF THE UNIT MATRIX ! THAT WAS ADDED TO G TO ACHIEVE POSITIVE DEFINITENESS. ! W(LW) : THE ELEMENTS OF W(.) ARE USED FOR WORKING SPACE. THE LENGTH ! OF W MUST NOT BE LESS THAN (1.5*NMAX*NMAX + 10*NMAX + M). ! WHEN INFO = 0 ON RETURN, THE LAGRANGE MULTIPLIERS OF THE ! FINAL ACTIVE CONSTRAINTS ARE HELD IN W(K), K=1,2,...,NACT. ! THE VALUES OF N, M, MEQ, MMAX, MN, MNN AND NMAX AND THE ELEMENTS OF ! A, B, GRAD AND G ARE NOT ALTERED. ! THE FOLLOWING INTEGERS ARE USED TO PARTITION W: ! THE FIRST N ELEMENTS OF W HOLD LAGRANGE MULTIPLIER ESTIMATES. ! W(IWZ+I+(N-1)*J) HOLDS THE MATRIX ELEMENT Z(I,J). ! W(IWR+I+0.5*J*(J-1)) HOLDS THE UPPER TRIANGULAR MATRIX ! ELEMENT R(I,J). THE SUBSEQUENT N COMPONENTS OF W MAY BE ! TREATED AS AN EXTRA COLUMN OF R(.,.). ! W(IWW-N+I) (I=1,2,...,N) ARE USED FOR TEMPORARY STORAGE. ! W(IWW+I) (I=1,2,...,N) ARE USED FOR TEMPORARY STORAGE. ! W(IWD+I) (I=1,2,...,N) HOLDS G(I,I) DURING THE CALCULATION. ! W(IWX+I) (I=1,2,...,N) HOLDS VARIABLES THAT WILL BE USED TO ! TEST THAT THE ITERATIONS INCREASE THE OBJECTIVE FUNCTION. ! W(IWA+K) (K=1,2,...,M) USUALLY HOLDS THE RECIPROCAL OF THE ! LENGTH OF THE K-TH CONSTRAINT, BUT ITS SIGN INDICATES ! WHETHER THE CONSTRAINT IS ACTIVE. ! AUTHOR: K. SCHITTKOWSKI, ! MATHEMATISCHES INSTITUT, ! UNIVERSITAET BAYREUTH, ! 8580 BAYREUTH, ! GERMANY, F.R. ! AUTHOR OF ORIGINAL VERSION: ! M.J.D. POWELL, DAMTP, ! UNIVERSITY OF CAMBRIDGE, SILVER STREET ! CAMBRIDGE, ! ENGLAND ! REFERENCE: M.J.D. POWELL: ZQPCVX, A FORTRAN SUBROUTINE FOR CONVEX ! PROGRAMMING, REPORT DAMTP/1983/NA17, UNIVERSITY OF ! CAMBRIDGE, ENGLAND, 1983. ! VERSION : 2.0 (MARCH, 1987) !************************************************************************* !subroutine ql0002(n,m,meq,mmax,mn,mnn,nmax,lql,a,b,grad,g, & ! xl,xu,x,nact,iact,maxit,vsmall,info,diag,w,lw) integer mmax,nmax,n,lw,nflag,iwwn,iact(n) _REAL_ a(mmax,n),b(mmax),grad(n),g(nmax,n),x(n), & w(lw),xl(n),xu(n),fmax,fmin integer m,meq,mn,mnn,nact,info,maxit _REAL_ cvmax,diag,diagr,fdiff,fdiffa,ga,gb,parinc,parnew & ,ratio,res,step,sum,sumx,sumy,suma,sumb,sumc,temp,tempa, & vsmall,xmag,xmagr,zero,one,two,onha,vfact ! INTRINSIC FUNCTIONS: DMAX1,sqrt,abs,fmin integer iwz,iwr,iww,iwd,iwa,ifinc,kfinc,k,i,ia,id,ii,ir,ira, & irb,j,nm,iz,iza,iterc,itref,jfinc,iflag,iws,is,k1,iw,kk,ip, & ipp,il,iu,ju,kflag,lflag,jflag,kdrop,nu,mflag,knext,ix,iwx, & iwy,iy,jl logical lql,lower ! INITIAL ADDRESSES iwz=nmax iwr=iwz+nmax*nmax iww=iwr+(nmax*(nmax+3))/2 iwd=iww+nmax iwx=iwd+nmax iwa=iwx+nmax ! SET SOME CONSTANTS. zero=0.d+0 one=1.d+0 two=2.d+0 onha=1.5d+0 vfact=1.d+0 ! SET SOME PARAMETERS. ! NUMBER LESS THAN VSMALL ARE ASSUMED TO BE NEGLIGIBLE. ! THE MULTIPLE OF I THAT IS ADDED TO G IS AT MOST DIAGR TIMES ! THE LEAST MULTIPLE OF I THAT GIVES POSITIVE DEFINITENESS. ! X IS RE-INITIALISED IF ITS MAGNITUDE IS REDUCED BY THE ! FACTOR XMAGR. ! A CHECK IS MADE FOR AN INCREASE IN F EVERY IFINC ITERATIONS, ! AFTER KFINC ITERATIONS ARE COMPLETED. diagr=two diag=zero xmagr=1.0d-2 ifinc=3 kfinc=max0(10,n) ! FIND THE RECIPROCALS OF THE LENGTHS OF THE CONSTRAINT NORMALS. ! RETURN IF A CONSTRAINT IS INFEASIBLE DUE TO A ZERO NORMAL. nact=0 if (m <= 0) goto 45 do k=1,m sum=zero do i=1,n sum=sum+a(k,i)**2 end do if (sum > zero) goto 20 if (b(k) == zero) goto 30 info=-k if (k <= meq) goto 730 if (b(k)) 30,30,730 20 sum=one/sqrt(sum) 30 ia=iwa+k w(ia)=sum end do 45 do k=1,n ia=iwa+m+k w(ia)=one end do ! IF NECESSARY INCREASE THE DIAGONAL ELEMENTS OF G. if (.not. lql) goto 165 do i=1,n id=iwd+i w(id)=g(i,i) diag=fmax(diag,vsmall-w(id)) if (i == n) cycle ii=i+1 do j=ii,n ga=-fmin(w(id),g(j,j)) gb=abs(w(id)-g(j,j))+abs(g(i,j)) if (gb > zero) ga=ga+g(i,j)**2/gb diag=fmax(diag,ga) end do end do if (diag <= zero) goto 90 70 diag=diagr*diag do i=1,n id=iwd+i g(i,i)=diag+w(id) end do ! FORM THE CHOLESKY FACTORISATION OF G. THE TRANSPOSE ! OF THE FACTOR WILL BE PLACED IN THE R-PARTITION OF W. 90 ir=iwr do j=1,n ira=iwr irb=ir+1 do i=1,j temp=g(i,j) if (i == 1) goto 110 do k=irb,ir ira=ira+1 temp=temp-w(k)*w(ira) end do 110 ir=ir+1 ira=ira+1 if (i < j) w(ir)=temp/w(ira) end do if (temp < vsmall) goto 140 w(ir)=sqrt(temp) end do goto 170 ! INCREASE FURTHER THE DIAGONAL ELEMENT OF G. 140 w(j)=one sumx=one k=j 150 sum=zero ira=ir-1 do i=k,j sum=sum-w(ira)*w(i) ira=ira+i end do ir=ir-k k=k-1 w(k)=sum/w(ir) sumx=sumx+w(k)**2 if (k >= 2) goto 150 diag=diag+vsmall-temp/sumx goto 70 ! STORE THE CHOLESKY FACTORISATION IN THE R-PARTITION ! OF W. 165 ir=iwr do i=1,n do j=1,i ir=ir+1 w(ir)=g(j,i) end do end do ! SET Z THE INVERSE OF THE MATRIX IN R. 170 nm=n-1 do i=1,n iz=iwz+i if (i == 1) goto 190 do j=2,i w(iz)=zero iz=iz+n end do 190 ir=iwr+(i+i*i)/2 w(iz)=one/w(ir) if (i == n) cycle iza=iz do j=i,nm ir=ir+i sum=zero do k=iza,iz,n sum=sum+w(k)*w(ir) ir=ir+1 end do iz=iz+n w(iz)=-sum/w(ir) end do end do ! SET THE INITIAL VALUES OF SOME VARIABLES. ! ITERC COUNTS THE NUMBER OF ITERATIONS. ! ITREF IS SET TO ONE WHEN ITERATIVE REFINEMENT IS REQUIRED. ! JFINC INDICATES WHEN TO TEST FOR AN INCREASE IN F. iterc=1 itref=0 jfinc=-kfinc ! SET X TO ZERO AND SET THE CORRESPONDING RESIDUALS OF THE ! KUHN-TUCKER CONDITIONS. 230 iflag=1 iws=iww-n do i=1,n x(i)=zero iw=iww+i w(iw)=grad(i) if (i > nact) cycle w(i)=zero is=iws+i k=iact(i) if (k <= m) goto 235 if (k > mn) goto 234 k1=k-m w(is)=xl(k1) cycle 234 k1=k-mn w(is)=-xu(k1) cycle 235 w(is)=b(k) end do xmag=zero vfact=1.d+0 if (nact) 340,340,280 ! SET THE RESIDUALS OF THE KUHN-TUCKER CONDITIONS FOR GENERAL X. 250 iflag=2 iws=iww-n do i=1,n iw=iww+i w(iw)=grad(i) if (lql) goto 259 id=iwd+i w(id)=zero do j=i,n w(id)=w(id)+g(i,j)*x(j) end do do j=1,i id=iwd+j w(iw)=w(iw)+g(j,i)*w(id) end do cycle 259 do j=1,n w(iw)=w(iw)+g(i,j)*x(j) end do end do if (nact == 0) goto 340 do k=1,nact kk=iact(k) is=iws+k if (kk > m) goto 265 w(is)=b(kk) do i=1,n iw=iww+i w(iw)=w(iw)-w(k)*a(kk,i) w(is)=w(is)-x(i)*a(kk,i) end do cycle 265 if (kk > mn) goto 266 k1=kk-m iw=iww+k1 w(iw)=w(iw)-w(k) w(is)=xl(k1)-x(k1) cycle 266 k1=kk-mn iw=iww+k1 w(iw)=w(iw)+w(k) w(is)=-xu(k1)+x(k1) end do ! PRE-MULTIPLY THE VECTOR IN THE S-PARTITION OF W BY THE ! INVERS OF R TRANSPOSE. 280 ir=iwr ip=iww+1 ipp=iww+n il=iws+1 iu=iws+nact do i=il,iu sum=zero if (i == il) goto 300 ju=i-1 do j=il,ju ir=ir+1 sum=sum+w(ir)*w(j) end do 300 ir=ir+1 w(i)=(w(i)-sum)/w(ir) end do ! SHIFT X TO SATISFY THE ACTIVE CONSTRAINTS AND MAKE THE ! CORRESPONDING CHANGE TO THE GRADIENT RESIDUALS. do i=1,n iz=iwz+i sum=zero do j=il,iu sum=sum+w(j)*w(iz) iz=iz+n end do x(i)=x(i)+sum if (lql) goto 329 id=iwd+i w(id)=zero do j=i,n w(id)=w(id)+g(i,j)*sum end do iw=iww+i do j=1,i id=iwd+j w(iw)=w(iw)+g(j,i)*w(id) end do cycle 329 do j=1,n iw=iww+j w(iw)=w(iw)+sum*g(i,j) end do end do ! FORM THE SCALAR PRODUCT OF THE CURRENT GRADIENT RESIDUALS ! WITH EACH COLUMN OF Z. 340 kflag=1 goto 930 350 if (nact == n) goto 380 ! SHIFT X SO THAT IT SATISFIES THE REMAINING KUHN-TUCKER ! CONDITIONS. il=iws+nact+1 iza=iwz+nact*n do i=1,n sum=zero iz=iza+i do j=il,iww sum=sum+w(iz)*w(j) iz=iz+n end do x(i)=x(i)-sum end do info=0 if (nact == 0) goto 410 ! UPDATE THE LAGRANGE MULTIPLIERS. 380 lflag=3 goto 740 390 do k=1,nact iw=iww+k w(k)=w(k)+w(iw) end do ! REVISE THE VALUES OF XMAG. ! BRANCH IF ITERATIVE REFINEMENT IS REQUIRED. 410 jflag=1 goto 910 420 if (iflag == itref) goto 250 ! DELETE A CONSTRAINT IF A LAGRANGE MULTIPLIER OF AN ! INEQUALITY CONSTRAINT IS NEGATIVE. kdrop=0 goto 440 430 kdrop=kdrop+1 if (w(kdrop) >= zero) goto 440 if (iact(kdrop) <= meq) goto 440 nu=nact mflag=1 goto 800 440 if (kdrop < nact) goto 430 ! SEEK THE GREATEAST NORMALISED CONSTRAINT VIOLATION, DISREGARDING ! ANY THAT MAY BE DUE TO COMPUTER ROUNDING ERRORS. 450 cvmax=zero if (m <= 0) goto 481 do k=1,m ia=iwa+k if (w(ia) <= zero) cycle sum=-b(k) do i=1,n sum=sum+x(i)*a(k,i) end do sumx=-sum*w(ia) if (k <= meq) sumx=abs(sumx) if (sumx <= cvmax) cycle temp=abs(b(k)) do i=1,n temp=temp+abs(x(i)*a(k,i)) end do tempa=temp+abs(sum) if (tempa <= temp) cycle temp=temp+onha*abs(sum) if (temp <= tempa) cycle cvmax=sumx res=sum knext=k end do 481 do k=1,n lower=.true. ia=iwa+m+k if (w(ia) <= zero) goto 485 sum=xl(k)-x(k) if (sum) 482,485,483 482 sum=x(k)-xu(k) lower=.false. 483 if (sum <= cvmax) goto 485 cvmax=sum res=-sum knext=k+m if (lower) goto 485 knext=k+mn end do 485 continue ! TEST FOR CONVERGENCE info=0 if (cvmax <= vsmall) goto 700 ! RETURN IF, DUE TO ROUNDING ERRORS, THE ACTUAL CHANGE IN ! X MAY NOT INCREASE THE OBJECTIVE FUNCTION jfinc=jfinc+1 if (jfinc == 0) goto 510 if (jfinc /= ifinc) goto 530 fdiff=zero fdiffa=zero do i=1,n sum=two*grad(i) sumx=abs(sum) if (lql) goto 489 id=iwd+i w(id)=zero do j=i,n ix=iwx+j w(id)=w(id)+g(i,j)*(w(ix)+x(j)) end do do j=1,i id=iwd+j temp=g(j,i)*w(id) sum=sum+temp sumx=sumx+abs(temp) end do goto 495 489 do j=1,n ix=iwx+j temp=g(i,j)*(w(ix)+x(j)) sum=sum+temp sumx=sumx+abs(temp) end do 495 ix=iwx+i fdiff=fdiff+sum*(x(i)-w(ix)) fdiffa=fdiffa+sumx*abs(x(i)-w(ix)) end do info=2 sum=fdiffa+fdiff if (sum <= fdiffa) goto 700 temp=fdiffa+onha*fdiff if (temp <= sum) goto 700 jfinc=0 info=0 510 do i=1,n ix=iwx+i w(ix)=x(i) end do ! FORM THE SCALAR PRODUCT OF THE NEW CONSTRAINT NORMAL WITH EACH ! COLUMN OF Z. PARNEW WILL BECOME THE LAGRANGE MULTIPLIER OF ! THE NEW CONSTRAINT. 530 iterc=iterc+1 if (iterc <= maxit) goto 531 info=1 goto 710 531 continue iws=iwr+(nact+nact*nact)/2 if (knext > m) goto 541 do i=1,n iw=iww+i w(iw)=a(knext,i) end do goto 549 541 do i=1,n iw=iww+i w(iw)=zero end do k1=knext-m if (k1 > n) goto 545 iw=iww+k1 w(iw)=one iz=iwz+k1 do i=1,n is=iws+i w(is)=w(iz) iz=iz+n end do goto 550 545 k1=knext-mn iw=iww+k1 w(iw)=-one iz=iwz+k1 do i=1,n is=iws+i w(is)=-w(iz) iz=iz+n end do goto 550 549 kflag=2 goto 930 550 parnew=zero ! APPLY GIVENS ROTATIONS TO MAKE THE LAST (N-NACT-2) SCALAR ! PRODUCTS EQUAL TO ZERO. if (nact == n) goto 570 nu=n nflag=1 goto 860 ! BRANCH IF THERE IS NO NEED TO DELETE A CONSTRAINT. 560 is=iws+nact if (nact == 0) goto 640 suma=zero sumb=zero sumc=zero iz=iwz+nact*n do i=1,n iz=iz+1 iw=iww+i suma=suma+w(iw)*w(iz) sumb=sumb+abs(w(iw)*w(iz)) sumc=sumc+w(iz)**2 end do temp=sumb+.1d+0*abs(suma) tempa=sumb+.2d+0*abs(suma) if (temp <= sumb) goto 570 if (tempa <= temp) goto 570 if (sumb > vsmall) goto 5 goto 570 5 sumc=sqrt(sumc) ia=iwa+knext if (knext <= m) sumc=sumc/w(ia) temp=sumc+.1d+0*abs(suma) tempa=sumc+.2d+0*abs(suma) if (temp <= sumc) goto 567 if (tempa <= temp) goto 567 goto 640 ! CALCULATE THE MULTIPLIERS FOR THE NEW CONSTRAINT NORMAL ! EXPRESSED IN TERMS OF THE ACTIVE CONSTRAINT NORMALS. ! THEN WORK OUT WHICH CONTRAINT TO DROP. 567 lflag=4 goto 740 570 lflag=1 goto 740 ! COMPLETE THE TEST FOR LINEARLY DEPENDENT CONSTRAINTS. 571 if (knext > m) goto 574 do i=1,n suma=a(knext,i) sumb=abs(suma) if (nact == 0) goto 581 do k=1,nact kk=iact(k) if (kk <= m) goto 568 kk=kk-m temp=zero if (kk == i) temp=w(iww+kk) kk=kk-n if (kk == i) temp=-w(iww+kk) goto 569 568 continue iw=iww+k temp=w(iw)*a(kk,i) 569 continue suma=suma-temp sumb=sumb+abs(temp) end do 581 if (suma <= vsmall) cycle temp=sumb+.1d+0*abs(suma) tempa=sumb+.2d+0*abs(suma) if (temp <= sumb) cycle if (tempa <= temp) cycle goto 630 end do lflag=1 goto 775 574 k1=knext-m if (k1 > n) k1=k1-n do i=1,n suma=zero if (i /= k1) goto 575 suma=one if (knext > mn) suma=-one 575 sumb=abs(suma) if (nact == 0) goto 582 do k=1,nact kk=iact(k) if (kk <= m) goto 579 kk=kk-m temp=zero if (kk == i) temp=w(iww+kk) kk=kk-n if (kk == i) temp=-w(iww+kk) goto 576 579 iw=iww+k temp=w(iw)*a(kk,i) 576 suma=suma-temp end do 577 sumb=sumb+abs(temp) 582 temp=sumb+.1d+0*abs(suma) tempa=sumb+.2d+0*abs(suma) if (temp <= sumb) cycle if (tempa <= temp) cycle goto 630 end do lflag=1 goto 775 ! BRANCH IF THE CONTRAINTS ARE INCONSISTENT. 580 info=-knext if (kdrop == 0) goto 700 parinc=ratio parnew=parinc ! REVISE THE LAGRANGE MULTIPLIERS OF THE ACTIVE CONSTRAINTS. 590 if (nact == 0) goto 601 do k=1,nact iw=iww+k w(k)=w(k)-parinc*w(iw) if (iact(k) > meq) w(k)=fmax(zero,w(k)) end do 601 if (kdrop == 0) goto 680 ! DELETE THE CONSTRAINT TO BE DROPPED. ! SHIFT THE VECTOR OF SCALAR PRODUCTS. ! THEN, IF APPROPRIATE, MAKE ONE MORE SCALAR PRODUCT ZERO. nu=nact+1 mflag=2 goto 800 610 iws=iws-nact-1 nu=min0(n,nu) do i=1,nu is=iws+i j=is+nact w(is)=w(j+1) end do nflag=2 goto 860 ! CALCULATE THE STEP TO THE VIOLATED CONSTRAINT. 630 is=iws+nact 640 sumy=w(is+1) step=-res/sumy parinc=step/sumy if (nact == 0) goto 660 ! CALCULATE THE CHANGES TO THE LAGRANGE MULTIPLIERS, AND REDUCE ! THE STEP ALONG THE NEW SEARCH DIRECTION IF NECESSARY. lflag=2 goto 740 650 if (kdrop == 0) goto 660 temp=one-ratio/parinc if (temp <= zero) kdrop=0 if (kdrop == 0) goto 660 step=ratio*sumy parinc=ratio res=temp*res ! UPDATE X AND THE LAGRANGE MULTIPIERS. ! DROP A CONSTRAINT IF THE FULL STEP IS NOT TAKEN. 660 iwy=iwz+nact*n do i=1,n iy=iwy+i x(i)=x(i)+step*w(iy) end do parnew=parnew+parinc if (nact >= 1) goto 590 ! ADD THE NEW CONSTRAINT TO THE ACTIVE SET. 680 nact=nact+1 w(nact)=parnew iact(nact)=knext ia=iwa+knext if (knext > mn) ia=ia-n w(ia)=-w(ia) ! ESTIMATE THE MAGNITUDE OF X. THEN BEGIN A NEW ITERATION, ! RE-INITILISING X IF THIS MAGNITUDE IS SMALL. jflag=2 goto 910 690 if (sum < (xmagr*xmag)) goto 230 if (itref) 450,450,250 ! INITIATE ITERATIVE REFINEMENT IF IT HAS NOT YET BEEN USED, ! OR RETURN AFTER RESTORING THE DIAGONAL ELEMENTS OF G. 700 if (iterc == 0) goto 710 itref=itref+1 jfinc=-1 if (itref == 1) goto 250 710 if (.not. lql) return do i=1,n id=iwd+i g(i,i)=w(id) end do 730 return ! THE REMAINING INSTRUCTIONS ARE USED AS SUBROUTINES. !******************************************************************** ! CALCULATE THE LAGRANGE MULTIPLIERS BY PRE-MULTIPLYING THE ! VECTOR IN THE S-PARTITION OF W BY THE INVERSE OF R. 740 ir=iwr+(nact+nact*nact)/2 i=nact sum=zero goto 770 750 ira=ir-1 sum=zero if (nact == 0) goto 761 do j=i,nact iw=iww+j sum=sum+w(ira)*w(iw) ira=ira+j end do 761 ir=ir-i i=i-1 770 iw=iww+i is=iws+i w(iw)=(w(is)-sum)/w(ir) if (i > 1) goto 750 if (lflag == 3) goto 390 if (lflag == 4) goto 571 ! CALCULATE THE NEXT CONSTRAINT TO DROP. 775 ip=iww+1 ipp=iww+nact kdrop=0 if (nact == 0) goto 791 do k=1,nact if (iact(k) <= meq) cycle iw=iww+k if ((res*w(iw)) >= zero) cycle temp=w(k)/w(iw) if (kdrop == 0) goto 780 if (abs(temp) >= abs(ratio)) cycle 780 kdrop=k ratio=temp end do 791 goto (580,650), lflag !******************************************************************** ! DROP THE CONSTRAINT IN POSITION KDROP IN THE ACTIVE SET. 800 ia=iwa+iact(kdrop) if (iact(kdrop) > mn) ia=ia-n w(ia)=-w(ia) if (kdrop == nact) goto 850 ! SET SOME INDICES AND CALCULATE THE ELEMENTS OF THE NEXT ! GIVENS ROTATION. iz=iwz+kdrop*n ir=iwr+(kdrop+kdrop*kdrop)/2 810 ira=ir ir=ir+kdrop+1 temp=fmax(abs(w(ir-1)),abs(w(ir))) sum=temp*sqrt((w(ir-1)/temp)**2+(w(ir)/temp)**2) ga=w(ir-1)/sum gb=w(ir)/sum ! EXCHANGE THE COLUMNS OF R. do i=1,kdrop ira=ira+1 j=ira-kdrop temp=w(ira) w(ira)=w(j) w(j)=temp end do w(ir)=zero ! APPLY THE ROTATION TO THE ROWS OF R. w(j)=sum kdrop=kdrop+1 do i=kdrop,nu temp=ga*w(ira)+gb*w(ira+1) w(ira+1)=ga*w(ira+1)-gb*w(ira) w(ira)=temp ira=ira+i end do ! APPLY THE ROTATION TO THE COLUMNS OF Z. do i=1,n iz=iz+1 j=iz-n temp=ga*w(j)+gb*w(iz) w(iz)=ga*w(iz)-gb*w(j) w(j)=temp end do ! REVISE IACT AND THE LAGRANGE MULTIPLIERS. iact(kdrop-1)=iact(kdrop) w(kdrop-1)=w(kdrop) if (kdrop < nact) goto 810 850 nact=nact-1 goto (250,610), mflag !******************************************************************** ! APPLY GIVENS ROTATION TO REDUCE SOME OF THE SCALAR ! PRODUCTS IN THE S-PARTITION OF W TO ZERO. 860 iz=iwz+nu*n 870 iz=iz-n 880 is=iws+nu nu=nu-1 if (nu == nact) goto 900 if (w(is) == zero) goto 870 temp=fmax(abs(w(is-1)),abs(w(is))) sum=temp*sqrt((w(is-1)/temp)**2+(w(is)/temp)**2) ga=w(is-1)/sum gb=w(is)/sum w(is-1)=sum do i=1,n k=iz+n temp=ga*w(iz)+gb*w(k) w(k)=ga*w(k)-gb*w(iz) w(iz)=temp iz=iz-1 end do goto 880 900 goto (560,630), nflag !******************************************************************** ! CALCULATE THE MAGNITUDE OF X AN REVISE XMAG. 910 sum=zero do i=1,n sum=sum+abs(x(i))*vfact*(abs(grad(i))+abs(g(i,i)*x(i))) if (lql) cycle if (sum < 1.d-30) cycle vfact=1.d-10*vfact sum=1.d-10*sum xmag=1.d-10*xmag end do 925 xmag=fmax(xmag,sum) goto (420,690), jflag !******************************************************************** ! PRE-MULTIPLY THE VECTOR IN THE W-PARTITION OF W BY Z TRANSPOSE. 930 jl=iww+1 iz=iwz do i=1,n is=iws+i w(is)=zero iwwn=iww+n do j=jl,iwwn iz=iz+1 w(is)=w(is)+w(iz)*w(j) end do end do goto (350,550), kflag return end subroutine ql0002