diff options
author | Gard Spreemann <gspreemann@gmail.com> | 2017-02-02 21:23:01 +0100 |
---|---|---|
committer | Gard Spreemann <gspreemann@gmail.com> | 2017-02-02 21:23:01 +0100 |
commit | 173802e9f8d98a01019d3f5aff055f1f04479974 (patch) | |
tree | c47b4169dee35def4cdeaaad4cdebaaae5c722da |
upstream/3.0+dfsg.3: Start git repository for package by importingdfsg/3.0+dfsg.3
upstream's 3.0 tarball DFSG-cleaned as currently used in published
Debian package.
-rw-r--r-- | License.txt | 71 | ||||
-rw-r--r-- | Makefile | 37 | ||||
-rw-r--r-- | OUTPUTS/output_77_1 | 85 | ||||
-rw-r--r-- | OUTPUTS/output_77_2 | 57 | ||||
-rw-r--r-- | OUTPUTS/output_77_3 | 222 | ||||
-rw-r--r-- | OUTPUTS/output_90_1 | 85 | ||||
-rw-r--r-- | OUTPUTS/output_90_2 | 57 | ||||
-rw-r--r-- | OUTPUTS/output_90_3 | 222 | ||||
-rw-r--r-- | README | 232 | ||||
-rwxr-xr-x | algorithm.pdf | bin | 0 -> 260255 bytes | |||
-rwxr-xr-x | code.pdf | bin | 0 -> 180639 bytes | |||
-rw-r--r-- | driver1.f | 321 | ||||
-rwxr-xr-x | driver1.f90 | 297 | ||||
-rw-r--r-- | driver2.f | 232 | ||||
-rwxr-xr-x | driver2.f90 | 220 | ||||
-rw-r--r-- | driver3.f | 273 | ||||
-rw-r--r-- | driver3.f90 | 256 | ||||
-rw-r--r-- | iterate.dat | 49 | ||||
-rw-r--r-- | lbfgsb.f | 3949 | ||||
-rw-r--r-- | timer.f | 32 |
20 files changed, 6697 insertions, 0 deletions
diff --git a/License.txt b/License.txt new file mode 100644 index 0000000..9ee303a --- /dev/null +++ b/License.txt @@ -0,0 +1,71 @@ +3-clause license ("New BSD License" or "Modified BSD License") +New BSD License +Author Regents of the University of California +Publisher Public Domain +Published July 22, 1999[8] +DFSG compatible Yes[7] +FSF approved Yes[1] +OSI approved Yes[3] +GPL compatible Yes[1] +Copyleft No[1] +Copyfree Yes +Linking from code with a different license Yes + +The advertising clause was removed from the license text in the official BSD on July 22, 1999 by William Hoskins, Director of +the Office of Technology Licensing for UC Berkeley.[8] Other BSD distributions removed the clause, but many similar clauses +remain in BSD-derived code from other sources, and unrelated code using a derived license. + +While the original license is sometimes referred to as "BSD-old", the resulting 3-clause version is sometimes referred to by +"BSD-new." Other names include "New BSD", "revised BSD", "BSD-3", or "3-clause BSD". This version has been vetted as +an Open source license by the OSI as the "The BSD License".[3] The Free Software Foundation, which refers to the license +as the "Modified BSD License", states that it is compatible with the GNU GPL. The FSF encourages users to be specific +when referring to the license by name (i.e. not simply referring to it as "a BSD license" or "BSD-style") to avoid confusion with +the original BSD license.[1] + +This version allows unlimited redistribution for any purpose as long as its copyright notices and the license's disclaimers of +warranty are maintained. The license also contains a clause restricting use of the names of contributors for endorsement of a +derived work without specific permission. + +Copyright (c) <year>, <copyright holder> +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the <organization> nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +References + +1. ^ a b c d e "Various Licenses and Comments about Them - GNU Project - Free Software Foundation (FSF): Modified BSD license". + Free Software Foundation. Retrieved 02 October 2010. +2. ^ a b c d e "Various Licenses and Comments about Them - GNU Project - Free Software Foundation (FSF): FreeBSD license". + Free Software Foundation. Retrieved 02 October 2010. +3. ^ a b c d e f "Open Source Initiative OSI - The BSD License:Licensing". Open Source Initiative. Retrieved 06 December 2009. +4. ^ a b c d "Various Licenses and Comments about Them - GNU Project - Free Software Foundation (FSF): Original BSD license". + Free Software Foundation. Retrieved 02 October 2010. +5. ^ The year given is the year 4.3BSD-Tahoe was released. Whether this is the first use of the license is not known. +6. ^ The year given is the year 4.3BSD-Reno was released. Whether this is the first use of the license is not known. +7. ^ a b "Debian -- License information". Debian. Retrieved 18 February 2010. +8. ^ a b c "To All Licensees, Distributors of Any Version of BSD". University of California, Berkeley. 1999-07-22. Retrieved + 2006-11-15. +9. ^ Richard Stallman. "The BSD License Problem". Free Software Foundation. Retrieved 2006-11-15. +10. ^ "The FreeBSD Copyright". The FreeBSD Project. Retrieved 6 December 2009. +11. ^ "NetBSD Licensing and Redistribution". The NetBSD Foundation. Retrieved 06 December 2009. diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..7006149 --- /dev/null +++ b/Makefile @@ -0,0 +1,37 @@ +FC = gfortran + +FFLAGS = -O -Wall -fbounds-check -g -Wno-uninitialized + +DRIVER1_77 = driver1.f +DRIVER2_77 = driver2.f +DRIVER3_77 = driver3.f + +DRIVER1_90 = driver1.f90 +DRIVER2_90 = driver2.f90 +DRIVER3_90 = driver3.f90 + +LBFGSB = lbfgsb.f +LINPACK = linpack.f +BLAS = blas.f +TIMER = timer.f + +all : lbfgsb_77_1 lbfgsb_77_2 lbfgsb_77_3 lbfgsb_90_1 lbfgsb_90_2 lbfgsb_90_3 + + +lbfgsb_77_1 : $(DRIVER1_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) + $(FC) $(FFLAGS) $(DRIVER1_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_77_1 + +lbfgsb_77_2 : $(DRIVER2_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) + $(FC) $(FFLAGS) $(DRIVER2_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_77_2 + +lbfgsb_77_3 : $(DRIVER3_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) + $(FC) $(FFLAGS) $(DRIVER3_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_77_3 + +lbfgsb_90_1 : $(DRIVER1_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) + $(FC) $(FFLAGS) $(DRIVER1_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_90_1 + +lbfgsb_90_2 : $(DRIVER2_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) + $(FC) $(FFLAGS) $(DRIVER2_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_90_2 + +lbfgsb_90_3 : $(DRIVER3_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) + $(FC) $(FFLAGS) $(DRIVER3_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_90_3 diff --git a/OUTPUTS/output_77_1 b/OUTPUTS/output_77_1 new file mode 100644 index 0000000..52dc375 --- /dev/null +++ b/OUTPUTS/output_77_1 @@ -0,0 +1,85 @@ + + Solving sample problem. + (f = 0.0 at the optimal solution.) + +RUNNING THE L-BFGS-B CODE + + * * * + +Machine precision = 2.220D-16 + N = 25 M = 5 + +At X0 0 variables are exactly at the bounds + +At iterate 0 f= 3.46000D+03 |proj g|= 1.03000D+02 + +At iterate 1 f= 2.39769D+03 |proj g|= 6.50700D+01 + +At iterate 2 f= 1.43806D+02 |proj g|= 3.64039D+01 + +At iterate 3 f= 7.28161D+01 |proj g|= 2.29042D+01 + +At iterate 4 f= 1.60308D+01 |proj g|= 6.95409D+00 + +At iterate 5 f= 5.18725D+00 |proj g|= 9.05481D+00 + +At iterate 6 f= 2.12716D+00 |proj g|= 1.96729D+01 + +At iterate 7 f= 2.07568D-01 |proj g|= 2.12849D+00 + +At iterate 8 f= 5.32739D-02 |proj g|= 8.32469D-01 + +At iterate 9 f= 1.30450D-02 |proj g|= 4.27926D-01 + +At iterate 10 f= 3.86031D-03 |proj g|= 2.00812D-01 + +At iterate 11 f= 7.45653D-04 |proj g|= 1.37723D-01 + +At iterate 12 f= 3.54016D-04 |proj g|= 1.21274D-01 + +At iterate 13 f= 7.42511D-05 |proj g|= 2.97814D-02 + +At iterate 14 f= 3.74062D-05 |proj g|= 1.72742D-02 + +At iterate 15 f= 1.09832D-05 |proj g|= 2.86815D-02 + +At iterate 16 f= 3.90757D-06 |proj g|= 8.08085D-03 + +At iterate 17 f= 1.99502D-06 |proj g|= 3.47648D-03 + +At iterate 18 f= 8.25796D-07 |proj g|= 2.25283D-03 + +At iterate 19 f= 1.99224D-07 |proj g|= 1.45773D-03 + +At iterate 20 f= 5.75772D-08 |proj g|= 1.48245D-03 + +At iterate 21 f= 1.46326D-08 |proj g|= 5.44304D-04 + +At iterate 22 f= 2.36329D-09 |proj g|= 2.25244D-04 + +At iterate 23 f= 1.08349D-09 |proj g|= 1.72052D-04 + + * * * + +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 + + * * * + + N Tit Tnf Tnint Skip Nact Projg F + 25 23 28 47 0 0 1.721D-04 1.083D-09 + F = 1.083490083518441E-009 + +CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH + + Cauchy time 1.000E-03 seconds. + Subspace minimization time 0.000E+00 seconds. + Line search time 0.000E+00 seconds. + + Total User time 2.000E-03 seconds. + diff --git a/OUTPUTS/output_77_2 b/OUTPUTS/output_77_2 new file mode 100644 index 0000000..f12b184 --- /dev/null +++ b/OUTPUTS/output_77_2 @@ -0,0 +1,57 @@ + + Solving sample problem. + (f = 0.0 at the optimal solution.) + +Iterate 1 nfg = 5 f = 2.39769D+03 |proj g| = 6.50700D+01 +Iterate 2 nfg = 6 f = 1.43806D+02 |proj g| = 3.64039D+01 +Iterate 3 nfg = 7 f = 7.28161D+01 |proj g| = 2.29042D+01 +Iterate 4 nfg = 8 f = 1.60308D+01 |proj g| = 6.95409D+00 +Iterate 5 nfg = 9 f = 5.18725D+00 |proj g| = 9.05481D+00 +Iterate 6 nfg = 10 f = 2.12716D+00 |proj g| = 1.96729D+01 +Iterate 7 nfg = 11 f = 2.07568D-01 |proj g| = 2.12849D+00 +Iterate 8 nfg = 12 f = 5.32739D-02 |proj g| = 8.32469D-01 +Iterate 9 nfg = 13 f = 1.30450D-02 |proj g| = 4.27926D-01 +Iterate 10 nfg = 14 f = 3.86031D-03 |proj g| = 2.00812D-01 +Iterate 11 nfg = 15 f = 7.45653D-04 |proj g| = 1.37723D-01 +Iterate 12 nfg = 16 f = 3.54016D-04 |proj g| = 1.21274D-01 +Iterate 13 nfg = 17 f = 7.42511D-05 |proj g| = 2.97814D-02 +Iterate 14 nfg = 18 f = 3.74062D-05 |proj g| = 1.72742D-02 +Iterate 15 nfg = 19 f = 1.09832D-05 |proj g| = 2.86815D-02 +Iterate 16 nfg = 21 f = 3.90757D-06 |proj g| = 8.08085D-03 +Iterate 17 nfg = 22 f = 1.99502D-06 |proj g| = 3.47648D-03 +Iterate 18 nfg = 23 f = 8.25796D-07 |proj g| = 2.25283D-03 +Iterate 19 nfg = 24 f = 1.99224D-07 |proj g| = 1.45773D-03 +Iterate 20 nfg = 25 f = 5.75772D-08 |proj g| = 1.48245D-03 +Iterate 21 nfg = 26 f = 1.46326D-08 |proj g| = 5.44304D-04 +Iterate 22 nfg = 27 f = 2.36329D-09 |proj g| = 2.25244D-04 +Iterate 23 nfg = 28 f = 1.08349D-09 |proj g| = 1.72052D-04 +Iterate 24 nfg = 29 f = 3.49043D-10 |proj g| = 5.78820D-05 +Iterate 25 nfg = 30 f = 8.43998D-11 |proj g| = 2.04504D-05 +Iterate 26 nfg = 31 f = 2.29812D-11 |proj g| = 1.74854D-05 +Iterate 27 nfg = 32 f = 5.56238D-12 |proj g| = 1.29349D-05 +Iterate 28 nfg = 33 f = 6.17073D-13 |proj g| = 3.52399D-06 +Iterate 29 nfg = 34 f = 2.06920D-13 |proj g| = 9.74653D-07 +Iterate 30 nfg = 35 f = 1.07949D-13 |proj g| = 6.84427D-07 +Iterate 31 nfg = 36 f = 3.27200D-14 |proj g| = 7.46523D-07 +Iterate 32 nfg = 38 f = 2.20380D-14 |proj g| = 6.32502D-07 +Iterate 33 nfg = 39 f = 1.37189D-14 |proj g| = 2.83064D-07 +Iterate 34 nfg = 40 f = 8.06341D-15 |proj g| = 1.25600D-07 +Iterate 35 nfg = 41 f = 6.30621D-15 |proj g| = 6.63279D-08 +Iterate 36 nfg = 42 f = 6.12906D-15 |proj g| = 9.93648D-08 +Iterate 37 nfg = 43 f = 5.83506D-15 |proj g| = 1.10012D-08 +Iterate 38 nfg = 44 f = 5.82125D-15 |proj g| = 8.09738D-09 +Iterate 39 nfg = 45 f = 5.81201D-15 |proj g| = 2.04600D-08 +Iterate 40 nfg = 47 f = 5.80895D-15 |proj g| = 4.71360D-09 +Iterate 41 nfg = 48 f = 5.80786D-15 |proj g| = 1.97091D-09 +Iterate 42 nfg = 49 f = 5.80716D-15 |proj g| = 1.10044D-09 +Iterate 43 nfg = 50 f = 5.80705D-15 |proj g| = 3.49651D-10 +Iterate 44 nfg = 51 f = 5.80703D-15 |proj g| = 2.29788D-10 +Iterate 45 nfg = 52 f = 5.80702D-15 |proj g| = 1.51178D-10 +Iterate 46 nfg = 53 f = 5.80702D-15 |proj g| = 6.61970D-11 + STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL + Final X= + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0001D+00 1.0002D+00 + 1.0003D+00 1.0006D+00 1.0013D+00 1.0026D+00 1.0052D+00 1.0105D+00 + 1.0210D+00 1.0425D+00 1.0867D+00 1.1810D+00 1.3947D+00 1.9452D+00 + 3.7837D+00 diff --git a/OUTPUTS/output_77_3 b/OUTPUTS/output_77_3 new file mode 100644 index 0000000..612f3d3 --- /dev/null +++ b/OUTPUTS/output_77_3 @@ -0,0 +1,222 @@ + + Solving sample problem. + (f = 0.0 at the optimal solution.) + +Iterate 1 nfg = 5 f = 9.98226D+04 |proj g| = 4.92220D+01 +Iterate 2 nfg = 6 f = 5.97672D+03 |proj g| = 7.29734D+00 +Iterate 3 nfg = 7 f = 3.02885D+03 |proj g| = 7.40956D+00 +Iterate 4 nfg = 8 f = 7.16214D+02 |proj g| = 6.00527D+00 +Iterate 5 nfg = 9 f = 2.20946D+02 |proj g| = 4.08022D+00 +Iterate 6 nfg = 10 f = 5.61473D+01 |proj g| = 5.17435D+00 +Iterate 7 nfg = 11 f = 3.39145D+01 |proj g| = 2.09576D+01 +Iterate 8 nfg = 12 f = 1.36320D+00 |proj g| = 9.93231D+00 +Iterate 9 nfg = 13 f = 3.50771D-01 |proj g| = 5.91099D+00 +Iterate 10 nfg = 14 f = 1.52466D-02 |proj g| = 6.13600D-01 +Iterate 11 nfg = 15 f = 7.36312D-03 |proj g| = 3.39884D-01 +Iterate 12 nfg = 16 f = 1.25124D-03 |proj g| = 1.17803D-01 +Iterate 13 nfg = 17 f = 4.68586D-04 |proj g| = 7.22477D-02 +Iterate 14 nfg = 18 f = 5.50022D-05 |proj g| = 4.40414D-02 +Iterate 15 nfg = 19 f = 1.82587D-05 |proj g| = 1.49985D-02 +Iterate 16 nfg = 20 f = 6.79533D-06 |proj g| = 7.05805D-03 +Iterate 17 nfg = 21 f = 1.59803D-06 |proj g| = 3.57847D-03 +Iterate 18 nfg = 22 f = 4.67168D-07 |proj g| = 2.96329D-03 +Iterate 19 nfg = 24 f = 1.61224D-07 |proj g| = 1.61125D-03 +Iterate 20 nfg = 25 f = 5.27054D-08 |proj g| = 6.79899D-04 +Iterate 21 nfg = 26 f = 2.21706D-08 |proj g| = 3.98553D-04 +Iterate 22 nfg = 27 f = 7.85699D-09 |proj g| = 8.29301D-04 +Iterate 23 nfg = 28 f = 2.57367D-09 |proj g| = 3.95843D-04 +Iterate 24 nfg = 29 f = 9.39513D-10 |proj g| = 8.04745D-05 +Iterate 25 nfg = 30 f = 5.42208D-10 |proj g| = 5.75104D-05 +Iterate 26 nfg = 31 f = 1.13904D-10 |proj g| = 3.88455D-05 +Iterate 27 nfg = 32 f = 2.19270D-11 |proj g| = 2.48929D-05 +Iterate 28 nfg = 34 f = 5.89247D-12 |proj g| = 5.98559D-06 +Iterate 29 nfg = 35 f = 3.05886D-12 |proj g| = 3.88360D-06 +Iterate 30 nfg = 36 f = 6.64922D-13 |proj g| = 2.57655D-06 +Iterate 31 nfg = 37 f = 1.82984D-13 |proj g| = 1.54163D-06 +Iterate 32 nfg = 38 f = 5.06020D-14 |proj g| = 1.46067D-06 +Iterate 33 nfg = 39 f = 2.03549D-14 |proj g| = 2.88905D-07 +Iterate 34 nfg = 40 f = 1.27697D-14 |proj g| = 2.09331D-07 +Iterate 35 nfg = 41 f = 2.59126D-15 |proj g| = 2.75406D-07 +Iterate 36 nfg = 43 f = 1.13639D-15 |proj g| = 1.13897D-07 +Iterate 37 nfg = 44 f = 5.29416D-16 |proj g| = 6.22050D-08 +Iterate 38 nfg = 45 f = 1.24096D-16 |proj g| = 2.30404D-08 +Iterate 39 nfg = 46 f = 3.14643D-17 |proj g| = 1.56779D-08 +Iterate 40 nfg = 47 f = 1.56638D-17 |proj g| = 1.19498D-08 +Iterate 41 nfg = 48 f = 2.80290D-18 |proj g| = 4.46387D-09 +Iterate 42 nfg = 49 f = 1.32837D-18 |proj g| = 2.55881D-09 +Iterate 43 nfg = 50 f = 2.54553D-19 |proj g| = 1.17505D-09 +Iterate 44 nfg = 52 f = 8.23941D-20 |proj g| = 9.01753D-10 +Iterate 45 nfg = 53 f = 2.79611D-20 |proj g| = 8.62939D-10 +Iterate 46 nfg = 54 f = 1.08724D-20 |proj g| = 2.44455D-10 +Iterate 47 nfg = 55 f = 5.32479D-21 |proj g| = 1.66311D-10 +Iterate 48 nfg = 56 f = 1.26618D-21 |proj g| = 1.90418D-10 +Iterate 49 nfg = 58 f = 5.35057D-22 |proj g| = 9.73814D-11 + STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL + Final X= + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0001D+00 1.0002D+00 1.0003D+00 1.0007D+00 1.0014D+00 + 1.0028D+00 1.0055D+00 1.0111D+00 1.0223D+00 1.0451D+00 1.0923D+00 + 1.1930D+00 1.4233D+00 2.0257D+00 4.1035D+00 diff --git a/OUTPUTS/output_90_1 b/OUTPUTS/output_90_1 new file mode 100644 index 0000000..008c044 --- /dev/null +++ b/OUTPUTS/output_90_1 @@ -0,0 +1,85 @@ + + Solving sample problem. + (f = 0.0 at the optimal solution.) + +RUNNING THE L-BFGS-B CODE + + * * * + +Machine precision = 2.220D-16 + N = 25 M = 5 + +At X0 0 variables are exactly at the bounds + +At iterate 0 f= 3.46000D+03 |proj g|= 1.03000D+02 + +At iterate 1 f= 2.39769D+03 |proj g|= 6.50700D+01 + +At iterate 2 f= 1.43806D+02 |proj g|= 3.64039D+01 + +At iterate 3 f= 7.28161D+01 |proj g|= 2.29042D+01 + +At iterate 4 f= 1.60308D+01 |proj g|= 6.95409D+00 + +At iterate 5 f= 5.18725D+00 |proj g|= 9.05481D+00 + +At iterate 6 f= 2.12716D+00 |proj g|= 1.96729D+01 + +At iterate 7 f= 2.07568D-01 |proj g|= 2.12849D+00 + +At iterate 8 f= 5.32739D-02 |proj g|= 8.32469D-01 + +At iterate 9 f= 1.30450D-02 |proj g|= 4.27926D-01 + +At iterate 10 f= 3.86031D-03 |proj g|= 2.00812D-01 + +At iterate 11 f= 7.45653D-04 |proj g|= 1.37723D-01 + +At iterate 12 f= 3.54016D-04 |proj g|= 1.21274D-01 + +At iterate 13 f= 7.42511D-05 |proj g|= 2.97814D-02 + +At iterate 14 f= 3.74062D-05 |proj g|= 1.72742D-02 + +At iterate 15 f= 1.09832D-05 |proj g|= 2.86815D-02 + +At iterate 16 f= 3.90757D-06 |proj g|= 8.08085D-03 + +At iterate 17 f= 1.99502D-06 |proj g|= 3.47648D-03 + +At iterate 18 f= 8.25796D-07 |proj g|= 2.25283D-03 + +At iterate 19 f= 1.99224D-07 |proj g|= 1.45773D-03 + +At iterate 20 f= 5.75772D-08 |proj g|= 1.48245D-03 + +At iterate 21 f= 1.46326D-08 |proj g|= 5.44304D-04 + +At iterate 22 f= 2.36329D-09 |proj g|= 2.25244D-04 + +At iterate 23 f= 1.08349D-09 |proj g|= 1.72052D-04 + + * * * + +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 + + * * * + + N Tit Tnf Tnint Skip Nact Projg F + 25 23 28 47 0 0 1.721D-04 1.083D-09 + F = 1.083490083461424E-009 + +CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH + + Cauchy time 0.000E+00 seconds. + Subspace minimization time 0.000E+00 seconds. + Line search time 0.000E+00 seconds. + + Total User time 1.000E-03 seconds. + diff --git a/OUTPUTS/output_90_2 b/OUTPUTS/output_90_2 new file mode 100644 index 0000000..159875c --- /dev/null +++ b/OUTPUTS/output_90_2 @@ -0,0 +1,57 @@ + + Solving sample problem. + (f = 0.0 at the optimal solution.) + +Iterate 1 nfg = 5 f = 2.39769D+03 |proj g| = 6.50700D+01 +Iterate 2 nfg = 6 f = 1.43806D+02 |proj g| = 3.64039D+01 +Iterate 3 nfg = 7 f = 7.28161D+01 |proj g| = 2.29042D+01 +Iterate 4 nfg = 8 f = 1.60308D+01 |proj g| = 6.95409D+00 +Iterate 5 nfg = 9 f = 5.18725D+00 |proj g| = 9.05481D+00 +Iterate 6 nfg = 10 f = 2.12716D+00 |proj g| = 1.96729D+01 +Iterate 7 nfg = 11 f = 2.07568D-01 |proj g| = 2.12849D+00 +Iterate 8 nfg = 12 f = 5.32739D-02 |proj g| = 8.32469D-01 +Iterate 9 nfg = 13 f = 1.30450D-02 |proj g| = 4.27926D-01 +Iterate 10 nfg = 14 f = 3.86031D-03 |proj g| = 2.00812D-01 +Iterate 11 nfg = 15 f = 7.45653D-04 |proj g| = 1.37723D-01 +Iterate 12 nfg = 16 f = 3.54016D-04 |proj g| = 1.21274D-01 +Iterate 13 nfg = 17 f = 7.42511D-05 |proj g| = 2.97814D-02 +Iterate 14 nfg = 18 f = 3.74062D-05 |proj g| = 1.72742D-02 +Iterate 15 nfg = 19 f = 1.09832D-05 |proj g| = 2.86815D-02 +Iterate 16 nfg = 21 f = 3.90757D-06 |proj g| = 8.08085D-03 +Iterate 17 nfg = 22 f = 1.99502D-06 |proj g| = 3.47648D-03 +Iterate 18 nfg = 23 f = 8.25796D-07 |proj g| = 2.25283D-03 +Iterate 19 nfg = 24 f = 1.99224D-07 |proj g| = 1.45773D-03 +Iterate 20 nfg = 25 f = 5.75772D-08 |proj g| = 1.48245D-03 +Iterate 21 nfg = 26 f = 1.46326D-08 |proj g| = 5.44304D-04 +Iterate 22 nfg = 27 f = 2.36329D-09 |proj g| = 2.25244D-04 +Iterate 23 nfg = 28 f = 1.08349D-09 |proj g| = 1.72052D-04 +Iterate 24 nfg = 29 f = 3.49043D-10 |proj g| = 5.78820D-05 +Iterate 25 nfg = 30 f = 8.43998D-11 |proj g| = 2.04504D-05 +Iterate 26 nfg = 31 f = 2.29812D-11 |proj g| = 1.74854D-05 +Iterate 27 nfg = 32 f = 5.56238D-12 |proj g| = 1.29349D-05 +Iterate 28 nfg = 33 f = 6.17073D-13 |proj g| = 3.52399D-06 +Iterate 29 nfg = 34 f = 2.06920D-13 |proj g| = 9.74653D-07 +Iterate 30 nfg = 35 f = 1.07949D-13 |proj g| = 6.84427D-07 +Iterate 31 nfg = 36 f = 3.27200D-14 |proj g| = 7.46523D-07 +Iterate 32 nfg = 38 f = 2.20380D-14 |proj g| = 6.32502D-07 +Iterate 33 nfg = 39 f = 1.37189D-14 |proj g| = 2.83064D-07 +Iterate 34 nfg = 40 f = 8.06341D-15 |proj g| = 1.25600D-07 +Iterate 35 nfg = 41 f = 6.30621D-15 |proj g| = 6.63280D-08 +Iterate 36 nfg = 42 f = 6.12906D-15 |proj g| = 9.93648D-08 +Iterate 37 nfg = 43 f = 5.83506D-15 |proj g| = 1.10011D-08 +Iterate 38 nfg = 44 f = 5.82125D-15 |proj g| = 8.09738D-09 +Iterate 39 nfg = 45 f = 5.81201D-15 |proj g| = 2.04599D-08 +Iterate 40 nfg = 47 f = 5.80895D-15 |proj g| = 4.71361D-09 +Iterate 41 nfg = 48 f = 5.80786D-15 |proj g| = 1.97091D-09 +Iterate 42 nfg = 49 f = 5.80716D-15 |proj g| = 1.10044D-09 +Iterate 43 nfg = 50 f = 5.80705D-15 |proj g| = 3.49639D-10 +Iterate 44 nfg = 51 f = 5.80703D-15 |proj g| = 2.29804D-10 +Iterate 45 nfg = 52 f = 5.80702D-15 |proj g| = 1.51167D-10 +Iterate 46 nfg = 53 f = 5.80702D-15 |proj g| = 6.62041D-11 + STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL + Final X= + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0001D+00 1.0002D+00 + 1.0003D+00 1.0006D+00 1.0013D+00 1.0026D+00 1.0052D+00 1.0105D+00 + 1.0210D+00 1.0425D+00 1.0867D+00 1.1810D+00 1.3947D+00 1.9452D+00 + 3.7837D+00 diff --git a/OUTPUTS/output_90_3 b/OUTPUTS/output_90_3 new file mode 100644 index 0000000..24e13da --- /dev/null +++ b/OUTPUTS/output_90_3 @@ -0,0 +1,222 @@ + + Solving sample problem. + (f = 0.0 at the optimal solution.) + +Iterate 1 nfg = 5 f = 9.98226D+04 |proj g| = 4.92220D+01 +Iterate 2 nfg = 6 f = 5.97672D+03 |proj g| = 7.29734D+00 +Iterate 3 nfg = 7 f = 3.02885D+03 |proj g| = 7.40956D+00 +Iterate 4 nfg = 8 f = 7.16214D+02 |proj g| = 6.00527D+00 +Iterate 5 nfg = 9 f = 2.20946D+02 |proj g| = 4.08022D+00 +Iterate 6 nfg = 10 f = 5.61473D+01 |proj g| = 5.17435D+00 +Iterate 7 nfg = 11 f = 3.39145D+01 |proj g| = 2.09576D+01 +Iterate 8 nfg = 12 f = 1.36320D+00 |proj g| = 9.93231D+00 +Iterate 9 nfg = 13 f = 3.50771D-01 |proj g| = 5.91099D+00 +Iterate 10 nfg = 14 f = 1.52466D-02 |proj g| = 6.13600D-01 +Iterate 11 nfg = 15 f = 7.36312D-03 |proj g| = 3.39884D-01 +Iterate 12 nfg = 16 f = 1.25124D-03 |proj g| = 1.17803D-01 +Iterate 13 nfg = 17 f = 4.68586D-04 |proj g| = 7.22477D-02 +Iterate 14 nfg = 18 f = 5.50022D-05 |proj g| = 4.40414D-02 +Iterate 15 nfg = 19 f = 1.82587D-05 |proj g| = 1.49985D-02 +Iterate 16 nfg = 20 f = 6.79533D-06 |proj g| = 7.05805D-03 +Iterate 17 nfg = 21 f = 1.59803D-06 |proj g| = 3.57847D-03 +Iterate 18 nfg = 22 f = 4.67168D-07 |proj g| = 2.96329D-03 +Iterate 19 nfg = 24 f = 1.61224D-07 |proj g| = 1.61125D-03 +Iterate 20 nfg = 25 f = 5.27054D-08 |proj g| = 6.79899D-04 +Iterate 21 nfg = 26 f = 2.21706D-08 |proj g| = 3.98553D-04 +Iterate 22 nfg = 27 f = 7.85699D-09 |proj g| = 8.29301D-04 +Iterate 23 nfg = 28 f = 2.57367D-09 |proj g| = 3.95843D-04 +Iterate 24 nfg = 29 f = 9.39513D-10 |proj g| = 8.04745D-05 +Iterate 25 nfg = 30 f = 5.42208D-10 |proj g| = 5.75104D-05 +Iterate 26 nfg = 31 f = 1.13904D-10 |proj g| = 3.88455D-05 +Iterate 27 nfg = 32 f = 2.19270D-11 |proj g| = 2.48929D-05 +Iterate 28 nfg = 34 f = 5.89247D-12 |proj g| = 5.98559D-06 +Iterate 29 nfg = 35 f = 3.05886D-12 |proj g| = 3.88360D-06 +Iterate 30 nfg = 36 f = 6.64922D-13 |proj g| = 2.57655D-06 +Iterate 31 nfg = 37 f = 1.82984D-13 |proj g| = 1.54163D-06 +Iterate 32 nfg = 38 f = 5.06020D-14 |proj g| = 1.46067D-06 +Iterate 33 nfg = 39 f = 2.03549D-14 |proj g| = 2.88905D-07 +Iterate 34 nfg = 40 f = 1.27697D-14 |proj g| = 2.09331D-07 +Iterate 35 nfg = 41 f = 2.59125D-15 |proj g| = 2.75406D-07 +Iterate 36 nfg = 43 f = 1.13639D-15 |proj g| = 1.13897D-07 +Iterate 37 nfg = 44 f = 5.29416D-16 |proj g| = 6.22050D-08 +Iterate 38 nfg = 45 f = 1.24096D-16 |proj g| = 2.30404D-08 +Iterate 39 nfg = 46 f = 3.14643D-17 |proj g| = 1.56779D-08 +Iterate 40 nfg = 47 f = 1.56638D-17 |proj g| = 1.19498D-08 +Iterate 41 nfg = 48 f = 2.80291D-18 |proj g| = 4.46388D-09 +Iterate 42 nfg = 49 f = 1.32837D-18 |proj g| = 2.55882D-09 +Iterate 43 nfg = 50 f = 2.54555D-19 |proj g| = 1.17507D-09 +Iterate 44 nfg = 52 f = 8.23938D-20 |proj g| = 9.01720D-10 +Iterate 45 nfg = 53 f = 2.79609D-20 |proj g| = 8.62899D-10 +Iterate 46 nfg = 54 f = 1.08725D-20 |proj g| = 2.44452D-10 +Iterate 47 nfg = 55 f = 5.32476D-21 |proj g| = 1.66320D-10 +Iterate 48 nfg = 56 f = 1.26624D-21 |proj g| = 1.90396D-10 +Iterate 49 nfg = 58 f = 5.35121D-22 |proj g| = 9.74083D-11 + STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL + Final X= + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 + 1.0000D+00 1.0001D+00 1.0002D+00 1.0003D+00 1.0007D+00 1.0014D+00 + 1.0028D+00 1.0055D+00 1.0111D+00 1.0223D+00 1.0451D+00 1.0923D+00 + 1.1930D+00 1.4233D+00 2.0257D+00 4.1035D+00 @@ -0,0 +1,232 @@ + +L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” or +“3-clause license”) +Please read attached file License.txt + + L-BFGS-B (version 3.0) march, 2011 + +This directory contains the modified/corrected limited memory +code for solving bound constrained optimization problems. + +_______________________________________________________________________ + + o f77 routines + + blas.f subset of blas + linpack.f subset of linpack + lbfgsb.f main routine + timer.f routine to compute CPU time + + o f77 and f90 drivers + + driver1.f driver1.f90 + driver2.f driver2.f90 + driver3.f driver3.f90 + + We recommend that the user read driver1.f/driver1.f90, which + give a good idea of how the code works + + o Makefile Linux/Unix users. + + To run the code on a Linux/Unix machine simply type + + make + + This will create the files + + x.lbfgsb_77_1 x.lbfgsb_90_1 + x.lbfgsb_77_2 x.lbfgsb_90_2 + x.lbfgsb_77_3 x.lbfgsb_90_3 + + which are the executable files for the drivers in the package. + One can execute the program by typing the name of the executable + file. + + o Output files are provided for each executable file in directory + OUPUTS + + output_77_1 output_90_1 + output_77_2 output_90_2 + output_77_3 output_90_3 + + o The following four articles are also included + + [1] algorithm.pdf - describes the algorithm + [2] compact.pdf - presents the compact form of matrices + [3] code.pdf - describes the code + [4] acm-remark.pdf - describes the modification/correction + + The most useful articles for those who only wish to use the code + are [3,4]. Users interested in understanding the algorithm should + read [1] and possibly also [2]. + +________________________________________________________________________ + +For questions and help contact + + Jorge Nocedal <nocedal@eecs.northwestern.edu> + Jose Luis Morales <jmorales@itam.mx> +________________________________________________________________________ + + + 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 + 5)nmax + 11mmax^2 + 8mmax 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 0<iprint<99 print also f and |proj g| every iprint iterations; +c iprint=99 print details of every iteration except n-vectors; +c iprint=100 print also the changes of active set and final x; +c iprint>100 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 +************************************************************************ + diff --git a/algorithm.pdf b/algorithm.pdf Binary files differnew file mode 100755 index 0000000..bfc9af6 --- /dev/null +++ b/algorithm.pdf diff --git a/code.pdf b/code.pdf Binary files differnew file mode 100755 index 0000000..b257b27 --- /dev/null +++ b/code.pdf diff --git a/driver1.f b/driver1.f new file mode 100644 index 0000000..b898e5d --- /dev/null +++ b/driver1.f @@ -0,0 +1,321 @@ +c +c L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” +c or “3-clause license”) +c Please read attached file License.txt +c +c DRIVER 1 in Fortran 77 +c -------------------------------------------------------------- +c SIMPLE DRIVER FOR L-BFGS-B (version 3.0) +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 March 2011 (latest revision) +c Optimization Center at Northwestern University +c Instituto Tecnologico Autonomo de Mexico +c +c Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778: +c L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained +c Optimization" (2011). To appear in ACM Transactions on +c Mathematical Software, + +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 + 5)nmax + 11mmax^2 + 8mmax 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 0<iprint<99 print also f and |proj g| every iprint iterations; +c iprint=99 print details of every iteration except n-vectors; +c iprint=100 print also the changes of active set and final x; +c iprint>100 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 -------------------------------------------------------------- + + 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 + 5*nmax + 11*mmax*mmax + 8*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 ============================ diff --git a/driver1.f90 b/driver1.f90 new file mode 100755 index 0000000..67fa540 --- /dev/null +++ b/driver1.f90 @@ -0,0 +1,297 @@ +! +! L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” +! or “3-clause license”) +! Please read attached file License.txt +! +! +! DRIVER1 in Fortran 90 +! -------------------------------------------------------------- +! +! L-BFGS-B is a code for solving large nonlinear optimization +! problems with simple bounds on the variables. +! +! The code can also be used for unconstrained problems and is +! as efficient for these problems as the earlier limited memory +! code L-BFGS. +! +! This is the simplest driver in the package. It uses all the +! default settings of the code. +! +! +! References: +! +! [1] 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, pp. 1190--1208. +! +! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN +! Subroutines for Large Scale Bound Constrained Optimization'' +! Tech. Report, NAM-11, 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.) +! +! * * * +! +! March 2011 (latest revision) +! Optimization Center at Northwestern University +! Instituto Tecnologico Autonomo de Mexico +! +! Jorge Nocedal and Jose Luis Morales +! +! -------------------------------------------------------------- +! DESCRIPTION OF THE VARIABLES IN L-BFGS-B +! -------------------------------------------------------------- +! +! n is an INTEGER variable that must be set by the user to the +! number of variables. It is not altered by the routine. +! +! m is an INTEGER variable that must be set by the user to the +! number of corrections used in the limited memory matrix. +! It is not altered by the routine. Values of m < 3 are +! not recommended, and large values of m can result in excessive +! computing time. The range 3 <= m <= 20 is recommended. +! +! x is a DOUBLE PRECISION array of length n. On initial entry +! it must be set by the user to the values of the initial +! estimate of the solution vector. Upon successful exit, it +! contains the values of the variables at the best point +! found (usually an approximate solution). +! +! l is a DOUBLE PRECISION array of length n that must be set by +! the user to the values of the lower bounds on the variables. If +! the i-th variable has no lower bound, l(i) need not be defined. +! +! u is a DOUBLE PRECISION array of length n that must be set by +! the user to the values of the upper bounds on the variables. If +! the i-th variable has no upper bound, u(i) need not be defined. +! +! nbd is an INTEGER array of dimension n that must be set by the +! user to the type of bounds imposed on the variables: +! nbd(i)=0 if x(i) is unbounded, +! 1 if x(i) has only a lower bound, +! 2 if x(i) has both lower and upper bounds, +! 3 if x(i) has only an upper bound. +! +! f is a DOUBLE PRECISION variable. If the routine setulb returns +! with task(1:2)= 'FG', then f must be set by the user to +! contain the value of the function at the point x. +! +! g is a DOUBLE PRECISION array of length n. If the routine setulb +! returns with taskb(1:2)= 'FG', then g must be set by the user to +! contain the components of the gradient at the point x. +! +! factr is a DOUBLE PRECISION variable that must be set by the user. +! It is a tolerance in the termination test for the algorithm. +! The iteration will stop when +! +! (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch +! +! where epsmch is the machine precision which is automatically +! generated by the code. Typical values for factr on a computer +! with 15 digits of accuracy in double precision are: +! factr=1.d+12 for low accuracy; +! 1.d+7 for moderate accuracy; +! 1.d+1 for extremely high accuracy. +! The user can suppress this termination test by setting factr=0. +! +! pgtol is a double precision variable. +! On entry pgtol >= 0 is specified by the user. The iteration +! will stop when +! +! max{|proj g_i | i = 1, ..., n} <= pgtol +! +! where pg_i is the ith component of the projected gradient. +! The user can suppress this termination test by setting pgtol=0. +! +! wa is a DOUBLE PRECISION array of length +! (2mmax + 5)nmax + 11mmax^2 + 8mmax used as workspace. +! This array must not be altered by the user. +! +! iwa is an INTEGER array of length 3nmax used as +! workspace. This array must not be altered by the user. +! +! task is a CHARACTER string of length 60. +! On first entry, it must be set to 'START'. +! On a return with task(1:2)='FG', the user must evaluate the +! function f and gradient g at the returned value of x. +! On a return with task(1:5)='NEW_X', an iteration of the +! algorithm has concluded, and f and g contain f(x) and g(x) +! respectively. The user can decide whether to continue or stop +! the iteration. +! When +! task(1:4)='CONV', the termination test in L-BFGS-B has been +! satisfied; +! task(1:4)='ABNO', the routine has terminated abnormally +! without being able to satisfy the termination conditions, +! x contains the best approximation found, +! f and g contain f(x) and g(x) respectively; +! task(1:5)='ERROR', the routine has detected an error in the +! input parameters; +! On exit with task = 'CONV', 'ABNO' or 'ERROR', the variable task +! contains additional information that the user can print. +! This array should not be altered unless the user wants to +! stop the run for some reason. See driver2 or driver3 +! for a detailed explanation on how to stop the run +! by assigning task(1:4)='STOP' in the driver. +! +! iprint is an INTEGER variable that must be set by the user. +! It controls the frequency and type of output generated: +! iprint<0 no output is generated; +! iprint=0 print only one line at the last iteration; +! 0<iprint<99 print also f and |proj g| every iprint iterations; +! iprint=99 print details of every iteration except n-vectors; +! iprint=100 print also the changes of active set and final x; +! iprint>100 print details of every iteration including x and g; +! When iprint > 0, the file iterate.dat will be created to +! summarize the iteration. +! +! csave is a CHARACTER working array of length 60. +! +! lsave is a LOGICAL working array of dimension 4. +! On exit with task = 'NEW_X', the following information is +! available: +! lsave(1) = .true. the initial x did not satisfy the bounds; +! lsave(2) = .true. the problem contains bounds; +! lsave(3) = .true. each variable has upper and lower bounds. +! +! isave is an INTEGER working array of dimension 44. +! On exit with task = 'NEW_X', it contains information that +! the user may want to access: +! isave(30) = the current iteration number; +! isave(34) = the total number of function and gradient +! evaluations; +! isave(36) = the number of function value or gradient +! evaluations in the current iteration; +! isave(38) = the number of free variables in the current +! iteration; +! isave(39) = the number of active constraints at the current +! iteration; +! +! see the subroutine setulb.f for a description of other +! information contained in isave +! +! dsave is a DOUBLE PRECISION working array of dimension 29. +! On exit with task = 'NEW_X', it contains information that +! the user may want to access: +! dsave(2) = the value of f at the previous iteration; +! dsave(5) = the machine precision epsmch generated by the code; +! dsave(13) = the infinity norm of the projected gradient; +! +! see the subroutine setulb.f for a description of other +! information contained in dsave +! +! -------------------------------------------------------------- +! END OF THE DESCRIPTION OF THE VARIABLES IN L-BFGS-B +! -------------------------------------------------------------- +! + program driver +! +! This simple driver demonstrates how to call the L-BFGS-B code to +! solve a sample problem (the extended Rosenbrock function +! subject to bounds on the variables). The dimension n of this +! problem is variable. + + implicit none +! +! Declare variables and parameters needed by the code. +! Note thar we wish to have output at every iteration. +! iprint=1 +! +! We also specify the tolerances in the stopping criteria. +! factr = 1.0d+7, pgtol = 1.0d-5 +! +! A description of all these variables is given at the beginning +! of the driver +! + integer, parameter :: n = 25, m = 5, iprint = 1 + integer, parameter :: dp = kind(1.0d0) + real(dp), parameter :: factr = 1.0d+7, pgtol = 1.0d-5 +! + character(len=60) :: task, csave + logical :: lsave(4) + integer :: isave(44) + real(dp) :: f + real(dp) :: dsave(29) + integer, allocatable :: nbd(:), iwa(:) + real(dp), allocatable :: x(:), l(:), u(:), g(:), wa(:) + +! Declare a few additional variables for this sample problem + + real(dp) :: t1, t2 + integer :: i + +! Allocate dynamic arrays + + allocate ( nbd(n), x(n), l(n), u(n), g(n) ) + allocate ( iwa(3*n) ) + allocate ( wa(2*m*n + 5*n + 11*m*m + 8*m) ) +! + do 10 i=1, n, 2 + nbd(i) = 2 + l(i) = 1.0d0 + u(i) = 1.0d2 + 10 continue + +! 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 + +! 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.)',/) + +! We start the iteration by initializing task. + + task = 'START' + +! The beginning of the loop + + do while(task(1:2).eq.'FG'.or.task.eq.'NEW_X'.or. & + task.eq.'START') + +! 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 + + 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 + +! 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 + + end if + end do + +! end of loop do while + + + end program driver + diff --git a/driver2.f b/driver2.f new file mode 100644 index 0000000..1ac91c9 --- /dev/null +++ b/driver2.f @@ -0,0 +1,232 @@ +c +c L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” +c or “3-clause license”) +c Please read attached file License.txt +c +c DRIVER 2 in Fortran 77 +c -------------------------------------------------------------- +c CUSTOMIZED DRIVER FOR L-BFGS-B (version 3.0) +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 February 2011 (latest revision) +c Optimization Center at Northwestern University +c Instituto Tecnologico Autonomo de Mexico +c +c Jorge Nocedal and Jose Luis Morales +c Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778: +c L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained +c Optimization" (2011). To appear in ACM Transactions on +c Mathematical Software, +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+5*nmax+11*mmax*mmax+8*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 ============================ + diff --git a/driver2.f90 b/driver2.f90 new file mode 100755 index 0000000..48e53dd --- /dev/null +++ b/driver2.f90 @@ -0,0 +1,220 @@ +! +! L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” +! or “3-clause license”) +! Please read attached file License.txt +! +! +! DRIVER 2 in Fortran 90 +! -------------------------------------------------------------- +! CUSTOMIZED DRIVER FOR L-BFGS-B +! -------------------------------------------------------------- +! +! L-BFGS-B is a code for solving large nonlinear optimization +! problems with simple bounds on the variables. +! +! The code can also be used for unconstrained problems and is +! as efficient for these problems as the earlier limited memory +! code L-BFGS. +! +! This driver illustrates how to control the termination of the +! run and how to design customized output. +! +! References: +! +! [1] 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, pp. 1190--1208. +! +! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN +! Subroutines for Large Scale Bound Constrained Optimization'' +! Tech. Report, NAM-11, 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.) +! +! * * * +! +! February 2011 (latest revision) +! Optimization Center at Northwestern University +! Instituto Tecnologico Autonomo de Mexico +! +! Jorge Nocedal and Jose Luis Morales +! +! ************** + program driver + +! This driver shows how to replace the default stopping test +! by other termination criteria. It also illustrates how to +! print the values of several parameters during the course of +! the iteration. The sample problem used here is the same as in +! DRIVER1 (the extended Rosenbrock function with bounds on the +! variables). + + implicit none + +! Declare variables and parameters needed by the code. +! +! Note that we suppress the default output (iprint = -1) +! We suppress both code-supplied stopping tests because the +! user is providing his/her own stopping criteria. + + integer, parameter :: n = 25, m = 5, iprint = -1 + integer, parameter :: dp = kind(1.0d0) + real(dp), parameter :: factr = 0.0d0, pgtol = 0.0d0 + + character(len=60) :: task, csave + logical :: lsave(4) + integer :: isave(44) + real(dp) :: f + real(dp) :: dsave(29) + integer, allocatable :: nbd(:), iwa(:) + real(dp), allocatable :: x(:), l(:), u(:), g(:), wa(:) +! + real(dp) :: t1, t2 + integer :: i + + allocate ( nbd(n), x(n), l(n), u(n), g(n) ) + allocate ( iwa(3*n) ) + allocate ( wa(2*m*n + 5*n + 11*m*m + 8*m) ) +! +! This driver shows how to replace the default stopping test +! by other termination criteria. It also illustrates how to +! print the values of several parameters during the course of +! the iteration. The sample problem used here is the same as in +! DRIVER1 (the extended Rosenbrock function with bounds on the +! variables). +! We now specify nbd which defines the bounds on the variables: +! l specifies the lower bounds, +! u specifies the upper bounds. + +! 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 + +! 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 + +! We now define the starting point. + + do 14 i=1, n + x(i)=3.0d0 + 14 continue + +! 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.)',/) + + +! We start the iteration by initializing task. +! + task = 'START' + +! ------- the beginning of the loop ---------- + + do while( task(1:2).eq.'FG'.or.task.eq.'NEW_X'.or. & + task.eq.'START') + +! 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 + +! the minimization routine has returned to request the +! function f and gradient g values at the current x. + +! 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 + +! 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 +! + else +! + if (task(1:5) .eq. 'NEW_X') then +! +! the minimization routine has returned with a new iterate. +! At this point have the opportunity of stopping the iteration +! or observing the values of certain parameters +! +! First are two examples of stopping tests. + +! Note: task(1:4) must be assigned the value 'STOP' to terminate +! the iteration and ensure that the final results are +! printed in the default format. The rest of the character +! string TASK may be used to store other information. + +! 1) Terminate if the total number of f and g evaluations +! exceeds 99. + + if (isave(34) .ge. 99) & + task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT' + +! 2) Terminate if |proj g|/(1+|f|) < 1.0d-10, where +! "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' + +! We now wish to print the following information at each +! iteration: +! +! 1) the current iteration number, isave(30), +! 2) the total number of f and g evaluations, isave(34), +! 3) the value of the objective function f, +! 4) the norm of the projected gradient, dsve(13) +! +! See the comments at the end of driver1 for a description +! 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) + +! If the run is to be terminated, we print also the information +! 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) + end if + + end if + end if + + end do +! ---------- the end of the loop ------------- + +! If task is neither FG nor NEW_X we terminate execution. + + end program driver + +!======================= The end of driver2 ============================ + diff --git a/driver3.f b/driver3.f new file mode 100644 index 0000000..8d1371d --- /dev/null +++ b/driver3.f @@ -0,0 +1,273 @@ +c +c L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” +c or “3-clause license”) +c Please read attached file License.txt +c +c DRIVER 3 in Fortran 77 +c -------------------------------------------------------------- +c TIME-CONTROLLED DRIVER FOR L-BFGS-B (version 3.0) +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 February 2011 (latest revision) +c Optimization Center at Northwestern University +c Instituto Tecnologico Autonomo de Mexico +c +c Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778: +c L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained +c Optimization" (2011). To appear in ACM Transactions on +c Mathematical Software, +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+5*nmax+11*mmax*mmax+8*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+11*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 ============================ + diff --git a/driver3.f90 b/driver3.f90 new file mode 100644 index 0000000..4cc9b36 --- /dev/null +++ b/driver3.f90 @@ -0,0 +1,256 @@ +! +! L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” +! or “3-clause license”) +! Please read attached file License.txt +! +! DRIVER 3 in Fortran 90 +! -------------------------------------------------------------- +! TIME-CONTROLLED DRIVER FOR L-BFGS-B +! -------------------------------------------------------------- +! +! L-BFGS-B is a code for solving large nonlinear optimization +! problems with simple bounds on the variables. +! +! The code can also be used for unconstrained problems and is +! as efficient for these problems as the earlier limited memory +! code L-BFGS. +! +! This driver shows how to terminate a run after some prescribed +! CPU time has elapsed, and how to print the desired information +! before exiting. +! +! References: +! +! [1] 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, pp. 1190--1208. +! +! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN +! Subroutines for Large Scale Bound Constrained Optimization'' +! Tech. Report, NAM-11, 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.) +! +! * * * +! +! February 2011 (latest revision) +! Optimization Center at Northwestern University +! Instituto Tecnologico Autonomo de Mexico +! +! Jorge Nocedal and Jose Luis Morales +! +! ************** + + program driver + +! This time-controlled driver shows that it is possible to terminate +! a run by elapsed CPU time, and yet be able to print all desired +! information. This driver also illustrates the use of two +! stopping criteria that may be used in conjunction with a limit +! on execution time. The sample problem used here is the same as in +! driver1 and driver2 (the extended Rosenbrock function with bounds +! on the variables). + + implicit none + +! We specify a limit on the CPU time (tlimit = 10 seconds) +! +! We suppress the default output (iprint = -1). The user could +! also elect to use the default output by choosing iprint >= 0.) +! We suppress the code-supplied stopping tests because we will +! provide our own termination conditions +! We specify the dimension n of the sample problem and the number +! m of limited memory corrections stored. + + integer, parameter :: n = 1000, m = 10, iprint = -1 + integer, parameter :: dp = kind(1.0d0) + real(dp), parameter :: factr = 0.0d0, pgtol = 0.0d0, & + tlimit = 10.0d0 +! + character(len=60) :: task, csave + logical :: lsave(4) + integer :: isave(44) + real(dp) :: f + real(dp) :: dsave(29) + integer, allocatable :: nbd(:), iwa(:) + real(dp), allocatable :: x(:), l(:), u(:), g(:), wa(:) +! + real(dp) :: t1, t2, time1, time2 + integer :: i, j + + allocate ( nbd(n), x(n), l(n), u(n), g(n) ) + allocate ( iwa(3*n) ) + allocate ( wa(2*m*n + 5*n + 11*m*m + 8*m) ) + +! This time-controlled driver shows that it is possible to terminate +! a run by elapsed CPU time, and yet be able to print all desired +! information. This driver also illustrates the use of two +! stopping criteria that may be used in conjunction with a limit +! on execution time. The sample problem used here is the same as in +! driver1 and driver2 (the extended Rosenbrock function with bounds +! on the variables). + +! We now specify nbd which defines the bounds on the variables: +! l specifies the lower bounds, +! u specifies the upper bounds. + +! 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 + +! 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 + +! We now define the starting point. + + do 14 i=1, n + x(i)=3.0d0 + 14 continue + +! 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.)',/) + +! We start the iteration by initializing task. + + task = 'START' + +! ------- the beginning of the loop ---------- + +! We begin counting the CPU time. + + call timer(time1) + + do while( task(1:2).eq.'FG'.or.task.eq.'NEW_X'.or. & + task.eq.'START') + +! 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 + +! the minimization routine has returned to request the +! function f and gradient g values at the current x. +! 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.' + +! Note: Assigning task(1:4)='STOP' will terminate the run; +! setting task(7:9)='CPU' will restore the information at +! the latest iterate generated by the code so that it can +! be correctly printed by the driver. + +! In this driver we have chosen to disable the +! printing options of the code (we set iprint=-1); +! instead we are using customized output: we print the +! latest value of x, the corresponding function value f and +! the norm of the projected gradient |proj g|. + +! We print out the information contained in task. + + write (6,*) task + +! We print the latest iterate contained in wa(j+1:j+n), where + + j = 3*n+2*m*n+11*m**2 + write (6,*) 'Latest iterate X =' + write (6,'((1x,1p, 6(1x,d11.4)))') (wa(i),i = j+1,j+n) + +! We print the function value f and the norm of the projected +! gradient |proj g| at the last iterate; they are stored in +! 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 + +! The time limit has not been reached and we compute +! 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 + +! 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 + +! go back to the minimization routine. + else + + if (task(1:5) .eq. 'NEW_X') then + +! the minimization routine has returned with a new iterate. +! The time limit has not been reached, and we test whether +! the following two stopping tests are satisfied: + +! 1) Terminate if the total number of f and g evaluations +! exceeds 900. + + if (isave(34) .ge. 900) & + task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT' + +! 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' + +! We wish to print the following information at each iteration: +! 1) the current iteration number, isave(30), +! 2) the total number of f and g evaluations, isave(34), +! 3) the value of the objective function f, +! 4) the norm of the projected gradient, dsve(13) +! +! See the comments at the end of driver1 for a description +! 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) + +! If the run is to be terminated, we print also the information +! 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 + + endif + end if + end do + +! If task is neither FG nor NEW_X we terminate execution. + + end program driver + +!======================= The end of driver3 ============================ + diff --git a/iterate.dat b/iterate.dat new file mode 100644 index 0000000..b9dbf62 --- /dev/null +++ b/iterate.dat @@ -0,0 +1,49 @@ +RUNNING THE L-BFGS-B CODE + +it = iteration number +nf = number of function evaluations +nseg = 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 = 2.220D-16 + N = 25 M = 5 + + it nf nseg nact sub itls stepl tstep projg f + 0 1 - - - - - - 1.030D+02 3.460D+03 + 1 5 25 24 --- 3 1.2D-02 4.2D+00 6.507D+01 2.398D+03 + 2 6 1 0 con 0 1.0D+00 4.2D+00 3.640D+01 1.438D+02 + 3 7 1 0 con 0 1.0D+00 6.8D-01 2.290D+01 7.282D+01 + 4 8 1 0 con 0 1.0D+00 1.1D+00 6.954D+00 1.603D+01 + 5 9 1 0 con 0 1.0D+00 6.5D-01 9.055D+00 5.187D+00 + 6 10 1 0 bnd 0 1.0D+00 6.1D-01 1.967D+01 2.127D+00 + 7 11 1 0 bnd 0 1.0D+00 1.8D-01 2.128D+00 2.076D-01 + 8 12 1 0 bnd 0 1.0D+00 9.0D-02 8.325D-01 5.327D-02 + 9 13 1 0 bnd 0 1.0D+00 5.3D-02 4.279D-01 1.305D-02 + 10 14 1 1 bnd 0 1.0D+00 2.5D-02 2.008D-01 3.860D-03 + 11 15 1 2 bnd 0 1.0D+00 2.0D-02 1.377D-01 7.457D-04 + 12 16 1 0 bnd 0 1.0D+00 7.8D-03 1.213D-01 3.540D-04 + 13 17 1 2 con 0 1.0D+00 3.1D-03 2.978D-02 7.425D-05 + 14 18 1 0 con 0 1.0D+00 1.1D-03 1.727D-02 3.741D-05 + 15 19 1 0 con 0 1.0D+00 2.1D-03 2.868D-02 1.098D-05 + 16 21 1 0 con 1 4.6D-01 4.8D-04 8.081D-03 3.908D-06 + 17 22 1 0 con 0 1.0D+00 2.2D-04 3.476D-03 1.995D-06 + 18 23 1 0 con 0 1.0D+00 2.2D-04 2.253D-03 8.258D-07 + 19 24 1 0 con 0 1.0D+00 2.4D-04 1.458D-03 1.992D-07 + 20 25 1 0 con 0 1.0D+00 1.8D-04 1.482D-03 5.758D-08 + 21 26 1 0 con 0 1.0D+00 6.8D-05 5.443D-04 1.463D-08 + 22 27 1 0 con 0 1.0D+00 2.3D-05 2.252D-04 2.363D-09 + 23 28 1 0 con 0 1.0D+00 8.0D-06 1.721D-04 1.083D-09 + +CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH + + Total User time 3.337E-03 seconds. + diff --git a/lbfgsb.f b/lbfgsb.f new file mode 100644 index 0000000..9c9e7d9 --- /dev/null +++ b/lbfgsb.f @@ -0,0 +1,3949 @@ +c +c L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” +c or “3-clause license”) +c Please read attached file License.txt +c +c=========== L-BFGS-B (version 3.0. April 25, 2011 =================== +c +c This is a modified version of L-BFGS-B. Minor changes in the updated +c code appear preceded by a line comment as follows +c +c c-jlm-jn +c +c Major changes are described in the accompanying paper: +c +c Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778: +c L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained +c Optimization" (2011). To appear in ACM Transactions on +c Mathematical Software, +c +c The paper describes an improvement and a correction to Algorithm 778. +c It is shown that the performance of the algorithm can be improved +c significantly by making a relatively simple modication to the subspace +c minimization phase. The correction concerns an error caused by the use +c of routine dpmeps to estimate machine precision. +c +c The total work space **wa** required by the new version is +c +c 2*m*n + 11m*m + 5*n + 8*m +c +c the old version required +c +c 2*m*n + 12m*m + 4*n + 12*m +c +c +c J. Nocedal Department of Electrical Engineering and +c Computer Science. +c Northwestern University. Evanston, IL. USA +c +c +c J.L Morales Departamento de Matematicas, +c Instituto Tecnologico Autonomo de Mexico +c Mexico D.F. Mexico. +c +c March 2011 +c +c============================================================================= + 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), +c +c-jlm-jn + + wa(2*m*n + 5*n + 11*m*m + 8*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 + 5)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 0<iprint<99 print also f and |proj g| every iprint iterations; +c iprint=99 print details of every iteration except n-vectors; +c iprint=100 print also the changes of active set and final x; +c iprint>100 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 ************ +c-jlm-jn + integer lws,lr,lz,lt,ld,lxp,lwa, + + lwy,lsy,lss,lwt,lwn,lsnd + + if (task .eq. 'START') then + isave(1) = m*n + isave(2) = m**2 + isave(3) = 4*m**2 + isave(4) = 1 ! ws m*n + isave(5) = isave(4) + isave(1) ! wy m*n + isave(6) = isave(5) + isave(1) ! wsy m**2 + isave(7) = isave(6) + isave(2) ! wss m**2 + isave(8) = isave(7) + isave(2) ! wt m**2 + isave(9) = isave(8) + isave(2) ! wn 4*m**2 + isave(10) = isave(9) + isave(3) ! wsnd 4*m**2 + isave(11) = isave(10) + isave(3) ! wz n + isave(12) = isave(11) + n ! wr n + isave(13) = isave(12) + n ! wd n + isave(14) = isave(13) + n ! wt n + isave(15) = isave(14) + n ! wxp n + isave(16) = isave(15) + n ! wa 8*m + endif + lws = isave(4) + lwy = isave(5) + lsy = isave(6) + lss = isave(7) + lwt = isave(8) + lwn = isave(9) + lsnd = isave(10) + lz = isave(11) + lr = isave(12) + ld = isave(13) + lt = isave(14) + lxp = isave(15) + lwa = isave(16) + + call mainlb(n,m,x,l,u,nbd,f,g,factr,pgtol, + + wa(lws),wa(lwy),wa(lsy),wa(lss), wa(lwt), + + wa(lwn),wa(lsnd),wa(lz),wa(lr),wa(ld),wa(lt),wa(lxp), + + wa(lwa), + + 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, wt, wn, snd, z, r, d, t, xp, wa, + + index, iwhere, indx2, task, + + iprint, csave, lsave, isave, dsave) + implicit none + 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), +c-jlm-jn + + xp(n), + + wa(8*m), + + ws(n, m), wy(n, m), sy(m, m), ss(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), xp(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 xp is used to safeguard the projected Newton direction +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 0<iprint<99 print also f and |proj g| every iprint iterations; +c iprint=99 print details of every iteration except n-vectors; +c iprint=100 print also the changes of active set and final x; +c iprint>100 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 +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, + + nseg,nfgv,info,ifun, + + iword,nfree,nact,ileave,nenter + double precision theta,fold,ddot,dr,rr,tol, + + 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 + + epsmch = epsilon(one) + + call timer(time1) + +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. + iback = 0 + itail = 0 + iword = 0 + nact = 0 + ileave = 0 + nenter = 0 + fold = zero + dnorm = zero + cpu1 = zero + gd = zero + stpmx = zero + sbgnrm = zero + stp = zero + gdold = zero + dtd = zero + +c for operation counts: + iter = 0 + nfgv = 0 + nseg = 0 + nintol = 0 + nskip = 0 + nfree = n + ifun = 0 +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 + + itfile = 8 + if (iprint .ge. 1) then +c open a summary file 'iterate.dat' + open (8, file = 'iterate.dat', status = 'unknown') + 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,nseg,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) + nseg = isave(12) + nfgv = isave(13) + info = isave(14) + ifun = isave(15) + iword = isave(16) + nfree = isave(17) + nact = isave(18) + ileave = isave(19) + nenter = isave(20) + + theta = dsave(1) + fold = dsave(2) + tol = dsave(3) + dnorm = dsave(4) + epsmch = dsave(5) + cpu1 = dsave(6) + cachyt = dsave(7) + sbtime = dsave(8) + lnscht = dsave(9) + time1 = dsave(10) + gd = dsave(11) + stpmx = dsave(12) + sbgnrm = dsave(13) + stp = dsave(14) + gdold = dsave(15) + dtd = dsave(16) + +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 + nseg = 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),nseg, + + iprint, sbgnrm, info, epsmch) + 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 + nseg + +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=theta*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-jlm-jn call the direct method. + + call subsm( n, m, nfree, index, l, u, nbd, z, r, xp, ws, wy, + + theta, x, g, col, head, iword, wa, wn, iprint, info) + 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,nseg,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,nseg,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) = nseg + isave(13) = nfgv + isave(14) = info + isave(15) = ifun + isave(16) = iword + isave(17) = nfree + isave(18) = nact + isave(19) = ileave + isave(20) = nenter + + dsave(1) = theta + dsave(2) = fold + dsave(3) = tol + dsave(4) = dnorm + dsave(5) = epsmch + dsave(6) = cpu1 + dsave(7) = cachyt + dsave(8) = sbtime + dsave(9) = lnscht + dsave(10) = time1 + dsave(11) = gd + dsave(12) = stpmx + dsave(13) = sbgnrm + dsave(14) = stp + dsave(15) = gdold + dsave(16) = dtd + + 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, nseg, iprint, sbgnrm, info, epsmch) + implicit none + integer n, m, head, col, nseg, iprint, info, + + nbd(n), iorder(n), iwhere(n) + double precision theta, epsmch, + + x(n), l(n), u(n), g(n), t(n), d(n), xcp(n), + + 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 nseg is an integer variable. +c On exit nseg 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 0<iprint<99 print also f and |proj g| every iprint iterations; +c iprint=99 print details of every iteration except n-vectors; +c iprint=100 print also the changes of active set and final x; +c iprint>100 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, + + f2_org + 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 + f2_org = f2 + 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 + nseg = 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) nseg,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. + + nseg = nseg + 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 + + f2 = max(epsmch*f2_org,f2) + 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) nseg,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) + +c write(6,*) ' k = index(i) ', k +c write(6,*) ' 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. + write(6,*)' ascent direction in projection gd = ', gd + 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',/, + + 'nseg = 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,'nseg',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, nseg, word, iword, iback, stp, xstep) + + character*3 word + integer n, iprint, itfile, iter, nfgv, nact, nseg, + + 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,nseg,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, nseg, word, iback, stp, xstep, k, + + cachyt, sbtime, lnscht) + + character*60 task + character*3 word + integer n, iprint, info, itfile, iter, nfgv, nintol, + + nskip, nact, nseg, 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,nseg,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',4x,'Tit',5x,'Tnf',2x,'Tnint',2x, + + 'Skip',2x,'Nact',5x,'Projg',8x,'F') + 3005 format (i5,2(1x,i6),(1x,i6),(2x,i4),(1x,i5),1p,2(2x,d10.3)) + 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, xp, ws, wy, + + theta, xx, gg, + + col, head, iword, wv, wn, iprint, info ) + implicit none + integer n, m, nsub, col, head, iword, iprint, info, + + ind(nsub), nbd(n) + double precision theta, + + l(n), u(n), x(n), d(n), xp(n), xx(n), gg(n), + + ws(n, m), wy(n, m), + + wv(2*m), wn(2*m, 2*m) + +c ********************************************************************** +c +c This routine contains the major changes in the updated version. +c The changes are described in the accompanying paper +c +c Jose Luis Morales, Jorge Nocedal +c "Remark On Algorithm 788: L-BFGS-B: Fortran Subroutines for Large-Scale +c Bound Constrained Optimization". Decemmber 27, 2010. +c +c J.L. Morales Departamento de Matematicas, +c Instituto Tecnologico Autonomo de Mexico +c Mexico D.F. +c +c J, Nocedal Department of Electrical Engineering and +c Computer Science. +c Northwestern University. Evanston, IL. USA +c +c January 17, 2011 +c +c ********************************************************************** +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 xp is a double precision array of dimension n. +c used to safeguard the projected Newton direction +c +c xx is a double precision array of dimension n +c On entry it holds the current iterate +c On output it is unchanged + +c gg is a double precision array of dimension n +c On entry it holds the gradient at the current iterate +c On output it is unchanged +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 0<iprint<99 print also f and |proj g| every iprint iterations; +c iprint=99 print details of every iteration except n-vectors; +c iprint=100 print also the changes of active set and final x; +c iprint>100 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, xk, dk, temp1, temp2 + double precision one,zero + parameter (one=1.0d0,zero=0.0d0) +c + double precision dd_p + + 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 + + call dscal( nsub, one/theta, d, 1 ) +c +c----------------------------------------------------------------- +c Let us try the projection, d is the Newton direction + + iword = 0 + + call dcopy ( n, x, 1, xp, 1 ) +c + do 50 i=1, nsub + k = ind(i) + dk = d(i) + xk = x(k) + if ( nbd(k) .ne. 0 ) then +c + if ( nbd(k).eq.1 ) then ! lower bounds only + x(k) = max( l(k), xk + dk ) + if ( x(k).eq.l(k) ) iword = 1 + else +c + if ( nbd(k).eq.2 ) then ! upper and lower bounds + xk = max( l(k), xk + dk ) + x(k) = min( u(k), xk ) + if ( x(k).eq.l(k) .or. x(k).eq.u(k) ) iword = 1 + else +c + if ( nbd(k).eq.3 ) then ! upper bounds only + x(k) = min( u(k), xk + dk ) + if ( x(k).eq.u(k) ) iword = 1 + end if + end if + end if +c + else ! free variables + x(k) = xk + dk + end if + 50 continue +c + if ( iword.eq.0 ) then + go to 911 + end if +c +c check sign of the directional derivative +c + dd_p = zero + do 55 i=1, n + dd_p = dd_p + (x(i) - xx(i))*gg(i) + 55 continue + if ( dd_p .gt.zero ) then + call dcopy( n, xp, 1, x, 1 ) + write(6,*) ' Positive dir derivative in projection ' + write(6,*) ' Using the backtracking step ' + else + go to 911 + endif +c +c----------------------------------------------------------------- +c + alpha = one + temp1 = alpha + ibd = 0 + 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 +cccccc + 911 continue + + if (iprint .ge. 99) write (6,1004) + + 1001 format (/,'----------------SUBSM entered-----------------',/) + 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 + + return + 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 + + return + end + @@ -0,0 +1,32 @@ +c +c L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” +c or “3-clause license”) +c Please read attached file License.txt +c + subroutine timer(ttime) + double precision ttime +c + real temp +c +c This routine computes cpu time in double precision; it makes use of +c the intrinsic f90 cpu_time therefore a conversion type is +c needed. +c +c J.L Morales Departamento de Matematicas, +c Instituto Tecnologico Autonomo de Mexico +c Mexico D.F. +c +c J.L Nocedal Department of Electrical Engineering and +c Computer Science. +c Northwestern University. Evanston, IL. USA +c +c January 21, 2011 +c + temp = sngl(ttime) + call cpu_time(temp) + ttime = dble(temp) + + return + + end + |