From 173802e9f8d98a01019d3f5aff055f1f04479974 Mon Sep 17 00:00:00 2001 From: Gard Spreemann Date: Thu, 2 Feb 2017 21:23:01 +0100 Subject: upstream/3.0+dfsg.3: Start git repository for package by importing upstream's 3.0 tarball DFSG-cleaned as currently used in published Debian package. --- License.txt | 71 + Makefile | 37 + OUTPUTS/output_77_1 | 85 ++ OUTPUTS/output_77_2 | 57 + OUTPUTS/output_77_3 | 222 +++ OUTPUTS/output_90_1 | 85 ++ OUTPUTS/output_90_2 | 57 + OUTPUTS/output_90_3 | 222 +++ README | 232 +++ algorithm.pdf | Bin 0 -> 260255 bytes code.pdf | Bin 0 -> 180639 bytes driver1.f | 321 +++++ driver1.f90 | 297 ++++ driver2.f | 232 +++ driver2.f90 | 220 +++ driver3.f | 273 ++++ driver3.f90 | 256 ++++ iterate.dat | 49 + lbfgsb.f | 3949 +++++++++++++++++++++++++++++++++++++++++++++++++++ timer.f | 32 + 20 files changed, 6697 insertions(+) create mode 100644 License.txt create mode 100644 Makefile create mode 100644 OUTPUTS/output_77_1 create mode 100644 OUTPUTS/output_77_2 create mode 100644 OUTPUTS/output_77_3 create mode 100644 OUTPUTS/output_90_1 create mode 100644 OUTPUTS/output_90_2 create mode 100644 OUTPUTS/output_90_3 create mode 100644 README create mode 100755 algorithm.pdf create mode 100755 code.pdf create mode 100644 driver1.f create mode 100755 driver1.f90 create mode 100644 driver2.f create mode 100755 driver2.f90 create mode 100644 driver3.f create mode 100644 driver3.f90 create mode 100644 iterate.dat create mode 100644 lbfgsb.f create mode 100644 timer.f 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) , +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 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 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 diff --git a/README b/README new file mode 100644 index 0000000..dc040e6 --- /dev/null +++ b/README @@ -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 + Jose Luis Morales +________________________________________________________________________ + + + 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 0100 print details of every iteration including x and g; +c When iprint > 0, the file iterate.dat will be created to +c summarize the iteration. +c +c csave is a CHARACTER working array of length 60. +c +c lsave is a LOGICAL working array of dimension 4. +c On exit with task = 'NEW_X', the following information is +c available: +c lsave(1) = .true. the initial x did not satisfy the bounds; +c lsave(2) = .true. the problem contains bounds; +c lsave(3) = .true. each variable has upper and lower bounds. +c +c isave is an INTEGER working array of dimension 44. +c On exit with task = 'NEW_X', it contains information that +c the user may want to access: +c isave(30) = the current iteration number; +c isave(34) = the total number of function and gradient +c evaluations; +c isave(36) = the number of function value or gradient +c evaluations in the current iteration; +c isave(38) = the number of free variables in the current +c iteration; +c isave(39) = the number of active constraints at the current +c iteration; +c +c see the subroutine setulb.f for a description of other +c information contained in isave +c +c dsave is a DOUBLE PRECISION working array of dimension 29. +c On exit with task = 'NEW_X', it contains information that +c the user may want to access: +c dsave(2) = the value of f at the previous iteration; +c dsave(5) = the machine precision epsmch generated by the code; +c dsave(13) = the infinity norm of the projected gradient; +c +c see the subroutine setulb.f for a description of other +c information contained in dsave +c +************************************************************************ + diff --git a/algorithm.pdf b/algorithm.pdf new file mode 100755 index 0000000..bfc9af6 Binary files /dev/null and b/algorithm.pdf differ diff --git a/code.pdf b/code.pdf new file mode 100755 index 0000000..b257b27 Binary files /dev/null and b/code.pdf differ 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 0100 print details of every iteration including x and g; +c When iprint > 0, the file iterate.dat will be created to +c summarize the iteration. +c +c csave is a CHARACTER working array of length 60. +c +c lsave is a LOGICAL working array of dimension 4. +c On exit with task = 'NEW_X', the following information is +c available: +c lsave(1) = .true. the initial x did not satisfy the bounds; +c lsave(2) = .true. the problem contains bounds; +c lsave(3) = .true. each variable has upper and lower bounds. +c +c isave is an INTEGER working array of dimension 44. +c On exit with task = 'NEW_X', it contains information that +c the user may want to access: +c isave(30) = the current iteration number; +c isave(34) = the total number of function and gradient +c evaluations; +c isave(36) = the number of function value or gradient +c evaluations in the current iteration; +c isave(38) = the number of free variables in the current +c iteration; +c isave(39) = the number of active constraints at the current +c iteration; +c +c see the subroutine setulb.f for a description of other +c information contained in isave +c +c dsave is a DOUBLE PRECISION working array of dimension 29. +c On exit with task = 'NEW_X', it contains information that +c the user may want to access: +c dsave(2) = the value of f at the previous iteration; +c dsave(5) = the machine precision epsmch generated by the code; +c dsave(13) = the infinity norm of the projected gradient; +c +c see the subroutine setulb.f for a description of other +c information contained in dsave +c +c -------------------------------------------------------------- +c END OF THE DESCRIPTION OF THE VARIABLES IN L-BFGS-B +c -------------------------------------------------------------- + + 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; +! 0100 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 0100 print details of every iteration including x and g; +c When iprint > 0, the file iterate.dat will be created to +c summarize the iteration. +c +c csave is a working string of characters of length 60. +c +c lsave is a logical working array of dimension 4. +c On exit with 'task' = NEW_X, the following information is +c available: +c If lsave(1) = .true. then the initial X has been replaced by +c its projection in the feasible set; +c If lsave(2) = .true. then the problem is constrained; +c If lsave(3) = .true. then each variable has upper and lower +c bounds; +c +c isave is an integer working array of dimension 44. +c On exit with 'task' = NEW_X, the following information is +c available: +c isave(22) = the total number of intervals explored in the +c search of Cauchy points; +c isave(26) = the total number of skipped BFGS updates before +c the current iteration; +c isave(30) = the number of current iteration; +c isave(31) = the total number of BFGS updates prior the current +c iteration; +c isave(33) = the number of intervals explored in the search of +c Cauchy point in the current iteration; +c isave(34) = the total number of function and gradient +c evaluations; +c isave(36) = the number of function value or gradient +c evaluations in the current iteration; +c if isave(37) = 0 then the subspace argmin is within the box; +c if isave(37) = 1 then the subspace argmin is beyond the box; +c isave(38) = the number of free variables in the current +c iteration; +c isave(39) = the number of active constraints in the current +c iteration; +c n + 1 - isave(40) = the number of variables leaving the set of +c active constraints in the current iteration; +c isave(41) = the number of variables entering the set of active +c constraints in the current iteration. +c +c dsave is a double precision working array of dimension 29. +c On exit with 'task' = NEW_X, the following information is +c available: +c dsave(1) = current 'theta' in the BFGS matrix; +c dsave(2) = f(x) in the previous iteration; +c dsave(3) = factr*epsmch; +c dsave(4) = 2-norm of the line search direction vector; +c dsave(5) = the machine precision epsmch generated by the code; +c dsave(7) = the accumulated time spent on searching for +c Cauchy points; +c dsave(8) = the accumulated time spent on +c subspace minimization; +c dsave(9) = the accumulated time spent on line search; +c dsave(11) = the slope of the line search function at +c the current point of line search; +c dsave(12) = the maximum relative step length imposed in +c line search; +c dsave(13) = the infinity norm of the projected gradient; +c dsave(14) = the relative step length in the line search; +c dsave(15) = the slope of the line search function at +c the starting point of the line search; +c dsave(16) = the square of the 2-norm of the line search +c direction vector. +c +c Subprograms called: +c +c L-BFGS-B Library ... mainlb. +c +c +c References: +c +c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited +c memory algorithm for bound constrained optimization'', +c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. +c +c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a +c limited memory FORTRAN code for solving bound constrained +c optimization problems'', Tech. Report, NAM-11, EECS Department, +c Northwestern University, 1994. +c +c (Postscript files of these papers are available via anonymous +c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) +c +c * * * +c +c NEOS, November 1994. (Latest revision June 1996.) +c Optimization Technology Center. +c Argonne National Laboratory and Northwestern University. +c Written by +c Ciyou Zhu +c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. +c +c +c ************ +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 0100 print details of every iteration including x and g; +c When iprint > 0, the file iterate.dat will be created to +c summarize the iteration. +c +c csave is a working string of characters of length 60. +c +c lsave is a logical working array of dimension 4. +c +c isave is an integer working array of dimension 23. +c +c dsave is a double precision working array of dimension 29. +c +c +c Subprograms called +c +c L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk, +c +c errclb, prn1lb, prn2lb, prn3lb, active, projgr, +c +c freev, cmprlb, matupd, formt. +c +c Minpack2 Library ... timer +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 0100 print details of every iteration including x and g; +c When iprint > 0, the file iterate.dat will be created to +c summarize the iteration. +c +c sbgnrm is a double precision variable. +c On entry sbgnrm is the norm of the projected gradient at x. +c On exit sbgnrm is unchanged. +c +c info is an integer variable. +c On entry info is 0. +c On exit info = 0 for normal return, +c = nonzero for abnormal return when the the system +c used in routine bmv is singular. +c +c Subprograms called: +c +c L-BFGS-B Library ... hpsolb, bmv. +c +c Linpack ... dscal dcopy, daxpy. +c +c +c References: +c +c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited +c memory algorithm for bound constrained optimization'', +c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. +c +c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN +c Subroutines for Large Scale Bound Constrained Optimization'' +c Tech. Report, NAM-11, EECS Department, Northwestern University, +c 1994. +c +c (Postscript files of these papers are available via anonymous +c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) +c +c * * * +c +c NEOS, November 1994. (Latest revision June 1996.) +c Optimization Technology Center. +c Argonne National Laboratory and Northwestern University. +c Written by +c Ciyou Zhu +c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. +c +c +c ************ + + logical xlower,xupper,bnded + integer i,j,col2,nfree,nbreak,pointr, + + ibp,nleft,ibkmin,iter + double precision f1,f2,dt,dtm,tsum,dibp,zibp,dibp2,bkmin, + + tu,tl,wmc,wmp,wmw,ddot,tj,tj0,neggi,sbgnrm, + + 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 0100 print details of every iteration including x and g; +c When iprint > 0, the file iterate.dat will be created to +c summarize the iteration. +c +c info is an integer variable. +c On entry info is unspecified. +c On exit info = 0 for normal return, +c = nonzero for abnormal return +c when the matrix K is ill-conditioned. +c +c Subprograms called: +c +c Linpack dtrsl. +c +c +c References: +c +c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited +c memory algorithm for bound constrained optimization'', +c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. +c +c +c +c * * * +c +c NEOS, November 1994. (Latest revision June 1996.) +c Optimization Technology Center. +c Argonne National Laboratory and Northwestern University. +c Written by +c Ciyou Zhu +c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. +c +c +c ************ + + integer pointr,m2,col2,ibd,jy,js,i,j,k + double precision alpha, 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 + diff --git a/timer.f b/timer.f new file mode 100644 index 0000000..5535ccf --- /dev/null +++ b/timer.f @@ -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 + -- cgit v1.2.3