This subroutine solves bound constrained optimization problems by using the compact formula of the limited memory BFGS updates.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n |
the number of variables. |
||
| integer, | intent(in) | :: | m |
the maximum number of variable metric corrections allowed in the limited memory matrix. |
||
| real(kind=wp), | intent(inout) | :: | x(n) |
On entry x is an approximation to the solution. On exit x is the current approximation. |
||
| real(kind=wp), | intent(in) | :: | l(n) |
the lower bound of x. |
||
| real(kind=wp), | intent(in) | :: | u(n) |
the upper bound of x. |
||
| integer, | intent(in) | :: | Nbd(n) |
the type of bounds imposed on the variables, and must be specified as follows:
|
||
| real(kind=wp), | intent(inout) | :: | f |
On first entry f is unspecified. On final exit f is the value of the function at x. |
||
| real(kind=wp), | intent(inout) | :: | g(n) |
On first entry g is unspecified. On final exit g is the value of the gradient at x. |
||
| real(kind=wp), | intent(in) | :: | Factr |
On entry
where epsmch is the machine precision, which is automatically generated by the code. |
||
| real(kind=wp), | intent(in) | :: | Pgtol |
On entry pgtol >= 0 is specified by the user. The iteration will stop when
where pg_i is the ith component of the projected gradient. |
||
| real(kind=wp) | :: | Ws(n,m) |
working array used to store information defining the limited memory BFGS matrix: stores S, the matrix of s-vectors; |
|||
| real(kind=wp) | :: | Wy(n,m) |
working array used to store information defining the limited memory BFGS matrix: stores Y, the matrix of y-vectors; |
|||
| real(kind=wp) | :: | Sy(m,m) |
working array used to store information defining the limited memory BFGS matrix: stores S'Y; |
|||
| real(kind=wp) | :: | Ss(m,m) |
working array used to store information defining the limited memory BFGS matrix: stores S'S; |
|||
| real(kind=wp) | :: | Wt(m,m) |
working array used to store information defining the limited memory BFGS matrix: stores the Cholesky factorization of (theta*S'S+LD^(-1)L'); see eq. (2.26) in [3]. |
|||
| real(kind=wp) | :: | Wn(2*m,2*m) |
working array used to store the LEL^T factorization of the indefinite matrix where |
|||
| real(kind=wp) | :: | Snd(2*m,2*m) |
working array used to store the lower triangular part of |
|||
| real(kind=wp) | :: | z(n) |
working array. used at different times to store the Cauchy point and the Newton point. |
|||
| real(kind=wp) | :: | r(n) |
working array. |
|||
| real(kind=wp) | :: | d(n) |
working array. |
|||
| real(kind=wp) | :: | t(n) |
working array. |
|||
| real(kind=wp) | :: | Xp(n) |
working array. used to safeguard the projected Newton direction |
|||
| real(kind=wp) | :: | Wa(8*m) |
working array. |
|||
| integer | :: | Index(n) |
working array. In subroutine freev, index is used to store the free and fixed variables at the Generalized Cauchy Point (GCP). |
|||
| integer | :: | Iwhere(n) |
working array used to record the status of the vector x for GCP computation:
|
|||
| integer | :: | Indx2(n) |
integer working array. Within subroutine cauchy, indx2 corresponds to the array iorder. In subroutine freev, a list of variables entering and leaving the free set is stored in indx2, and it is passed on to subroutine formk with this information. |
|||
| character(len=60), | intent(inout) | :: | Task |
indicates the current job when entering and leaving this subroutine. |
||
| integer, | intent(in) | :: | Iprint |
Controls the frequency and type of output generated:
When |
||
| character(len=60) | :: | Csave |
working string |
|||
| logical | :: | Lsave(4) |
working array |
|||
| integer | :: | Isave(23) |
working array |
|||
| real(kind=wp) | :: | Dsave(29) |
working array |
|||
| character(len=*), | intent(in), | optional | :: | iteration_file |
The name of the iteration file if used ( |
subroutine mainlb(n,m,x,l,u,Nbd,f,g,Factr,Pgtol,Ws,Wy,Sy,Ss,Wt,Wn, & Snd,z,r,d,t,Xp,Wa,Index,Iwhere,Indx2,Task, & Iprint,Csave,Lsave,Isave,Dsave,iteration_file) implicit none character(len=60),intent(inout) :: Task !! indicates the current job when !! entering and leaving this subroutine. character(len=60) :: Csave !! working string logical :: Lsave(4) !! working array integer :: Isave(23) !! working array real(wp) :: Dsave(29) !! working array integer,intent(in) :: n !! the number of variables. integer,intent(in) :: m !! the maximum number of variable metric !! corrections allowed in the limited memory matrix. integer,intent(in) :: Iprint !! Controls the frequency and type of output generated: !! !! * `iprint<0 ` no output is generated; !! * `iprint=0 ` print only one line at the last iteration; !! * `0<iprint<99` print also f and |proj g| every iprint iterations; !! * `iprint=99 ` print details of every iteration except n-vectors; !! * `iprint=100 ` print also the changes of active set and final x; !! * `iprint>100 ` print details of every iteration including x and g; !! !! When `iprint > 0`, the file `iterate.dat` will be created to !! summarize the iteration. integer,intent(in) :: Nbd(n) !! the type of bounds imposed on the !! variables, and must be specified as follows: !! !! * `nbd(i)=0` if `x(i)` is unbounded, !! * `nbd(i)=1` if `x(i)` has only a lower bound, !! * `nbd(i)=2` if `x(i)` has both lower and upper bounds, !! * `nbd(i)=3` if `x(i)` has only an upper bound. integer :: Index(n) !! working array. !! In subroutine [[freev]], index is used to store the free and fixed !! variables at the Generalized Cauchy Point (GCP). integer :: Iwhere(n) !! working array used to record !! the status of the vector x for GCP computation: !! !! * `iwhere(i)= 0 or -3` if `x(i)` is free and has bounds, !! * `iwhere(i)= 1 ` if `x(i)` is fixed at l(i), and l(i) /= u(i) !! * `iwhere(i)= 2 ` if `x(i)` is fixed at u(i), and u(i) /= l(i) !! * `iwhere(i)= 3 ` if `x(i)` is always fixed, i.e., u(i)=x(i)=l(i) !! * `iwhere(i)=-1 ` if `x(i)` is always free, i.e., no bounds on it. integer :: Indx2(n) !! integer working array. !! Within subroutine [[cauchy]], indx2 corresponds to the array iorder. !! In subroutine [[freev]], a list of variables entering and leaving !! the free set is stored in indx2, and it is passed on to !! subroutine [[formk]] with this information. real(wp),intent(inout) :: f !! On first entry f is unspecified. !! On final exit f is the value of the function at x. real(wp),intent(in) :: Factr !! On entry `factr >= 0` is specified by the user. The iteration !! will stop when !! !! `(f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch` !! !! where epsmch is the machine precision, which is automatically !! generated by the code. real(wp),intent(in) :: Pgtol !! On entry pgtol >= 0 is specified by the user. The iteration !! will stop when !! !! `max{|proj g_i | i = 1, ..., n} <= pgtol` !! !! where pg_i is the ith component of the projected gradient. real(wp),intent(inout) :: x(n) !! On entry x is an approximation to the solution. !! On exit x is the current approximation. real(wp),intent(in) :: l(n) !! the lower bound of x. real(wp),intent(in) :: u(n) !! the upper bound of x. real(wp),intent(inout) :: g(n) !! On first entry g is unspecified. !! On final exit g is the value of the gradient at x. real(wp) :: z(n) !! working array. !! used at different times to store the Cauchy point and the Newton point. real(wp) :: r(n) !! working array. real(wp) :: d(n) !! working array. real(wp) :: t(n) !! working array. real(wp) :: Xp(n) !! working array. !! used to safeguard the projected Newton direction real(wp) :: Wa(8*m) !! working array. real(wp) :: Ws(n,m) !! working array used to store information defining the limited memory BFGS matrix: stores S, the matrix of s-vectors; real(wp) :: Wy(n,m) !! working array used to store information defining the limited memory BFGS matrix: stores Y, the matrix of y-vectors; real(wp) :: Sy(m,m) !! working array used to store information defining the limited memory BFGS matrix: stores S'Y; real(wp) :: Ss(m,m) !! working array used to store information defining the limited memory BFGS matrix: stores S'S; real(wp) :: Wt(m,m) !! working array used to store information defining the limited memory BFGS matrix: stores the Cholesky factorization !! of (theta*S'S+LD^(-1)L'); see eq. (2.26) in [3]. real(wp) :: Wn(2*m,2*m) !! working array !! used to store the LEL^T factorization of the indefinite matrix !!``` !! K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] !! [L_a -R_z theta*S'AA'S ] !!``` !! where !!``` !! E = [-I 0] !! [ 0 I] !!``` real(wp) :: Snd(2*m,2*m) !! working array !! used to store the lower triangular part of !!``` !! N = [Y' ZZ'Y L_a'+R_z'] !! [L_a +R_z S'AA'S ] !!``` character(len=*),intent(in),optional :: iteration_file !! The name of the iteration file if used (`Iprint>=1`). !! If not specified, then 'iterate.dat' is used. logical :: prjctd , cnstnd , boxed , updatd , wrk character(len=3) :: word integer :: i , k , nintol , itfile , iback , nskip , head , col , & iter , itail , iupdat , nseg , nfgv , info , ifun , & iword , nfree , nact , ileave , nenter real(wp) :: theta , fold , dr , rr , tol , xstep , & sbgnrm , ddum , dnorm , dtd , epsmch , cpu1 , & cpu2 , cachyt , sbtime , lnscht , time1 , time2 , & gd , gdold , stp , stpmx , time logical :: compute_infinity_norm_of_projected_gradient logical :: prelims logical :: linesearch if ( Task=='START' ) then epsmch = epsilon(one) call cpu_time(time1) ! Initialize counters and scalars when task='START'. ! for the limited memory BFGS matrices: col = 0 head = 1 theta = one iupdat = 0 updatd = .false. iback = 0 itail = 0 iword = 0 nact = 0 ileave = 0 nenter = 0 fold = zero dnorm = zero cpu1 = zero gd = zero stpmx = zero sbgnrm = zero stp = zero gdold = zero dtd = zero ! for operation counts: iter = 0 nfgv = 0 nseg = 0 nintol = 0 nskip = 0 nfree = n ifun = 0 ! for stopping tolerance: tol = Factr*epsmch ! for measuring running time: cachyt = 0 sbtime = 0 lnscht = 0 ! 'word' records the status of subspace solutions. word = '---' ! 'info' records the termination information. info = 0 ! open a summary file 'iterate.dat' if ( Iprint>=1 ) then if (present(iteration_file)) then open (newunit=itfile,file=trim(iteration_file),status='unknown') else ! use default name if not specified open (newunit=itfile,file='iterate.dat',status='unknown') end if end if ! Check the input arguments for errors. call errclb(n,m,Factr,l,u,Nbd,Task,info,k) if ( Task(1:5)=='ERROR' ) then call prn3lb(n,x,f,Task,Iprint,info,itfile,iter,nfgv,nintol, & nskip,nact,sbgnrm,zero,nseg,word,iback,stp, & xstep,k,cachyt,sbtime,lnscht) return endif call prn1lb(n,m,l,u,x,Iprint,itfile,epsmch) ! Initialize iwhere & project x onto the feasible set. call active(n,l,u,Nbd,x,Iwhere,Iprint,prjctd,cnstnd,boxed) ! The end of the initialization. call start() return end if ! restore local variables. prjctd = Lsave(1) cnstnd = Lsave(2) boxed = Lsave(3) updatd = Lsave(4) nintol = Isave(1) itfile = Isave(3) iback = Isave(4) nskip = Isave(5) head = Isave(6) col = Isave(7) itail = Isave(8) iter = Isave(9) iupdat = Isave(10) nseg = Isave(12) nfgv = Isave(13) info = Isave(14) ifun = Isave(15) iword = Isave(16) nfree = Isave(17) nact = Isave(18) ileave = Isave(19) nenter = Isave(20) theta = Dsave(1) fold = Dsave(2) tol = Dsave(3) dnorm = Dsave(4) epsmch = Dsave(5) cpu1 = Dsave(6) cachyt = Dsave(7) sbtime = Dsave(8) lnscht = Dsave(9) time1 = Dsave(10) gd = Dsave(11) stpmx = Dsave(12) sbgnrm = Dsave(13) stp = Dsave(14) gdold = Dsave(15) dtd = Dsave(16) ! After returning from the driver go to the point ! where execution is to resume. compute_infinity_norm_of_projected_gradient = .true. prelims = .true. linesearch = .true. if ( Task(1:5)=='FG_LN' ) then compute_infinity_norm_of_projected_gradient = .false. prelims = .false. else if ( Task(1:5)=='NEW_X' ) then compute_infinity_norm_of_projected_gradient = .false. prelims = .false. linesearch = .false. else if ( Task(1:5)/='FG_ST' ) then if ( Task(1:4)=='STOP' ) then if ( Task(7:9)=='CPU' ) then ! restore the previous iterate. call dcopy(n,t,1,x,1) call dcopy(n,r,1,g,1) f = fold endif call finish() else call start() endif return end if if (compute_infinity_norm_of_projected_gradient) then ! Compute the infinity norm of the (-) projected gradient. nfgv = 1 call projgr(n,l,u,Nbd,x,g,sbgnrm) if ( Iprint>=1 ) then write (output_unit,'(/,a,i5,4x,a,1p,d12.5,4x,a,1p,d12.5)') & 'At iterate', iter , 'f= ', f , '|proj g|= ', sbgnrm write (itfile,'(2(1x,i4),5x,a,5x,a,3x,a,5x,a,5x,a,8x,a,3x,1p,2(1x,d10.3))') & iter , nfgv , '-', '-', '-', '-', '-', '-', sbgnrm , f endif if ( sbgnrm<=Pgtol ) then ! terminate the algorithm. Task = 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL' call finish() return endif end if ! ----------------- the beginning of the loop -------------------------- main_loop : do if (prelims) then if ( Iprint>=99 ) write (output_unit,'(//,A,i5)') 'ITERATION ', iter + 1 iword = -1 if ( .not.cnstnd .and. col>0 ) then ! skip the search for GCP. call dcopy(n,x,1,z,1) wrk = updatd nseg = 0 else ! Compute the Generalized Cauchy Point (GCP). call cpu_time(cpu1) call cauchy(n,x,l,u,Nbd,g,Indx2,Iwhere,t,d,z,m,Wy,Ws,Sy,Wt,theta, & col,head,Wa(1),Wa(2*m+1),Wa(4*m+1),Wa(6*m+1),nseg, & Iprint,sbgnrm,info,epsmch) if ( info/=0 ) then ! singular triangular system detected; refresh the lbfgs memory. if ( Iprint>=1 ) write (output_unit,'(/,A,/,A)') & ' Singular triangular system detected;',& ' refresh the lbfgs memory and restart the iteration.' info = 0 col = 0 head = 1 theta = one iupdat = 0 updatd = .false. call cpu_time(cpu2) cachyt = cachyt + cpu2 - cpu1 call continue_loop() cycle main_loop endif call cpu_time(cpu2) cachyt = cachyt + cpu2 - cpu1 nintol = nintol + nseg ! Count the entering and leaving variables for iter > 0; ! find the index set of free and active variables at the GCP. call freev(n,nfree,Index,nenter,ileave,Indx2,Iwhere,wrk,updatd, & cnstnd,Iprint,iter) nact = n - nfree end if if ( nfree==0 .or. col==0 ) then ! If there are no free variables or B=theta*I, then ! skip the subspace minimization. else ! Subspace minimization. call cpu_time(cpu1) ! Form the LEL^T factorization of the indefinite ! matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] ! [L_a -R_z theta*S'AA'S ] ! where E = [-I 0] ! [ 0 I] if ( wrk ) call formk(n,nfree,Index,nenter,ileave,Indx2,iupdat, & updatd,Wn,Snd,m,Ws,Wy,Sy,theta,col,head, & info) if ( info/=0 ) then ! nonpositive definiteness in Cholesky factorization; ! refresh the lbfgs memory and restart the iteration. if ( Iprint>=1 ) write (output_unit,'(/,a,/,a)') & ' Nonpositive definiteness in Cholesky factorization in formk;',& ' refresh the lbfgs memory and restart the iteration.' info = 0 col = 0 head = 1 theta = one iupdat = 0 updatd = .false. call cpu_time(cpu2) sbtime = sbtime + cpu2 - cpu1 call continue_loop() cycle main_loop endif ! compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) ! from 'cauchy'). call cmprlb(n,m,x,g,Ws,Wy,Sy,Wt,z,r,Wa,Index,theta,col,head,nfree,& cnstnd,info) if ( info==0 ) then !-jlm-jn call the direct method. call subsm(n,m,nfree,Index,l,u,Nbd,z,r,Xp,Ws,Wy,theta,x,g,col, & head,iword,Wa,Wn,Iprint,info) end if if ( info/=0 ) then ! singular triangular system detected; ! refresh the lbfgs memory and restart the iteration. if ( Iprint>=1 ) write (output_unit,'(/,A,/,A)') & ' Singular triangular system detected;',& ' refresh the lbfgs memory and restart the iteration.' info = 0 col = 0 head = 1 theta = one iupdat = 0 updatd = .false. call cpu_time(cpu2) sbtime = sbtime + cpu2 - cpu1 call continue_loop() cycle main_loop endif call cpu_time(cpu2) sbtime = sbtime + cpu2 - cpu1 end if ! Line search and optimality tests. ! Generate the search direction d:=z-x. do i = 1 , n d(i) = z(i) - x(i) enddo call cpu_time(cpu1) end if ! ------------------------------------ if (linesearch) then call lnsrlb(n,l,u,Nbd,x,f,fold,gd,gdold,g,d,r,t,z,stp,dnorm,dtd, & xstep,stpmx,iter,ifun,iback,nfgv,info,Task,boxed, & cnstnd,Csave,Isave(22),Dsave(17)) if ( info/=0 .or. iback>=20 ) then ! restore the previous iterate. call dcopy(n,t,1,x,1) call dcopy(n,r,1,g,1) f = fold if ( col==0 ) then ! abnormal termination. if ( info==0 ) then info = -9 ! restore the actual number of f and g evaluations etc. nfgv = nfgv - 1 ifun = ifun - 1 iback = iback - 1 endif Task = 'ABNORMAL_TERMINATION_IN_LNSRCH' iter = iter + 1 call finish() return else ! refresh the lbfgs memory and restart the iteration. if ( Iprint>=1 ) write (output_unit,'(/,a,/,a)') & ' Bad direction in the line search;',& ' refresh the lbfgs memory and restart the iteration.' if ( info==0 ) nfgv = nfgv - 1 info = 0 col = 0 head = 1 theta = one iupdat = 0 updatd = .false. Task = 'RESTART_FROM_LNSRCH' call cpu_time(cpu2) lnscht = lnscht + cpu2 - cpu1 call continue_loop() cycle main_loop endif else if ( Task(1:5)=='FG_LN' ) then ! return to the driver for calculating f and g; reenter at 666. call save_locals() return else ! calculate and print out the quantities related to the new X. call cpu_time(cpu2) lnscht = lnscht + cpu2 - cpu1 iter = iter + 1 ! Compute the infinity norm of the projected (-)gradient. call projgr(n,l,u,Nbd,x,g,sbgnrm) ! Print iteration information. call prn2lb(n,x,f,g,Iprint,itfile,iter,nfgv,nact,sbgnrm,nseg, & word,iword,iback,stp,xstep) call save_locals() return endif end if ! ------------------------------------ ! Test for termination. if ( sbgnrm<=Pgtol ) then ! terminate the algorithm. Task = 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL' call finish() return endif ddum = max(abs(fold),abs(f),one) if ( (fold-f)<=tol*ddum ) then ! terminate the algorithm. Task = 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH' if ( iback>=10 ) info = -5 ! i.e., to issue a warning if iback>10 in the line search. call finish() return endif ! Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. do i = 1 , n r(i) = g(i) - r(i) enddo rr = ddot(n,r,1,r,1) if ( stp==one ) then dr = gd - gdold ddum = -gdold else dr = (gd-gdold)*stp call dscal(n,stp,d,1) ddum = -gdold*stp endif if ( dr<=epsmch*ddum ) then ! skip the L-BFGS update. nskip = nskip + 1 updatd = .false. if ( Iprint>=1 ) write (output_unit,'(a,1p,e10.3,a,1p,e10.3,a)') & ' ys=', dr,' -gs=' , ddum, ' BFGS update SKIPPED' call continue_loop() cycle main_loop endif ! Update the L-BFGS matrix. updatd = .true. iupdat = iupdat + 1 ! Update matrices WS and WY and form the middle matrix in B. call matupd(n,m,Ws,Wy,Sy,Ss,d,r,itail,iupdat,col,head,theta,rr,dr,& stp,dtd) ! Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; ! Store T in the upper triangular of the array wt; ! Cholesky factorize T to J*J' with ! J' stored in the upper triangular of wt. call formt(m,Wt,Sy,Ss,col,theta,info) if ( info/=0 ) then ! nonpositive definiteness in Cholesky factorization; ! refresh the lbfgs memory and restart the iteration. if ( Iprint>=1 ) write (output_unit,'(/,a,/,a)') & ' Nonpositive definiteness in Cholesky factorization in formt;',& ' refresh the lbfgs memory and restart the iteration.' info = 0 col = 0 head = 1 theta = one iupdat = 0 updatd = .false. endif ! Now the inverse of the middle matrix in B is ! [ D^(1/2) O ] [ -D^(1/2) D^(-1/2)*L' ] ! [ -L*D^(-1/2) J ] [ 0 J' ] call continue_loop() end do main_loop contains subroutine continue_loop() !! prepare for next loop iteration prelims = .true. linesearch = .true. end subroutine continue_loop subroutine start() !! return to the driver to calculate f and g Task = 'FG_START' call save_locals() end subroutine start subroutine finish() !! before returning call cpu_time(time2) time = time2 - time1 call prn3lb(n,x,f,Task,Iprint,info,itfile,iter,nfgv,nintol,nskip, & nact,sbgnrm,time,nseg,word,iback,stp,xstep,k,cachyt, & sbtime,lnscht) call save_locals() end subroutine finish subroutine save_locals() !! Save local variables. Lsave(1) = prjctd Lsave(2) = cnstnd Lsave(3) = boxed Lsave(4) = updatd Isave(1) = nintol Isave(3) = itfile Isave(4) = iback Isave(5) = nskip Isave(6) = head Isave(7) = col Isave(8) = itail Isave(9) = iter Isave(10) = iupdat Isave(12) = nseg Isave(13) = nfgv Isave(14) = info Isave(15) = ifun Isave(16) = iword Isave(17) = nfree Isave(18) = nact Isave(19) = ileave Isave(20) = nenter Dsave(1) = theta Dsave(2) = fold Dsave(3) = tol Dsave(4) = dnorm Dsave(5) = epsmch Dsave(6) = cpu1 Dsave(7) = cachyt Dsave(8) = sbtime Dsave(9) = lnscht Dsave(10) = time1 Dsave(11) = gd Dsave(12) = stpmx Dsave(13) = sbgnrm Dsave(14) = stp Dsave(15) = gdold Dsave(16) = dtd end subroutine save_locals end subroutine mainlb