# To unbundle, sh this file echo Index 1>&2 sed 's/.//' >Index <<'EOF Index' -file lbfgs_bcm -for bound constrained optimization problems -alg limited memory BFGS method -by Ciyou Zhu in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal -contact ciyou@eecs.nwu.edu or nocedal@eecs.nwu.edu -ref R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited - memory algorithm for bound constrained optimization'', - SIAM J. Scientific Computing 16 (1995), no. 5. - C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a limited - memory FORTRAN code for solving bound constrained optimization - problems'', Tech. Report, EECS Department, - Northwestern University, 1994. - (Postscript files of these papers are available via anonymous - ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) -size 224kB -prec double - EOF Index echo Readme 1>&2 sed 's/.//' >Readme <<'EOF Readme' - - L-BFGS-B - -This directory contains the limited memory -code for solving bound constrained -optimization problems. - -To use the code copy the following files: - -Makefile (for Unix users) -driver1.f -driver2.f -driver3.f -routines.f - -We recommend that the user read driver1.f, which -gives a good idea of how the code works as well as -directions on how to adapt the timing routine to particular -systems. - -The following three articles are also included - -[1] algorithm.ps - describes the algorithm -[2] compact.ps - presents the compact form of matrices -[3] code.ps - describes the code - -The most useful article for those who only wish to use the code is [3]. -Users interested in understanding the algorithm should read [1] and -possibly also [2]. - -_______________________________________________________________________ - -For questions and help contact - - Ciyou Zhu or Jorge Nocedal - ciyou@eecs.nwu.edu nocedal@eecs.nwu.edu - -_______________________________________________________________________ - - -To run the code on a Unix machine simply -type - make - -This will create (among others) the file -x.lbfgsb1* which is the executable file for the -simplest driver in the package. One can execute -the program by typing - - x.lbfgsb1 - - -************************************************************************ - How to use L-BFGS-B -************************************************************************ - - The simplest way to use the code is to modify one of the -drivers provided in the package. Most users will only need to make -a few changes to the drivers to run their applications. - - L-BFGS-B is written in FORTRAN 77, in double precision. The -user is required to calculate the function value f and its gradient g. -In order to allow the user complete control over these computations, -reverse communication is used. The routine setulb.f must be called -repeatedly under the control of the variable task. The calling -statement of L-BFGS-B is - - call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint, - + csave,lsave,isave,dsave) - - - Following is a description of all the parameters used in this call. - -c n is an INTEGER variable that must be set by the user to the -c number of variables. It is not altered by the routine. -c -c m is an INTEGER variable that must be set by the user to the -c number of corrections used in the limited memory matrix. -c It is not altered by the routine. Values of m < 3 are -c not recommended, and large values of m can result in excessive -c computing time. The range 3 <= m <= 20 is recommended. -c -c x is a DOUBLE PRECISION array of length n. On initial entry -c it must be set by the user to the values of the initial -c estimate of the solution vector. Upon successful exit, it -c contains the values of the variables at the best point -c found (usually an approximate solution). -c -c l is a DOUBLE PRECISION array of length n that must be set by -c the user to the values of the lower bounds on the variables. If -c the i-th variable has no lower bound, l(i) need not be defined. -c -c u is a DOUBLE PRECISION array of length n that must be set by -c the user to the values of the upper bounds on the variables. If -c the i-th variable has no upper bound, u(i) need not be defined. -c -c nbd is an INTEGER array of dimension n that must be set by the -c user to the type of bounds imposed on the variables: -c nbd(i)=0 if x(i) is unbounded, -c 1 if x(i) has only a lower bound, -c 2 if x(i) has both lower and upper bounds, -c 3 if x(i) has only an upper bound. -c -c f is a DOUBLE PRECISION variable. If the routine setulb returns -c with task(1:2)= 'FG', then f must be set by the user to -c contain the value of the function at the point x. -c -c g is a DOUBLE PRECISION array of length n. If the routine setulb -c returns with taskb(1:2)= 'FG', then g must be set by the user to -c contain the components of the gradient at the point x. -c -c factr is a DOUBLE PRECISION variable that must be set by the user. -c It is a tolerance in the termination test for the algorithm. -c The iteration will stop when -c -c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch -c -c where epsmch is the machine precision which is automatically -c generated by the code. Typical values for factr on a computer -c with 15 digits of accuracy in double precision are: -c factr=1.d+12 for low accuracy; -c 1.d+7 for moderate accuracy; -c 1.d+1 for extremely high accuracy. -c The user can suppress this termination test by setting factr=0. -c -c pgtol is a double precision variable. -c On entry pgtol >= 0 is specified by the user. The iteration -c will stop when -c -c max{|proj g_i | i = 1, ..., n} <= pgtol -c -c where pg_i is the ith component of the projected gradient. -c The user can suppress this termination test by setting pgtol=0. -c -c wa is a DOUBLE PRECISION array of length -c (2mmax + 4)nmax + 12mmax^2 + 12mmax used as workspace. -c This array must not be altered by the user. -c -c iwa is an INTEGER array of length 3nmax used as -c workspace. This array must not be altered by the user. -c -c task is a CHARACTER string of length 60. -c On first entry, it must be set to 'START'. -c On a return with task(1:2)='FG', the user must evaluate the -c function f and gradient g at the returned value of x. -c On a return with task(1:5)='NEW_X', an iteration of the -c algorithm has concluded, and f and g contain f(x) and g(x) -c respectively. The user can decide whether to continue or stop -c the iteration. -c When -c task(1:4)='CONV', the termination test in L-BFGS-B has been -c satisfied; -c task(1:4)='ABNO', the routine has terminated abnormally -c without being able to satisfy the termination conditions, -c x contains the best approximation found, -c f and g contain f(x) and g(x) respectively; -c task(1:5)='ERROR', the routine has detected an error in the -c input parameters; -c On exit with task = 'CONV', 'ABNO' or 'ERROR', the variable task -c contains additional information that the user can print. -c This array should not be altered unless the user wants to -c stop the run for some reason. See driver2 or driver3 -c for a detailed explanation on how to stop the run -c by assigning task(1:4)='STOP' in the driver. -c -c iprint is an INTEGER variable that must be set by the user. -c It controls the frequency and type of output generated: -c iprint<0 no output is generated; -c iprint=0 print only one line at the last iteration; -c 0100 print details of every iteration including x and g; -c When iprint > 0, the file iterate.dat will be created to -c summarize the iteration. -c -c csave is a CHARACTER working array of length 60. -c -c lsave is a LOGICAL working array of dimension 4. -c On exit with task = 'NEW_X', the following information is -c available: -c lsave(1) = .true. the initial x did not satisfy the bounds; -c lsave(2) = .true. the problem contains bounds; -c lsave(3) = .true. each variable has upper and lower bounds. -c -c isave is an INTEGER working array of dimension 44. -c On exit with task = 'NEW_X', it contains information that -c the user may want to access: -c isave(30) = the current iteration number; -c isave(34) = the total number of function and gradient -c evaluations; -c isave(36) = the number of function value or gradient -c evaluations in the current iteration; -c isave(38) = the number of free variables in the current -c iteration; -c isave(39) = the number of active constraints at the current -c iteration; -c -c see the subroutine setulb.f for a description of other -c information contained in isave -c -c dsave is a DOUBLE PRECISION working array of dimension 29. -c On exit with task = 'NEW_X', it contains information that -c the user may want to access: -c dsave(2) = the value of f at the previous iteration; -c dsave(5) = the machine precision epsmch generated by the code; -c dsave(13) = the infinity norm of the projected gradient; -c -c see the subroutine setulb.f for a description of other -c information contained in dsave - - EOF Readme echo Makefile 1>&2 sed 's/.//' >Makefile <<'EOF Makefile' -FC = f77 - -FFLAGS = -O - -DRIVER1 = driver1.o -DRIVER2 = driver2.o -DRIVER3 = driver3.o - -ROUTINES = routines.o - -all : lbfgsb1 lbfgsb2 lbfgsb3 - -lbfgsb1 : $(DRIVER1) $(ROUTINES) - $(FC) $(FFLAGS) $(DRIVER1) $(ROUTINES) -o x.lbfgsb1 - -lbfgsb2 : $(DRIVER2) $(ROUTINES) - $(FC) $(FFLAGS) $(DRIVER2) $(ROUTINES) -o x.lbfgsb2 - -lbfgsb3 : $(DRIVER3) $(ROUTINES) - $(FC) $(FFLAGS) $(DRIVER3) $(ROUTINES) -o x.lbfgsb3 - EOF Makefile echo driver1.f 1>&2 sed 's/.//' >driver1.f <<'EOF driver1.f' -c DRIVER 1 -c -------------------------------------------------------------- -c SIMPLE DRIVER FOR L-BFGS-B -c -------------------------------------------------------------- -c -c L-BFGS-B is a code for solving large nonlinear optimization -c problems with simple bounds on the variables. -c -c The code can also be used for unconstrained problems and is -c as efficient for these problems as the earlier limited memory -c code L-BFGS. -c -c This is the simplest driver in the package. It uses all the -c default settings of the code. -c -c -c References: -c -c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited -c memory algorithm for bound constrained optimization'', -c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. -c -c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN -c Subroutines for Large Scale Bound Constrained Optimization'' -c Tech. Report, NAM-11, EECS Department, Northwestern University, -c 1994. -c -c -c (Postscript files of these papers are available via anonymous -c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c NOTE: The user should adapt the subroutine 'timer' if 'etime' is -c not available on the system. An example for system -c AIX Version 3.2 is available at the end of this driver. -c -c ************** - - program driver - -c This simple driver demonstrates how to call the L-BFGS-B code to -c solve a sample problem (the extended Rosenbrock function -c subject to bounds on the variables). The dimension n of this -c problem is variable. - - integer nmax, mmax - parameter (nmax=1024, mmax=17) -c nmax is the dimension of the largest problem to be solved. -c mmax is the maximum number of limited memory corrections. - -c Declare the variables needed by the code. -c A description of all these variables is given at the end of -c the driver. - - character*60 task, csave - logical lsave(4) - integer n, m, iprint, - + nbd(nmax), iwa(3*nmax), isave(44) - double precision f, factr, pgtol, - + x(nmax), l(nmax), u(nmax), g(nmax), dsave(29), - + wa(2*mmax*nmax+4*nmax+12*mmax*mmax+12*mmax) - -c Declare a few additional variables for this sample problem. - - double precision t1, t2 - integer i - -c We wish to have output at every iteration. - - iprint = 1 - -c We specify the tolerances in the stopping criteria. - - factr=1.0d+7 - pgtol=1.0d-5 - -c We specify the dimension n of the sample problem and the number -c m of limited memory corrections stored. (n and m should not -c exceed the limits nmax and mmax respectively.) - - n=25 - m=5 - -c We now provide nbd which defines the bounds on the variables: -c l specifies the lower bounds, -c u specifies the upper bounds. - -c First set bounds on the odd-numbered variables. - - do 10 i=1,n,2 - nbd(i)=2 - l(i)=1.0d0 - u(i)=1.0d2 - 10 continue - -c Next set bounds on the even-numbered variables. - - do 12 i=2,n,2 - nbd(i)=2 - l(i)=-1.0d2 - u(i)=1.0d2 - 12 continue - -c We now define the starting point. - - do 14 i=1,n - x(i)=3.0d0 - 14 continue - - write (6,16) - 16 format(/,5x, 'Solving sample problem.', - + /,5x, ' (f = 0.0 at the optimal solution.)',/) - -c We start the iteration by initializing task. -c - task = 'START' - -c ------- the beginning of the loop ---------- - - 111 continue - -c This is the call to the L-BFGS-B code. - - call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint, - + csave,lsave,isave,dsave) - - if (task(1:2) .eq. 'FG') then -c the minimization routine has returned to request the -c function f and gradient g values at the current x. - -c Compute function value f for the sample problem. - - f=.25d0*(x(1)-1.d0)**2 - do 20 i=2,n - f=f+(x(i)-x(i-1)**2)**2 - 20 continue - f=4.d0*f - -c Compute gradient g for the sample problem. - - t1=x(2)-x(1)**2 - g(1)=2.d0*(x(1)-1.d0)-1.6d1*x(1)*t1 - do 22 i=2,n-1 - t2=t1 - t1=x(i+1)-x(i)**2 - g(i)=8.d0*t2-1.6d1*x(i)*t1 - 22 continue - g(n)=8.d0*t1 - -c go back to the minimization routine. - goto 111 - endif -c - if (task(1:5) .eq. 'NEW_X') goto 111 -c the minimization routine has returned with a new iterate, -c and we have opted to continue the iteration. - -c ---------- the end of the loop ------------- - -c If task is neither FG nor NEW_X we terminate execution. - - stop - - end - -c======================= The end of driver1 ============================ - -c -------------------------------------------------------------- -c DESCRIPTION OF THE VARIABLES IN L-BFGS-B -c -------------------------------------------------------------- -c -c n is an INTEGER variable that must be set by the user to the -c number of variables. It is not altered by the routine. -c -c m is an INTEGER variable that must be set by the user to the -c number of corrections used in the limited memory matrix. -c It is not altered by the routine. Values of m < 3 are -c not recommended, and large values of m can result in excessive -c computing time. The range 3 <= m <= 20 is recommended. -c -c x is a DOUBLE PRECISION array of length n. On initial entry -c it must be set by the user to the values of the initial -c estimate of the solution vector. Upon successful exit, it -c contains the values of the variables at the best point -c found (usually an approximate solution). -c -c l is a DOUBLE PRECISION array of length n that must be set by -c the user to the values of the lower bounds on the variables. If -c the i-th variable has no lower bound, l(i) need not be defined. -c -c u is a DOUBLE PRECISION array of length n that must be set by -c the user to the values of the upper bounds on the variables. If -c the i-th variable has no upper bound, u(i) need not be defined. -c -c nbd is an INTEGER array of dimension n that must be set by the -c user to the type of bounds imposed on the variables: -c nbd(i)=0 if x(i) is unbounded, -c 1 if x(i) has only a lower bound, -c 2 if x(i) has both lower and upper bounds, -c 3 if x(i) has only an upper bound. -c -c f is a DOUBLE PRECISION variable. If the routine setulb returns -c with task(1:2)= 'FG', then f must be set by the user to -c contain the value of the function at the point x. -c -c g is a DOUBLE PRECISION array of length n. If the routine setulb -c returns with taskb(1:2)= 'FG', then g must be set by the user to -c contain the components of the gradient at the point x. -c -c factr is a DOUBLE PRECISION variable that must be set by the user. -c It is a tolerance in the termination test for the algorithm. -c The iteration will stop when -c -c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch -c -c where epsmch is the machine precision which is automatically -c generated by the code. Typical values for factr on a computer -c with 15 digits of accuracy in double precision are: -c factr=1.d+12 for low accuracy; -c 1.d+7 for moderate accuracy; -c 1.d+1 for extremely high accuracy. -c The user can suppress this termination test by setting factr=0. -c -c pgtol is a double precision variable. -c On entry pgtol >= 0 is specified by the user. The iteration -c will stop when -c -c max{|proj g_i | i = 1, ..., n} <= pgtol -c -c where pg_i is the ith component of the projected gradient. -c The user can suppress this termination test by setting pgtol=0. -c -c wa is a DOUBLE PRECISION array of length -c (2mmax + 4)nmax + 12mmax^2 + 12mmax used as workspace. -c This array must not be altered by the user. -c -c iwa is an INTEGER array of length 3nmax used as -c workspace. This array must not be altered by the user. -c -c task is a CHARACTER string of length 60. -c On first entry, it must be set to 'START'. -c On a return with task(1:2)='FG', the user must evaluate the -c function f and gradient g at the returned value of x. -c On a return with task(1:5)='NEW_X', an iteration of the -c algorithm has concluded, and f and g contain f(x) and g(x) -c respectively. The user can decide whether to continue or stop -c the iteration. -c When -c task(1:4)='CONV', the termination test in L-BFGS-B has been -c satisfied; -c task(1:4)='ABNO', the routine has terminated abnormally -c without being able to satisfy the termination conditions, -c x contains the best approximation found, -c f and g contain f(x) and g(x) respectively; -c task(1:5)='ERROR', the routine has detected an error in the -c input parameters; -c On exit with task = 'CONV', 'ABNO' or 'ERROR', the variable task -c contains additional information that the user can print. -c This array should not be altered unless the user wants to -c stop the run for some reason. See driver2 or driver3 -c for a detailed explanation on how to stop the run -c by assigning task(1:4)='STOP' in the driver. -c -c iprint is an INTEGER variable that must be set by the user. -c It controls the frequency and type of output generated: -c iprint<0 no output is generated; -c iprint=0 print only one line at the last iteration; -c 0100 print details of every iteration including x and g; -c When iprint > 0, the file iterate.dat will be created to -c summarize the iteration. -c -c csave is a CHARACTER working array of length 60. -c -c lsave is a LOGICAL working array of dimension 4. -c On exit with task = 'NEW_X', the following information is -c available: -c lsave(1) = .true. the initial x did not satisfy the bounds; -c lsave(2) = .true. the problem contains bounds; -c lsave(3) = .true. each variable has upper and lower bounds. -c -c isave is an INTEGER working array of dimension 44. -c On exit with task = 'NEW_X', it contains information that -c the user may want to access: -c isave(30) = the current iteration number; -c isave(34) = the total number of function and gradient -c evaluations; -c isave(36) = the number of function value or gradient -c evaluations in the current iteration; -c isave(38) = the number of free variables in the current -c iteration; -c isave(39) = the number of active constraints at the current -c iteration; -c -c see the subroutine setulb.f for a description of other -c information contained in isave -c -c dsave is a DOUBLE PRECISION working array of dimension 29. -c On exit with task = 'NEW_X', it contains information that -c the user may want to access: -c dsave(2) = the value of f at the previous iteration; -c dsave(5) = the machine precision epsmch generated by the code; -c dsave(13) = the infinity norm of the projected gradient; -c -c see the subroutine setulb.f for a description of other -c information contained in dsave -c -c -------------------------------------------------------------- -c END OF THE DESCRIPTION OF THE VARIABLES IN L-BFGS-B -c -------------------------------------------------------------- -c -c << An example of subroutine 'timer' for AIX Version 3.2 >> -c -c subroutine timer(ttime) -c double precision ttime -c integer itemp, integer mclock -c itemp = mclock() -c ttime = dble(itemp)*1.0d-2 -c return -c end -c----------------------------------------------------------------------- EOF driver1.f echo driver2.f 1>&2 sed 's/.//' >driver2.f <<'EOF driver2.f' -c DRIVER 2 -c -------------------------------------------------------------- -c CUSTOMIZED DRIVER FOR L-BFGS-B -c -------------------------------------------------------------- -c -c L-BFGS-B is a code for solving large nonlinear optimization -c problems with simple bounds on the variables. -c -c The code can also be used for unconstrained problems and is -c as efficient for these problems as the earlier limited memory -c code L-BFGS. -c -c This driver illustrates how to control the termination of the -c run and how to design customized output. -c -c References: -c -c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited -c memory algorithm for bound constrained optimization'', -c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. -c -c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN -c Subroutines for Large Scale Bound Constrained Optimization'' -c Tech. Report, NAM-11, EECS Department, Northwestern University, -c 1994. -c -c -c (Postscript files of these papers are available via anonymous -c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************** - - program driver - -c This driver shows how to replace the default stopping test -c by other termination criteria. It also illustrates how to -c print the values of several parameters during the course of -c the iteration. The sample problem used here is the same as in -c DRIVER1 (the extended Rosenbrock function with bounds on the -c variables). - - integer nmax, mmax - parameter (nmax=1024, mmax=17) -c nmax is the dimension of the largest problem to be solved. -c mmax is the maximum number of limited memory corrections. - -c Declare the variables needed by the code. -c A description of all these variables is given at the end of -c driver1. - - character*60 task, csave - logical lsave(4) - integer n, m, iprint, - + nbd(nmax), iwa(3*nmax), isave(44) - double precision f, factr, pgtol, - + x(nmax), l(nmax), u(nmax), g(nmax), dsave(29), - + wa(2*mmax*nmax+4*nmax+12*mmax*mmax+12*mmax) - -c Declare a few additional variables for the sample problem. - - double precision t1, t2 - integer i - -c We suppress the default output. - - iprint = -1 - -c We suppress both code-supplied stopping tests because the -c user is providing his own stopping criteria. - - factr=0.0d0 - pgtol=0.0d0 - -c We specify the dimension n of the sample problem and the number -c m of limited memory corrections stored. (n and m should not -c exceed the limits nmax and mmax respectively.) - - n=25 - m=5 - -c We now specify nbd which defines the bounds on the variables: -c l specifies the lower bounds, -c u specifies the upper bounds. - -c First set bounds on the odd numbered variables. - - do 10 i=1,n,2 - nbd(i)=2 - l(i)=1.0d0 - u(i)=1.0d2 - 10 continue - -c Next set bounds on the even numbered variables. - - do 12 i=2,n,2 - nbd(i)=2 - l(i)=-1.0d2 - u(i)=1.0d2 - 12 continue - -c We now define the starting point. - - do 14 i=1,n - x(i)=3.0d0 - 14 continue - -c We now write the heading of the output. - - write (6,16) - 16 format(/,5x, 'Solving sample problem.', - + /,5x, ' (f = 0.0 at the optimal solution.)',/) - - -c We start the iteration by initializing task. -c - task = 'START' - -c ------- the beginning of the loop ---------- - - 111 continue - -c This is the call to the L-BFGS-B code. - - call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint, - + csave,lsave,isave,dsave) - - if (task(1:2) .eq. 'FG') then -c the minimization routine has returned to request the -c function f and gradient g values at the current x. - -c Compute function value f for the sample problem. - - f=.25d0*(x(1)-1.d0)**2 - do 20 i=2,n - f=f+(x(i)-x(i-1)**2)**2 - 20 continue - f=4.d0*f - -c Compute gradient g for the sample problem. - - t1=x(2)-x(1)**2 - g(1)=2.d0*(x(1)-1.d0)-1.6d1*x(1)*t1 - do 22 i=2,n-1 - t2=t1 - t1=x(i+1)-x(i)**2 - g(i)=8.d0*t2-1.6d1*x(i)*t1 - 22 continue - g(n)=8.d0*t1 - -c go back to the minimization routine. - goto 111 - endif -c - if (task(1:5) .eq. 'NEW_X') then -c -c the minimization routine has returned with a new iterate. -c At this point have the opportunity of stopping the iteration -c or observing the values of certain parameters -c -c First are two examples of stopping tests. - -c Note: task(1:4) must be assigned the value 'STOP' to terminate -c the iteration and ensure that the final results are -c printed in the default format. The rest of the character -c string TASK may be used to store other information. - -c 1) Terminate if the total number of f and g evaluations -c exceeds 99. - - if (isave(34) .ge. 99) - + task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT' - -c 2) Terminate if |proj g|/(1+|f|) < 1.0d-10, where -c "proj g" denoted the projected gradient - - if (dsave(13) .le. 1.d-10*(1.0d0 + abs(f))) - + task='STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL' - -c We now wish to print the following information at each -c iteration: -c -c 1) the current iteration number, isave(30), -c 2) the total number of f and g evaluations, isave(34), -c 3) the value of the objective function f, -c 4) the norm of the projected gradient, dsve(13) -c -c See the comments at the end of driver1 for a description -c of the variables isave and dsave. - - write (6,'(2(a,i5,4x),a,1p,d12.5,4x,a,1p,d12.5)') 'Iterate' - + ,isave(30),'nfg =',isave(34),'f =',f,'|proj g| =',dsave(13) - -c If the run is to be terminated, we print also the information -c contained in task as well as the final value of x. - - if (task(1:4) .eq. 'STOP') then - write (6,*) task - write (6,*) 'Final X=' - write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n) - endif - -c go back to the minimization routine. - goto 111 - - endif - -c ---------- the end of the loop ------------- - -c If task is neither FG nor NEW_X we terminate execution. - - stop - - end - -c======================= The end of driver2 ============================ - EOF driver2.f echo driver3.f 1>&2 sed 's/.//' >driver3.f <<'EOF driver3.f' -c DRIVER 3 -c -------------------------------------------------------------- -c TIME-CONTROLLED DRIVER FOR L-BFGS-B -c -------------------------------------------------------------- -c -c L-BFGS-B is a code for solving large nonlinear optimization -c problems with simple bounds on the variables. -c -c The code can also be used for unconstrained problems and is -c as efficient for these problems as the earlier limited memory -c code L-BFGS. -c -c This driver shows how to terminate a run after some prescribed -c CPU time has elapsed, and how to print the desired information -c before exiting. -c -c References: -c -c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited -c memory algorithm for bound constrained optimization'', -c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. -c -c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN -c Subroutines for Large Scale Bound Constrained Optimization'' -c Tech. Report, NAM-11, EECS Department, Northwestern University, -c 1994. -c -c -c (Postscript files of these papers are available via anonymous -c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************** - - program driver - -c This time-controlled driver shows that it is possible to terminate -c a run by elapsed CPU time, and yet be able to print all desired -c information. This driver also illustrates the use of two -c stopping criteria that may be used in conjunction with a limit -c on execution time. The sample problem used here is the same as in -c driver1 and driver2 (the extended Rosenbrock function with bounds -c on the variables). - - integer nmax, mmax - parameter (nmax=1024,mmax=17) -c nmax is the dimension of the largest problem to be solved. -c mmax is the maximum number of limited memory corrections. - -c Declare the variables needed by the code. -c A description of all these variables is given at the end of -c driver1. - - character*60 task, csave - logical lsave(4) - integer n, m, iprint, - + nbd(nmax), iwa(3*nmax), isave(44) - double precision f, factr, pgtol, - + x(nmax), l(nmax), u(nmax), g(nmax), dsave(29), - + wa(2*mmax*nmax+4*nmax+12*mmax*mmax+12*mmax) - -c Declare a few additional variables for the sample problem -c and for keeping track of time. - - double precision t1, t2, time1, time2, tlimit - integer i, j - -c We specify a limite on the CPU time (in seconds). - - tlimit = 0.2 - -c We suppress the default output. (The user could also elect to -c use the default output by choosing iprint >= 0.) - - iprint = -1 - -c We suppress the code-supplied stopping tests because we will -c provide our own termination conditions - - factr=0.0d0 - pgtol=0.0d0 - -c We specify the dimension n of the sample problem and the number -c m of limited memory corrections stored. (n and m should not -c exceed the limits nmax and mmax respectively.) - - n=1000 - m=10 - -c We now specify nbd which defines the bounds on the variables: -c l specifies the lower bounds, -c u specifies the upper bounds. - -c First set bounds on the odd-numbered variables. - - do 10 i=1,n,2 - nbd(i)=2 - l(i)=1.0d0 - u(i)=1.0d2 - 10 continue - -c Next set bounds on the even-numbered variables. - - do 12 i=2,n,2 - nbd(i)=2 - l(i)=-1.0d2 - u(i)=1.0d2 - 12 continue - -c We now define the starting point. - - do 14 i=1,n - x(i)=3.0d0 - 14 continue - -c We now write the heading of the output. - - write (6,16) - 16 format(/,5x, 'Solving sample problem.', - + /,5x, ' (f = 0.0 at the optimal solution.)',/) - -c We start the iteration by initializing task. -c - task = 'START' - -c ------- the beginning of the loop ---------- - -c We begin counting the CPU time. - - call timer(time1) - - 111 continue - -c This is the call to the L-BFGS-B code. - - call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint, - + csave,lsave,isave,dsave) - - if (task(1:2) .eq. 'FG') then -c the minimization routine has returned to request the -c function f and gradient g values at the current x. -c Before evaluating f and g we check the CPU time spent. - - call timer(time2) - if (time2-time1 .gt. tlimit) then - task='STOP: CPU EXCEEDING THE TIME LIMIT.' - -c Note: Assigning task(1:4)='STOP' will terminate the run; -c setting task(7:9)='CPU' will restore the information at -c the latest iterate generated by the code so that it can -c be correctly printed by the driver. - -c In this driver we have chosen to disable the -c printing options of the code (we set iprint=-1); -c instead we are using customized output: we print the -c latest value of x, the corresponding function value f and -c the norm of the projected gradient |proj g|. - -c We print out the information contained in task. - - write (6,*) task - -c We print the latest iterate contained in wa(j+1:j+n), where -c - j = 3*n+2*m*n+12*m**2 - write (6,*) 'Latest iterate X =' - write (6,'((1x,1p, 6(1x,d11.4)))') (wa(i),i = j+1,j+n) - -c We print the function value f and the norm of the projected -c gradient |proj g| at the last iterate; they are stored in -c dsave(2) and dsave(13) respectively. - - write (6,'(a,1p,d12.5,4x,a,1p,d12.5)') - + 'At latest iterate f =',dsave(2),'|proj g| =',dsave(13) - - else - -c The time limit has not been reached and we compute -c the function value f for the sample problem. - - f=.25d0*(x(1)-1.d0)**2 - do 20 i=2,n - f=f+(x(i)-x(i-1)**2)**2 - 20 continue - f=4.d0*f - -c Compute gradient g for the sample problem. - - t1=x(2)-x(1)**2 - g(1)=2.d0*(x(1)-1.d0)-1.6d1*x(1)*t1 - do 22 i=2,n-1 - t2=t1 - t1=x(i+1)-x(i)**2 - g(i)=8.d0*t2-1.6d1*x(i)*t1 - 22 continue - g(n)=8.d0*t1 - - endif - -c go back to the minimization routine. - goto 111 - endif -c - if (task(1:5) .eq. 'NEW_X') then -c the minimization routine has returned with a new iterate. -c The time limit has not been reached, and we test whether -c the following two stopping tests are satisfied: - -c 1) Terminate if the total number of f and g evaluations -c exceeds 900. - - if (isave(34) .ge. 900) - + task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT' - -c 2) Terminate if |proj g|/(1+|f|) < 1.0d-10. - - if (dsave(13) .le. 1.d-10*(1.0d0 + abs(f))) - + task='STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL' - -c We wish to print the following information at each iteration: -c 1) the current iteration number, isave(30), -c 2) the total number of f and g evaluations, isave(34), -c 3) the value of the objective function f, -c 4) the norm of the projected gradient, dsve(13) -c -c See the comments at the end of driver1 for a description -c of the variables isave and dsave. - - - write (6,'(2(a,i5,4x),a,1p,d12.5,4x,a,1p,d12.5)') 'Iterate' - + ,isave(30),'nfg =',isave(34),'f =',f,'|proj g| =',dsave(13) - -c If the run is to be terminated, we print also the information -c contained in task as well as the final value of x. - - - if (task(1:4) .eq. 'STOP') then - write (6,*) task - write (6,*) 'Final X=' - write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n) - endif - -c go back to the minimization routine. - goto 111 - - endif - -c ---------- the end of the loop ------------- - -c If task is neither FG nor NEW_X we terminate execution. - - stop - - end - -c======================= The end of driver3 ============================ - EOF driver3.f echo routines.f 1>&2 sed 's/.//' >routines.f <<'EOF routines.f' - - subroutine setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, - + task, iprint, csave, lsave, isave, dsave) - - character*60 task, csave - logical lsave(4) - integer n, m, iprint, - + nbd(n), iwa(3*n), isave(44) - double precision f, factr, pgtol, x(n), l(n), u(n), g(n), - + wa(2*m*n+4*n+12*m*m+12*m), dsave(29) - -c ************ -c -c Subroutine setulb -c -c This subroutine partitions the working arrays wa and iwa, and -c then uses the limited memory BFGS method to solve the bound -c constrained optimization problem by calling mainlb. -c (The direct method will be used in the subspace minimization.) -c -c n is an integer variable. -c On entry n is the dimension of the problem. -c On exit n is unchanged. -c -c m is an integer variable. -c On entry m is the maximum number of variable metric corrections -c used to define the limited memory matrix. -c On exit m is unchanged. -c -c x is a double precision array of dimension n. -c On entry x is an approximation to the solution. -c On exit x is the current approximation. -c -c l is a double precision array of dimension n. -c On entry l is the lower bound on x. -c On exit l is unchanged. -c -c u is a double precision array of dimension n. -c On entry u is the upper bound on x. -c On exit u is unchanged. -c -c nbd is an integer array of dimension n. -c On entry nbd represents the type of bounds imposed on the -c variables, and must be specified as follows: -c nbd(i)=0 if x(i) is unbounded, -c 1 if x(i) has only a lower bound, -c 2 if x(i) has both lower and upper bounds, and -c 3 if x(i) has only an upper bound. -c On exit nbd is unchanged. -c -c f is a double precision variable. -c On first entry f is unspecified. -c On final exit f is the value of the function at x. -c -c g is a double precision array of dimension n. -c On first entry g is unspecified. -c On final exit g is the value of the gradient at x. -c -c factr is a double precision variable. -c On entry factr >= 0 is specified by the user. The iteration -c will stop when -c -c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch -c -c where epsmch is the machine precision, which is automatically -c generated by the code. Typical values for factr: 1.d+12 for -c low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely -c high accuracy. -c On exit factr is unchanged. -c -c pgtol is a double precision variable. -c On entry pgtol >= 0 is specified by the user. The iteration -c will stop when -c -c max{|proj g_i | i = 1, ..., n} <= pgtol -c -c where pg_i is the ith component of the projected gradient. -c On exit pgtol is unchanged. -c -c wa is a double precision working array of length -c (2mmax + 4)nmax + 12mmax^2 + 12mmax. -c -c iwa is an integer working array of length 3nmax. -c -c task is a working string of characters of length 60 indicating -c the current job when entering and quitting this subroutine. -c -c iprint is an integer variable that must be set by the user. -c It controls the frequency and type of output generated: -c iprint<0 no output is generated; -c iprint=0 print only one line at the last iteration; -c 0100 print details of every iteration including x and g; -c When iprint > 0, the file iterate.dat will be created to -c summarize the iteration. -c -c csave is a working string of characters of length 60. -c -c lsave is a logical working array of dimension 4. -c On exit with 'task' = NEW_X, the following information is -c available: -c If lsave(1) = .true. then the initial X has been replaced by -c its projection in the feasible set; -c If lsave(2) = .true. then the problem is constrained; -c If lsave(3) = .true. then each variable has upper and lower -c bounds; -c -c isave is an integer working array of dimension 44. -c On exit with 'task' = NEW_X, the following information is -c available: -c isave(22) = the total number of intervals explored in the -c search of Cauchy points; -c isave(26) = the total number of skipped BFGS updates before -c the current iteration; -c isave(30) = the number of current iteration; -c isave(31) = the total number of BFGS updates prior the current -c iteration; -c isave(33) = the number of intervals explored in the search of -c Cauchy point in the current iteration; -c isave(34) = the total number of function and gradient -c evaluations; -c isave(36) = the number of function value or gradient -c evaluations in the current iteration; -c if isave(37) = 0 then the subspace argmin is within the box; -c if isave(37) = 1 then the subspace argmin is beyond the box; -c isave(38) = the number of free variables in the current -c iteration; -c isave(39) = the number of active constraints in the current -c iteration; -c n + 1 - isave(40) = the number of variables leaving the set of -c active constraints in the current iteration; -c isave(41) = the number of variables entering the set of active -c constraints in the current iteration. -c -c dsave is a double precision working array of dimension 29. -c On exit with 'task' = NEW_X, the following information is -c available: -c dsave(1) = current 'theta' in the BFGS matrix; -c dsave(2) = f(x) in the previous iteration; -c dsave(3) = factr*epsmch; -c dsave(4) = 2-norm of the line search direction vector; -c dsave(5) = the machine precision epsmch generated by the code; -c dsave(7) = the accumulated time spent on searching for -c Cauchy points; -c dsave(8) = the accumulated time spent on -c subspace minimization; -c dsave(9) = the accumulated time spent on line search; -c dsave(11) = the slope of the line search function at -c the current point of line search; -c dsave(12) = the maximum relative step length imposed in -c line search; -c dsave(13) = the infinity norm of the projected gradient; -c dsave(14) = the relative step length in the line search; -c dsave(15) = the slope of the line search function at -c the starting point of the line search; -c dsave(16) = the square of the 2-norm of the line search -c direction vector. -c -c Subprograms called: -c -c L-BFGS-B Library ... mainlb. -c -c -c References: -c -c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited -c memory algorithm for bound constrained optimization'', -c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. -c -c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a -c limited memory FORTRAN code for solving bound constrained -c optimization problems'', Tech. Report, NAM-11, EECS Department, -c Northwestern University, 1994. -c -c (Postscript files of these papers are available via anonymous -c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer l1,l2,l3,lws,lr,lz,lt,ld,lsg,lwa,lyg, - + lsgo,lwy,lsy,lss,lyy,lwt,lwn,lsnd,lygo - - if (task .eq. 'START') then - isave(1) = m*n - isave(2) = m**2 - isave(3) = 4*m**2 - isave(4) = 1 - isave(5) = isave(4) + isave(1) - isave(6) = isave(5) + isave(1) - isave(7) = isave(6) + isave(2) - isave(8) = isave(7) + isave(2) - isave(9) = isave(8) + isave(2) - isave(10) = isave(9) + isave(2) - isave(11) = isave(10) + isave(3) - isave(12) = isave(11) + isave(3) - isave(13) = isave(12) + n - isave(14) = isave(13) + n - isave(15) = isave(14) + n - isave(16) = isave(15) + n - isave(17) = isave(16) + 8*m - isave(18) = isave(17) + m - isave(19) = isave(18) + m - isave(20) = isave(19) + m - endif - l1 = isave(1) - l2 = isave(2) - l3 = isave(3) - lws = isave(4) - lwy = isave(5) - lsy = isave(6) - lss = isave(7) - lyy = isave(8) - lwt = isave(9) - lwn = isave(10) - lsnd = isave(11) - lz = isave(12) - lr = isave(13) - ld = isave(14) - lt = isave(15) - lwa = isave(16) - lsg = isave(17) - lsgo = isave(18) - lyg = isave(19) - lygo = isave(20) - - call mainlb(n,m,x,l,u,nbd,f,g,factr,pgtol, - + wa(lws),wa(lwy),wa(lsy),wa(lss),wa(lyy),wa(lwt), - + wa(lwn),wa(lsnd),wa(lz),wa(lr),wa(ld),wa(lt), - + wa(lwa),wa(lsg),wa(lsgo),wa(lyg),wa(lygo), - + iwa(1),iwa(n+1),iwa(2*n+1),task,iprint, - + csave,lsave,isave(22),dsave) - - return - - end - -c======================= The end of setulb ============================= - - subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy, - + sy, ss, yy, wt, wn, snd, z, r, d, t, wa, sg, - + sgo, yg, ygo, index, iwhere, indx2, task, - + iprint, csave, lsave, isave, dsave) - - character*60 task, csave - logical lsave(4) - integer n, m, iprint, nbd(n), index(n), - + iwhere(n), indx2(n), isave(23) - double precision f, factr, pgtol, - + x(n), l(n), u(n), g(n), z(n), r(n), d(n), t(n), - + wa(8*m), sg(m), sgo(m), yg(m), ygo(m), - + ws(n, m), wy(n, m), sy(m, m), ss(m, m), yy(m, m), - + wt(m, m), wn(2*m, 2*m), snd(2*m, 2*m), dsave(29) - -c ************ -c -c Subroutine mainlb -c -c This subroutine solves bound constrained optimization problems by -c using the compact formula of the limited memory BFGS updates. -c -c n is an integer variable. -c On entry n is the number of variables. -c On exit n is unchanged. -c -c m is an integer variable. -c On entry m is the maximum number of variable metric -c corrections allowed in the limited memory matrix. -c On exit m is unchanged. -c -c x is a double precision array of dimension n. -c On entry x is an approximation to the solution. -c On exit x is the current approximation. -c -c l is a double precision array of dimension n. -c On entry l is the lower bound of x. -c On exit l is unchanged. -c -c u is a double precision array of dimension n. -c On entry u is the upper bound of x. -c On exit u is unchanged. -c -c nbd is an integer array of dimension n. -c On entry nbd represents the type of bounds imposed on the -c variables, and must be specified as follows: -c nbd(i)=0 if x(i) is unbounded, -c 1 if x(i) has only a lower bound, -c 2 if x(i) has both lower and upper bounds, -c 3 if x(i) has only an upper bound. -c On exit nbd is unchanged. -c -c f is a double precision variable. -c On first entry f is unspecified. -c On final exit f is the value of the function at x. -c -c g is a double precision array of dimension n. -c On first entry g is unspecified. -c On final exit g is the value of the gradient at x. -c -c factr is a double precision variable. -c On entry factr >= 0 is specified by the user. The iteration -c will stop when -c -c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch -c -c where epsmch is the machine precision, which is automatically -c generated by the code. -c On exit factr is unchanged. -c -c pgtol is a double precision variable. -c On entry pgtol >= 0 is specified by the user. The iteration -c will stop when -c -c max{|proj g_i | i = 1, ..., n} <= pgtol -c -c where pg_i is the ith component of the projected gradient. -c On exit pgtol is unchanged. -c -c ws, wy, sy, and wt are double precision working arrays used to -c store the following information defining the limited memory -c BFGS matrix: -c ws, of dimension n x m, stores S, the matrix of s-vectors; -c wy, of dimension n x m, stores Y, the matrix of y-vectors; -c sy, of dimension m x m, stores S'Y; -c ss, of dimension m x m, stores S'S; -c yy, of dimension m x m, stores Y'Y; -c wt, of dimension m x m, stores the Cholesky factorization -c of (theta*S'S+LD^(-1)L'); see eq. -c (2.26) in [3]. -c -c wn is a double precision working array of dimension 2m x 2m -c used to store the LEL^T factorization of the indefinite matrix -c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] -c [L_a -R_z theta*S'AA'S ] -c -c where E = [-I 0] -c [ 0 I] -c -c snd is a double precision working array of dimension 2m x 2m -c used to store the lower triangular part of -c N = [Y' ZZ'Y L_a'+R_z'] -c [L_a +R_z S'AA'S ] -c -c z(n),r(n),d(n),t(n),wa(8*m) are double precision working arrays. -c z is used at different times to store the Cauchy point and -c the Newton point. -c -c sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays. -c -c index is an integer working array of dimension n. -c In subroutine freev, index is used to store the free and fixed -c variables at the Generalized Cauchy Point (GCP). -c -c iwhere is an integer working array of dimension n used to record -c the status of the vector x for GCP computation. -c iwhere(i)=0 or -3 if x(i) is free and has bounds, -c 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) -c 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) -c 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) -c -1 if x(i) is always free, i.e., no bounds on it. -c -c indx2 is an integer working array of dimension n. -c Within subroutine cauchy, indx2 corresponds to the array iorder. -c In subroutine freev, a list of variables entering and leaving -c the free set is stored in indx2, and it is passed on to -c subroutine formk with this information. -c -c task is a working string of characters of length 60 indicating -c the current job when entering and leaving this subroutine. -c -c iprint is an INTEGER variable that must be set by the user. -c It controls the frequency and type of output generated: -c iprint<0 no output is generated; -c iprint=0 print only one line at the last iteration; -c 0100 print details of every iteration including x and g; -c When iprint > 0, the file iterate.dat will be created to -c summarize the iteration. -c -c csave is a working string of characters of length 60. -c -c lsave is a logical working array of dimension 4. -c -c isave is an integer working array of dimension 23. -c -c dsave is a double precision working array of dimension 29. -c -c -c Subprograms called -c -c L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk, -c -c errclb, prn1lb, prn2lb, prn3lb, active, projgr, -c -c freev, cmprlb, matupd, formt. -c -c Minpack2 Library ... timer, dpmeps. -c -c Linpack Library ... dcopy, ddot. -c -c -c References: -c -c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited -c memory algorithm for bound constrained optimization'', -c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. -c -c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN -c Subroutines for Large Scale Bound Constrained Optimization'' -c Tech. Report, NAM-11, EECS Department, Northwestern University, -c 1994. -c -c [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of -c Quasi-Newton Matrices and their use in Limited Memory Methods'', -c Mathematical Programming 63 (1994), no. 4, pp. 129-156. -c -c (Postscript files of these papers are available via anonymous -c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - logical prjctd,cnstnd,boxed,updatd,wrk - character*3 word - integer i,k,nintol,itfile,iback,nskip, - + head,col,iter,itail,iupdat, - + nint,nfgv,info,ifun, - + iword,nfree,nact,ileave,nenter - double precision theta,fold,ddot,dr,rr,tol,dpmeps, - + xstep,sbgnrm,ddum,dnorm,dtd,epsmch, - + cpu1,cpu2,cachyt,sbtime,lnscht,time1,time2, - + gd,gdold,stp,stpmx,time - double precision one,zero - parameter (one=1.0d0,zero=0.0d0) - - if (task .eq. 'START') then - - call timer(time1) - -c Generate the current machine precision. - - epsmch = dpmeps() - -c Initialize counters and scalars when task='START'. - -c for the limited memory BFGS matrices: - col = 0 - head = 1 - theta = one - iupdat = 0 - updatd = .false. - -c for operation counts: - iter = 0 - nfgv = 0 - nint = 0 - nintol = 0 - nskip = 0 - nfree = n - -c for stopping tolerance: - tol = factr*epsmch - -c for measuring running time: - cachyt = 0 - sbtime = 0 - lnscht = 0 - -c 'word' records the status of subspace solutions. - word = '---' - -c 'info' records the termination information. - info = 0 - - if (iprint .ge. 1) then -c open a summary file 'iterate.dat' - open (8, file = 'iterate.dat', status = 'unknown') - itfile = 8 - endif - -c Check the input arguments for errors. - - call errclb(n,m,factr,l,u,nbd,task,info,k) - if (task(1:5) .eq. 'ERROR') then - call prn3lb(n,x,f,task,iprint,info,itfile, - + iter,nfgv,nintol,nskip,nact,sbgnrm, - + zero,nint,word,iback,stp,xstep,k, - + cachyt,sbtime,lnscht) - return - endif - - call prn1lb(n,m,l,u,x,iprint,itfile,epsmch) - -c Initialize iwhere & project x onto the feasible set. - - call active(n,l,u,nbd,x,iwhere,iprint,prjctd,cnstnd,boxed) - -c The end of the initialization. - - else -c 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) - nint = 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) - -c After returning from the driver go to the point where execution -c is to resume. - - if (task(1:5) .eq. 'FG_LN') goto 666 - if (task(1:5) .eq. 'NEW_X') goto 777 - if (task(1:5) .eq. 'FG_ST') goto 111 - if (task(1:4) .eq. 'STOP') then - if (task(7:9) .eq. 'CPU') then -c restore the previous iterate. - call dcopy(n,t,1,x,1) - call dcopy(n,r,1,g,1) - f = fold - endif - goto 999 - endif - endif - -c Compute f0 and g0. - - task = 'FG_START' -c return to the driver to calculate f and g; reenter at 111. - goto 1000 - 111 continue - nfgv = 1 - -c Compute the infinity norm of the (-) projected gradient. - - call projgr(n,l,u,nbd,x,g,sbgnrm) - - if (iprint .ge. 1) then - write (6,1002) iter,f,sbgnrm - write (itfile,1003) iter,nfgv,sbgnrm,f - endif - if (sbgnrm .le. pgtol) then -c terminate the algorithm. - task = 'CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL' - goto 999 - endif - -c ----------------- the beginning of the loop -------------------------- - - 222 continue - if (iprint .ge. 99) write (6,1001) iter + 1 - iword = -1 -c - if (.not. cnstnd .and. col .gt. 0) then -c skip the search for GCP. - call dcopy(n,x,1,z,1) - wrk = updatd - nint = 0 - goto 333 - endif - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c -c Compute the Generalized Cauchy Point (GCP). -c -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - call timer(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),nint, - + sg,yg,iprint,sbgnrm,info) - if (info .ne. 0) then -c singular triangular system detected; refresh the lbfgs memory. - if(iprint .ge. 1) write (6, 1005) - info = 0 - col = 0 - head = 1 - theta = one - iupdat = 0 - updatd = .false. - call timer(cpu2) - cachyt = cachyt + cpu2 - cpu1 - goto 222 - endif - call timer(cpu2) - cachyt = cachyt + cpu2 - cpu1 - nintol = nintol + nint - -c Count the entering and leaving variables for iter > 0; -c 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 - - 333 continue - -c If there are no free variables or B=I, then -c skip the subspace minimization. - - if (nfree .eq. 0 .or. col .eq. 0) goto 555 - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c -c Subspace minimization. -c -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - call timer(cpu1) - -c Form the LEL^T factorization of the indefinite -c matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] -c [L_a -R_z theta*S'AA'S ] -c where E = [-I 0] -c [ 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 .ne. 0) then -c nonpositive definiteness in Cholesky factorization; -c refresh the lbfgs memory and restart the iteration. - if(iprint .ge. 1) write (6, 1006) - info = 0 - col = 0 - head = 1 - theta = one - iupdat = 0 - updatd = .false. - call timer(cpu2) - sbtime = sbtime + cpu2 - cpu1 - goto 222 - endif - -c compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) -c from 'cauchy'). - call cmprlb(n,m,x,g,ws,wy,sy,wt,z,r,wa,index, - + theta,col,head,nfree,cnstnd,info) - if (info .ne. 0) goto 444 -c call the direct method. - call subsm(n,m,nfree,index,l,u,nbd,z,r,ws,wy,theta, - + col,head,iword,wa,wn,iprint,info) - 444 continue - if (info .ne. 0) then -c singular triangular system detected; -c refresh the lbfgs memory and restart the iteration. - if(iprint .ge. 1) write (6, 1005) - info = 0 - col = 0 - head = 1 - theta = one - iupdat = 0 - updatd = .false. - call timer(cpu2) - sbtime = sbtime + cpu2 - cpu1 - goto 222 - endif - - call timer(cpu2) - sbtime = sbtime + cpu2 - cpu1 - 555 continue - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c -c Line search and optimality tests. -c -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c Generate the search direction d:=z-x. - - do 40 i = 1, n - d(i) = z(i) - x(i) - 40 continue - call timer(cpu1) - 666 continue - 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 .ne. 0 .or. iback .ge. 20) then -c restore the previous iterate. - call dcopy(n,t,1,x,1) - call dcopy(n,r,1,g,1) - f = fold - if (col .eq. 0) then -c abnormal termination. - if (info .eq. 0) then - info = -9 -c 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 - goto 999 - else -c refresh the lbfgs memory and restart the iteration. - if(iprint .ge. 1) write (6, 1008) - if (info .eq. 0) nfgv = nfgv - 1 - info = 0 - col = 0 - head = 1 - theta = one - iupdat = 0 - updatd = .false. - task = 'RESTART_FROM_LNSRCH' - call timer(cpu2) - lnscht = lnscht + cpu2 - cpu1 - goto 222 - endif - else if (task(1:5) .eq. 'FG_LN') then -c return to the driver for calculating f and g; reenter at 666. - goto 1000 - else -c calculate and print out the quantities related to the new X. - call timer(cpu2) - lnscht = lnscht + cpu2 - cpu1 - iter = iter + 1 - -c Compute the infinity norm of the projected (-)gradient. - - call projgr(n,l,u,nbd,x,g,sbgnrm) - -c Print iteration information. - - call prn2lb(n,x,f,g,iprint,itfile,iter,nfgv,nact, - + sbgnrm,nint,word,iword,iback,stp,xstep) - goto 1000 - endif - 777 continue - -c Test for termination. - - if (sbgnrm .le. pgtol) then -c terminate the algorithm. - task = 'CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL' - goto 999 - endif - - ddum = max(abs(fold), abs(f), one) - if ((fold - f) .le. tol*ddum) then -c terminate the algorithm. - task = 'CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH' - if (iback .ge. 10) info = -5 -c i.e., to issue a warning if iback>10 in the line search. - goto 999 - endif - -c Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. - - do 42 i = 1, n - r(i) = g(i) - r(i) - 42 continue - rr = ddot(n,r,1,r,1) - if (stp .eq. one) then - dr = gd - gdold - ddum = -gdold - else - dr = (gd - gdold)*stp - call dscal(n,stp,d,1) - ddum = -gdold*stp - endif - - if (dr .le. epsmch*ddum) then -c skip the L-BFGS update. - nskip = nskip + 1 - updatd = .false. - if (iprint .ge. 1) write (6,1004) dr, ddum - goto 888 - endif - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c -c Update the L-BFGS matrix. -c -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - updatd = .true. - iupdat = iupdat + 1 - -c 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) - -c Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; -c Store T in the upper triangular of the array wt; -c Cholesky factorize T to J*J' with -c J' stored in the upper triangular of wt. - - call formt(m,wt,sy,ss,col,theta,info) - - if (info .ne. 0) then -c nonpositive definiteness in Cholesky factorization; -c refresh the lbfgs memory and restart the iteration. - if(iprint .ge. 1) write (6, 1007) - info = 0 - col = 0 - head = 1 - theta = one - iupdat = 0 - updatd = .false. - goto 222 - endif - -c Now the inverse of the middle matrix in B is - -c [ D^(1/2) O ] [ -D^(1/2) D^(-1/2)*L' ] -c [ -L*D^(-1/2) J ] [ 0 J' ] - - 888 continue - -c -------------------- the end of the loop ----------------------------- - - goto 222 - 999 continue - call timer(time2) - time = time2 - time1 - call prn3lb(n,x,f,task,iprint,info,itfile, - + iter,nfgv,nintol,nskip,nact,sbgnrm, - + time,nint,word,iback,stp,xstep,k, - + cachyt,sbtime,lnscht) - 1000 continue - -c 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) = nint - 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 - - 1001 format (//,'ITERATION ',i5) - 1002 format - + (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5) - 1003 format (2(1x,i4),5x,'-',5x,'-',3x,'-',5x,'-',5x,'-',8x,'-',3x, - + 1p,2(1x,d10.3)) - 1004 format (' ys=',1p,e10.3,' -gs=',1p,e10.3,' BFGS update SKIPPED') - 1005 format (/, - +' Singular triangular system detected;',/, - +' refresh the lbfgs memory and restart the iteration.') - 1006 format (/, - +' Nonpositive definiteness in Cholesky factorization in formk;',/, - +' refresh the lbfgs memory and restart the iteration.') - 1007 format (/, - +' Nonpositive definiteness in Cholesky factorization in formt;',/, - +' refresh the lbfgs memory and restart the iteration.') - 1008 format (/, - +' Bad direction in the line search;',/, - +' refresh the lbfgs memory and restart the iteration.') - - return - - end - -c======================= The end of mainlb ============================= - - subroutine active(n, l, u, nbd, x, iwhere, iprint, - + prjctd, cnstnd, boxed) - - logical prjctd, cnstnd, boxed - integer n, iprint, nbd(n), iwhere(n) - double precision x(n), l(n), u(n) - -c ************ -c -c Subroutine active -c -c This subroutine initializes iwhere and projects the initial x to -c the feasible set if necessary. -c -c iwhere is an integer array of dimension n. -c On entry iwhere is unspecified. -c On exit iwhere(i)=-1 if x(i) has no bounds -c 3 if l(i)=u(i) -c 0 otherwise. -c In cauchy, iwhere is given finer gradations. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer nbdd,i - double precision zero - parameter (zero=0.0d0) - -c Initialize nbdd, prjctd, cnstnd and boxed. - - nbdd = 0 - prjctd = .false. - cnstnd = .false. - boxed = .true. - -c Project the initial x to the easible set if necessary. - - do 10 i = 1, n - if (nbd(i) .gt. 0) then - if (nbd(i) .le. 2 .and. x(i) .le. l(i)) then - if (x(i) .lt. l(i)) then - prjctd = .true. - x(i) = l(i) - endif - nbdd = nbdd + 1 - else if (nbd(i) .ge. 2 .and. x(i) .ge. u(i)) then - if (x(i) .gt. u(i)) then - prjctd = .true. - x(i) = u(i) - endif - nbdd = nbdd + 1 - endif - endif - 10 continue - -c Initialize iwhere and assign values to cnstnd and boxed. - - do 20 i = 1, n - if (nbd(i) .ne. 2) boxed = .false. - if (nbd(i) .eq. 0) then -c this variable is always free - iwhere(i) = -1 - -c otherwise set x(i)=mid(x(i), u(i), l(i)). - else - cnstnd = .true. - if (nbd(i) .eq. 2 .and. u(i) - l(i) .le. zero) then -c this variable is always fixed - iwhere(i) = 3 - else - iwhere(i) = 0 - endif - endif - 20 continue - - if (iprint .ge. 0) then - if (prjctd) write (6,*) - + 'The initial X is infeasible. Restart with its projection.' - if (.not. cnstnd) - + write (6,*) 'This problem is unconstrained.' - endif - - if (iprint .gt. 0) write (6,1001) nbdd - - 1001 format (/,'At X0 ',i9,' variables are exactly at the bounds') - - return - - end - -c======================= The end of active ============================= - - subroutine bmv(m, sy, wt, col, v, p, info) - - integer m, col, info - double precision sy(m, m), wt(m, m), v(2*col), p(2*col) - -c ************ -c -c Subroutine bmv -c -c This subroutine computes the product of the 2m x 2m middle matrix -c in the compact L-BFGS formula of B and a 2m vector v; -c it returns the product in p. -c -c m is an integer variable. -c On entry m is the maximum number of variable metric corrections -c used to define the limited memory matrix. -c On exit m is unchanged. -c -c sy is a double precision array of dimension m x m. -c On entry sy specifies the matrix S'Y. -c On exit sy is unchanged. -c -c wt is a double precision array of dimension m x m. -c On entry wt specifies the upper triangular matrix J' which is -c the Cholesky factor of (thetaS'S+LD^(-1)L'). -c On exit wt is unchanged. -c -c col is an integer variable. -c On entry col specifies the number of s-vectors (or y-vectors) -c stored in the compact L-BFGS formula. -c On exit col is unchanged. -c -c v is a double precision array of dimension 2col. -c On entry v specifies vector v. -c On exit v is unchanged. -c -c p is a double precision array of dimension 2col. -c On entry p is unspecified. -c On exit p is the product Mv. -c -c info is an integer variable. -c On entry info is unspecified. -c On exit info = 0 for normal return, -c = nonzero for abnormal return when the system -c to be solved by dtrsl is singular. -c -c Subprograms called: -c -c Linpack ... dtrsl. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer i,k,i2 - double precision sum - - if (col .eq. 0) return - -c PART I: solve [ D^(1/2) O ] [ p1 ] = [ v1 ] -c [ -L*D^(-1/2) J ] [ p2 ] [ v2 ]. - -c solve Jp2=v2+LD^(-1)v1. - p(col + 1) = v(col + 1) - do 20 i = 2, col - i2 = col + i - sum = 0.0d0 - do 10 k = 1, i - 1 - sum = sum + sy(i,k)*v(k)/sy(k,k) - 10 continue - p(i2) = v(i2) + sum - 20 continue -c Solve the triangular system - call dtrsl(wt,m,col,p(col+1),11,info) - if (info .ne. 0) return - -c solve D^(1/2)p1=v1. - do 30 i = 1, col - p(i) = v(i)/sqrt(sy(i,i)) - 30 continue - -c PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ] -c [ 0 J' ] [ p2 ] [ p2 ]. - -c solve J^Tp2=p2. - call dtrsl(wt,m,col,p(col+1),01,info) - if (info .ne. 0) return - -c compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2) -c =-D^(-1/2)p1+D^(-1)L'p2. - do 40 i = 1, col - p(i) = -p(i)/sqrt(sy(i,i)) - 40 continue - do 60 i = 1, col - sum = 0.d0 - do 50 k = i + 1, col - sum = sum + sy(k,i)*p(col+k)/sy(i,i) - 50 continue - p(i) = p(i) + sum - 60 continue - - return - - end - -c======================== The end of bmv =============================== - - subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp, - + m, wy, ws, sy, wt, theta, col, head, p, c, wbp, - + v, nint, sg, yg, iprint, sbgnrm, info) - - integer n, m, head, col, nint, iprint, info, - + nbd(n), iorder(n), iwhere(n) - double precision theta, - + x(n), l(n), u(n), g(n), t(n), d(n), xcp(n), - + sg(m), yg(m), wy(n, col), ws(n, col), sy(m, m), - + wt(m, m), p(2*m), c(2*m), wbp(2*m), v(2*m) - -c ************ -c -c Subroutine cauchy -c -c For given x, l, u, g (with sbgnrm > 0), and a limited memory -c BFGS matrix B defined in terms of matrices WY, WS, WT, and -c scalars head, col, and theta, this subroutine computes the -c generalized Cauchy point (GCP), defined as the first local -c minimizer of the quadratic -c -c Q(x + s) = g's + 1/2 s'Bs -c -c along the projected gradient direction P(x-tg,l,u). -c The routine returns the GCP in xcp. -c -c n is an integer variable. -c On entry n is the dimension of the problem. -c On exit n is unchanged. -c -c x is a double precision array of dimension n. -c On entry x is the starting point for the GCP computation. -c On exit x is unchanged. -c -c l is a double precision array of dimension n. -c On entry l is the lower bound of x. -c On exit l is unchanged. -c -c u is a double precision array of dimension n. -c On entry u is the upper bound of x. -c On exit u is unchanged. -c -c nbd is an integer array of dimension n. -c On entry nbd represents the type of bounds imposed on the -c variables, and must be specified as follows: -c nbd(i)=0 if x(i) is unbounded, -c 1 if x(i) has only a lower bound, -c 2 if x(i) has both lower and upper bounds, and -c 3 if x(i) has only an upper bound. -c On exit nbd is unchanged. -c -c g is a double precision array of dimension n. -c On entry g is the gradient of f(x). g must be a nonzero vector. -c On exit g is unchanged. -c -c iorder is an integer working array of dimension n. -c iorder will be used to store the breakpoints in the piecewise -c linear path and free variables encountered. On exit, -c iorder(1),...,iorder(nleft) are indices of breakpoints -c which have not been encountered; -c iorder(nleft+1),...,iorder(nbreak) are indices of -c encountered breakpoints; and -c iorder(nfree),...,iorder(n) are indices of variables which -c have no bound constraits along the search direction. -c -c iwhere is an integer array of dimension n. -c On entry iwhere indicates only the permanently fixed (iwhere=3) -c or free (iwhere= -1) components of x. -c On exit iwhere records the status of the current x variables. -c iwhere(i)=-3 if x(i) is free and has bounds, but is not moved -c 0 if x(i) is free and has bounds, and is moved -c 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) -c 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) -c 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) -c -1 if x(i) is always free, i.e., it has no bounds. -c -c t is a double precision working array of dimension n. -c t will be used to store the break points. -c -c d is a double precision array of dimension n used to store -c the Cauchy direction P(x-tg)-x. -c -c xcp is a double precision array of dimension n used to return the -c GCP on exit. -c -c m is an integer variable. -c On entry m is the maximum number of variable metric corrections -c used to define the limited memory matrix. -c On exit m is unchanged. -c -c ws, wy, sy, and wt are double precision arrays. -c On entry they store information that defines the -c limited memory BFGS matrix: -c ws(n,m) stores S, a set of s-vectors; -c wy(n,m) stores Y, a set of y-vectors; -c sy(m,m) stores S'Y; -c wt(m,m) stores the -c Cholesky factorization of (theta*S'S+LD^(-1)L'). -c On exit these arrays are unchanged. -c -c theta is a double precision variable. -c On entry theta is the scaling factor specifying B_0 = theta I. -c On exit theta is unchanged. -c -c col is an integer variable. -c On entry col is the actual number of variable metric -c corrections stored so far. -c On exit col is unchanged. -c -c head is an integer variable. -c On entry head is the location of the first s-vector (or y-vector) -c in S (or Y). -c On exit col is unchanged. -c -c p is a double precision working array of dimension 2m. -c p will be used to store the vector p = W^(T)d. -c -c c is a double precision working array of dimension 2m. -c c will be used to store the vector c = W^(T)(xcp-x). -c -c wbp is a double precision working array of dimension 2m. -c wbp will be used to store the row of W corresponding -c to a breakpoint. -c -c v is a double precision working array of dimension 2m. -c -c nint is an integer variable. -c On exit nint records the number of quadratic segments explored -c in searching for the GCP. -c -c sg and yg are double precision arrays of dimension m. -c On entry sg and yg store S'g and Y'g correspondingly. -c On exit they are unchanged. -c -c iprint is an INTEGER variable that must be set by the user. -c It controls the frequency and type of output generated: -c iprint<0 no output is generated; -c iprint=0 print only one line at the last iteration; -c 0100 print details of every iteration including x and g; -c When iprint > 0, the file iterate.dat will be created to -c summarize the iteration. -c -c sbgnrm is a double precision variable. -c On entry sbgnrm is the norm of the projected gradient at x. -c On exit sbgnrm is unchanged. -c -c info is an integer variable. -c On entry info is 0. -c On exit info = 0 for normal return, -c = nonzero for abnormal return when the the system -c used in routine bmv is singular. -c -c Subprograms called: -c -c L-BFGS-B Library ... hpsolb, bmv. -c -c Linpack ... dscal dcopy, daxpy. -c -c -c References: -c -c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited -c memory algorithm for bound constrained optimization'', -c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. -c -c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN -c Subroutines for Large Scale Bound Constrained Optimization'' -c Tech. Report, NAM-11, EECS Department, Northwestern University, -c 1994. -c -c (Postscript files of these papers are available via anonymous -c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - logical xlower,xupper,bnded - integer i,j,col2,nfree,nbreak,pointr, - + ibp,nleft,ibkmin,iter - double precision f1,f2,dt,dtm,tsum,dibp,zibp,dibp2,bkmin, - + tu,tl,wmc,wmp,wmw,ddot,tj,tj0,neggi,sbgnrm - double precision one,zero - parameter (one=1.0d0,zero=0.0d0) - -c Check the status of the variables, reset iwhere(i) if necessary; -c compute the Cauchy direction d and the breakpoints t; initialize -c the derivative f1 and the vector p = W'd (for theta = 1). - - if (sbgnrm .le. zero) then - if (iprint .ge. 0) write (6,*) 'Subgnorm = 0. GCP = X.' - call dcopy(n,x,1,xcp,1) - return - endif - bnded = .true. - nfree = n + 1 - nbreak = 0 - ibkmin = 0 - bkmin = zero - col2 = 2*col - f1 = zero - if (iprint .ge. 99) write (6,3010) - -c We set p to zero and build it up as we determine d. - - do 20 i = 1, col2 - p(i) = zero - 20 continue - -c In the following loop we determine for each variable its bound -c status and its breakpoint, and update p accordingly. -c Smallest breakpoint is identified. - - do 50 i = 1, n - neggi = -g(i) - if (iwhere(i) .ne. 3 .and. iwhere(i) .ne. -1) then -c if x(i) is not a constant and has bounds, -c compute the difference between x(i) and its bounds. - if (nbd(i) .le. 2) tl = x(i) - l(i) - if (nbd(i) .ge. 2) tu = u(i) - x(i) - -c If a variable is close enough to a bound -c we treat it as at bound. - xlower = nbd(i) .le. 2 .and. tl .le. zero - xupper = nbd(i) .ge. 2 .and. tu .le. zero - -c reset iwhere(i). - iwhere(i) = 0 - if (xlower) then - if (neggi .le. zero) iwhere(i) = 1 - else if (xupper) then - if (neggi .ge. zero) iwhere(i) = 2 - else - if (abs(neggi) .le. zero) iwhere(i) = -3 - endif - endif - pointr = head - if (iwhere(i) .ne. 0 .and. iwhere(i) .ne. -1) then - d(i) = zero - else - d(i) = neggi - f1 = f1 - neggi*neggi -c calculate p := p - W'e_i* (g_i). - do 40 j = 1, col - p(j) = p(j) + wy(i,pointr)* neggi - p(col + j) = p(col + j) + ws(i,pointr)*neggi - pointr = mod(pointr,m) + 1 - 40 continue - if (nbd(i) .le. 2 .and. nbd(i) .ne. 0 - + .and. neggi .lt. zero) then -c x(i) + d(i) is bounded; compute t(i). - nbreak = nbreak + 1 - iorder(nbreak) = i - t(nbreak) = tl/(-neggi) - if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then - bkmin = t(nbreak) - ibkmin = nbreak - endif - else if (nbd(i) .ge. 2 .and. neggi .gt. zero) then -c x(i) + d(i) is bounded; compute t(i). - nbreak = nbreak + 1 - iorder(nbreak) = i - t(nbreak) = tu/neggi - if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then - bkmin = t(nbreak) - ibkmin = nbreak - endif - else -c x(i) + d(i) is not bounded. - nfree = nfree - 1 - iorder(nfree) = i - if (abs(neggi) .gt. zero) bnded = .false. - endif - endif - 50 continue - -c The indices of the nonzero components of d are now stored -c in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n). -c The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin. - - if (theta .ne. one) then -c complete the initialization of p for theta not= one. - call dscal(col,theta,p(col+1),1) - endif - -c Initialize GCP xcp = x. - - call dcopy(n,x,1,xcp,1) - - if (nbreak .eq. 0 .and. nfree .eq. n + 1) then -c is a zero vector, return with the initial xcp as GCP. - if (iprint .gt. 100) write (6,1010) (xcp(i), i = 1, n) - return - endif - -c Initialize c = W'(xcp - x) = 0. - - do 60 j = 1, col2 - c(j) = zero - 60 continue - -c Initialize derivative f2. - - f2 = -theta*f1 - if (col .gt. 0) then - call bmv(m,sy,wt,col,p,v,info) - if (info .ne. 0) return - f2 = f2 - ddot(col2,v,1,p,1) - endif - dtm = -f1/f2 - tsum = zero - nint = 1 - if (iprint .ge. 99) - + write (6,*) 'There are ',nbreak,' breakpoints ' - -c If there are no breakpoints, locate the GCP and return. - - if (nbreak .eq. 0) goto 888 - - nleft = nbreak - iter = 1 - - - tj = zero - -c------------------- the beginning of the loop ------------------------- - - 777 continue - -c Find the next smallest breakpoint; -c compute dt = t(nleft) - t(nleft + 1). - - tj0 = tj - if (iter .eq. 1) then -c Since we already have the smallest breakpoint we need not do -c heapsort yet. Often only one breakpoint is used and the -c cost of heapsort is avoided. - tj = bkmin - ibp = iorder(ibkmin) - else - if (iter .eq. 2) then -c Replace the already used smallest breakpoint with the -c breakpoint numbered nbreak > nlast, before heapsort call. - if (ibkmin .ne. nbreak) then - t(ibkmin) = t(nbreak) - iorder(ibkmin) = iorder(nbreak) - endif -c Update heap structure of breakpoints -c (if iter=2, initialize heap). - endif - call hpsolb(nleft,t,iorder,iter-2) - tj = t(nleft) - ibp = iorder(nleft) - endif - - dt = tj - tj0 - - if (dt .ne. zero .and. iprint .ge. 100) then - write (6,4011) nint,f1,f2 - write (6,5010) dt - write (6,6010) dtm - endif - -c If a minimizer is within this interval, locate the GCP and return. - - if (dtm .lt. dt) goto 888 - -c Otherwise fix one variable and -c reset the corresponding component of d to zero. - - tsum = tsum + dt - nleft = nleft - 1 - iter = iter + 1 - dibp = d(ibp) - d(ibp) = zero - if (dibp .gt. zero) then - zibp = u(ibp) - x(ibp) - xcp(ibp) = u(ibp) - iwhere(ibp) = 2 - else - zibp = l(ibp) - x(ibp) - xcp(ibp) = l(ibp) - iwhere(ibp) = 1 - endif - if (iprint .ge. 100) write (6,*) 'Variable ',ibp,' is fixed.' - if (nleft .eq. 0 .and. nbreak .eq. n) then -c all n variables are fixed, -c return with xcp as GCP. - dtm = dt - goto 999 - endif - -c Update the derivative information. - - nint = nint + 1 - dibp2 = dibp**2 - -c Update f1 and f2. - -c temporarily set f1 and f2 for col=0. - f1 = f1 + dt*f2 + dibp2 - theta*dibp*zibp - f2 = f2 - theta*dibp2 - - if (col .gt. 0) then -c update c = c + dt*p. - call daxpy(col2,dt,p,1,c,1) - -c choose wbp, -c the row of W corresponding to the breakpoint encountered. - pointr = head - do 70 j = 1,col - wbp(j) = wy(ibp,pointr) - wbp(col + j) = theta*ws(ibp,pointr) - pointr = mod(pointr,m) + 1 - 70 continue - -c compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'. - call bmv(m,sy,wt,col,wbp,v,info) - if (info .ne. 0) return - wmc = ddot(col2,c,1,v,1) - wmp = ddot(col2,p,1,v,1) - wmw = ddot(col2,wbp,1,v,1) - -c update p = p - dibp*wbp. - call daxpy(col2,-dibp,wbp,1,p,1) - -c complete updating f1 and f2 while col > 0. - f1 = f1 + dibp*wmc - f2 = f2 + 2.0d0*dibp*wmp - dibp2*wmw - endif - - if (nleft .gt. 0) then - dtm = -f1/f2 - goto 777 -c to repeat the loop for unsearched intervals. - else if(bnded) then - f1 = zero - f2 = zero - dtm = zero - else - dtm = -f1/f2 - endif - -c------------------- the end of the loop ------------------------------- - - 888 continue - if (iprint .ge. 99) then - write (6,*) - write (6,*) 'GCP found in this segment' - write (6,4010) nint,f1,f2 - write (6,6010) dtm - endif - if (dtm .le. zero) dtm = zero - tsum = tsum + dtm - -c Move free variables (i.e., the ones w/o breakpoints) and -c the variables whose breakpoints haven't been reached. - - call daxpy(n,tsum,d,1,xcp,1) - - 999 continue - -c Update c = c + dtm*p = W'(x^c - x) -c which will be used in computing r = Z'(B(x^c - x) + g). - - if (col .gt. 0) call daxpy(col2,dtm,p,1,c,1) - if (iprint .gt. 100) write (6,1010) (xcp(i),i = 1,n) - if (iprint .ge. 99) write (6,2010) - - 1010 format ('Cauchy X = ',/,(4x,1p,6(1x,d11.4))) - 2010 format (/,'---------------- exit CAUCHY----------------------',/) - 3010 format (/,'---------------- CAUCHY entered-------------------') - 4010 format ('Piece ',i3,' --f1, f2 at start point ',1p,2(1x,d11.4)) - 4011 format (/,'Piece ',i3,' --f1, f2 at start point ', - + 1p,2(1x,d11.4)) - 5010 format ('Distance to the next break point = ',1p,d11.4) - 6010 format ('Distance to the stationary point = ',1p,d11.4) - - return - - end - -c====================== The end of cauchy ============================== - - subroutine cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index, - + theta, col, head, nfree, cnstnd, info) - - logical cnstnd - integer n, m, col, head, nfree, info, index(n) - double precision theta, - + x(n), g(n), z(n), r(n), wa(4*m), - + ws(n, m), wy(n, m), sy(m, m), wt(m, m) - -c ************ -c -c Subroutine cmprlb -c -c This subroutine computes r=-Z'B(xcp-xk)-Z'g by using -c wa(2m+1)=W'(xcp-x) from subroutine cauchy. -c -c Subprograms called: -c -c L-BFGS-B Library ... bmv. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer i,j,k,pointr - double precision a1,a2 - - if (.not. cnstnd .and. col .gt. 0) then - do 26 i = 1, n - r(i) = -g(i) - 26 continue - else - do 30 i = 1, nfree - k = index(i) - r(i) = -theta*(z(k) - x(k)) - g(k) - 30 continue - call bmv(m,sy,wt,col,wa(2*m+1),wa(1),info) - if (info .ne. 0) then - info = -8 - return - endif - pointr = head - do 34 j = 1, col - a1 = wa(j) - a2 = theta*wa(col + j) - do 32 i = 1, nfree - k = index(i) - r(i) = r(i) + wy(k,pointr)*a1 + ws(k,pointr)*a2 - 32 continue - pointr = mod(pointr,m) + 1 - 34 continue - endif - - return - - end - -c======================= The end of cmprlb ============================= - - subroutine errclb(n, m, factr, l, u, nbd, task, info, k) - - character*60 task - integer n, m, info, k, nbd(n) - double precision factr, l(n), u(n) - -c ************ -c -c Subroutine errclb -c -c This subroutine checks the validity of the input data. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer i - double precision one,zero - parameter (one=1.0d0,zero=0.0d0) - -c Check the input arguments for errors. - - if (n .le. 0) task = 'ERROR: N .LE. 0' - if (m .le. 0) task = 'ERROR: M .LE. 0' - if (factr .lt. zero) task = 'ERROR: FACTR .LT. 0' - -c Check the validity of the arrays nbd(i), u(i), and l(i). - - do 10 i = 1, n - if (nbd(i) .lt. 0 .or. nbd(i) .gt. 3) then -c return - task = 'ERROR: INVALID NBD' - info = -6 - k = i - endif - if (nbd(i) .eq. 2) then - if (l(i) .gt. u(i)) then -c return - task = 'ERROR: NO FEASIBLE SOLUTION' - info = -7 - k = i - endif - endif - 10 continue - - return - - end - -c======================= The end of errclb ============================= - - subroutine formk(n, nsub, ind, nenter, ileave, indx2, iupdat, - + updatd, wn, wn1, m, ws, wy, sy, theta, col, - + head, info) - - integer n, nsub, m, col, head, nenter, ileave, iupdat, - + info, ind(n), indx2(n) - double precision theta, wn(2*m, 2*m), wn1(2*m, 2*m), - + ws(n, m), wy(n, m), sy(m, m) - logical updatd - -c ************ -c -c Subroutine formk -c -c This subroutine forms the LEL^T factorization of the indefinite -c -c matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] -c [L_a -R_z theta*S'AA'S ] -c where E = [-I 0] -c [ 0 I] -c The matrix K can be shown to be equal to the matrix M^[-1]N -c occurring in section 5.1 of [1], as well as to the matrix -c Mbar^[-1] Nbar in section 5.3. -c -c n is an integer variable. -c On entry n is the dimension of the problem. -c On exit n is unchanged. -c -c nsub is an integer variable -c On entry nsub is the number of subspace variables in free set. -c On exit nsub is not changed. -c -c ind is an integer array of dimension nsub. -c On entry ind specifies the indices of subspace variables. -c On exit ind is unchanged. -c -c nenter is an integer variable. -c On entry nenter is the number of variables entering the -c free set. -c On exit nenter is unchanged. -c -c ileave is an integer variable. -c On entry indx2(ileave),...,indx2(n) are the variables leaving -c the free set. -c On exit ileave is unchanged. -c -c indx2 is an integer array of dimension n. -c On entry indx2(1),...,indx2(nenter) are the variables entering -c the free set, while indx2(ileave),...,indx2(n) are the -c variables leaving the free set. -c On exit indx2 is unchanged. -c -c iupdat is an integer variable. -c On entry iupdat is the total number of BFGS updates made so far. -c On exit iupdat is unchanged. -c -c updatd is a logical variable. -c On entry 'updatd' is true if the L-BFGS matrix is updatd. -c On exit 'updatd' is unchanged. -c -c wn is a double precision array of dimension 2m x 2m. -c On entry wn is unspecified. -c On exit the upper triangle of wn stores the LEL^T factorization -c of the 2*col x 2*col indefinite matrix -c [-D -Y'ZZ'Y/theta L_a'-R_z' ] -c [L_a -R_z theta*S'AA'S ] -c -c wn1 is a double precision array of dimension 2m x 2m. -c On entry wn1 stores the lower triangular part of -c [Y' ZZ'Y L_a'+R_z'] -c [L_a+R_z S'AA'S ] -c in the previous iteration. -c On exit wn1 stores the corresponding updated matrices. -c The purpose of wn1 is just to store these inner products -c so they can be easily updated and inserted into wn. -c -c m is an integer variable. -c On entry m is the maximum number of variable metric corrections -c used to define the limited memory matrix. -c On exit m is unchanged. -c -c ws, wy, sy, and wtyy are double precision arrays; -c theta is a double precision variable; -c col is an integer variable; -c head is an integer variable. -c On entry they store the information defining the -c limited memory BFGS matrix: -c ws(n,m) stores S, a set of s-vectors; -c wy(n,m) stores Y, a set of y-vectors; -c sy(m,m) stores S'Y; -c wtyy(m,m) stores the Cholesky factorization -c of (theta*S'S+LD^(-1)L') -c theta is the scaling factor specifying B_0 = theta I; -c col is the number of variable metric corrections stored; -c head is the location of the 1st s- (or y-) vector in S (or Y). -c On exit they are unchanged. -c -c info is an integer variable. -c On entry info is unspecified. -c On exit info = 0 for normal return; -c = -1 when the 1st Cholesky factorization failed; -c = -2 when the 2st Cholesky factorization failed. -c -c Subprograms called: -c -c Linpack ... dcopy, dpofa, dtrsl. -c -c -c References: -c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited -c memory algorithm for bound constrained optimization'', -c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. -c -c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a -c limited memory FORTRAN code for solving bound constrained -c optimization problems'', Tech. Report, NAM-11, EECS Department, -c Northwestern University, 1994. -c -c (Postscript files of these papers are available via anonymous -c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer m2,ipntr,jpntr,iy,is,jy,js,is1,js1,k1,i,k, - + col2,pbegin,pend,dbegin,dend,upcl - double precision ddot,temp1,temp2,temp3,temp4 - double precision one,zero - parameter (one=1.0d0,zero=0.0d0) - -c Form the lower triangular part of -c WN1 = [Y' ZZ'Y L_a'+R_z'] -c [L_a+R_z S'AA'S ] -c where L_a is the strictly lower triangular part of S'AA'Y -c R_z is the upper triangular part of S'ZZ'Y. - - if (updatd) then - if (iupdat .gt. m) then -c shift old part of WN1. - do 10 jy = 1, m - 1 - js = m + jy - call dcopy(m-jy,wn1(jy+1,jy+1),1,wn1(jy,jy),1) - call dcopy(m-jy,wn1(js+1,js+1),1,wn1(js,js),1) - call dcopy(m-1,wn1(m+2,jy+1),1,wn1(m+1,jy),1) - 10 continue - endif - -c put new rows in blocks (1,1), (2,1) and (2,2). - pbegin = 1 - pend = nsub - dbegin = nsub + 1 - dend = n - iy = col - is = m + col - ipntr = head + col - 1 - if (ipntr .gt. m) ipntr = ipntr - m - jpntr = head - do 20 jy = 1, col - js = m + jy - temp1 = zero - temp2 = zero - temp3 = zero -c compute element jy of row 'col' of Y'ZZ'Y - do 15 k = pbegin, pend - k1 = ind(k) - temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr) - 15 continue -c compute elements jy of row 'col' of L_a and S'AA'S - do 16 k = dbegin, dend - k1 = ind(k) - temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr) - temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr) - 16 continue - wn1(iy,jy) = temp1 - wn1(is,js) = temp2 - wn1(is,jy) = temp3 - jpntr = mod(jpntr,m) + 1 - 20 continue - -c put new column in block (2,1). - jy = col - jpntr = head + col - 1 - if (jpntr .gt. m) jpntr = jpntr - m - ipntr = head - do 30 i = 1, col - is = m + i - temp3 = zero -c compute element i of column 'col' of R_z - do 25 k = pbegin, pend - k1 = ind(k) - temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr) - 25 continue - ipntr = mod(ipntr,m) + 1 - wn1(is,jy) = temp3 - 30 continue - upcl = col - 1 - else - upcl = col - endif - -c modify the old parts in blocks (1,1) and (2,2) due to changes -c in the set of free variables. - ipntr = head - do 45 iy = 1, upcl - is = m + iy - jpntr = head - do 40 jy = 1, iy - js = m + jy - temp1 = zero - temp2 = zero - temp3 = zero - temp4 = zero - do 35 k = 1, nenter - k1 = indx2(k) - temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr) - temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr) - 35 continue - do 36 k = ileave, n - k1 = indx2(k) - temp3 = temp3 + wy(k1,ipntr)*wy(k1,jpntr) - temp4 = temp4 + ws(k1,ipntr)*ws(k1,jpntr) - 36 continue - wn1(iy,jy) = wn1(iy,jy) + temp1 - temp3 - wn1(is,js) = wn1(is,js) - temp2 + temp4 - jpntr = mod(jpntr,m) + 1 - 40 continue - ipntr = mod(ipntr,m) + 1 - 45 continue - -c modify the old parts in block (2,1). - ipntr = head - do 60 is = m + 1, m + upcl - jpntr = head - do 55 jy = 1, upcl - temp1 = zero - temp3 = zero - do 50 k = 1, nenter - k1 = indx2(k) - temp1 = temp1 + ws(k1,ipntr)*wy(k1,jpntr) - 50 continue - do 51 k = ileave, n - k1 = indx2(k) - temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr) - 51 continue - if (is .le. jy + m) then - wn1(is,jy) = wn1(is,jy) + temp1 - temp3 - else - wn1(is,jy) = wn1(is,jy) - temp1 + temp3 - endif - jpntr = mod(jpntr,m) + 1 - 55 continue - ipntr = mod(ipntr,m) + 1 - 60 continue - -c Form the upper triangle of WN = [D+Y' ZZ'Y/theta -L_a'+R_z' ] -c [-L_a +R_z S'AA'S*theta] - - m2 = 2*m - do 70 iy = 1, col - is = col + iy - is1 = m + iy - do 65 jy = 1, iy - js = col + jy - js1 = m + jy - wn(jy,iy) = wn1(iy,jy)/theta - wn(js,is) = wn1(is1,js1)*theta - 65 continue - do 66 jy = 1, iy - 1 - wn(jy,is) = -wn1(is1,jy) - 66 continue - do 67 jy = iy, col - wn(jy,is) = wn1(is1,jy) - 67 continue - wn(iy,iy) = wn(iy,iy) + sy(iy,iy) - 70 continue - -c Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')] -c [(-L_a +R_z)L'^-1 S'AA'S*theta ] - -c first Cholesky factor (1,1) block of wn to get LL' -c with L' stored in the upper triangle of wn. - call dpofa(wn,m2,col,info) - if (info .ne. 0) then - info = -1 - return - endif -c then form L^-1(-L_a'+R_z') in the (1,2) block. - col2 = 2*col - do 71 js = col+1 ,col2 - call dtrsl(wn,m2,col,wn(1,js),11,info) - 71 continue - -c Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the -c upper triangle of (2,2) block of wn. - - - do 72 is = col+1, col2 - do 74 js = is, col2 - wn(is,js) = wn(is,js) + ddot(col,wn(1,is),1,wn(1,js),1) - 74 continue - 72 continue - -c Cholesky factorization of (2,2) block of wn. - - call dpofa(wn(col+1,col+1),m2,col,info) - if (info .ne. 0) then - info = -2 - return - endif - - return - - end - -c======================= The end of formk ============================== - - subroutine formt(m, wt, sy, ss, col, theta, info) - - integer m, col, info - double precision theta, wt(m, m), sy(m, m), ss(m, m) - -c ************ -c -c Subroutine formt -c -c This subroutine forms the upper half of the pos. def. and symm. -c T = theta*SS + L*D^(-1)*L', stores T in the upper triangle -c of the array wt, and performs the Cholesky factorization of T -c to produce J*J', with J' stored in the upper triangle of wt. -c -c Subprograms called: -c -c Linpack ... dpofa. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer i,j,k,k1 - double precision ddum - double precision zero - parameter (zero=0.0d0) - - -c Form the upper half of T = theta*SS + L*D^(-1)*L', -c store T in the upper triangle of the array wt. - - do 52 j = 1, col - wt(1,j) = theta*ss(1,j) - 52 continue - do 55 i = 2, col - do 54 j = i, col - k1 = min(i,j) - 1 - ddum = zero - do 53 k = 1, k1 - ddum = ddum + sy(i,k)*sy(j,k)/sy(k,k) - 53 continue - wt(i,j) = ddum + theta*ss(i,j) - 54 continue - 55 continue - -c Cholesky factorize T to J*J' with -c J' stored in the upper triangle of wt. - - call dpofa(wt,m,col,info) - if (info .ne. 0) then - info = -3 - endif - - return - - end - -c======================= The end of formt ============================== - - subroutine freev(n, nfree, index, nenter, ileave, indx2, - + iwhere, wrk, updatd, cnstnd, iprint, iter) - - integer n, nfree, nenter, ileave, iprint, iter, - + index(n), indx2(n), iwhere(n) - logical wrk, updatd, cnstnd - -c ************ -c -c Subroutine freev -c -c This subroutine counts the entering and leaving variables when -c iter > 0, and finds the index set of free and active variables -c at the GCP. -c -c cnstnd is a logical variable indicating whether bounds are present -c -c index is an integer array of dimension n -c for i=1,...,nfree, index(i) are the indices of free variables -c for i=nfree+1,...,n, index(i) are the indices of bound variables -c On entry after the first iteration, index gives -c the free variables at the previous iteration. -c On exit it gives the free variables based on the determination -c in cauchy using the array iwhere. -c -c indx2 is an integer array of dimension n -c On entry indx2 is unspecified. -c On exit with iter>0, indx2 indicates which variables -c have changed status since the previous iteration. -c For i= 1,...,nenter, indx2(i) have changed from bound to free. -c For i= ileave+1,...,n, indx2(i) have changed from free to bound. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer iact,i,k - - nenter = 0 - ileave = n + 1 - if (iter .gt. 0 .and. cnstnd) then -c count the entering and leaving variables. - do 20 i = 1, nfree - k = index(i) - if (iwhere(k) .gt. 0) then - ileave = ileave - 1 - indx2(ileave) = k - if (iprint .ge. 100) write (6,*) - + 'Variable ',k,' leaves the set of free variables' - endif - 20 continue - do 22 i = 1 + nfree, n - k = index(i) - if (iwhere(k) .le. 0) then - nenter = nenter + 1 - indx2(nenter) = k - if (iprint .ge. 100) write (6,*) - + 'Variable ',k,' enters the set of free variables' - endif - 22 continue - if (iprint .ge. 99) write (6,*) - + n+1-ileave,' variables leave; ',nenter,' variables enter' - endif - wrk = (ileave .lt. n+1) .or. (nenter .gt. 0) .or. updatd - -c Find the index set of free and active variables at the GCP. - - nfree = 0 - iact = n + 1 - do 24 i = 1, n - if (iwhere(i) .le. 0) then - nfree = nfree + 1 - index(nfree) = i - else - iact = iact - 1 - index(iact) = i - endif - 24 continue - if (iprint .ge. 99) write (6,*) - + nfree,' variables are free at GCP ',iter + 1 - - return - - end - -c======================= The end of freev ============================== - - subroutine hpsolb(n, t, iorder, iheap) - integer iheap, n, iorder(n) - double precision t(n) - -c ************ -c -c Subroutine hpsolb -c -c This subroutine sorts out the least element of t, and puts the -c remaining elements of t in a heap. -c -c n is an integer variable. -c On entry n is the dimension of the arrays t and iorder. -c On exit n is unchanged. -c -c t is a double precision array of dimension n. -c On entry t stores the elements to be sorted, -c On exit t(n) stores the least elements of t, and t(1) to t(n-1) -c stores the remaining elements in the form of a heap. -c -c iorder is an integer array of dimension n. -c On entry iorder(i) is the index of t(i). -c On exit iorder(i) is still the index of t(i), but iorder may be -c permuted in accordance with t. -c -c iheap is an integer variable specifying the task. -c On entry iheap should be set as follows: -c iheap .eq. 0 if t(1) to t(n) is not in the form of a heap, -c iheap .ne. 0 if otherwise. -c On exit iheap is unchanged. -c -c -c References: -c Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT. -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c ************ - - integer i,j,k,indxin,indxou - double precision ddum,out - - if (iheap .eq. 0) then - -c Rearrange the elements t(1) to t(n) to form a heap. - - do 20 k = 2, n - ddum = t(k) - indxin = iorder(k) - -c Add ddum to the heap. - i = k - 10 continue - if (i.gt.1) then - j = i/2 - if (ddum .lt. t(j)) then - t(i) = t(j) - iorder(i) = iorder(j) - i = j - goto 10 - endif - endif - t(i) = ddum - iorder(i) = indxin - 20 continue - endif - -c Assign to 'out' the value of t(1), the least member of the heap, -c and rearrange the remaining members to form a heap as -c elements 1 to n-1 of t. - - if (n .gt. 1) then - i = 1 - out = t(1) - indxou = iorder(1) - ddum = t(n) - indxin = iorder(n) - -c Restore the heap - 30 continue - j = i+i - if (j .le. n-1) then - if (t(j+1) .lt. t(j)) j = j+1 - if (t(j) .lt. ddum ) then - t(i) = t(j) - iorder(i) = iorder(j) - i = j - goto 30 - endif - endif - t(i) = ddum - iorder(i) = indxin - -c Put the least member in t(n). - - t(n) = out - iorder(n) = indxou - endif - - return - - end - -c====================== The end of hpsolb ============================== - - subroutine 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, dsave) - - character*60 task, csave - logical boxed, cnstnd - integer n, iter, ifun, iback, nfgv, info, - + nbd(n), isave(2) - double precision f, fold, gd, gdold, stp, dnorm, dtd, xstep, - + stpmx, x(n), l(n), u(n), g(n), d(n), r(n), t(n), - + z(n), dsave(13) -c ********** -c -c Subroutine lnsrlb -c -c This subroutine calls subroutine dcsrch from the Minpack2 library -c to perform the line search. Subroutine dscrch is safeguarded so -c that all trial points lie within the feasible region. -c -c Subprograms called: -c -c Minpack2 Library ... dcsrch. -c -c Linpack ... dtrsl, ddot. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ********** - - integer i - double precision ddot,a1,a2 - double precision one,zero,big - parameter (one=1.0d0,zero=0.0d0,big=1.0d+10) - double precision ftol,gtol,xtol - parameter (ftol=1.0d-3,gtol=0.9d0,xtol=0.1d0) - - if (task(1:5) .eq. 'FG_LN') goto 556 - - dtd = ddot(n,d,1,d,1) - dnorm = sqrt(dtd) - -c Determine the maximum step length. - - stpmx = big - if (cnstnd) then - if (iter .eq. 0) then - stpmx = one - else - do 43 i = 1, n - a1 = d(i) - if (nbd(i) .ne. 0) then - if (a1 .lt. zero .and. nbd(i) .le. 2) then - a2 = l(i) - x(i) - if (a2 .ge. zero) then - stpmx = zero - else if (a1*stpmx .lt. a2) then - stpmx = a2/a1 - endif - else if (a1 .gt. zero .and. nbd(i) .ge. 2) then - a2 = u(i) - x(i) - if (a2 .le. zero) then - stpmx = zero - else if (a1*stpmx .gt. a2) then - stpmx = a2/a1 - endif - endif - endif - 43 continue - endif - endif - - if (iter .eq. 0 .and. .not. boxed) then - stp = min(one/dnorm, stpmx) - else - stp = one - endif - - call dcopy(n,x,1,t,1) - call dcopy(n,g,1,r,1) - fold = f - ifun = 0 - iback = 0 - csave = 'START' - 556 continue - gd = ddot(n,g,1,d,1) - if (ifun .eq. 0) then - gdold=gd - if (gd .ge. zero) then -c the directional derivative >=0. -c Line search is impossible. - info = -4 - return - endif - endif - - call dcsrch(f,gd,stp,ftol,gtol,xtol,zero,stpmx,csave,isave,dsave) - - xstep = stp*dnorm - if (csave(1:4) .ne. 'CONV' .and. csave(1:4) .ne. 'WARN') then - task = 'FG_LNSRCH' - ifun = ifun + 1 - nfgv = nfgv + 1 - iback = ifun - 1 - if (stp .eq. one) then - call dcopy(n,z,1,x,1) - else - do 41 i = 1, n - x(i) = stp*d(i) + t(i) - 41 continue - endif - else - task = 'NEW_X' - endif - - return - - end - -c======================= The end of lnsrlb ============================= - - subroutine matupd(n, m, ws, wy, sy, ss, d, r, itail, - + iupdat, col, head, theta, rr, dr, stp, dtd) - - integer n, m, itail, iupdat, col, head - double precision theta, rr, dr, stp, dtd, d(n), r(n), - + ws(n, m), wy(n, m), sy(m, m), ss(m, m) - -c ************ -c -c Subroutine matupd -c -c This subroutine updates matrices WS and WY, and forms the -c middle matrix in B. -c -c Subprograms called: -c -c Linpack ... dcopy, ddot. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer j,pointr - double precision ddot - double precision one - parameter (one=1.0d0) - -c Set pointers for matrices WS and WY. - - if (iupdat .le. m) then - col = iupdat - itail = mod(head+iupdat-2,m) + 1 - else - itail = mod(itail,m) + 1 - head = mod(head,m) + 1 - endif - -c Update matrices WS and WY. - - call dcopy(n,d,1,ws(1,itail),1) - call dcopy(n,r,1,wy(1,itail),1) - -c Set theta=yy/ys. - - theta = rr/dr - -c Form the middle matrix in B. - -c update the upper triangle of SS, -c and the lower triangle of SY: - if (iupdat .gt. m) then -c move old information - do 50 j = 1, col - 1 - call dcopy(j,ss(2,j+1),1,ss(1,j),1) - call dcopy(col-j,sy(j+1,j+1),1,sy(j,j),1) - 50 continue - endif -c add new information: the last row of SY -c and the last column of SS: - pointr = head - do 51 j = 1, col - 1 - sy(col,j) = ddot(n,d,1,wy(1,pointr),1) - ss(j,col) = ddot(n,ws(1,pointr),1,d,1) - pointr = mod(pointr,m) + 1 - 51 continue - if (stp .eq. one) then - ss(col,col) = dtd - else - ss(col,col) = stp*stp*dtd - endif - sy(col,col) = dr - - return - - end - -c======================= The end of matupd ============================= - - subroutine prn1lb(n, m, l, u, x, iprint, itfile, epsmch) - - integer n, m, iprint, itfile - double precision epsmch, x(n), l(n), u(n) - -c ************ -c -c Subroutine prn1lb -c -c This subroutine prints the input data, initial point, upper and -c lower bounds of each variable, machine precision, as well as -c the headings of the output. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer i - - if (iprint .ge. 0) then - write (6,7001) epsmch - write (6,*) 'N = ',n,' M = ',m - if (iprint .ge. 1) then - write (itfile,2001) epsmch - write (itfile,*)'N = ',n,' M = ',m - write (itfile,9001) - if (iprint .gt. 100) then - write (6,1004) 'L =',(l(i),i = 1,n) - write (6,1004) 'X0 =',(x(i),i = 1,n) - write (6,1004) 'U =',(u(i),i = 1,n) - endif - endif - endif - - 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4))) - 2001 format ('RUNNING THE L-BFGS-B CODE',/,/, - + 'it = iteration number',/, - + 'nf = number of function evaluations',/, - + 'nint = number of segments explored during the Cauchy search',/, - + 'nact = number of active bounds at the generalized Cauchy point' - + ,/, - + 'sub = manner in which the subspace minimization terminated:' - + ,/,' con = converged, bnd = a bound was reached',/, - + 'itls = number of iterations performed in the line search',/, - + 'stepl = step length used',/, - + 'tstep = norm of the displacement (total step)',/, - + 'projg = norm of the projected gradient',/, - + 'f = function value',/,/, - + ' * * *',/,/, - + 'Machine precision =',1p,d10.3) - 7001 format ('RUNNING THE L-BFGS-B CODE',/,/, - + ' * * *',/,/, - + 'Machine precision =',1p,d10.3) - 9001 format (/,3x,'it',3x,'nf',2x,'nint',2x,'nact',2x,'sub',2x,'itls', - + 2x,'stepl',4x,'tstep',5x,'projg',8x,'f') - - return - - end - -c======================= The end of prn1lb ============================= - - subroutine prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact, - + sbgnrm, nint, word, iword, iback, stp, xstep) - - character*3 word - integer n, iprint, itfile, iter, nfgv, nact, nint, - + iword, iback - double precision f, sbgnrm, stp, xstep, x(n), g(n) - -c ************ -c -c Subroutine prn2lb -c -c This subroutine prints out new information after a successful -c line search. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer i,imod - -c 'word' records the status of subspace solutions. - if (iword .eq. 0) then -c the subspace minimization converged. - word = 'con' - else if (iword .eq. 1) then -c the subspace minimization stopped at a bound. - word = 'bnd' - else if (iword .eq. 5) then -c the truncated Newton step has been used. - word = 'TNT' - else - word = '---' - endif - if (iprint .ge. 99) then - write (6,*) 'LINE SEARCH',iback,' times; norm of step = ',xstep - write (6,2001) iter,f,sbgnrm - if (iprint .gt. 100) then - write (6,1004) 'X =',(x(i), i = 1, n) - write (6,1004) 'G =',(g(i), i = 1, n) - endif - else if (iprint .gt. 0) then - imod = mod(iter,iprint) - if (imod .eq. 0) write (6,2001) iter,f,sbgnrm - endif - if (iprint .ge. 1) write (itfile,3001) - + iter,nfgv,nint,nact,word,iback,stp,xstep,sbgnrm,f - - 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4))) - 2001 format - + (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5) - 3001 format(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),1p,2(1x,d10.3)) - - return - - end - -c======================= The end of prn2lb ============================= - - subroutine prn3lb(n, x, f, task, iprint, info, itfile, - + iter, nfgv, nintol, nskip, nact, sbgnrm, - + time, nint, word, iback, stp, xstep, k, - + cachyt, sbtime, lnscht) - - character*60 task - character*3 word - integer n, iprint, info, itfile, iter, nfgv, nintol, - + nskip, nact, nint, iback, k - double precision f, sbgnrm, time, stp, xstep, cachyt, sbtime, - + lnscht, x(n) - -c ************ -c -c Subroutine prn3lb -c -c This subroutine prints out information when either a built-in -c convergence test is satisfied or when an error message is -c generated. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer i - - if (task(1:5) .eq. 'ERROR') goto 999 - - if (iprint .ge. 0) then - write (6,3003) - write (6,3004) - write(6,3005) n,iter,nfgv,nintol,nskip,nact,sbgnrm,f - if (iprint .ge. 100) then - write (6,1004) 'X =',(x(i),i = 1,n) - endif - if (iprint .ge. 1) write (6,*) ' F =',f - endif - 999 continue - if (iprint .ge. 0) then - write (6,3009) task - if (info .ne. 0) then - if (info .eq. -1) write (6,9011) - if (info .eq. -2) write (6,9012) - if (info .eq. -3) write (6,9013) - if (info .eq. -4) write (6,9014) - if (info .eq. -5) write (6,9015) - if (info .eq. -6) write (6,*)' Input nbd(',k,') is invalid.' - if (info .eq. -7) - + write (6,*)' l(',k,') > u(',k,'). No feasible solution.' - if (info .eq. -8) write (6,9018) - if (info .eq. -9) write (6,9019) - endif - if (iprint .ge. 1) write (6,3007) cachyt,sbtime,lnscht - write (6,3008) time - if (iprint .ge. 1) then - if (info .eq. -4 .or. info .eq. -9) then - write (itfile,3002) - + iter,nfgv,nint,nact,word,iback,stp,xstep - endif - write (itfile,3009) task - if (info .ne. 0) then - if (info .eq. -1) write (itfile,9011) - if (info .eq. -2) write (itfile,9012) - if (info .eq. -3) write (itfile,9013) - if (info .eq. -4) write (itfile,9014) - if (info .eq. -5) write (itfile,9015) - if (info .eq. -8) write (itfile,9018) - if (info .eq. -9) write (itfile,9019) - endif - write (itfile,3008) time - endif - endif - - 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4))) - 3002 format(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),6x,'-',10x,'-') - 3003 format (/, - + ' * * *',/,/, - + 'Tit = total number of iterations',/, - + 'Tnf = total number of function evaluations',/, - + 'Tnint = total number of segments explored during', - + ' Cauchy searches',/, - + 'Skip = number of BFGS updates skipped',/, - + 'Nact = number of active bounds at final generalized', - + ' Cauchy point',/, - + 'Projg = norm of the final projected gradient',/, - + 'F = final function value',/,/, - + ' * * *') - 3004 format (/,3x,'N',3x,'Tit',2x,'Tnf',2x,'Tnint',2x, - + 'Skip',2x,'Nact',5x,'Projg',8x,'F') - 3005 format (i5,2(1x,i4),(1x,i6),(2x,i4),(1x,i5),1p,2(2x,d10.3)) - 3006 format (i5,2(1x,i4),2(1x,i6),(1x,i4),(1x,i5),7x,'-',10x,'-') - 3007 format (/,' Cauchy time',1p,e10.3,' seconds.',/ - + ' Subspace minimization time',1p,e10.3,' seconds.',/ - + ' Line search time',1p,e10.3,' seconds.') - 3008 format (/,' Total User time',1p,e10.3,' seconds.',/) - 3009 format (/,a60) - 9011 format (/, - +' Matrix in 1st Cholesky factorization in formk is not Pos. Def.') - 9012 format (/, - +' Matrix in 2st Cholesky factorization in formk is not Pos. Def.') - 9013 format (/, - +' Matrix in the Cholesky factorization in formt is not Pos. Def.') - 9014 format (/, - +' Derivative >= 0, backtracking line search impossible.',/, - +' Previous x, f and g restored.',/, - +' Possible causes: 1 error in function or gradient evaluation;',/, - +' 2 rounding errors dominate computation.') - 9015 format (/, - +' Warning: more than 10 function and gradient',/, - +' evaluations in the last line search. Termination',/, - +' may possibly be caused by a bad search direction.') - 9018 format (/,' The triangular system is singular.') - 9019 format (/, - +' Line search cannot locate an adequate point after 20 function',/ - +,' and gradient evaluations. Previous x, f and g restored.',/, - +' Possible causes: 1 error in function or gradient evaluation;',/, - +' 2 rounding error dominate computation.') - - return - - end - -c======================= The end of prn3lb ============================= - - subroutine projgr(n, l, u, nbd, x, g, sbgnrm) - - integer n, nbd(n) - double precision sbgnrm, x(n), l(n), u(n), g(n) - -c ************ -c -c Subroutine projgr -c -c This subroutine computes the infinity norm of the projected -c gradient. -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer i - double precision gi - double precision one,zero - parameter (one=1.0d0,zero=0.0d0) - - sbgnrm = zero - do 15 i = 1, n - gi = g(i) - if (nbd(i) .ne. 0) then - if (gi .lt. zero) then - if (nbd(i) .ge. 2) gi = max((x(i)-u(i)),gi) - else - if (nbd(i) .le. 2) gi = min((x(i)-l(i)),gi) - endif - endif - sbgnrm = max(sbgnrm,abs(gi)) - 15 continue - - return - - end - -c======================= The end of projgr ============================= - - subroutine subsm(n, m, nsub, ind, l, u, nbd, x, d, ws, wy, theta, - + col, head, iword, wv, wn, iprint, info) - - integer n, m, nsub, col, head, iword, iprint, info, - + ind(nsub), nbd(n) - double precision theta, - + l(n), u(n), x(n), d(n), - + ws(n, m), wy(n, m), - + wv(2*m), wn(2*m, 2*m) - -c ************ -c -c Subroutine subsm -c -c Given xcp, l, u, r, an index set that specifies -c the active set at xcp, and an l-BFGS matrix B -c (in terms of WY, WS, SY, WT, head, col, and theta), -c this subroutine computes an approximate solution -c of the subspace problem -c -c (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp) -c -c subject to l<=x<=u -c x_i=xcp_i for all i in A(xcp) -c -c along the subspace unconstrained Newton direction -c -c d = -(Z'BZ)^(-1) r. -c -c The formula for the Newton direction, given the L-BFGS matrix -c and the Sherman-Morrison formula, is -c -c d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r. -c -c where -c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] -c [L_a -R_z theta*S'AA'S ] -c -c Note that this procedure for computing d differs -c from that described in [1]. One can show that the matrix K is -c equal to the matrix M^[-1]N in that paper. -c -c n is an integer variable. -c On entry n is the dimension of the problem. -c On exit n is unchanged. -c -c m is an integer variable. -c On entry m is the maximum number of variable metric corrections -c used to define the limited memory matrix. -c On exit m is unchanged. -c -c nsub is an integer variable. -c On entry nsub is the number of free variables. -c On exit nsub is unchanged. -c -c ind is an integer array of dimension nsub. -c On entry ind specifies the coordinate indices of free variables. -c On exit ind is unchanged. -c -c l is a double precision array of dimension n. -c On entry l is the lower bound of x. -c On exit l is unchanged. -c -c u is a double precision array of dimension n. -c On entry u is the upper bound of x. -c On exit u is unchanged. -c -c nbd is a integer array of dimension n. -c On entry nbd represents the type of bounds imposed on the -c variables, and must be specified as follows: -c nbd(i)=0 if x(i) is unbounded, -c 1 if x(i) has only a lower bound, -c 2 if x(i) has both lower and upper bounds, and -c 3 if x(i) has only an upper bound. -c On exit nbd is unchanged. -c -c x is a double precision array of dimension n. -c On entry x specifies the Cauchy point xcp. -c On exit x(i) is the minimizer of Q over the subspace of -c free variables. -c -c d is a double precision array of dimension n. -c On entry d is the reduced gradient of Q at xcp. -c On exit d is the Newton direction of Q. -c -c ws and wy are double precision arrays; -c theta is a double precision variable; -c col is an integer variable; -c head is an integer variable. -c On entry they store the information defining the -c limited memory BFGS matrix: -c ws(n,m) stores S, a set of s-vectors; -c wy(n,m) stores Y, a set of y-vectors; -c theta is the scaling factor specifying B_0 = theta I; -c col is the number of variable metric corrections stored; -c head is the location of the 1st s- (or y-) vector in S (or Y). -c On exit they are unchanged. -c -c iword is an integer variable. -c On entry iword is unspecified. -c On exit iword specifies the status of the subspace solution. -c iword = 0 if the solution is in the box, -c 1 if some bound is encountered. -c -c wv is a double precision working array of dimension 2m. -c -c wn is a double precision array of dimension 2m x 2m. -c On entry the upper triangle of wn stores the LEL^T factorization -c of the indefinite matrix -c -c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] -c [L_a -R_z theta*S'AA'S ] -c where E = [-I 0] -c [ 0 I] -c On exit wn is unchanged. -c -c iprint is an INTEGER variable that must be set by the user. -c It controls the frequency and type of output generated: -c iprint<0 no output is generated; -c iprint=0 print only one line at the last iteration; -c 0100 print details of every iteration including x and g; -c When iprint > 0, the file iterate.dat will be created to -c summarize the iteration. -c -c info is an integer variable. -c On entry info is unspecified. -c On exit info = 0 for normal return, -c = nonzero for abnormal return -c when the matrix K is ill-conditioned. -c -c Subprograms called: -c -c Linpack dtrsl. -c -c -c References: -c -c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited -c memory algorithm for bound constrained optimization'', -c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. -c -c -c -c * * * -c -c NEOS, November 1994. (Latest revision June 1996.) -c Optimization Technology Center. -c Argonne National Laboratory and Northwestern University. -c Written by -c Ciyou Zhu -c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. -c -c -c ************ - - integer pointr,m2,col2,ibd,jy,js,i,j,k - double precision alpha,dk,temp1,temp2 - double precision one,zero - parameter (one=1.0d0,zero=0.0d0) - - if (nsub .le. 0) return - if (iprint .ge. 99) write (6,1001) - -c Compute wv = W'Zd. - - pointr = head - do 20 i = 1, col - temp1 = zero - temp2 = zero - do 10 j = 1, nsub - k = ind(j) - temp1 = temp1 + wy(k,pointr)*d(j) - temp2 = temp2 + ws(k,pointr)*d(j) - 10 continue - wv(i) = temp1 - wv(col + i) = theta*temp2 - pointr = mod(pointr,m) + 1 - 20 continue - -c Compute wv:=K^(-1)wv. - - m2 = 2*m - col2 = 2*col - call dtrsl(wn,m2,col2,wv,11,info) - if (info .ne. 0) return - do 25 i = 1, col - wv(i) = -wv(i) - 25 continue - call dtrsl(wn,m2,col2,wv,01,info) - if (info .ne. 0) return - -c Compute d = (1/theta)d + (1/theta**2)Z'W wv. - - pointr = head - do 40 jy = 1, col - js = col + jy - do 30 i = 1, nsub - k = ind(i) - d(i) = d(i) + wy(k,pointr)*wv(jy)/theta - + + ws(k,pointr)*wv(js) - 30 continue - pointr = mod(pointr,m) + 1 - 40 continue - do 50 i = 1, nsub - d(i) = d(i)/theta - 50 continue - -c Backtrack to the feasible region. - - alpha = one - temp1 = alpha - do 60 i = 1, nsub - k = ind(i) - dk = d(i) - if (nbd(k) .ne. 0) then - if (dk .lt. zero .and. nbd(k) .le. 2) then - temp2 = l(k) - x(k) - if (temp2 .ge. zero) then - temp1 = zero - else if (dk*alpha .lt. temp2) then - temp1 = temp2/dk - endif - else if (dk .gt. zero .and. nbd(k) .ge. 2) then - temp2 = u(k) - x(k) - if (temp2 .le. zero) then - temp1 = zero - else if (dk*alpha .gt. temp2) then - temp1 = temp2/dk - endif - endif - if (temp1 .lt. alpha) then - alpha = temp1 - ibd = i - endif - endif - 60 continue - - if (alpha .lt. one) then - dk = d(ibd) - k = ind(ibd) - if (dk .gt. zero) then - x(k) = u(k) - d(ibd) = zero - else if (dk .lt. zero) then - x(k) = l(k) - d(ibd) = zero - endif - endif - do 70 i = 1, nsub - k = ind(i) - x(k) = x(k) + alpha*d(i) - 70 continue - - if (iprint .ge. 99) then - if (alpha .lt. one) then - write (6,1002) alpha - else - write (6,*) 'SM solution inside the box' - end if - if (iprint .gt.100) write (6,1003) (x(i),i=1,n) - endif - - if (alpha .lt. one) then - iword = 1 - else - iword = 0 - endif - if (iprint .ge. 99) write (6,1004) - - 1001 format (/,'----------------SUBSM entered-----------------',/) - 1002 format ( 'ALPHA = ',f7.5,' backtrack to the BOX') - 1003 format ('Subspace solution X = ',/,(4x,1p,6(1x,d11.4))) - 1004 format (/,'----------------exit SUBSM --------------------',/) - - return - - end - -c====================== The end of subsm =============================== - - subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax, - + task,isave,dsave) - character*(*) task - integer isave(2) - double precision f,g,stp,ftol,gtol,xtol,stpmin,stpmax - double precision dsave(13) -c ********** -c -c Subroutine dcsrch -c -c This subroutine finds a step that satisfies a sufficient -c decrease condition and a curvature condition. -c -c Each call of the subroutine updates an interval with -c endpoints stx and sty. The interval is initially chosen -c so that it contains a minimizer of the modified function -c -c psi(stp) = f(stp) - f(0) - ftol*stp*f'(0). -c -c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the -c interval is chosen so that it contains a minimizer of f. -c -c The algorithm is designed to find a step that satisfies -c the sufficient decrease condition -c -c f(stp) <= f(0) + ftol*stp*f'(0), -c -c and the curvature condition -c -c abs(f'(stp)) <= gtol*abs(f'(0)). -c -c If ftol is less than gtol and if, for example, the function -c is bounded below, then there is always a step which satisfies -c both conditions. -c -c If no step can be found that satisfies both conditions, then -c the algorithm stops with a warning. In this case stp only -c satisfies the sufficient decrease condition. -c -c A typical invocation of dcsrch has the following outline: -c -c task = 'START' -c 10 continue -c call dcsrch( ... ) -c if (task .eq. 'FG') then -c Evaluate the function and the gradient at stp -c goto 10 -c end if -c -c NOTE: The user must no alter work arrays between calls. -c -c The subroutine statement is -c -c subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax, -c task,isave,dsave) -c where -c -c f is a double precision variable. -c On initial entry f is the value of the function at 0. -c On subsequent entries f is the value of the -c function at stp. -c On exit f is the value of the function at stp. -c -c g is a double precision variable. -c On initial entry g is the derivative of the function at 0. -c On subsequent entries g is the derivative of the -c function at stp. -c On exit g is the derivative of the function at stp. -c -c stp is a double precision variable. -c On entry stp is the current estimate of a satisfactory -c step. On initial entry, a positive initial estimate -c must be provided. -c On exit stp is the current estimate of a satisfactory step -c if task = 'FG'. If task = 'CONV' then stp satisfies -c the sufficient decrease and curvature condition. -c -c ftol is a double precision variable. -c On entry ftol specifies a nonnegative tolerance for the -c sufficient decrease condition. -c On exit ftol is unchanged. -c -c gtol is a double precision variable. -c On entry gtol specifies a nonnegative tolerance for the -c curvature condition. -c On exit gtol is unchanged. -c -c xtol is a double precision variable. -c On entry xtol specifies a nonnegative relative tolerance -c for an acceptable step. The subroutine exits with a -c warning if the relative difference between sty and stx -c is less than xtol. -c On exit xtol is unchanged. -c -c stpmin is a double precision variable. -c On entry stpmin is a nonnegative lower bound for the step. -c On exit stpmin is unchanged. -c -c stpmax is a double precision variable. -c On entry stpmax is a nonnegative upper bound for the step. -c On exit stpmax is unchanged. -c -c task is a character variable of length at least 60. -c On initial entry task must be set to 'START'. -c On exit task indicates the required action: -c -c If task(1:2) = 'FG' then evaluate the function and -c derivative at stp and call dcsrch again. -c -c If task(1:4) = 'CONV' then the search is successful. -c -c If task(1:4) = 'WARN' then the subroutine is not able -c to satisfy the convergence conditions. The exit value of -c stp contains the best point found during the search. -c -c If task(1:5) = 'ERROR' then there is an error in the -c input arguments. -c -c On exit with convergence, a warning or an error, the -c variable task contains additional information. -c -c isave is an integer work array of dimension 2. -c -c dsave is a double precision work array of dimension 13. -c -c Subprograms called -c -c MINPACK-2 ... dcstep -c -c MINPACK-1 Project. June 1983. -c Argonne National Laboratory. -c Jorge J. More' and David J. Thuente. -c -c MINPACK-2 Project. October 1993. -c Argonne National Laboratory and University of Minnesota. -c Brett M. Averick, Richard G. Carter, and Jorge J. More'. -c -c ********** - double precision zero,p5,p66 - parameter(zero=0.0d0,p5=0.5d0,p66=0.66d0) - double precision xtrapl,xtrapu - parameter(xtrapl=1.1d0,xtrapu=4.0d0) - - logical brackt - integer stage - double precision finit,ftest,fm,fx,fxm,fy,fym,ginit,gtest, - + gm,gx,gxm,gy,gym,stx,sty,stmin,stmax,width,width1 - -c Initialization block. - - if (task(1:5) .eq. 'START') then - -c Check the input arguments for errors. - - if (stp .lt. stpmin) task = 'ERROR: STP .LT. STPMIN' - if (stp .gt. stpmax) task = 'ERROR: STP .GT. STPMAX' - if (g .ge. zero) task = 'ERROR: INITIAL G .GE. ZERO' - if (ftol .lt. zero) task = 'ERROR: FTOL .LT. ZERO' - if (gtol .lt. zero) task = 'ERROR: GTOL .LT. ZERO' - if (xtol .lt. zero) task = 'ERROR: XTOL .LT. ZERO' - if (stpmin .lt. zero) task = 'ERROR: STPMIN .LT. ZERO' - if (stpmax .lt. stpmin) task = 'ERROR: STPMAX .LT. STPMIN' - -c Exit if there are errors on input. - - if (task(1:5) .eq. 'ERROR') return - -c Initialize local variables. - - brackt = .false. - stage = 1 - finit = f - ginit = g - gtest = ftol*ginit - width = stpmax - stpmin - width1 = width/p5 - -c The variables stx, fx, gx contain the values of the step, -c function, and derivative at the best step. -c The variables sty, fy, gy contain the value of the step, -c function, and derivative at sty. -c The variables stp, f, g contain the values of the step, -c function, and derivative at stp. - - stx = zero - fx = finit - gx = ginit - sty = zero - fy = finit - gy = ginit - stmin = zero - stmax = stp + xtrapu*stp - task = 'FG' - - goto 1000 - - else - -c Restore local variables. - - if (isave(1) .eq. 1) then - brackt = .true. - else - brackt = .false. - endif - stage = isave(2) - ginit = dsave(1) - gtest = dsave(2) - gx = dsave(3) - gy = dsave(4) - finit = dsave(5) - fx = dsave(6) - fy = dsave(7) - stx = dsave(8) - sty = dsave(9) - stmin = dsave(10) - stmax = dsave(11) - width = dsave(12) - width1 = dsave(13) - - endif - -c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the -c algorithm enters the second stage. - - ftest = finit + stp*gtest - if (stage .eq. 1 .and. f .le. ftest .and. g .ge. zero) - + stage = 2 - -c Test for warnings. - - if (brackt .and. (stp .le. stmin .or. stp .ge. stmax)) - + task = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS' - if (brackt .and. stmax - stmin .le. xtol*stmax) - + task = 'WARNING: XTOL TEST SATISFIED' - if (stp .eq. stpmax .and. f .le. ftest .and. g .le. gtest) - + task = 'WARNING: STP = STPMAX' - if (stp .eq. stpmin .and. (f .gt. ftest .or. g .ge. gtest)) - + task = 'WARNING: STP = STPMIN' - -c Test for convergence. - - if (f .le. ftest .and. abs(g) .le. gtol*(-ginit)) - + task = 'CONVERGENCE' - -c Test for termination. - - if (task(1:4) .eq. 'WARN' .or. task(1:4) .eq. 'CONV') goto 1000 - -c A modified function is used to predict the step during the -c first stage if a lower function value has been obtained but -c the decrease is not sufficient. - - if (stage .eq. 1 .and. f .le. fx .and. f .gt. ftest) then - -c Define the modified function and derivative values. - - fm = f - stp*gtest - fxm = fx - stx*gtest - fym = fy - sty*gtest - gm = g - gtest - gxm = gx - gtest - gym = gy - gtest - -c Call dcstep to update stx, sty, and to compute the new step. - - call dcstep(stx,fxm,gxm,sty,fym,gym,stp,fm,gm, - + brackt,stmin,stmax) - -c Reset the function and derivative values for f. - - fx = fxm + stx*gtest - fy = fym + sty*gtest - gx = gxm + gtest - gy = gym + gtest - - else - -c Call dcstep to update stx, sty, and to compute the new step. - - call dcstep(stx,fx,gx,sty,fy,gy,stp,f,g, - + brackt,stmin,stmax) - - endif - -c Decide if a bisection step is needed. - - if (brackt) then - if (abs(sty-stx) .ge. p66*width1) stp = stx + p5*(sty - stx) - width1 = width - width = abs(sty-stx) - endif - -c Set the minimum and maximum steps allowed for stp. - - if (brackt) then - stmin = min(stx,sty) - stmax = max(stx,sty) - else - stmin = stp + xtrapl*(stp - stx) - stmax = stp + xtrapu*(stp - stx) - endif - -c Force the step to be within the bounds stpmax and stpmin. - - stp = max(stp,stpmin) - stp = min(stp,stpmax) - -c If further progress is not possible, let stp be the best -c point obtained during the search. - - if (brackt .and. (stp .le. stmin .or. stp .ge. stmax) - + .or. (brackt .and. stmax-stmin .le. xtol*stmax)) stp = stx - -c Obtain another function and derivative. - - task = 'FG' - - 1000 continue - -c Save local variables. - - if (brackt) then - isave(1) = 1 - else - isave(1) = 0 - endif - isave(2) = stage - dsave(1) = ginit - dsave(2) = gtest - dsave(3) = gx - dsave(4) = gy - dsave(5) = finit - dsave(6) = fx - dsave(7) = fy - dsave(8) = stx - dsave(9) = sty - dsave(10) = stmin - dsave(11) = stmax - dsave(12) = width - dsave(13) = width1 - - end - -c====================== The end of dcsrch ============================== - - subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, - + stpmin,stpmax) - logical brackt - double precision stx,fx,dx,sty,fy,dy,stp,fp,dp,stpmin,stpmax -c ********** -c -c Subroutine dcstep -c -c This subroutine computes a safeguarded step for a search -c procedure and updates an interval that contains a step that -c satisfies a sufficient decrease and a curvature condition. -c -c The parameter stx contains the step with the least function -c value. If brackt is set to .true. then a minimizer has -c been bracketed in an interval with endpoints stx and sty. -c The parameter stp contains the current step. -c The subroutine assumes that if brackt is set to .true. then -c -c min(stx,sty) < stp < max(stx,sty), -c -c and that the derivative at stx is negative in the direction -c of the step. -c -c The subroutine statement is -c -c subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, -c stpmin,stpmax) -c -c where -c -c stx is a double precision variable. -c On entry stx is the best step obtained so far and is an -c endpoint of the interval that contains the minimizer. -c On exit stx is the updated best step. -c -c fx is a double precision variable. -c On entry fx is the function at stx. -c On exit fx is the function at stx. -c -c dx is a double precision variable. -c On entry dx is the derivative of the function at -c stx. The derivative must be negative in the direction of -c the step, that is, dx and stp - stx must have opposite -c signs. -c On exit dx is the derivative of the function at stx. -c -c sty is a double precision variable. -c On entry sty is the second endpoint of the interval that -c contains the minimizer. -c On exit sty is the updated endpoint of the interval that -c contains the minimizer. -c -c fy is a double precision variable. -c On entry fy is the function at sty. -c On exit fy is the function at sty. -c -c dy is a double precision variable. -c On entry dy is the derivative of the function at sty. -c On exit dy is the derivative of the function at the exit sty. -c -c stp is a double precision variable. -c On entry stp is the current step. If brackt is set to .true. -c then on input stp must be between stx and sty. -c On exit stp is a new trial step. -c -c fp is a double precision variable. -c On entry fp is the function at stp -c On exit fp is unchanged. -c -c dp is a double precision variable. -c On entry dp is the the derivative of the function at stp. -c On exit dp is unchanged. -c -c brackt is an logical variable. -c On entry brackt specifies if a minimizer has been bracketed. -c Initially brackt must be set to .false. -c On exit brackt specifies if a minimizer has been bracketed. -c When a minimizer is bracketed brackt is set to .true. -c -c stpmin is a double precision variable. -c On entry stpmin is a lower bound for the step. -c On exit stpmin is unchanged. -c -c stpmax is a double precision variable. -c On entry stpmax is an upper bound for the step. -c On exit stpmax is unchanged. -c -c MINPACK-1 Project. June 1983 -c Argonne National Laboratory. -c Jorge J. More' and David J. Thuente. -c -c MINPACK-2 Project. October 1993. -c Argonne National Laboratory and University of Minnesota. -c Brett M. Averick and Jorge J. More'. -c -c ********** - double precision zero,p66,two,three - parameter(zero=0.0d0,p66=0.66d0,two=2.0d0,three=3.0d0) - - double precision gamma,p,q,r,s,sgnd,stpc,stpf,stpq,theta - - sgnd = dp*(dx/abs(dx)) - -c First case: A higher function value. The minimum is bracketed. -c If the cubic step is closer to stx than the quadratic step, the -c cubic step is taken, otherwise the average of the cubic and -c quadratic steps is taken. - - if (fp .gt. fx) then - theta = three*(fx - fp)/(stp - stx) + dx + dp - s = max(abs(theta),abs(dx),abs(dp)) - gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp/s)) - if (stp .lt. stx) gamma = -gamma - p = (gamma - dx) + theta - q = ((gamma - dx) + gamma) + dp - r = p/q - stpc = stx + r*(stp - stx) - stpq = stx + ((dx/((fx - fp)/(stp - stx) + dx))/two)* - + (stp - stx) - if (abs(stpc-stx) .lt. abs(stpq-stx)) then - stpf = stpc - else - stpf = stpc + (stpq - stpc)/two - endif - brackt = .true. - -c Second case: A lower function value and derivatives of opposite -c sign. The minimum is bracketed. If the cubic step is farther from -c stp than the secant step, the cubic step is taken, otherwise the -c secant step is taken. - - else if (sgnd .lt. zero) then - theta = three*(fx - fp)/(stp - stx) + dx + dp - s = max(abs(theta),abs(dx),abs(dp)) - gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp/s)) - if (stp .gt. stx) gamma = -gamma - p = (gamma - dp) + theta - q = ((gamma - dp) + gamma) + dx - r = p/q - stpc = stp + r*(stx - stp) - stpq = stp + (dp/(dp - dx))*(stx - stp) - if (abs(stpc-stp) .gt. abs(stpq-stp)) then - stpf = stpc - else - stpf = stpq - endif - brackt = .true. - -c Third case: A lower function value, derivatives of the same sign, -c and the magnitude of the derivative decreases. - - else if (abs(dp) .lt. abs(dx)) then - -c The cubic step is computed only if the cubic tends to infinity -c in the direction of the step or if the minimum of the cubic -c is beyond stp. Otherwise the cubic step is defined to be the -c secant step. - - theta = three*(fx - fp)/(stp - stx) + dx + dp - s = max(abs(theta),abs(dx),abs(dp)) - -c The case gamma = 0 only arises if the cubic does not tend -c to infinity in the direction of the step. - - gamma = s*sqrt(max(zero,(theta/s)**2-(dx/s)*(dp/s))) - if (stp .gt. stx) gamma = -gamma - p = (gamma - dp) + theta - q = (gamma + (dx - dp)) + gamma - r = p/q - if (r .lt. zero .and. gamma .ne. zero) then - stpc = stp + r*(stx - stp) - else if (stp .gt. stx) then - stpc = stpmax - else - stpc = stpmin - endif - stpq = stp + (dp/(dp - dx))*(stx - stp) - - if (brackt) then - -c A minimizer has been bracketed. If the cubic step is -c closer to stp than the secant step, the cubic step is -c taken, otherwise the secant step is taken. - - if (abs(stpc-stp) .lt. abs(stpq-stp)) then - stpf = stpc - else - stpf = stpq - endif - if (stp .gt. stx) then - stpf = min(stp+p66*(sty-stp),stpf) - else - stpf = max(stp+p66*(sty-stp),stpf) - endif - else - -c A minimizer has not been bracketed. If the cubic step is -c farther from stp than the secant step, the cubic step is -c taken, otherwise the secant step is taken. - - if (abs(stpc-stp) .gt. abs(stpq-stp)) then - stpf = stpc - else - stpf = stpq - endif - stpf = min(stpmax,stpf) - stpf = max(stpmin,stpf) - endif - -c Fourth case: A lower function value, derivatives of the same sign, -c and the magnitude of the derivative does not decrease. If the -c minimum is not bracketed, the step is either stpmin or stpmax, -c otherwise the cubic step is taken. - - else - if (brackt) then - theta = three*(fp - fy)/(sty - stp) + dy + dp - s = max(abs(theta),abs(dy),abs(dp)) - gamma = s*sqrt((theta/s)**2 - (dy/s)*(dp/s)) - if (stp .gt. sty) gamma = -gamma - p = (gamma - dp) + theta - q = ((gamma - dp) + gamma) + dy - r = p/q - stpc = stp + r*(sty - stp) - stpf = stpc - else if (stp .gt. stx) then - stpf = stpmax - else - stpf = stpmin - endif - endif - -c Update the interval which contains a minimizer. - - if (fp .gt. fx) then - sty = stp - fy = fp - dy = dp - else - if (sgnd .lt. zero) then - sty = stx - fy = fx - dy = dx - endif - stx = stp - fx = fp - dx = dp - endif - -c Compute the new step. - - stp = stpf - - end - -c====================== The end of dcstep ============================== - - subroutine timer(ttime) - double precision ttime -c ********* -c -c Subroutine timer -c -c This subroutine is used to determine user time. In a typical -c application, the user time for a code segment requires calls -c to subroutine timer to determine the initial and final time. -c -c The subroutine statement is -c -c subroutine timer(ttime) -c -c where -c -c ttime is an output variable which specifies the user time. -c -c Argonne National Laboratory and University of Minnesota. -c MINPACK-2 Project. -c -c Modified October 1990 by Brett M. Averick. -c -c ********** - real temp - real tarray(2) - real etime - -c The first element of the array tarray specifies user time - - temp = etime(tarray) - - ttime = dble(tarray(1)) - - return - - end - -c====================== The end of timer =============================== - - double precision function dnrm2(n,x,incx) - integer n,incx - double precision x(n) -c ********** -c -c Function dnrm2 -c -c Given a vector x of length n, this function calculates the -c Euclidean norm of x with stride incx. -c -c The function statement is -c -c double precision function dnrm2(n,x,incx) -c -c where -c -c n is a positive integer input variable. -c -c x is an input array of length n. -c -c incx is a positive integer variable that specifies the -c stride of the vector. -c -c Subprograms called -c -c FORTRAN-supplied ... abs, max, sqrt -c -c MINPACK-2 Project. February 1991. -c Argonne National Laboratory. -c Brett M. Averick. -c -c ********** - integer i - double precision scale - - dnrm2 = 0.0d0 - scale = 0.0d0 - - do 10 i = 1, n, incx - scale = max(scale, abs(x(i))) - 10 continue - - if (scale .eq. 0.0d0) return - - do 20 i = 1, n, incx - dnrm2 = dnrm2 + (x(i)/scale)**2 - 20 continue - - dnrm2 = scale*sqrt(dnrm2) - - - return - - end - -c====================== The end of dnrm2 =============================== - - double precision function dpmeps() -c ********** -c -c Subroutine dpeps -c -c This subroutine computes the machine precision parameter -c dpmeps as the smallest floating point number such that -c 1 + dpmeps differs from 1. -c -c This subroutine is based on the subroutine machar described in -c -c W. J. Cody, -c MACHAR: A subroutine to dynamically determine machine parameters, -c ACM Transactions on Mathematical Software, 14, 1988, pages 303-311. -c -c The subroutine statement is: -c -c subroutine dpeps(dpmeps) -c -c where -c -c dpmeps is a double precision variable. -c On entry dpmeps need not be specified. -c On exit dpmeps is the machine precision. -c -c MINPACK-2 Project. February 1991. -c Argonne National Laboratory and University of Minnesota. -c Brett M. Averick. -c -c ******* - integer i,ibeta,irnd,it,itemp,negep - double precision a,b,beta,betain,betah,temp,tempa,temp1, - + zero,one,two - data zero,one,two /0.0d0,1.0d0,2.0d0/ - -c determine ibeta, beta ala malcolm. - - a = one - b = one - 10 continue - a = a + a - temp = a + one - temp1 = temp - a - if (temp1 - one .eq. zero) go to 10 - 20 continue - b = b + b - temp = a + b - itemp = int(temp - a) - if (itemp .eq. 0) go to 20 - ibeta = itemp - beta = dble(ibeta) - -c determine it, irnd. - - it = 0 - b = one - 30 continue - it = it + 1 - b = b * beta - temp = b + one - temp1 = temp - b - if (temp1 - one .eq. zero) go to 30 - irnd = 0 - betah = beta/two - temp = a + betah - if (temp - a .ne. zero) irnd = 1 - tempa = a + beta - temp = tempa + betah - if ((irnd .eq. 0) .and. (temp - tempa .ne. zero)) irnd = 2 - -c determine dpmeps. - - negep = it + 3 - betain = one/beta - a = one - do 40 i = 1, negep - a = a*betain - 40 continue - 50 continue - temp = one + a - if (temp - one .ne. zero) go to 60 - a = a*beta - go to 50 - 60 continue - dpmeps = a - if ((ibeta .eq. 2) .or. (irnd .eq. 0)) go to 70 - a = (a*(one + a))/two - temp = one + a - if (temp - one .ne. zero) dpmeps = a - - 70 return - - end - -c====================== The end of dpmeps ============================== - - subroutine daxpy(n,da,dx,incx,dy,incy) -c -c constant times a vector plus a vector. -c uses unrolled loops for increments equal to one. -c jack dongarra, linpack, 3/11/78. -c - double precision dx(1),dy(1),da - integer i,incx,incy,ix,iy,m,mp1,n -c - if(n.le.0)return - if (da .eq. 0.0d0) return - if(incx.eq.1.and.incy.eq.1)go to 20 -c -c code for unequal increments or equal increments -c not equal to 1 -c - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - dy(iy) = dy(iy) + da*dx(ix) - ix = ix + incx - iy = iy + incy - 10 continue - return -c -c code for both increments equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,4) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dy(i) = dy(i) + da*dx(i) - 30 continue - if( n .lt. 4 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,4 - dy(i) = dy(i) + da*dx(i) - dy(i + 1) = dy(i + 1) + da*dx(i + 1) - dy(i + 2) = dy(i + 2) + da*dx(i + 2) - dy(i + 3) = dy(i + 3) + da*dx(i + 3) - 50 continue - return - end - -c====================== The end of daxpy =============================== - - subroutine dcopy(n,dx,incx,dy,incy) -c -c copies a vector, x, to a vector, y. -c uses unrolled loops for increments equal to one. -c jack dongarra, linpack, 3/11/78. -c - double precision dx(1),dy(1) - integer i,incx,incy,ix,iy,m,mp1,n -c - if(n.le.0)return - if(incx.eq.1.and.incy.eq.1)go to 20 -c -c code for unequal increments or equal increments -c not equal to 1 -c - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - dy(iy) = dx(ix) - ix = ix + incx - iy = iy + incy - 10 continue - return -c -c code for both increments equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,7) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dy(i) = dx(i) - 30 continue - if( n .lt. 7 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,7 - dy(i) = dx(i) - dy(i + 1) = dx(i + 1) - dy(i + 2) = dx(i + 2) - dy(i + 3) = dx(i + 3) - dy(i + 4) = dx(i + 4) - dy(i + 5) = dx(i + 5) - dy(i + 6) = dx(i + 6) - 50 continue - return - end - -c====================== The end of dcopy =============================== - - double precision function ddot(n,dx,incx,dy,incy) -c -c forms the dot product of two vectors. -c uses unrolled loops for increments equal to one. -c jack dongarra, linpack, 3/11/78. -c - double precision dx(1),dy(1),dtemp - integer i,incx,incy,ix,iy,m,mp1,n -c - ddot = 0.0d0 - dtemp = 0.0d0 - if(n.le.0)return - if(incx.eq.1.and.incy.eq.1)go to 20 -c -c code for unequal increments or equal increments -c not equal to 1 -c - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - dtemp = dtemp + dx(ix)*dy(iy) - ix = ix + incx - iy = iy + incy - 10 continue - ddot = dtemp - return -c -c code for both increments equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,5) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dtemp = dtemp + dx(i)*dy(i) - 30 continue - if( n .lt. 5 ) go to 60 - 40 mp1 = m + 1 - do 50 i = mp1,n,5 - dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + - * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) - 50 continue - 60 ddot = dtemp - return - end - -c====================== The end of ddot ================================ - - subroutine dpofa(a,lda,n,info) - integer lda,n,info - double precision a(lda,1) -c -c dpofa factors a double precision symmetric positive definite -c matrix. -c -c dpofa is usually called by dpoco, but it can be called -c directly with a saving in time if rcond is not needed. -c (time for dpoco) = (1 + 18/n)*(time for dpofa) . -c -c on entry -c -c a double precision(lda, n) -c the symmetric matrix to be factored. only the -c diagonal and upper triangle are used. -c -c lda integer -c the leading dimension of the array a . -c -c n integer -c the order of the matrix a . -c -c on return -c -c a an upper triangular matrix r so that a = trans(r)*r -c where trans(r) is the transpose. -c the strict lower triangle is unaltered. -c if info .ne. 0 , the factorization is not complete. -c -c info integer -c = 0 for normal return. -c = k signals an error condition. the leading minor -c of order k is not positive definite. -c -c linpack. this version dated 08/14/78 . -c cleve moler, university of new mexico, argonne national lab. -c -c subroutines and functions -c -c blas ddot -c fortran sqrt -c -c internal variables -c - double precision ddot,t - double precision s - integer j,jm1,k -c begin block with ...exits to 40 -c -c - do 30 j = 1, n - info = j - s = 0.0d0 - jm1 = j - 1 - if (jm1 .lt. 1) go to 20 - do 10 k = 1, jm1 - t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1) - t = t/a(k,k) - a(k,j) = t - s = s + t*t - 10 continue - 20 continue - s = a(j,j) - s -c ......exit - if (s .le. 0.0d0) go to 40 - a(j,j) = sqrt(s) - 30 continue - info = 0 - 40 continue - return - end - -c====================== The end of dpofa =============================== - - subroutine dscal(n,da,dx,incx) -c -c scales a vector by a constant. -c uses unrolled loops for increment equal to one. -c jack dongarra, linpack, 3/11/78. -c modified 3/93 to return if incx .le. 0. -c - double precision da,dx(1) - integer i,incx,m,mp1,n,nincx -c - if( n.le.0 .or. incx.le.0 )return - if(incx.eq.1)go to 20 -c -c code for increment not equal to 1 -c - nincx = n*incx - do 10 i = 1,nincx,incx - dx(i) = da*dx(i) - 10 continue - return -c -c code for increment equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,5) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dx(i) = da*dx(i) - 30 continue - if( n .lt. 5 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,5 - dx(i) = da*dx(i) - dx(i + 1) = da*dx(i + 1) - dx(i + 2) = da*dx(i + 2) - dx(i + 3) = da*dx(i + 3) - dx(i + 4) = da*dx(i + 4) - 50 continue - return - end - -c====================== The end of dscal =============================== - - subroutine dtrsl(t,ldt,n,b,job,info) - integer ldt,n,job,info - double precision t(ldt,1),b(1) -c -c -c dtrsl solves systems of the form -c -c t * x = b -c or -c trans(t) * x = b -c -c where t is a triangular matrix of order n. here trans(t) -c denotes the transpose of the matrix t. -c -c on entry -c -c t double precision(ldt,n) -c t contains the matrix of the system. the zero -c elements of the matrix are not referenced, and -c the corresponding elements of the array can be -c used to store other information. -c -c ldt integer -c ldt is the leading dimension of the array t. -c -c n integer -c n is the order of the system. -c -c b double precision(n). -c b contains the right hand side of the system. -c -c job integer -c job specifies what kind of system is to be solved. -c if job is -c -c 00 solve t*x=b, t lower triangular, -c 01 solve t*x=b, t upper triangular, -c 10 solve trans(t)*x=b, t lower triangular, -c 11 solve trans(t)*x=b, t upper triangular. -c -c on return -c -c b b contains the solution, if info .eq. 0. -c otherwise b is unaltered. -c -c info integer -c info contains zero if the system is nonsingular. -c otherwise info contains the index of -c the first zero diagonal element of t. -c -c linpack. this version dated 08/14/78 . -c g. w. stewart, university of maryland, argonne national lab. -c -c subroutines and functions -c -c blas daxpy,ddot -c fortran mod -c -c internal variables -c - double precision ddot,temp - integer case,j,jj -c -c begin block permitting ...exits to 150 -c -c check for zero diagonal elements. -c - do 10 info = 1, n -c ......exit - if (t(info,info) .eq. 0.0d0) go to 150 - 10 continue - info = 0 -c -c determine the task and go to it. -c - case = 1 - if (mod(job,10) .ne. 0) case = 2 - if (mod(job,100)/10 .ne. 0) case = case + 2 - go to (20,50,80,110), case -c -c solve t*x=b for t lower triangular -c - 20 continue - b(1) = b(1)/t(1,1) - if (n .lt. 2) go to 40 - do 30 j = 2, n - temp = -b(j-1) - call daxpy(n-j+1,temp,t(j,j-1),1,b(j),1) - b(j) = b(j)/t(j,j) - 30 continue - 40 continue - go to 140 -c -c solve t*x=b for t upper triangular. -c - 50 continue - b(n) = b(n)/t(n,n) - if (n .lt. 2) go to 70 - do 60 jj = 2, n - j = n - jj + 1 - temp = -b(j+1) - call daxpy(j,temp,t(1,j+1),1,b(1),1) - b(j) = b(j)/t(j,j) - 60 continue - 70 continue - go to 140 -c -c solve trans(t)*x=b for t lower triangular. -c - 80 continue - b(n) = b(n)/t(n,n) - if (n .lt. 2) go to 100 - do 90 jj = 2, n - j = n - jj + 1 - b(j) = b(j) - ddot(jj-1,t(j+1,j),1,b(j+1),1) - b(j) = b(j)/t(j,j) - 90 continue - 100 continue - go to 140 -c -c solve trans(t)*x=b for t upper triangular. -c - 110 continue - b(1) = b(1)/t(1,1) - if (n .lt. 2) go to 130 - do 120 j = 2, n - b(j) = b(j) - ddot(j-1,t(1,j),1,b(1),1) - b(j) = b(j)/t(j,j) - 120 continue - 130 continue - 140 continue - 150 continue - return - end - -c====================== The end of dtrsl =============================== - EOF routines.f