summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorGard Spreemann <gspreemann@gmail.com>2017-02-02 21:23:01 +0100
committerGard Spreemann <gspreemann@gmail.com>2017-02-02 21:23:01 +0100
commit173802e9f8d98a01019d3f5aff055f1f04479974 (patch)
treec47b4169dee35def4cdeaaad4cdebaaae5c722da
upstream/3.0+dfsg.3: Start git repository for package by importingdfsg/3.0+dfsg.3dfsg/latest
upstream's 3.0 tarball DFSG-cleaned as currently used in published Debian package.
-rw-r--r--License.txt71
-rw-r--r--Makefile37
-rw-r--r--OUTPUTS/output_77_185
-rw-r--r--OUTPUTS/output_77_257
-rw-r--r--OUTPUTS/output_77_3222
-rw-r--r--OUTPUTS/output_90_185
-rw-r--r--OUTPUTS/output_90_257
-rw-r--r--OUTPUTS/output_90_3222
-rw-r--r--README232
-rwxr-xr-xalgorithm.pdfbin0 -> 260255 bytes
-rwxr-xr-xcode.pdfbin0 -> 180639 bytes
-rw-r--r--driver1.f321
-rwxr-xr-xdriver1.f90297
-rw-r--r--driver2.f232
-rwxr-xr-xdriver2.f90220
-rw-r--r--driver3.f273
-rw-r--r--driver3.f90256
-rw-r--r--iterate.dat49
-rw-r--r--lbfgsb.f3949
-rw-r--r--timer.f32
20 files changed, 6697 insertions, 0 deletions
diff --git a/License.txt b/License.txt
new file mode 100644
index 0000000..9ee303a
--- /dev/null
+++ b/License.txt
@@ -0,0 +1,71 @@
+3-clause license ("New BSD License" or "Modified BSD License")
+New BSD License
+Author Regents of the University of California
+Publisher Public Domain
+Published July 22, 1999[8]
+DFSG compatible Yes[7]
+FSF approved Yes[1]
+OSI approved Yes[3]
+GPL compatible Yes[1]
+Copyleft No[1]
+Copyfree Yes
+Linking from code with a different license Yes
+
+The advertising clause was removed from the license text in the official BSD on July 22, 1999 by William Hoskins, Director of
+the Office of Technology Licensing for UC Berkeley.[8] Other BSD distributions removed the clause, but many similar clauses
+remain in BSD-derived code from other sources, and unrelated code using a derived license.
+
+While the original license is sometimes referred to as "BSD-old", the resulting 3-clause version is sometimes referred to by
+"BSD-new." Other names include "New BSD", "revised BSD", "BSD-3", or "3-clause BSD". This version has been vetted as
+an Open source license by the OSI as the "The BSD License".[3] The Free Software Foundation, which refers to the license
+as the "Modified BSD License", states that it is compatible with the GNU GPL. The FSF encourages users to be specific
+when referring to the license by name (i.e. not simply referring to it as "a BSD license" or "BSD-style") to avoid confusion with
+the original BSD license.[1]
+
+This version allows unlimited redistribution for any purpose as long as its copyright notices and the license's disclaimers of
+warranty are maintained. The license also contains a clause restricting use of the names of contributors for endorsement of a
+derived work without specific permission.
+
+Copyright (c) <year>, <copyright holder>
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+ * Neither the name of the <organization> nor the
+ names of its contributors may be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+References
+
+1. ^ a b c d e "Various Licenses and Comments about Them - GNU Project - Free Software Foundation (FSF): Modified BSD license".
+ Free Software Foundation. Retrieved 02 October 2010.
+2. ^ a b c d e "Various Licenses and Comments about Them - GNU Project - Free Software Foundation (FSF): FreeBSD license".
+ Free Software Foundation. Retrieved 02 October 2010.
+3. ^ a b c d e f "Open Source Initiative OSI - The BSD License:Licensing". Open Source Initiative. Retrieved 06 December 2009.
+4. ^ a b c d "Various Licenses and Comments about Them - GNU Project - Free Software Foundation (FSF): Original BSD license".
+ Free Software Foundation. Retrieved 02 October 2010.
+5. ^ The year given is the year 4.3BSD-Tahoe was released. Whether this is the first use of the license is not known.
+6. ^ The year given is the year 4.3BSD-Reno was released. Whether this is the first use of the license is not known.
+7. ^ a b "Debian -- License information". Debian. Retrieved 18 February 2010.
+8. ^ a b c "To All Licensees, Distributors of Any Version of BSD". University of California, Berkeley. 1999-07-22. Retrieved
+ 2006-11-15.
+9. ^ Richard Stallman. "The BSD License Problem". Free Software Foundation. Retrieved 2006-11-15.
+10. ^ "The FreeBSD Copyright". The FreeBSD Project. Retrieved 6 December 2009.
+11. ^ "NetBSD Licensing and Redistribution". The NetBSD Foundation. Retrieved 06 December 2009.
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..7006149
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,37 @@
+FC = gfortran
+
+FFLAGS = -O -Wall -fbounds-check -g -Wno-uninitialized
+
+DRIVER1_77 = driver1.f
+DRIVER2_77 = driver2.f
+DRIVER3_77 = driver3.f
+
+DRIVER1_90 = driver1.f90
+DRIVER2_90 = driver2.f90
+DRIVER3_90 = driver3.f90
+
+LBFGSB = lbfgsb.f
+LINPACK = linpack.f
+BLAS = blas.f
+TIMER = timer.f
+
+all : lbfgsb_77_1 lbfgsb_77_2 lbfgsb_77_3 lbfgsb_90_1 lbfgsb_90_2 lbfgsb_90_3
+
+
+lbfgsb_77_1 : $(DRIVER1_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER)
+ $(FC) $(FFLAGS) $(DRIVER1_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_77_1
+
+lbfgsb_77_2 : $(DRIVER2_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER)
+ $(FC) $(FFLAGS) $(DRIVER2_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_77_2
+
+lbfgsb_77_3 : $(DRIVER3_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER)
+ $(FC) $(FFLAGS) $(DRIVER3_77) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_77_3
+
+lbfgsb_90_1 : $(DRIVER1_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER)
+ $(FC) $(FFLAGS) $(DRIVER1_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_90_1
+
+lbfgsb_90_2 : $(DRIVER2_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER)
+ $(FC) $(FFLAGS) $(DRIVER2_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_90_2
+
+lbfgsb_90_3 : $(DRIVER3_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER)
+ $(FC) $(FFLAGS) $(DRIVER3_90) $(LBFGSB) $(LINPACK) $(BLAS) $(TIMER) -o x.lbfgsb_90_3
diff --git a/OUTPUTS/output_77_1 b/OUTPUTS/output_77_1
new file mode 100644
index 0000000..52dc375
--- /dev/null
+++ b/OUTPUTS/output_77_1
@@ -0,0 +1,85 @@
+
+ Solving sample problem.
+ (f = 0.0 at the optimal solution.)
+
+RUNNING THE L-BFGS-B CODE
+
+ * * *
+
+Machine precision = 2.220D-16
+ N = 25 M = 5
+
+At X0 0 variables are exactly at the bounds
+
+At iterate 0 f= 3.46000D+03 |proj g|= 1.03000D+02
+
+At iterate 1 f= 2.39769D+03 |proj g|= 6.50700D+01
+
+At iterate 2 f= 1.43806D+02 |proj g|= 3.64039D+01
+
+At iterate 3 f= 7.28161D+01 |proj g|= 2.29042D+01
+
+At iterate 4 f= 1.60308D+01 |proj g|= 6.95409D+00
+
+At iterate 5 f= 5.18725D+00 |proj g|= 9.05481D+00
+
+At iterate 6 f= 2.12716D+00 |proj g|= 1.96729D+01
+
+At iterate 7 f= 2.07568D-01 |proj g|= 2.12849D+00
+
+At iterate 8 f= 5.32739D-02 |proj g|= 8.32469D-01
+
+At iterate 9 f= 1.30450D-02 |proj g|= 4.27926D-01
+
+At iterate 10 f= 3.86031D-03 |proj g|= 2.00812D-01
+
+At iterate 11 f= 7.45653D-04 |proj g|= 1.37723D-01
+
+At iterate 12 f= 3.54016D-04 |proj g|= 1.21274D-01
+
+At iterate 13 f= 7.42511D-05 |proj g|= 2.97814D-02
+
+At iterate 14 f= 3.74062D-05 |proj g|= 1.72742D-02
+
+At iterate 15 f= 1.09832D-05 |proj g|= 2.86815D-02
+
+At iterate 16 f= 3.90757D-06 |proj g|= 8.08085D-03
+
+At iterate 17 f= 1.99502D-06 |proj g|= 3.47648D-03
+
+At iterate 18 f= 8.25796D-07 |proj g|= 2.25283D-03
+
+At iterate 19 f= 1.99224D-07 |proj g|= 1.45773D-03
+
+At iterate 20 f= 5.75772D-08 |proj g|= 1.48245D-03
+
+At iterate 21 f= 1.46326D-08 |proj g|= 5.44304D-04
+
+At iterate 22 f= 2.36329D-09 |proj g|= 2.25244D-04
+
+At iterate 23 f= 1.08349D-09 |proj g|= 1.72052D-04
+
+ * * *
+
+Tit = total number of iterations
+Tnf = total number of function evaluations
+Tnint = total number of segments explored during Cauchy searches
+Skip = number of BFGS updates skipped
+Nact = number of active bounds at final generalized Cauchy point
+Projg = norm of the final projected gradient
+F = final function value
+
+ * * *
+
+ N Tit Tnf Tnint Skip Nact Projg F
+ 25 23 28 47 0 0 1.721D-04 1.083D-09
+ F = 1.083490083518441E-009
+
+CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
+
+ Cauchy time 1.000E-03 seconds.
+ Subspace minimization time 0.000E+00 seconds.
+ Line search time 0.000E+00 seconds.
+
+ Total User time 2.000E-03 seconds.
+
diff --git a/OUTPUTS/output_77_2 b/OUTPUTS/output_77_2
new file mode 100644
index 0000000..f12b184
--- /dev/null
+++ b/OUTPUTS/output_77_2
@@ -0,0 +1,57 @@
+
+ Solving sample problem.
+ (f = 0.0 at the optimal solution.)
+
+Iterate 1 nfg = 5 f = 2.39769D+03 |proj g| = 6.50700D+01
+Iterate 2 nfg = 6 f = 1.43806D+02 |proj g| = 3.64039D+01
+Iterate 3 nfg = 7 f = 7.28161D+01 |proj g| = 2.29042D+01
+Iterate 4 nfg = 8 f = 1.60308D+01 |proj g| = 6.95409D+00
+Iterate 5 nfg = 9 f = 5.18725D+00 |proj g| = 9.05481D+00
+Iterate 6 nfg = 10 f = 2.12716D+00 |proj g| = 1.96729D+01
+Iterate 7 nfg = 11 f = 2.07568D-01 |proj g| = 2.12849D+00
+Iterate 8 nfg = 12 f = 5.32739D-02 |proj g| = 8.32469D-01
+Iterate 9 nfg = 13 f = 1.30450D-02 |proj g| = 4.27926D-01
+Iterate 10 nfg = 14 f = 3.86031D-03 |proj g| = 2.00812D-01
+Iterate 11 nfg = 15 f = 7.45653D-04 |proj g| = 1.37723D-01
+Iterate 12 nfg = 16 f = 3.54016D-04 |proj g| = 1.21274D-01
+Iterate 13 nfg = 17 f = 7.42511D-05 |proj g| = 2.97814D-02
+Iterate 14 nfg = 18 f = 3.74062D-05 |proj g| = 1.72742D-02
+Iterate 15 nfg = 19 f = 1.09832D-05 |proj g| = 2.86815D-02
+Iterate 16 nfg = 21 f = 3.90757D-06 |proj g| = 8.08085D-03
+Iterate 17 nfg = 22 f = 1.99502D-06 |proj g| = 3.47648D-03
+Iterate 18 nfg = 23 f = 8.25796D-07 |proj g| = 2.25283D-03
+Iterate 19 nfg = 24 f = 1.99224D-07 |proj g| = 1.45773D-03
+Iterate 20 nfg = 25 f = 5.75772D-08 |proj g| = 1.48245D-03
+Iterate 21 nfg = 26 f = 1.46326D-08 |proj g| = 5.44304D-04
+Iterate 22 nfg = 27 f = 2.36329D-09 |proj g| = 2.25244D-04
+Iterate 23 nfg = 28 f = 1.08349D-09 |proj g| = 1.72052D-04
+Iterate 24 nfg = 29 f = 3.49043D-10 |proj g| = 5.78820D-05
+Iterate 25 nfg = 30 f = 8.43998D-11 |proj g| = 2.04504D-05
+Iterate 26 nfg = 31 f = 2.29812D-11 |proj g| = 1.74854D-05
+Iterate 27 nfg = 32 f = 5.56238D-12 |proj g| = 1.29349D-05
+Iterate 28 nfg = 33 f = 6.17073D-13 |proj g| = 3.52399D-06
+Iterate 29 nfg = 34 f = 2.06920D-13 |proj g| = 9.74653D-07
+Iterate 30 nfg = 35 f = 1.07949D-13 |proj g| = 6.84427D-07
+Iterate 31 nfg = 36 f = 3.27200D-14 |proj g| = 7.46523D-07
+Iterate 32 nfg = 38 f = 2.20380D-14 |proj g| = 6.32502D-07
+Iterate 33 nfg = 39 f = 1.37189D-14 |proj g| = 2.83064D-07
+Iterate 34 nfg = 40 f = 8.06341D-15 |proj g| = 1.25600D-07
+Iterate 35 nfg = 41 f = 6.30621D-15 |proj g| = 6.63279D-08
+Iterate 36 nfg = 42 f = 6.12906D-15 |proj g| = 9.93648D-08
+Iterate 37 nfg = 43 f = 5.83506D-15 |proj g| = 1.10012D-08
+Iterate 38 nfg = 44 f = 5.82125D-15 |proj g| = 8.09738D-09
+Iterate 39 nfg = 45 f = 5.81201D-15 |proj g| = 2.04600D-08
+Iterate 40 nfg = 47 f = 5.80895D-15 |proj g| = 4.71360D-09
+Iterate 41 nfg = 48 f = 5.80786D-15 |proj g| = 1.97091D-09
+Iterate 42 nfg = 49 f = 5.80716D-15 |proj g| = 1.10044D-09
+Iterate 43 nfg = 50 f = 5.80705D-15 |proj g| = 3.49651D-10
+Iterate 44 nfg = 51 f = 5.80703D-15 |proj g| = 2.29788D-10
+Iterate 45 nfg = 52 f = 5.80702D-15 |proj g| = 1.51178D-10
+Iterate 46 nfg = 53 f = 5.80702D-15 |proj g| = 6.61970D-11
+ STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL
+ Final X=
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0001D+00 1.0002D+00
+ 1.0003D+00 1.0006D+00 1.0013D+00 1.0026D+00 1.0052D+00 1.0105D+00
+ 1.0210D+00 1.0425D+00 1.0867D+00 1.1810D+00 1.3947D+00 1.9452D+00
+ 3.7837D+00
diff --git a/OUTPUTS/output_77_3 b/OUTPUTS/output_77_3
new file mode 100644
index 0000000..612f3d3
--- /dev/null
+++ b/OUTPUTS/output_77_3
@@ -0,0 +1,222 @@
+
+ Solving sample problem.
+ (f = 0.0 at the optimal solution.)
+
+Iterate 1 nfg = 5 f = 9.98226D+04 |proj g| = 4.92220D+01
+Iterate 2 nfg = 6 f = 5.97672D+03 |proj g| = 7.29734D+00
+Iterate 3 nfg = 7 f = 3.02885D+03 |proj g| = 7.40956D+00
+Iterate 4 nfg = 8 f = 7.16214D+02 |proj g| = 6.00527D+00
+Iterate 5 nfg = 9 f = 2.20946D+02 |proj g| = 4.08022D+00
+Iterate 6 nfg = 10 f = 5.61473D+01 |proj g| = 5.17435D+00
+Iterate 7 nfg = 11 f = 3.39145D+01 |proj g| = 2.09576D+01
+Iterate 8 nfg = 12 f = 1.36320D+00 |proj g| = 9.93231D+00
+Iterate 9 nfg = 13 f = 3.50771D-01 |proj g| = 5.91099D+00
+Iterate 10 nfg = 14 f = 1.52466D-02 |proj g| = 6.13600D-01
+Iterate 11 nfg = 15 f = 7.36312D-03 |proj g| = 3.39884D-01
+Iterate 12 nfg = 16 f = 1.25124D-03 |proj g| = 1.17803D-01
+Iterate 13 nfg = 17 f = 4.68586D-04 |proj g| = 7.22477D-02
+Iterate 14 nfg = 18 f = 5.50022D-05 |proj g| = 4.40414D-02
+Iterate 15 nfg = 19 f = 1.82587D-05 |proj g| = 1.49985D-02
+Iterate 16 nfg = 20 f = 6.79533D-06 |proj g| = 7.05805D-03
+Iterate 17 nfg = 21 f = 1.59803D-06 |proj g| = 3.57847D-03
+Iterate 18 nfg = 22 f = 4.67168D-07 |proj g| = 2.96329D-03
+Iterate 19 nfg = 24 f = 1.61224D-07 |proj g| = 1.61125D-03
+Iterate 20 nfg = 25 f = 5.27054D-08 |proj g| = 6.79899D-04
+Iterate 21 nfg = 26 f = 2.21706D-08 |proj g| = 3.98553D-04
+Iterate 22 nfg = 27 f = 7.85699D-09 |proj g| = 8.29301D-04
+Iterate 23 nfg = 28 f = 2.57367D-09 |proj g| = 3.95843D-04
+Iterate 24 nfg = 29 f = 9.39513D-10 |proj g| = 8.04745D-05
+Iterate 25 nfg = 30 f = 5.42208D-10 |proj g| = 5.75104D-05
+Iterate 26 nfg = 31 f = 1.13904D-10 |proj g| = 3.88455D-05
+Iterate 27 nfg = 32 f = 2.19270D-11 |proj g| = 2.48929D-05
+Iterate 28 nfg = 34 f = 5.89247D-12 |proj g| = 5.98559D-06
+Iterate 29 nfg = 35 f = 3.05886D-12 |proj g| = 3.88360D-06
+Iterate 30 nfg = 36 f = 6.64922D-13 |proj g| = 2.57655D-06
+Iterate 31 nfg = 37 f = 1.82984D-13 |proj g| = 1.54163D-06
+Iterate 32 nfg = 38 f = 5.06020D-14 |proj g| = 1.46067D-06
+Iterate 33 nfg = 39 f = 2.03549D-14 |proj g| = 2.88905D-07
+Iterate 34 nfg = 40 f = 1.27697D-14 |proj g| = 2.09331D-07
+Iterate 35 nfg = 41 f = 2.59126D-15 |proj g| = 2.75406D-07
+Iterate 36 nfg = 43 f = 1.13639D-15 |proj g| = 1.13897D-07
+Iterate 37 nfg = 44 f = 5.29416D-16 |proj g| = 6.22050D-08
+Iterate 38 nfg = 45 f = 1.24096D-16 |proj g| = 2.30404D-08
+Iterate 39 nfg = 46 f = 3.14643D-17 |proj g| = 1.56779D-08
+Iterate 40 nfg = 47 f = 1.56638D-17 |proj g| = 1.19498D-08
+Iterate 41 nfg = 48 f = 2.80290D-18 |proj g| = 4.46387D-09
+Iterate 42 nfg = 49 f = 1.32837D-18 |proj g| = 2.55881D-09
+Iterate 43 nfg = 50 f = 2.54553D-19 |proj g| = 1.17505D-09
+Iterate 44 nfg = 52 f = 8.23941D-20 |proj g| = 9.01753D-10
+Iterate 45 nfg = 53 f = 2.79611D-20 |proj g| = 8.62939D-10
+Iterate 46 nfg = 54 f = 1.08724D-20 |proj g| = 2.44455D-10
+Iterate 47 nfg = 55 f = 5.32479D-21 |proj g| = 1.66311D-10
+Iterate 48 nfg = 56 f = 1.26618D-21 |proj g| = 1.90418D-10
+Iterate 49 nfg = 58 f = 5.35057D-22 |proj g| = 9.73814D-11
+ STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL
+ Final X=
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0001D+00 1.0002D+00 1.0003D+00 1.0007D+00 1.0014D+00
+ 1.0028D+00 1.0055D+00 1.0111D+00 1.0223D+00 1.0451D+00 1.0923D+00
+ 1.1930D+00 1.4233D+00 2.0257D+00 4.1035D+00
diff --git a/OUTPUTS/output_90_1 b/OUTPUTS/output_90_1
new file mode 100644
index 0000000..008c044
--- /dev/null
+++ b/OUTPUTS/output_90_1
@@ -0,0 +1,85 @@
+
+ Solving sample problem.
+ (f = 0.0 at the optimal solution.)
+
+RUNNING THE L-BFGS-B CODE
+
+ * * *
+
+Machine precision = 2.220D-16
+ N = 25 M = 5
+
+At X0 0 variables are exactly at the bounds
+
+At iterate 0 f= 3.46000D+03 |proj g|= 1.03000D+02
+
+At iterate 1 f= 2.39769D+03 |proj g|= 6.50700D+01
+
+At iterate 2 f= 1.43806D+02 |proj g|= 3.64039D+01
+
+At iterate 3 f= 7.28161D+01 |proj g|= 2.29042D+01
+
+At iterate 4 f= 1.60308D+01 |proj g|= 6.95409D+00
+
+At iterate 5 f= 5.18725D+00 |proj g|= 9.05481D+00
+
+At iterate 6 f= 2.12716D+00 |proj g|= 1.96729D+01
+
+At iterate 7 f= 2.07568D-01 |proj g|= 2.12849D+00
+
+At iterate 8 f= 5.32739D-02 |proj g|= 8.32469D-01
+
+At iterate 9 f= 1.30450D-02 |proj g|= 4.27926D-01
+
+At iterate 10 f= 3.86031D-03 |proj g|= 2.00812D-01
+
+At iterate 11 f= 7.45653D-04 |proj g|= 1.37723D-01
+
+At iterate 12 f= 3.54016D-04 |proj g|= 1.21274D-01
+
+At iterate 13 f= 7.42511D-05 |proj g|= 2.97814D-02
+
+At iterate 14 f= 3.74062D-05 |proj g|= 1.72742D-02
+
+At iterate 15 f= 1.09832D-05 |proj g|= 2.86815D-02
+
+At iterate 16 f= 3.90757D-06 |proj g|= 8.08085D-03
+
+At iterate 17 f= 1.99502D-06 |proj g|= 3.47648D-03
+
+At iterate 18 f= 8.25796D-07 |proj g|= 2.25283D-03
+
+At iterate 19 f= 1.99224D-07 |proj g|= 1.45773D-03
+
+At iterate 20 f= 5.75772D-08 |proj g|= 1.48245D-03
+
+At iterate 21 f= 1.46326D-08 |proj g|= 5.44304D-04
+
+At iterate 22 f= 2.36329D-09 |proj g|= 2.25244D-04
+
+At iterate 23 f= 1.08349D-09 |proj g|= 1.72052D-04
+
+ * * *
+
+Tit = total number of iterations
+Tnf = total number of function evaluations
+Tnint = total number of segments explored during Cauchy searches
+Skip = number of BFGS updates skipped
+Nact = number of active bounds at final generalized Cauchy point
+Projg = norm of the final projected gradient
+F = final function value
+
+ * * *
+
+ N Tit Tnf Tnint Skip Nact Projg F
+ 25 23 28 47 0 0 1.721D-04 1.083D-09
+ F = 1.083490083461424E-009
+
+CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
+
+ Cauchy time 0.000E+00 seconds.
+ Subspace minimization time 0.000E+00 seconds.
+ Line search time 0.000E+00 seconds.
+
+ Total User time 1.000E-03 seconds.
+
diff --git a/OUTPUTS/output_90_2 b/OUTPUTS/output_90_2
new file mode 100644
index 0000000..159875c
--- /dev/null
+++ b/OUTPUTS/output_90_2
@@ -0,0 +1,57 @@
+
+ Solving sample problem.
+ (f = 0.0 at the optimal solution.)
+
+Iterate 1 nfg = 5 f = 2.39769D+03 |proj g| = 6.50700D+01
+Iterate 2 nfg = 6 f = 1.43806D+02 |proj g| = 3.64039D+01
+Iterate 3 nfg = 7 f = 7.28161D+01 |proj g| = 2.29042D+01
+Iterate 4 nfg = 8 f = 1.60308D+01 |proj g| = 6.95409D+00
+Iterate 5 nfg = 9 f = 5.18725D+00 |proj g| = 9.05481D+00
+Iterate 6 nfg = 10 f = 2.12716D+00 |proj g| = 1.96729D+01
+Iterate 7 nfg = 11 f = 2.07568D-01 |proj g| = 2.12849D+00
+Iterate 8 nfg = 12 f = 5.32739D-02 |proj g| = 8.32469D-01
+Iterate 9 nfg = 13 f = 1.30450D-02 |proj g| = 4.27926D-01
+Iterate 10 nfg = 14 f = 3.86031D-03 |proj g| = 2.00812D-01
+Iterate 11 nfg = 15 f = 7.45653D-04 |proj g| = 1.37723D-01
+Iterate 12 nfg = 16 f = 3.54016D-04 |proj g| = 1.21274D-01
+Iterate 13 nfg = 17 f = 7.42511D-05 |proj g| = 2.97814D-02
+Iterate 14 nfg = 18 f = 3.74062D-05 |proj g| = 1.72742D-02
+Iterate 15 nfg = 19 f = 1.09832D-05 |proj g| = 2.86815D-02
+Iterate 16 nfg = 21 f = 3.90757D-06 |proj g| = 8.08085D-03
+Iterate 17 nfg = 22 f = 1.99502D-06 |proj g| = 3.47648D-03
+Iterate 18 nfg = 23 f = 8.25796D-07 |proj g| = 2.25283D-03
+Iterate 19 nfg = 24 f = 1.99224D-07 |proj g| = 1.45773D-03
+Iterate 20 nfg = 25 f = 5.75772D-08 |proj g| = 1.48245D-03
+Iterate 21 nfg = 26 f = 1.46326D-08 |proj g| = 5.44304D-04
+Iterate 22 nfg = 27 f = 2.36329D-09 |proj g| = 2.25244D-04
+Iterate 23 nfg = 28 f = 1.08349D-09 |proj g| = 1.72052D-04
+Iterate 24 nfg = 29 f = 3.49043D-10 |proj g| = 5.78820D-05
+Iterate 25 nfg = 30 f = 8.43998D-11 |proj g| = 2.04504D-05
+Iterate 26 nfg = 31 f = 2.29812D-11 |proj g| = 1.74854D-05
+Iterate 27 nfg = 32 f = 5.56238D-12 |proj g| = 1.29349D-05
+Iterate 28 nfg = 33 f = 6.17073D-13 |proj g| = 3.52399D-06
+Iterate 29 nfg = 34 f = 2.06920D-13 |proj g| = 9.74653D-07
+Iterate 30 nfg = 35 f = 1.07949D-13 |proj g| = 6.84427D-07
+Iterate 31 nfg = 36 f = 3.27200D-14 |proj g| = 7.46523D-07
+Iterate 32 nfg = 38 f = 2.20380D-14 |proj g| = 6.32502D-07
+Iterate 33 nfg = 39 f = 1.37189D-14 |proj g| = 2.83064D-07
+Iterate 34 nfg = 40 f = 8.06341D-15 |proj g| = 1.25600D-07
+Iterate 35 nfg = 41 f = 6.30621D-15 |proj g| = 6.63280D-08
+Iterate 36 nfg = 42 f = 6.12906D-15 |proj g| = 9.93648D-08
+Iterate 37 nfg = 43 f = 5.83506D-15 |proj g| = 1.10011D-08
+Iterate 38 nfg = 44 f = 5.82125D-15 |proj g| = 8.09738D-09
+Iterate 39 nfg = 45 f = 5.81201D-15 |proj g| = 2.04599D-08
+Iterate 40 nfg = 47 f = 5.80895D-15 |proj g| = 4.71361D-09
+Iterate 41 nfg = 48 f = 5.80786D-15 |proj g| = 1.97091D-09
+Iterate 42 nfg = 49 f = 5.80716D-15 |proj g| = 1.10044D-09
+Iterate 43 nfg = 50 f = 5.80705D-15 |proj g| = 3.49639D-10
+Iterate 44 nfg = 51 f = 5.80703D-15 |proj g| = 2.29804D-10
+Iterate 45 nfg = 52 f = 5.80702D-15 |proj g| = 1.51167D-10
+Iterate 46 nfg = 53 f = 5.80702D-15 |proj g| = 6.62041D-11
+ STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL
+ Final X=
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0001D+00 1.0002D+00
+ 1.0003D+00 1.0006D+00 1.0013D+00 1.0026D+00 1.0052D+00 1.0105D+00
+ 1.0210D+00 1.0425D+00 1.0867D+00 1.1810D+00 1.3947D+00 1.9452D+00
+ 3.7837D+00
diff --git a/OUTPUTS/output_90_3 b/OUTPUTS/output_90_3
new file mode 100644
index 0000000..24e13da
--- /dev/null
+++ b/OUTPUTS/output_90_3
@@ -0,0 +1,222 @@
+
+ Solving sample problem.
+ (f = 0.0 at the optimal solution.)
+
+Iterate 1 nfg = 5 f = 9.98226D+04 |proj g| = 4.92220D+01
+Iterate 2 nfg = 6 f = 5.97672D+03 |proj g| = 7.29734D+00
+Iterate 3 nfg = 7 f = 3.02885D+03 |proj g| = 7.40956D+00
+Iterate 4 nfg = 8 f = 7.16214D+02 |proj g| = 6.00527D+00
+Iterate 5 nfg = 9 f = 2.20946D+02 |proj g| = 4.08022D+00
+Iterate 6 nfg = 10 f = 5.61473D+01 |proj g| = 5.17435D+00
+Iterate 7 nfg = 11 f = 3.39145D+01 |proj g| = 2.09576D+01
+Iterate 8 nfg = 12 f = 1.36320D+00 |proj g| = 9.93231D+00
+Iterate 9 nfg = 13 f = 3.50771D-01 |proj g| = 5.91099D+00
+Iterate 10 nfg = 14 f = 1.52466D-02 |proj g| = 6.13600D-01
+Iterate 11 nfg = 15 f = 7.36312D-03 |proj g| = 3.39884D-01
+Iterate 12 nfg = 16 f = 1.25124D-03 |proj g| = 1.17803D-01
+Iterate 13 nfg = 17 f = 4.68586D-04 |proj g| = 7.22477D-02
+Iterate 14 nfg = 18 f = 5.50022D-05 |proj g| = 4.40414D-02
+Iterate 15 nfg = 19 f = 1.82587D-05 |proj g| = 1.49985D-02
+Iterate 16 nfg = 20 f = 6.79533D-06 |proj g| = 7.05805D-03
+Iterate 17 nfg = 21 f = 1.59803D-06 |proj g| = 3.57847D-03
+Iterate 18 nfg = 22 f = 4.67168D-07 |proj g| = 2.96329D-03
+Iterate 19 nfg = 24 f = 1.61224D-07 |proj g| = 1.61125D-03
+Iterate 20 nfg = 25 f = 5.27054D-08 |proj g| = 6.79899D-04
+Iterate 21 nfg = 26 f = 2.21706D-08 |proj g| = 3.98553D-04
+Iterate 22 nfg = 27 f = 7.85699D-09 |proj g| = 8.29301D-04
+Iterate 23 nfg = 28 f = 2.57367D-09 |proj g| = 3.95843D-04
+Iterate 24 nfg = 29 f = 9.39513D-10 |proj g| = 8.04745D-05
+Iterate 25 nfg = 30 f = 5.42208D-10 |proj g| = 5.75104D-05
+Iterate 26 nfg = 31 f = 1.13904D-10 |proj g| = 3.88455D-05
+Iterate 27 nfg = 32 f = 2.19270D-11 |proj g| = 2.48929D-05
+Iterate 28 nfg = 34 f = 5.89247D-12 |proj g| = 5.98559D-06
+Iterate 29 nfg = 35 f = 3.05886D-12 |proj g| = 3.88360D-06
+Iterate 30 nfg = 36 f = 6.64922D-13 |proj g| = 2.57655D-06
+Iterate 31 nfg = 37 f = 1.82984D-13 |proj g| = 1.54163D-06
+Iterate 32 nfg = 38 f = 5.06020D-14 |proj g| = 1.46067D-06
+Iterate 33 nfg = 39 f = 2.03549D-14 |proj g| = 2.88905D-07
+Iterate 34 nfg = 40 f = 1.27697D-14 |proj g| = 2.09331D-07
+Iterate 35 nfg = 41 f = 2.59125D-15 |proj g| = 2.75406D-07
+Iterate 36 nfg = 43 f = 1.13639D-15 |proj g| = 1.13897D-07
+Iterate 37 nfg = 44 f = 5.29416D-16 |proj g| = 6.22050D-08
+Iterate 38 nfg = 45 f = 1.24096D-16 |proj g| = 2.30404D-08
+Iterate 39 nfg = 46 f = 3.14643D-17 |proj g| = 1.56779D-08
+Iterate 40 nfg = 47 f = 1.56638D-17 |proj g| = 1.19498D-08
+Iterate 41 nfg = 48 f = 2.80291D-18 |proj g| = 4.46388D-09
+Iterate 42 nfg = 49 f = 1.32837D-18 |proj g| = 2.55882D-09
+Iterate 43 nfg = 50 f = 2.54555D-19 |proj g| = 1.17507D-09
+Iterate 44 nfg = 52 f = 8.23938D-20 |proj g| = 9.01720D-10
+Iterate 45 nfg = 53 f = 2.79609D-20 |proj g| = 8.62899D-10
+Iterate 46 nfg = 54 f = 1.08725D-20 |proj g| = 2.44452D-10
+Iterate 47 nfg = 55 f = 5.32476D-21 |proj g| = 1.66320D-10
+Iterate 48 nfg = 56 f = 1.26624D-21 |proj g| = 1.90396D-10
+Iterate 49 nfg = 58 f = 5.35121D-22 |proj g| = 9.74083D-11
+ STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL
+ Final X=
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00 1.0000D+00
+ 1.0000D+00 1.0001D+00 1.0002D+00 1.0003D+00 1.0007D+00 1.0014D+00
+ 1.0028D+00 1.0055D+00 1.0111D+00 1.0223D+00 1.0451D+00 1.0923D+00
+ 1.1930D+00 1.4233D+00 2.0257D+00 4.1035D+00
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 <nocedal@eecs.northwestern.edu>
+ Jose Luis Morales <jmorales@itam.mx>
+________________________________________________________________________
+
+
+ How to use L-BFGS-B
+
+************************************************************************
+
+ The simplest way to use the code is to modify one of the
+drivers provided in the package. Most users will only need to make
+a few changes to the drivers to run their applications.
+
+ L-BFGS-B is written in FORTRAN 77, in double precision. The
+user is required to calculate the function value f and its gradient g.
+In order to allow the user complete control over these computations,
+reverse communication is used. The routine setulb.f must be called
+repeatedly under the control of the variable task. The calling
+statement of L-BFGS-B is
+
+ call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint,
+ + csave,lsave,isave,dsave)
+
+
+ Following is a description of all the parameters used in this call.
+
+c n is an INTEGER variable that must be set by the user to the
+c number of variables. It is not altered by the routine.
+c
+c m is an INTEGER variable that must be set by the user to the
+c number of corrections used in the limited memory matrix.
+c It is not altered by the routine. Values of m < 3 are
+c not recommended, and large values of m can result in excessive
+c computing time. The range 3 <= m <= 20 is recommended.
+c
+c x is a DOUBLE PRECISION array of length n. On initial entry
+c it must be set by the user to the values of the initial
+c estimate of the solution vector. Upon successful exit, it
+c contains the values of the variables at the best point
+c found (usually an approximate solution).
+c
+c l is a DOUBLE PRECISION array of length n that must be set by
+c the user to the values of the lower bounds on the variables. If
+c the i-th variable has no lower bound, l(i) need not be defined.
+c
+c u is a DOUBLE PRECISION array of length n that must be set by
+c the user to the values of the upper bounds on the variables. If
+c the i-th variable has no upper bound, u(i) need not be defined.
+c
+c nbd is an INTEGER array of dimension n that must be set by the
+c user to the type of bounds imposed on the variables:
+c nbd(i)=0 if x(i) is unbounded,
+c 1 if x(i) has only a lower bound,
+c 2 if x(i) has both lower and upper bounds,
+c 3 if x(i) has only an upper bound.
+c
+c f is a DOUBLE PRECISION variable. If the routine setulb returns
+c with task(1:2)= 'FG', then f must be set by the user to
+c contain the value of the function at the point x.
+c
+c g is a DOUBLE PRECISION array of length n. If the routine setulb
+c returns with taskb(1:2)= 'FG', then g must be set by the user to
+c contain the components of the gradient at the point x.
+c
+c factr is a DOUBLE PRECISION variable that must be set by the user.
+c It is a tolerance in the termination test for the algorithm.
+c The iteration will stop when
+c
+c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
+c
+c where epsmch is the machine precision which is automatically
+c generated by the code. Typical values for factr on a computer
+c with 15 digits of accuracy in double precision are:
+c factr=1.d+12 for low accuracy;
+c 1.d+7 for moderate accuracy;
+c 1.d+1 for extremely high accuracy.
+c The user can suppress this termination test by setting factr=0.
+c
+c pgtol is a double precision variable.
+c On entry pgtol >= 0 is specified by the user. The iteration
+c will stop when
+c
+c max{|proj g_i | i = 1, ..., n} <= pgtol
+c
+c where pg_i is the ith component of the projected gradient.
+c The user can suppress this termination test by setting pgtol=0.
+c
+c wa is a DOUBLE PRECISION array of length
+c (2mmax + 5)nmax + 11mmax^2 + 8mmax used as workspace.
+c This array must not be altered by the user.
+c
+c iwa is an INTEGER array of length 3nmax used as
+c workspace. This array must not be altered by the user.
+c
+c task is a CHARACTER string of length 60.
+c On first entry, it must be set to 'START'.
+c On a return with task(1:2)='FG', the user must evaluate the
+c function f and gradient g at the returned value of x.
+c On a return with task(1:5)='NEW_X', an iteration of the
+c algorithm has concluded, and f and g contain f(x) and g(x)
+c respectively. The user can decide whether to continue or stop
+c the iteration.
+c When
+c task(1:4)='CONV', the termination test in L-BFGS-B has been
+c satisfied;
+c task(1:4)='ABNO', the routine has terminated abnormally
+c without being able to satisfy the termination conditions,
+c x contains the best approximation found,
+c f and g contain f(x) and g(x) respectively;
+c task(1:5)='ERROR', the routine has detected an error in the
+c input parameters;
+c On exit with task = 'CONV', 'ABNO' or 'ERROR', the variable task
+c contains additional information that the user can print.
+c This array should not be altered unless the user wants to
+c stop the run for some reason. See driver2 or driver3
+c for a detailed explanation on how to stop the run
+c by assigning task(1:4)='STOP' in the driver.
+c
+c iprint is an INTEGER variable that must be set by the user.
+c It controls the frequency and type of output generated:
+c iprint<0 no output is generated;
+c iprint=0 print only one line at the last iteration;
+c 0<iprint<99 print also f and |proj g| every iprint iterations;
+c iprint=99 print details of every iteration except n-vectors;
+c iprint=100 print also the changes of active set and final x;
+c iprint>100 print details of every iteration including x and g;
+c When iprint > 0, the file iterate.dat will be created to
+c summarize the iteration.
+c
+c csave is a CHARACTER working array of length 60.
+c
+c lsave is a LOGICAL working array of dimension 4.
+c On exit with task = 'NEW_X', the following information is
+c available:
+c lsave(1) = .true. the initial x did not satisfy the bounds;
+c lsave(2) = .true. the problem contains bounds;
+c lsave(3) = .true. each variable has upper and lower bounds.
+c
+c isave is an INTEGER working array of dimension 44.
+c On exit with task = 'NEW_X', it contains information that
+c the user may want to access:
+c isave(30) = the current iteration number;
+c isave(34) = the total number of function and gradient
+c evaluations;
+c isave(36) = the number of function value or gradient
+c evaluations in the current iteration;
+c isave(38) = the number of free variables in the current
+c iteration;
+c isave(39) = the number of active constraints at the current
+c iteration;
+c
+c see the subroutine setulb.f for a description of other
+c information contained in isave
+c
+c dsave is a DOUBLE PRECISION working array of dimension 29.
+c On exit with task = 'NEW_X', it contains information that
+c the user may want to access:
+c dsave(2) = the value of f at the previous iteration;
+c dsave(5) = the machine precision epsmch generated by the code;
+c dsave(13) = the infinity norm of the projected gradient;
+c
+c see the subroutine setulb.f for a description of other
+c information contained in dsave
+c
+************************************************************************
+
diff --git a/algorithm.pdf b/algorithm.pdf
new file mode 100755
index 0000000..bfc9af6
--- /dev/null
+++ b/algorithm.pdf
Binary files differ
diff --git a/code.pdf b/code.pdf
new file mode 100755
index 0000000..b257b27
--- /dev/null
+++ b/code.pdf
Binary files 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 0<iprint<99 print also f and |proj g| every iprint iterations;
+c iprint=99 print details of every iteration except n-vectors;
+c iprint=100 print also the changes of active set and final x;
+c iprint>100 print details of every iteration including x and g;
+c When iprint > 0, the file iterate.dat will be created to
+c summarize the iteration.
+c
+c csave is a CHARACTER working array of length 60.
+c
+c lsave is a LOGICAL working array of dimension 4.
+c On exit with task = 'NEW_X', the following information is
+c available:
+c lsave(1) = .true. the initial x did not satisfy the bounds;
+c lsave(2) = .true. the problem contains bounds;
+c lsave(3) = .true. each variable has upper and lower bounds.
+c
+c isave is an INTEGER working array of dimension 44.
+c On exit with task = 'NEW_X', it contains information that
+c the user may want to access:
+c isave(30) = the current iteration number;
+c isave(34) = the total number of function and gradient
+c evaluations;
+c isave(36) = the number of function value or gradient
+c evaluations in the current iteration;
+c isave(38) = the number of free variables in the current
+c iteration;
+c isave(39) = the number of active constraints at the current
+c iteration;
+c
+c see the subroutine setulb.f for a description of other
+c information contained in isave
+c
+c dsave is a DOUBLE PRECISION working array of dimension 29.
+c On exit with task = 'NEW_X', it contains information that
+c the user may want to access:
+c dsave(2) = the value of f at the previous iteration;
+c dsave(5) = the machine precision epsmch generated by the code;
+c dsave(13) = the infinity norm of the projected gradient;
+c
+c see the subroutine setulb.f for a description of other
+c information contained in dsave
+c
+c --------------------------------------------------------------
+c END OF THE DESCRIPTION OF THE VARIABLES IN L-BFGS-B
+c --------------------------------------------------------------
+
+ program driver
+
+c This simple driver demonstrates how to call the L-BFGS-B code to
+c solve a sample problem (the extended Rosenbrock function
+c subject to bounds on the variables). The dimension n of this
+c problem is variable.
+
+ integer nmax, mmax
+ parameter (nmax=1024, mmax=17)
+c nmax is the dimension of the largest problem to be solved.
+c mmax is the maximum number of limited memory corrections.
+
+c Declare the variables needed by the code.
+c A description of all these variables is given at the end of
+c the driver.
+
+ character*60 task, csave
+ logical lsave(4)
+ integer n, m, iprint,
+ + nbd(nmax), iwa(3*nmax), isave(44)
+ double precision f, factr, pgtol,
+ + x(nmax), l(nmax), u(nmax), g(nmax), dsave(29),
+ + wa(2*mmax*nmax + 5*nmax + 11*mmax*mmax + 8*mmax)
+
+c Declare a few additional variables for this sample problem.
+
+ double precision t1, t2
+ integer i
+
+c We wish to have output at every iteration.
+
+ iprint = 1
+
+c We specify the tolerances in the stopping criteria.
+
+ factr=1.0d+7
+ pgtol=1.0d-5
+
+c We specify the dimension n of the sample problem and the number
+c m of limited memory corrections stored. (n and m should not
+c exceed the limits nmax and mmax respectively.)
+
+ n=25
+ m=5
+
+c We now provide nbd which defines the bounds on the variables:
+c l specifies the lower bounds,
+c u specifies the upper bounds.
+
+c First set bounds on the odd-numbered variables.
+
+ do 10 i=1,n,2
+ nbd(i)=2
+ l(i)=1.0d0
+ u(i)=1.0d2
+ 10 continue
+
+c Next set bounds on the even-numbered variables.
+
+ do 12 i=2,n,2
+ nbd(i)=2
+ l(i)=-1.0d2
+ u(i)=1.0d2
+ 12 continue
+
+c We now define the starting point.
+
+ do 14 i=1,n
+ x(i)=3.0d0
+ 14 continue
+
+ write (6,16)
+ 16 format(/,5x, 'Solving sample problem.',
+ + /,5x, ' (f = 0.0 at the optimal solution.)',/)
+
+c We start the iteration by initializing task.
+c
+ task = 'START'
+
+c ------- the beginning of the loop ----------
+
+ 111 continue
+
+c This is the call to the L-BFGS-B code.
+
+ call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint,
+ + csave,lsave,isave,dsave)
+
+ if (task(1:2) .eq. 'FG') then
+c the minimization routine has returned to request the
+c function f and gradient g values at the current x.
+
+c Compute function value f for the sample problem.
+
+ f=.25d0*(x(1)-1.d0)**2
+ do 20 i=2,n
+ f=f+(x(i)-x(i-1)**2)**2
+ 20 continue
+ f=4.d0*f
+
+c Compute gradient g for the sample problem.
+
+ t1=x(2)-x(1)**2
+ g(1)=2.d0*(x(1)-1.d0)-1.6d1*x(1)*t1
+ do 22 i=2,n-1
+ t2=t1
+ t1=x(i+1)-x(i)**2
+ g(i)=8.d0*t2-1.6d1*x(i)*t1
+ 22 continue
+ g(n)=8.d0*t1
+
+c go back to the minimization routine.
+ goto 111
+ endif
+c
+ if (task(1:5) .eq. 'NEW_X') goto 111
+c the minimization routine has returned with a new iterate,
+c and we have opted to continue the iteration.
+
+c ---------- the end of the loop -------------
+
+c If task is neither FG nor NEW_X we terminate execution.
+
+ stop
+
+ end
+
+c======================= The end of driver1 ============================
diff --git a/driver1.f90 b/driver1.f90
new file mode 100755
index 0000000..67fa540
--- /dev/null
+++ b/driver1.f90
@@ -0,0 +1,297 @@
+!
+! L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License”
+! or “3-clause license”)
+! Please read attached file License.txt
+!
+!
+! DRIVER1 in Fortran 90
+! --------------------------------------------------------------
+!
+! L-BFGS-B is a code for solving large nonlinear optimization
+! problems with simple bounds on the variables.
+!
+! The code can also be used for unconstrained problems and is
+! as efficient for these problems as the earlier limited memory
+! code L-BFGS.
+!
+! This is the simplest driver in the package. It uses all the
+! default settings of the code.
+!
+!
+! References:
+!
+! [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+! memory algorithm for bound constrained optimization'',
+! SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+!
+! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
+! Subroutines for Large Scale Bound Constrained Optimization''
+! Tech. Report, NAM-11, EECS Department, Northwestern University,
+! 1994.
+!
+!
+! (Postscript files of these papers are available via anonymous
+! ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
+!
+! * * *
+!
+! March 2011 (latest revision)
+! Optimization Center at Northwestern University
+! Instituto Tecnologico Autonomo de Mexico
+!
+! Jorge Nocedal and Jose Luis Morales
+!
+! --------------------------------------------------------------
+! DESCRIPTION OF THE VARIABLES IN L-BFGS-B
+! --------------------------------------------------------------
+!
+! n is an INTEGER variable that must be set by the user to the
+! number of variables. It is not altered by the routine.
+!
+! m is an INTEGER variable that must be set by the user to the
+! number of corrections used in the limited memory matrix.
+! It is not altered by the routine. Values of m < 3 are
+! not recommended, and large values of m can result in excessive
+! computing time. The range 3 <= m <= 20 is recommended.
+!
+! x is a DOUBLE PRECISION array of length n. On initial entry
+! it must be set by the user to the values of the initial
+! estimate of the solution vector. Upon successful exit, it
+! contains the values of the variables at the best point
+! found (usually an approximate solution).
+!
+! l is a DOUBLE PRECISION array of length n that must be set by
+! the user to the values of the lower bounds on the variables. If
+! the i-th variable has no lower bound, l(i) need not be defined.
+!
+! u is a DOUBLE PRECISION array of length n that must be set by
+! the user to the values of the upper bounds on the variables. If
+! the i-th variable has no upper bound, u(i) need not be defined.
+!
+! nbd is an INTEGER array of dimension n that must be set by the
+! user to the type of bounds imposed on the variables:
+! nbd(i)=0 if x(i) is unbounded,
+! 1 if x(i) has only a lower bound,
+! 2 if x(i) has both lower and upper bounds,
+! 3 if x(i) has only an upper bound.
+!
+! f is a DOUBLE PRECISION variable. If the routine setulb returns
+! with task(1:2)= 'FG', then f must be set by the user to
+! contain the value of the function at the point x.
+!
+! g is a DOUBLE PRECISION array of length n. If the routine setulb
+! returns with taskb(1:2)= 'FG', then g must be set by the user to
+! contain the components of the gradient at the point x.
+!
+! factr is a DOUBLE PRECISION variable that must be set by the user.
+! It is a tolerance in the termination test for the algorithm.
+! The iteration will stop when
+!
+! (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
+!
+! where epsmch is the machine precision which is automatically
+! generated by the code. Typical values for factr on a computer
+! with 15 digits of accuracy in double precision are:
+! factr=1.d+12 for low accuracy;
+! 1.d+7 for moderate accuracy;
+! 1.d+1 for extremely high accuracy.
+! The user can suppress this termination test by setting factr=0.
+!
+! pgtol is a double precision variable.
+! On entry pgtol >= 0 is specified by the user. The iteration
+! will stop when
+!
+! max{|proj g_i | i = 1, ..., n} <= pgtol
+!
+! where pg_i is the ith component of the projected gradient.
+! The user can suppress this termination test by setting pgtol=0.
+!
+! wa is a DOUBLE PRECISION array of length
+! (2mmax + 5)nmax + 11mmax^2 + 8mmax used as workspace.
+! This array must not be altered by the user.
+!
+! iwa is an INTEGER array of length 3nmax used as
+! workspace. This array must not be altered by the user.
+!
+! task is a CHARACTER string of length 60.
+! On first entry, it must be set to 'START'.
+! On a return with task(1:2)='FG', the user must evaluate the
+! function f and gradient g at the returned value of x.
+! On a return with task(1:5)='NEW_X', an iteration of the
+! algorithm has concluded, and f and g contain f(x) and g(x)
+! respectively. The user can decide whether to continue or stop
+! the iteration.
+! When
+! task(1:4)='CONV', the termination test in L-BFGS-B has been
+! satisfied;
+! task(1:4)='ABNO', the routine has terminated abnormally
+! without being able to satisfy the termination conditions,
+! x contains the best approximation found,
+! f and g contain f(x) and g(x) respectively;
+! task(1:5)='ERROR', the routine has detected an error in the
+! input parameters;
+! On exit with task = 'CONV', 'ABNO' or 'ERROR', the variable task
+! contains additional information that the user can print.
+! This array should not be altered unless the user wants to
+! stop the run for some reason. See driver2 or driver3
+! for a detailed explanation on how to stop the run
+! by assigning task(1:4)='STOP' in the driver.
+!
+! iprint is an INTEGER variable that must be set by the user.
+! It controls the frequency and type of output generated:
+! iprint<0 no output is generated;
+! iprint=0 print only one line at the last iteration;
+! 0<iprint<99 print also f and |proj g| every iprint iterations;
+! iprint=99 print details of every iteration except n-vectors;
+! iprint=100 print also the changes of active set and final x;
+! iprint>100 print details of every iteration including x and g;
+! When iprint > 0, the file iterate.dat will be created to
+! summarize the iteration.
+!
+! csave is a CHARACTER working array of length 60.
+!
+! lsave is a LOGICAL working array of dimension 4.
+! On exit with task = 'NEW_X', the following information is
+! available:
+! lsave(1) = .true. the initial x did not satisfy the bounds;
+! lsave(2) = .true. the problem contains bounds;
+! lsave(3) = .true. each variable has upper and lower bounds.
+!
+! isave is an INTEGER working array of dimension 44.
+! On exit with task = 'NEW_X', it contains information that
+! the user may want to access:
+! isave(30) = the current iteration number;
+! isave(34) = the total number of function and gradient
+! evaluations;
+! isave(36) = the number of function value or gradient
+! evaluations in the current iteration;
+! isave(38) = the number of free variables in the current
+! iteration;
+! isave(39) = the number of active constraints at the current
+! iteration;
+!
+! see the subroutine setulb.f for a description of other
+! information contained in isave
+!
+! dsave is a DOUBLE PRECISION working array of dimension 29.
+! On exit with task = 'NEW_X', it contains information that
+! the user may want to access:
+! dsave(2) = the value of f at the previous iteration;
+! dsave(5) = the machine precision epsmch generated by the code;
+! dsave(13) = the infinity norm of the projected gradient;
+!
+! see the subroutine setulb.f for a description of other
+! information contained in dsave
+!
+! --------------------------------------------------------------
+! END OF THE DESCRIPTION OF THE VARIABLES IN L-BFGS-B
+! --------------------------------------------------------------
+!
+ program driver
+!
+! This simple driver demonstrates how to call the L-BFGS-B code to
+! solve a sample problem (the extended Rosenbrock function
+! subject to bounds on the variables). The dimension n of this
+! problem is variable.
+
+ implicit none
+!
+! Declare variables and parameters needed by the code.
+! Note thar we wish to have output at every iteration.
+! iprint=1
+!
+! We also specify the tolerances in the stopping criteria.
+! factr = 1.0d+7, pgtol = 1.0d-5
+!
+! A description of all these variables is given at the beginning
+! of the driver
+!
+ integer, parameter :: n = 25, m = 5, iprint = 1
+ integer, parameter :: dp = kind(1.0d0)
+ real(dp), parameter :: factr = 1.0d+7, pgtol = 1.0d-5
+!
+ character(len=60) :: task, csave
+ logical :: lsave(4)
+ integer :: isave(44)
+ real(dp) :: f
+ real(dp) :: dsave(29)
+ integer, allocatable :: nbd(:), iwa(:)
+ real(dp), allocatable :: x(:), l(:), u(:), g(:), wa(:)
+
+! Declare a few additional variables for this sample problem
+
+ real(dp) :: t1, t2
+ integer :: i
+
+! Allocate dynamic arrays
+
+ allocate ( nbd(n), x(n), l(n), u(n), g(n) )
+ allocate ( iwa(3*n) )
+ allocate ( wa(2*m*n + 5*n + 11*m*m + 8*m) )
+!
+ do 10 i=1, n, 2
+ nbd(i) = 2
+ l(i) = 1.0d0
+ u(i) = 1.0d2
+ 10 continue
+
+! Next set bounds on the even-numbered variables.
+
+ do 12 i=2, n, 2
+ nbd(i) = 2
+ l(i) = -1.0d2
+ u(i) = 1.0d2
+ 12 continue
+
+! We now define the starting point.
+
+ do 14 i=1, n
+ x(i) = 3.0d0
+ 14 continue
+
+ write (6,16)
+ 16 format(/,5x, 'Solving sample problem.', &
+ /,5x, ' (f = 0.0 at the optimal solution.)',/)
+
+! We start the iteration by initializing task.
+
+ task = 'START'
+
+! The beginning of the loop
+
+ do while(task(1:2).eq.'FG'.or.task.eq.'NEW_X'.or. &
+ task.eq.'START')
+
+! This is the call to the L-BFGS-B code.
+
+ call setulb ( n, m, x, l, u, nbd, f, g, factr, pgtol, &
+ wa, iwa, task, iprint,&
+ csave, lsave, isave, dsave )
+
+ if (task(1:2) .eq. 'FG') then
+
+ f=.25d0*( x(1)-1.d0 )**2
+ do 20 i=2, n
+ f = f + ( x(i)-x(i-1 )**2 )**2
+ 20 continue
+ f = 4.d0*f
+
+! Compute gradient g for the sample problem.
+
+ t1 = x(2) - x(1)**2
+ g(1) = 2.d0*(x(1) - 1.d0) - 1.6d1*x(1)*t1
+ do 22 i=2, n-1
+ t2 = t1
+ t1 = x(i+1) - x(i)**2
+ g(i) = 8.d0*t2 - 1.6d1*x(i)*t1
+ 22 continue
+ g(n) = 8.d0*t1
+
+ end if
+ end do
+
+! end of loop do while
+
+
+ end program driver
+
diff --git a/driver2.f b/driver2.f
new file mode 100644
index 0000000..1ac91c9
--- /dev/null
+++ b/driver2.f
@@ -0,0 +1,232 @@
+c
+c L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License”
+c or “3-clause license”)
+c Please read attached file License.txt
+c
+c DRIVER 2 in Fortran 77
+c --------------------------------------------------------------
+c CUSTOMIZED DRIVER FOR L-BFGS-B (version 3.0)
+c --------------------------------------------------------------
+c
+c L-BFGS-B is a code for solving large nonlinear optimization
+c problems with simple bounds on the variables.
+c
+c The code can also be used for unconstrained problems and is
+c as efficient for these problems as the earlier limited memory
+c code L-BFGS.
+c
+c This driver illustrates how to control the termination of the
+c run and how to design customized output.
+c
+c References:
+c
+c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+c memory algorithm for bound constrained optimization'',
+c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+c
+c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
+c Subroutines for Large Scale Bound Constrained Optimization''
+c Tech. Report, NAM-11, EECS Department, Northwestern University,
+c 1994.
+c
+c
+c (Postscript files of these papers are available via anonymous
+c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
+c
+c * * *
+c
+c February 2011 (latest revision)
+c Optimization Center at Northwestern University
+c Instituto Tecnologico Autonomo de Mexico
+c
+c Jorge Nocedal and Jose Luis Morales
+c Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778:
+c L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained
+c Optimization" (2011). To appear in ACM Transactions on
+c Mathematical Software,
+c
+c **************
+
+ program driver
+
+c This driver shows how to replace the default stopping test
+c by other termination criteria. It also illustrates how to
+c print the values of several parameters during the course of
+c the iteration. The sample problem used here is the same as in
+c DRIVER1 (the extended Rosenbrock function with bounds on the
+c variables).
+
+ integer nmax, mmax
+ parameter (nmax=1024, mmax=17)
+c nmax is the dimension of the largest problem to be solved.
+c mmax is the maximum number of limited memory corrections.
+
+c Declare the variables needed by the code.
+c A description of all these variables is given at the end of
+c driver1.
+
+ character*60 task, csave
+ logical lsave(4)
+ integer n, m, iprint,
+ + nbd(nmax), iwa(3*nmax), isave(44)
+ double precision f, factr, pgtol,
+ + x(nmax), l(nmax), u(nmax), g(nmax), dsave(29),
+ + wa(2*mmax*nmax+5*nmax+11*mmax*mmax+8*mmax)
+
+c Declare a few additional variables for the sample problem.
+
+ double precision t1, t2
+ integer i
+
+c We suppress the default output.
+
+ iprint = -1
+
+c We suppress both code-supplied stopping tests because the
+c user is providing his own stopping criteria.
+
+ factr=0.0d0
+ pgtol=0.0d0
+
+c We specify the dimension n of the sample problem and the number
+c m of limited memory corrections stored. (n and m should not
+c exceed the limits nmax and mmax respectively.)
+
+ n=25
+ m=5
+
+c We now specify nbd which defines the bounds on the variables:
+c l specifies the lower bounds,
+c u specifies the upper bounds.
+
+c First set bounds on the odd numbered variables.
+
+ do 10 i=1,n,2
+ nbd(i)=2
+ l(i)=1.0d0
+ u(i)=1.0d2
+ 10 continue
+
+c Next set bounds on the even numbered variables.
+
+ do 12 i=2,n,2
+ nbd(i)=2
+ l(i)=-1.0d2
+ u(i)=1.0d2
+ 12 continue
+
+c We now define the starting point.
+
+ do 14 i=1,n
+ x(i)=3.0d0
+ 14 continue
+
+c We now write the heading of the output.
+
+ write (6,16)
+ 16 format(/,5x, 'Solving sample problem.',
+ + /,5x, ' (f = 0.0 at the optimal solution.)',/)
+
+c We start the iteration by initializing task.
+c
+ task = 'START'
+
+c ------- the beginning of the loop ----------
+
+ 111 continue
+
+c This is the call to the L-BFGS-B code.
+
+ call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint,
+ + csave,lsave,isave,dsave)
+
+ if (task(1:2) .eq. 'FG') then
+c the minimization routine has returned to request the
+c function f and gradient g values at the current x.
+
+c Compute function value f for the sample problem.
+
+ f=.25d0*(x(1)-1.d0)**2
+ do 20 i=2,n
+ f=f+(x(i)-x(i-1)**2)**2
+ 20 continue
+ f=4.d0*f
+
+c Compute gradient g for the sample problem.
+
+ t1=x(2)-x(1)**2
+ g(1)=2.d0*(x(1)-1.d0)-1.6d1*x(1)*t1
+ do 22 i=2,n-1
+ t2=t1
+ t1=x(i+1)-x(i)**2
+ g(i)=8.d0*t2-1.6d1*x(i)*t1
+ 22 continue
+ g(n)=8.d0*t1
+
+c go back to the minimization routine.
+ goto 111
+ endif
+c
+ if (task(1:5) .eq. 'NEW_X') then
+c
+c the minimization routine has returned with a new iterate.
+c At this point have the opportunity of stopping the iteration
+c or observing the values of certain parameters
+c
+c First are two examples of stopping tests.
+
+c Note: task(1:4) must be assigned the value 'STOP' to terminate
+c the iteration and ensure that the final results are
+c printed in the default format. The rest of the character
+c string TASK may be used to store other information.
+
+c 1) Terminate if the total number of f and g evaluations
+c exceeds 99.
+
+ if (isave(34) .ge. 99)
+ + task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT'
+
+c 2) Terminate if |proj g|/(1+|f|) < 1.0d-10, where
+c "proj g" denoted the projected gradient
+
+ if (dsave(13) .le. 1.d-10*(1.0d0 + abs(f)))
+ + task='STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL'
+
+c We now wish to print the following information at each
+c iteration:
+c
+c 1) the current iteration number, isave(30),
+c 2) the total number of f and g evaluations, isave(34),
+c 3) the value of the objective function f,
+c 4) the norm of the projected gradient, dsve(13)
+c
+c See the comments at the end of driver1 for a description
+c of the variables isave and dsave.
+
+ write (6,'(2(a,i5,4x),a,1p,d12.5,4x,a,1p,d12.5)') 'Iterate'
+ + ,isave(30),'nfg =',isave(34),'f =',f,'|proj g| =',dsave(13)
+
+c If the run is to be terminated, we print also the information
+c contained in task as well as the final value of x.
+
+ if (task(1:4) .eq. 'STOP') then
+ write (6,*) task
+ write (6,*) 'Final X='
+ write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n)
+ endif
+
+c go back to the minimization routine.
+ goto 111
+
+ endif
+
+c ---------- the end of the loop -------------
+
+c If task is neither FG nor NEW_X we terminate execution.
+
+ stop
+
+ end
+
+c======================= The end of driver2 ============================
+
diff --git a/driver2.f90 b/driver2.f90
new file mode 100755
index 0000000..48e53dd
--- /dev/null
+++ b/driver2.f90
@@ -0,0 +1,220 @@
+!
+! L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License”
+! or “3-clause license”)
+! Please read attached file License.txt
+!
+!
+! DRIVER 2 in Fortran 90
+! --------------------------------------------------------------
+! CUSTOMIZED DRIVER FOR L-BFGS-B
+! --------------------------------------------------------------
+!
+! L-BFGS-B is a code for solving large nonlinear optimization
+! problems with simple bounds on the variables.
+!
+! The code can also be used for unconstrained problems and is
+! as efficient for these problems as the earlier limited memory
+! code L-BFGS.
+!
+! This driver illustrates how to control the termination of the
+! run and how to design customized output.
+!
+! References:
+!
+! [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+! memory algorithm for bound constrained optimization'',
+! SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+!
+! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
+! Subroutines for Large Scale Bound Constrained Optimization''
+! Tech. Report, NAM-11, EECS Department, Northwestern University,
+! 1994.
+!
+!
+! (Postscript files of these papers are available via anonymous
+! ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
+!
+! * * *
+!
+! February 2011 (latest revision)
+! Optimization Center at Northwestern University
+! Instituto Tecnologico Autonomo de Mexico
+!
+! Jorge Nocedal and Jose Luis Morales
+!
+! **************
+ program driver
+
+! This driver shows how to replace the default stopping test
+! by other termination criteria. It also illustrates how to
+! print the values of several parameters during the course of
+! the iteration. The sample problem used here is the same as in
+! DRIVER1 (the extended Rosenbrock function with bounds on the
+! variables).
+
+ implicit none
+
+! Declare variables and parameters needed by the code.
+!
+! Note that we suppress the default output (iprint = -1)
+! We suppress both code-supplied stopping tests because the
+! user is providing his/her own stopping criteria.
+
+ integer, parameter :: n = 25, m = 5, iprint = -1
+ integer, parameter :: dp = kind(1.0d0)
+ real(dp), parameter :: factr = 0.0d0, pgtol = 0.0d0
+
+ character(len=60) :: task, csave
+ logical :: lsave(4)
+ integer :: isave(44)
+ real(dp) :: f
+ real(dp) :: dsave(29)
+ integer, allocatable :: nbd(:), iwa(:)
+ real(dp), allocatable :: x(:), l(:), u(:), g(:), wa(:)
+!
+ real(dp) :: t1, t2
+ integer :: i
+
+ allocate ( nbd(n), x(n), l(n), u(n), g(n) )
+ allocate ( iwa(3*n) )
+ allocate ( wa(2*m*n + 5*n + 11*m*m + 8*m) )
+!
+! This driver shows how to replace the default stopping test
+! by other termination criteria. It also illustrates how to
+! print the values of several parameters during the course of
+! the iteration. The sample problem used here is the same as in
+! DRIVER1 (the extended Rosenbrock function with bounds on the
+! variables).
+! We now specify nbd which defines the bounds on the variables:
+! l specifies the lower bounds,
+! u specifies the upper bounds.
+
+! First set bounds on the odd numbered variables.
+
+ do 10 i=1, n,2
+ nbd(i)=2
+ l(i)=1.0d0
+ u(i)=1.0d2
+ 10 continue
+
+! Next set bounds on the even numbered variables.
+
+ do 12 i=2, n,2
+ nbd(i)=2
+ l(i)=-1.0d2
+ u(i)=1.0d2
+ 12 continue
+
+! We now define the starting point.
+
+ do 14 i=1, n
+ x(i)=3.0d0
+ 14 continue
+
+! We now write the heading of the output.
+
+ write (6,16)
+ 16 format(/,5x, 'Solving sample problem.', &
+ /,5x, ' (f = 0.0 at the optimal solution.)',/)
+
+
+! We start the iteration by initializing task.
+!
+ task = 'START'
+
+! ------- the beginning of the loop ----------
+
+ do while( task(1:2).eq.'FG'.or.task.eq.'NEW_X'.or. &
+ task.eq.'START')
+
+! This is the call to the L-BFGS-B code.
+
+ call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint, &
+ csave,lsave,isave,dsave)
+
+ if (task(1:2) .eq. 'FG') then
+
+! the minimization routine has returned to request the
+! function f and gradient g values at the current x.
+
+! Compute function value f for the sample problem.
+
+ f =.25d0*(x(1) - 1.d0)**2
+ do 20 i=2,n
+ f = f + (x(i) - x(i-1)**2)**2
+ 20 continue
+ f = 4.d0*f
+
+! Compute gradient g for the sample problem.
+
+ t1 = x(2) - x(1)**2
+ g(1) = 2.d0*(x(1) - 1.d0) - 1.6d1*x(1)*t1
+ do 22 i= 2,n-1
+ t2 = t1
+ t1 = x(i+1) - x(i)**2
+ g(i) = 8.d0*t2 - 1.6d1*x(i)*t1
+ 22 continue
+ g(n)=8.d0*t1
+!
+ else
+!
+ if (task(1:5) .eq. 'NEW_X') then
+!
+! the minimization routine has returned with a new iterate.
+! At this point have the opportunity of stopping the iteration
+! or observing the values of certain parameters
+!
+! First are two examples of stopping tests.
+
+! Note: task(1:4) must be assigned the value 'STOP' to terminate
+! the iteration and ensure that the final results are
+! printed in the default format. The rest of the character
+! string TASK may be used to store other information.
+
+! 1) Terminate if the total number of f and g evaluations
+! exceeds 99.
+
+ if (isave(34) .ge. 99) &
+ task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT'
+
+! 2) Terminate if |proj g|/(1+|f|) < 1.0d-10, where
+! "proj g" denoted the projected gradient
+
+ if (dsave(13) .le. 1.d-10*(1.0d0 + abs(f))) &
+ task='STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL'
+
+! We now wish to print the following information at each
+! iteration:
+!
+! 1) the current iteration number, isave(30),
+! 2) the total number of f and g evaluations, isave(34),
+! 3) the value of the objective function f,
+! 4) the norm of the projected gradient, dsve(13)
+!
+! See the comments at the end of driver1 for a description
+! of the variables isave and dsave.
+
+ write (6,'(2(a,i5,4x),a,1p,d12.5,4x,a,1p,d12.5)') 'Iterate' &
+ , isave(30),'nfg =',isave(34),'f =',f,'|proj g| =',dsave(13)
+
+! If the run is to be terminated, we print also the information
+! contained in task as well as the final value of x.
+
+ if (task(1:4) .eq. 'STOP') then
+ write (6,*) task
+ write (6,*) 'Final X='
+ write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n)
+ end if
+
+ end if
+ end if
+
+ end do
+! ---------- the end of the loop -------------
+
+! If task is neither FG nor NEW_X we terminate execution.
+
+ end program driver
+
+!======================= The end of driver2 ============================
+
diff --git a/driver3.f b/driver3.f
new file mode 100644
index 0000000..8d1371d
--- /dev/null
+++ b/driver3.f
@@ -0,0 +1,273 @@
+c
+c L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License”
+c or “3-clause license”)
+c Please read attached file License.txt
+c
+c DRIVER 3 in Fortran 77
+c --------------------------------------------------------------
+c TIME-CONTROLLED DRIVER FOR L-BFGS-B (version 3.0)
+c --------------------------------------------------------------
+c
+c L-BFGS-B is a code for solving large nonlinear optimization
+c problems with simple bounds on the variables.
+c
+c The code can also be used for unconstrained problems and is
+c as efficient for these problems as the earlier limited memory
+c code L-BFGS.
+c
+c This driver shows how to terminate a run after some prescribed
+c CPU time has elapsed, and how to print the desired information
+c before exiting.
+c
+c References:
+c
+c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+c memory algorithm for bound constrained optimization'',
+c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+c
+c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
+c Subroutines for Large Scale Bound Constrained Optimization''
+c Tech. Report, NAM-11, EECS Department, Northwestern University,
+c 1994.
+c
+c
+c (Postscript files of these papers are available via anonymous
+c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
+c
+c * * *
+c
+c February 2011 (latest revision)
+c Optimization Center at Northwestern University
+c Instituto Tecnologico Autonomo de Mexico
+c
+c Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778:
+c L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained
+c Optimization" (2011). To appear in ACM Transactions on
+c Mathematical Software,
+c
+c
+c **************
+
+ program driver
+
+c This time-controlled driver shows that it is possible to terminate
+c a run by elapsed CPU time, and yet be able to print all desired
+c information. This driver also illustrates the use of two
+c stopping criteria that may be used in conjunction with a limit
+c on execution time. The sample problem used here is the same as in
+c driver1 and driver2 (the extended Rosenbrock function with bounds
+c on the variables).
+
+ integer nmax, mmax
+ parameter (nmax=1024,mmax=17)
+c nmax is the dimension of the largest problem to be solved.
+c mmax is the maximum number of limited memory corrections.
+
+c Declare the variables needed by the code.
+c A description of all these variables is given at the end of
+c driver1.
+
+ character*60 task, csave
+ logical lsave(4)
+ integer n, m, iprint,
+ + nbd(nmax), iwa(3*nmax), isave(44)
+ double precision f, factr, pgtol,
+ + x(nmax), l(nmax), u(nmax), g(nmax), dsave(29),
+ + wa(2*mmax*nmax+5*nmax+11*mmax*mmax+8*mmax)
+
+c Declare a few additional variables for the sample problem
+c and for keeping track of time.
+
+ double precision t1, t2, time1, time2, tlimit
+ integer i, j
+
+c We specify a limite on the CPU time (in seconds).
+
+ tlimit = 0.2
+
+c We suppress the default output. (The user could also elect to
+c use the default output by choosing iprint >= 0.)
+
+ iprint = -1
+
+c We suppress the code-supplied stopping tests because we will
+c provide our own termination conditions
+
+ factr=0.0d0
+ pgtol=0.0d0
+
+c We specify the dimension n of the sample problem and the number
+c m of limited memory corrections stored. (n and m should not
+c exceed the limits nmax and mmax respectively.)
+
+ n=1000
+ m=10
+
+c We now specify nbd which defines the bounds on the variables:
+c l specifies the lower bounds,
+c u specifies the upper bounds.
+
+c First set bounds on the odd-numbered variables.
+
+ do 10 i=1,n,2
+ nbd(i)=2
+ l(i)=1.0d0
+ u(i)=1.0d2
+ 10 continue
+
+c Next set bounds on the even-numbered variables.
+
+ do 12 i=2,n,2
+ nbd(i)=2
+ l(i)=-1.0d2
+ u(i)=1.0d2
+ 12 continue
+
+c We now define the starting point.
+
+ do 14 i=1,n
+ x(i)=3.0d0
+ 14 continue
+
+c We now write the heading of the output.
+
+ write (6,16)
+ 16 format(/,5x, 'Solving sample problem.',
+ + /,5x, ' (f = 0.0 at the optimal solution.)',/)
+
+c We start the iteration by initializing task.
+c
+ task = 'START'
+
+c ------- the beginning of the loop ----------
+
+c We begin counting the CPU time.
+
+ call timer(time1)
+
+ 111 continue
+
+c This is the call to the L-BFGS-B code.
+
+ call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint,
+ + csave,lsave,isave,dsave)
+
+ if (task(1:2) .eq. 'FG') then
+c the minimization routine has returned to request the
+c function f and gradient g values at the current x.
+c Before evaluating f and g we check the CPU time spent.
+
+ call timer(time2)
+ if (time2-time1 .gt. tlimit) then
+ task='STOP: CPU EXCEEDING THE TIME LIMIT.'
+
+c Note: Assigning task(1:4)='STOP' will terminate the run;
+c setting task(7:9)='CPU' will restore the information at
+c the latest iterate generated by the code so that it can
+c be correctly printed by the driver.
+
+c In this driver we have chosen to disable the
+c printing options of the code (we set iprint=-1);
+c instead we are using customized output: we print the
+c latest value of x, the corresponding function value f and
+c the norm of the projected gradient |proj g|.
+
+c We print out the information contained in task.
+
+ write (6,*) task
+
+c We print the latest iterate contained in wa(j+1:j+n), where
+c
+ j = 3*n+2*m*n+11*m**2
+ write (6,*) 'Latest iterate X ='
+ write (6,'((1x,1p, 6(1x,d11.4)))') (wa(i),i = j+1,j+n)
+
+c We print the function value f and the norm of the projected
+c gradient |proj g| at the last iterate; they are stored in
+c dsave(2) and dsave(13) respectively.
+
+ write (6,'(a,1p,d12.5,4x,a,1p,d12.5)')
+ + 'At latest iterate f =',dsave(2),'|proj g| =',dsave(13)
+
+ else
+
+c The time limit has not been reached and we compute
+c the function value f for the sample problem.
+
+ f=.25d0*(x(1)-1.d0)**2
+ do 20 i=2,n
+ f=f+(x(i)-x(i-1)**2)**2
+ 20 continue
+ f=4.d0*f
+
+c Compute gradient g for the sample problem.
+
+ t1=x(2)-x(1)**2
+ g(1)=2.d0*(x(1)-1.d0)-1.6d1*x(1)*t1
+ do 22 i=2,n-1
+ t2=t1
+ t1=x(i+1)-x(i)**2
+ g(i)=8.d0*t2-1.6d1*x(i)*t1
+ 22 continue
+ g(n)=8.d0*t1
+
+ endif
+
+c go back to the minimization routine.
+ goto 111
+ endif
+c
+ if (task(1:5) .eq. 'NEW_X') then
+c the minimization routine has returned with a new iterate.
+c The time limit has not been reached, and we test whether
+c the following two stopping tests are satisfied:
+
+c 1) Terminate if the total number of f and g evaluations
+c exceeds 900.
+
+ if (isave(34) .ge. 900)
+ + task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT'
+
+c 2) Terminate if |proj g|/(1+|f|) < 1.0d-10.
+
+ if (dsave(13) .le. 1.d-10*(1.0d0 + abs(f)))
+ + task='STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL'
+
+c We wish to print the following information at each iteration:
+c 1) the current iteration number, isave(30),
+c 2) the total number of f and g evaluations, isave(34),
+c 3) the value of the objective function f,
+c 4) the norm of the projected gradient, dsve(13)
+c
+c See the comments at the end of driver1 for a description
+c of the variables isave and dsave.
+
+
+ write (6,'(2(a,i5,4x),a,1p,d12.5,4x,a,1p,d12.5)') 'Iterate'
+ + ,isave(30),'nfg =',isave(34),'f =',f,'|proj g| =',dsave(13)
+
+c If the run is to be terminated, we print also the information
+c contained in task as well as the final value of x.
+
+
+ if (task(1:4) .eq. 'STOP') then
+ write (6,*) task
+ write (6,*) 'Final X='
+ write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n)
+ endif
+
+c go back to the minimization routine.
+ goto 111
+
+ endif
+
+c ---------- the end of the loop -------------
+
+c If task is neither FG nor NEW_X we terminate execution.
+
+ stop
+
+ end
+
+c======================= The end of driver3 ============================
+
diff --git a/driver3.f90 b/driver3.f90
new file mode 100644
index 0000000..4cc9b36
--- /dev/null
+++ b/driver3.f90
@@ -0,0 +1,256 @@
+!
+! L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License”
+! or “3-clause license”)
+! Please read attached file License.txt
+!
+! DRIVER 3 in Fortran 90
+! --------------------------------------------------------------
+! TIME-CONTROLLED DRIVER FOR L-BFGS-B
+! --------------------------------------------------------------
+!
+! L-BFGS-B is a code for solving large nonlinear optimization
+! problems with simple bounds on the variables.
+!
+! The code can also be used for unconstrained problems and is
+! as efficient for these problems as the earlier limited memory
+! code L-BFGS.
+!
+! This driver shows how to terminate a run after some prescribed
+! CPU time has elapsed, and how to print the desired information
+! before exiting.
+!
+! References:
+!
+! [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+! memory algorithm for bound constrained optimization'',
+! SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+!
+! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
+! Subroutines for Large Scale Bound Constrained Optimization''
+! Tech. Report, NAM-11, EECS Department, Northwestern University,
+! 1994.
+!
+!
+! (Postscript files of these papers are available via anonymous
+! ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
+!
+! * * *
+!
+! February 2011 (latest revision)
+! Optimization Center at Northwestern University
+! Instituto Tecnologico Autonomo de Mexico
+!
+! Jorge Nocedal and Jose Luis Morales
+!
+! **************
+
+ program driver
+
+! This time-controlled driver shows that it is possible to terminate
+! a run by elapsed CPU time, and yet be able to print all desired
+! information. This driver also illustrates the use of two
+! stopping criteria that may be used in conjunction with a limit
+! on execution time. The sample problem used here is the same as in
+! driver1 and driver2 (the extended Rosenbrock function with bounds
+! on the variables).
+
+ implicit none
+
+! We specify a limit on the CPU time (tlimit = 10 seconds)
+!
+! We suppress the default output (iprint = -1). The user could
+! also elect to use the default output by choosing iprint >= 0.)
+! We suppress the code-supplied stopping tests because we will
+! provide our own termination conditions
+! We specify the dimension n of the sample problem and the number
+! m of limited memory corrections stored.
+
+ integer, parameter :: n = 1000, m = 10, iprint = -1
+ integer, parameter :: dp = kind(1.0d0)
+ real(dp), parameter :: factr = 0.0d0, pgtol = 0.0d0, &
+ tlimit = 10.0d0
+!
+ character(len=60) :: task, csave
+ logical :: lsave(4)
+ integer :: isave(44)
+ real(dp) :: f
+ real(dp) :: dsave(29)
+ integer, allocatable :: nbd(:), iwa(:)
+ real(dp), allocatable :: x(:), l(:), u(:), g(:), wa(:)
+!
+ real(dp) :: t1, t2, time1, time2
+ integer :: i, j
+
+ allocate ( nbd(n), x(n), l(n), u(n), g(n) )
+ allocate ( iwa(3*n) )
+ allocate ( wa(2*m*n + 5*n + 11*m*m + 8*m) )
+
+! This time-controlled driver shows that it is possible to terminate
+! a run by elapsed CPU time, and yet be able to print all desired
+! information. This driver also illustrates the use of two
+! stopping criteria that may be used in conjunction with a limit
+! on execution time. The sample problem used here is the same as in
+! driver1 and driver2 (the extended Rosenbrock function with bounds
+! on the variables).
+
+! We now specify nbd which defines the bounds on the variables:
+! l specifies the lower bounds,
+! u specifies the upper bounds.
+
+! First set bounds on the odd-numbered variables.
+
+ do 10 i=1, n,2
+ nbd(i)=2
+ l(i)=1.0d0
+ u(i)=1.0d2
+ 10 continue
+
+! Next set bounds on the even-numbered variables.
+
+ do 12 i=2, n,2
+ nbd(i)=2
+ l(i)=-1.0d2
+ u(i)=1.0d2
+ 12 continue
+
+! We now define the starting point.
+
+ do 14 i=1, n
+ x(i)=3.0d0
+ 14 continue
+
+! We now write the heading of the output.
+
+ write (6,16)
+ 16 format(/,5x, 'Solving sample problem.',&
+ /,5x, ' (f = 0.0 at the optimal solution.)',/)
+
+! We start the iteration by initializing task.
+
+ task = 'START'
+
+! ------- the beginning of the loop ----------
+
+! We begin counting the CPU time.
+
+ call timer(time1)
+
+ do while( task(1:2).eq.'FG'.or.task.eq.'NEW_X'.or. &
+ task.eq.'START')
+
+! This is the call to the L-BFGS-B code.
+
+ call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa, &
+ task,iprint, csave,lsave,isave,dsave)
+
+ if (task(1:2) .eq. 'FG') then
+
+! the minimization routine has returned to request the
+! function f and gradient g values at the current x.
+! Before evaluating f and g we check the CPU time spent.
+
+ call timer(time2)
+ if (time2-time1 .gt. tlimit) then
+ task='STOP: CPU EXCEEDING THE TIME LIMIT.'
+
+! Note: Assigning task(1:4)='STOP' will terminate the run;
+! setting task(7:9)='CPU' will restore the information at
+! the latest iterate generated by the code so that it can
+! be correctly printed by the driver.
+
+! In this driver we have chosen to disable the
+! printing options of the code (we set iprint=-1);
+! instead we are using customized output: we print the
+! latest value of x, the corresponding function value f and
+! the norm of the projected gradient |proj g|.
+
+! We print out the information contained in task.
+
+ write (6,*) task
+
+! We print the latest iterate contained in wa(j+1:j+n), where
+
+ j = 3*n+2*m*n+11*m**2
+ write (6,*) 'Latest iterate X ='
+ write (6,'((1x,1p, 6(1x,d11.4)))') (wa(i),i = j+1,j+n)
+
+! We print the function value f and the norm of the projected
+! gradient |proj g| at the last iterate; they are stored in
+! dsave(2) and dsave(13) respectively.
+
+ write (6,'(a,1p,d12.5,4x,a,1p,d12.5)') &
+ 'At latest iterate f =',dsave(2),'|proj g| =',dsave(13)
+ else
+
+! The time limit has not been reached and we compute
+! the function value f for the sample problem.
+
+ f=.25d0*(x(1)-1.d0)**2
+ do 20 i=2, n
+ f=f+(x(i)-x(i-1)**2)**2
+ 20 continue
+ f=4.d0*f
+
+! Compute gradient g for the sample problem.
+
+ t1 = x(2) - x(1)**2
+ g(1) = 2.d0*(x(1)-1.d0)-1.6d1*x(1)*t1
+ do 22 i=2,n-1
+ t2=t1
+ t1=x(i+1)-x(i)**2
+ g(i)=8.d0*t2-1.6d1*x(i)*t1
+ 22 continue
+ g(n)=8.d0*t1
+ endif
+
+! go back to the minimization routine.
+ else
+
+ if (task(1:5) .eq. 'NEW_X') then
+
+! the minimization routine has returned with a new iterate.
+! The time limit has not been reached, and we test whether
+! the following two stopping tests are satisfied:
+
+! 1) Terminate if the total number of f and g evaluations
+! exceeds 900.
+
+ if (isave(34) .ge. 900) &
+ task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT'
+
+! 2) Terminate if |proj g|/(1+|f|) < 1.0d-10.
+
+ if (dsave(13) .le. 1.d-10*(1.0d0 + abs(f))) &
+ task='STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL'
+
+! We wish to print the following information at each iteration:
+! 1) the current iteration number, isave(30),
+! 2) the total number of f and g evaluations, isave(34),
+! 3) the value of the objective function f,
+! 4) the norm of the projected gradient, dsve(13)
+!
+! See the comments at the end of driver1 for a description
+! of the variables isave and dsave.
+
+ write (6,'(2(a,i5,4x),a,1p,d12.5,4x,a,1p,d12.5)') 'Iterate' &
+ ,isave(30),'nfg =',isave(34),'f =',f,'|proj g| =',dsave(13)
+
+! If the run is to be terminated, we print also the information
+! contained in task as well as the final value of x.
+
+ if (task(1:4) .eq. 'STOP') then
+ write (6,*) task
+ write (6,*) 'Final X='
+ write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n)
+ endif
+
+ endif
+ end if
+ end do
+
+! If task is neither FG nor NEW_X we terminate execution.
+
+ end program driver
+
+!======================= The end of driver3 ============================
+
diff --git a/iterate.dat b/iterate.dat
new file mode 100644
index 0000000..b9dbf62
--- /dev/null
+++ b/iterate.dat
@@ -0,0 +1,49 @@
+RUNNING THE L-BFGS-B CODE
+
+it = iteration number
+nf = number of function evaluations
+nseg = number of segments explored during the Cauchy search
+nact = number of active bounds at the generalized Cauchy point
+sub = manner in which the subspace minimization terminated:
+ con = converged, bnd = a bound was reached
+itls = number of iterations performed in the line search
+stepl = step length used
+tstep = norm of the displacement (total step)
+projg = norm of the projected gradient
+f = function value
+
+ * * *
+
+Machine precision = 2.220D-16
+ N = 25 M = 5
+
+ it nf nseg nact sub itls stepl tstep projg f
+ 0 1 - - - - - - 1.030D+02 3.460D+03
+ 1 5 25 24 --- 3 1.2D-02 4.2D+00 6.507D+01 2.398D+03
+ 2 6 1 0 con 0 1.0D+00 4.2D+00 3.640D+01 1.438D+02
+ 3 7 1 0 con 0 1.0D+00 6.8D-01 2.290D+01 7.282D+01
+ 4 8 1 0 con 0 1.0D+00 1.1D+00 6.954D+00 1.603D+01
+ 5 9 1 0 con 0 1.0D+00 6.5D-01 9.055D+00 5.187D+00
+ 6 10 1 0 bnd 0 1.0D+00 6.1D-01 1.967D+01 2.127D+00
+ 7 11 1 0 bnd 0 1.0D+00 1.8D-01 2.128D+00 2.076D-01
+ 8 12 1 0 bnd 0 1.0D+00 9.0D-02 8.325D-01 5.327D-02
+ 9 13 1 0 bnd 0 1.0D+00 5.3D-02 4.279D-01 1.305D-02
+ 10 14 1 1 bnd 0 1.0D+00 2.5D-02 2.008D-01 3.860D-03
+ 11 15 1 2 bnd 0 1.0D+00 2.0D-02 1.377D-01 7.457D-04
+ 12 16 1 0 bnd 0 1.0D+00 7.8D-03 1.213D-01 3.540D-04
+ 13 17 1 2 con 0 1.0D+00 3.1D-03 2.978D-02 7.425D-05
+ 14 18 1 0 con 0 1.0D+00 1.1D-03 1.727D-02 3.741D-05
+ 15 19 1 0 con 0 1.0D+00 2.1D-03 2.868D-02 1.098D-05
+ 16 21 1 0 con 1 4.6D-01 4.8D-04 8.081D-03 3.908D-06
+ 17 22 1 0 con 0 1.0D+00 2.2D-04 3.476D-03 1.995D-06
+ 18 23 1 0 con 0 1.0D+00 2.2D-04 2.253D-03 8.258D-07
+ 19 24 1 0 con 0 1.0D+00 2.4D-04 1.458D-03 1.992D-07
+ 20 25 1 0 con 0 1.0D+00 1.8D-04 1.482D-03 5.758D-08
+ 21 26 1 0 con 0 1.0D+00 6.8D-05 5.443D-04 1.463D-08
+ 22 27 1 0 con 0 1.0D+00 2.3D-05 2.252D-04 2.363D-09
+ 23 28 1 0 con 0 1.0D+00 8.0D-06 1.721D-04 1.083D-09
+
+CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
+
+ Total User time 3.337E-03 seconds.
+
diff --git a/lbfgsb.f b/lbfgsb.f
new file mode 100644
index 0000000..9c9e7d9
--- /dev/null
+++ b/lbfgsb.f
@@ -0,0 +1,3949 @@
+c
+c L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License”
+c or “3-clause license”)
+c Please read attached file License.txt
+c
+c=========== L-BFGS-B (version 3.0. April 25, 2011 ===================
+c
+c This is a modified version of L-BFGS-B. Minor changes in the updated
+c code appear preceded by a line comment as follows
+c
+c c-jlm-jn
+c
+c Major changes are described in the accompanying paper:
+c
+c Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778:
+c L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained
+c Optimization" (2011). To appear in ACM Transactions on
+c Mathematical Software,
+c
+c The paper describes an improvement and a correction to Algorithm 778.
+c It is shown that the performance of the algorithm can be improved
+c significantly by making a relatively simple modication to the subspace
+c minimization phase. The correction concerns an error caused by the use
+c of routine dpmeps to estimate machine precision.
+c
+c The total work space **wa** required by the new version is
+c
+c 2*m*n + 11m*m + 5*n + 8*m
+c
+c the old version required
+c
+c 2*m*n + 12m*m + 4*n + 12*m
+c
+c
+c J. Nocedal Department of Electrical Engineering and
+c Computer Science.
+c Northwestern University. Evanston, IL. USA
+c
+c
+c J.L Morales Departamento de Matematicas,
+c Instituto Tecnologico Autonomo de Mexico
+c Mexico D.F. Mexico.
+c
+c March 2011
+c
+c=============================================================================
+ subroutine setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa,
+ + task, iprint, csave, lsave, isave, dsave)
+
+ character*60 task, csave
+ logical lsave(4)
+ integer n, m, iprint,
+ + nbd(n), iwa(3*n), isave(44)
+ double precision f, factr, pgtol, x(n), l(n), u(n), g(n),
+c
+c-jlm-jn
+ + wa(2*m*n + 5*n + 11*m*m + 8*m), dsave(29)
+
+c ************
+c
+c Subroutine setulb
+c
+c This subroutine partitions the working arrays wa and iwa, and
+c then uses the limited memory BFGS method to solve the bound
+c constrained optimization problem by calling mainlb.
+c (The direct method will be used in the subspace minimization.)
+c
+c n is an integer variable.
+c On entry n is the dimension of the problem.
+c On exit n is unchanged.
+c
+c m is an integer variable.
+c On entry m is the maximum number of variable metric corrections
+c used to define the limited memory matrix.
+c On exit m is unchanged.
+c
+c x is a double precision array of dimension n.
+c On entry x is an approximation to the solution.
+c On exit x is the current approximation.
+c
+c l is a double precision array of dimension n.
+c On entry l is the lower bound on x.
+c On exit l is unchanged.
+c
+c u is a double precision array of dimension n.
+c On entry u is the upper bound on x.
+c On exit u is unchanged.
+c
+c nbd is an integer array of dimension n.
+c On entry nbd represents the type of bounds imposed on the
+c variables, and must be specified as follows:
+c nbd(i)=0 if x(i) is unbounded,
+c 1 if x(i) has only a lower bound,
+c 2 if x(i) has both lower and upper bounds, and
+c 3 if x(i) has only an upper bound.
+c On exit nbd is unchanged.
+c
+c f is a double precision variable.
+c On first entry f is unspecified.
+c On final exit f is the value of the function at x.
+c
+c g is a double precision array of dimension n.
+c On first entry g is unspecified.
+c On final exit g is the value of the gradient at x.
+c
+c factr is a double precision variable.
+c On entry factr >= 0 is specified by the user. The iteration
+c will stop when
+c
+c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
+c
+c where epsmch is the machine precision, which is automatically
+c generated by the code. Typical values for factr: 1.d+12 for
+c low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely
+c high accuracy.
+c On exit factr is unchanged.
+c
+c pgtol is a double precision variable.
+c On entry pgtol >= 0 is specified by the user. The iteration
+c will stop when
+c
+c max{|proj g_i | i = 1, ..., n} <= pgtol
+c
+c where pg_i is the ith component of the projected gradient.
+c On exit pgtol is unchanged.
+c
+c wa is a double precision working array of length
+c (2mmax + 5)nmax + 12mmax^2 + 12mmax.
+c
+c iwa is an integer working array of length 3nmax.
+c
+c task is a working string of characters of length 60 indicating
+c the current job when entering and quitting this subroutine.
+c
+c iprint is an integer variable that must be set by the user.
+c It controls the frequency and type of output generated:
+c iprint<0 no output is generated;
+c iprint=0 print only one line at the last iteration;
+c 0<iprint<99 print also f and |proj g| every iprint iterations;
+c iprint=99 print details of every iteration except n-vectors;
+c iprint=100 print also the changes of active set and final x;
+c iprint>100 print details of every iteration including x and g;
+c When iprint > 0, the file iterate.dat will be created to
+c summarize the iteration.
+c
+c csave is a working string of characters of length 60.
+c
+c lsave is a logical working array of dimension 4.
+c On exit with 'task' = NEW_X, the following information is
+c available:
+c If lsave(1) = .true. then the initial X has been replaced by
+c its projection in the feasible set;
+c If lsave(2) = .true. then the problem is constrained;
+c If lsave(3) = .true. then each variable has upper and lower
+c bounds;
+c
+c isave is an integer working array of dimension 44.
+c On exit with 'task' = NEW_X, the following information is
+c available:
+c isave(22) = the total number of intervals explored in the
+c search of Cauchy points;
+c isave(26) = the total number of skipped BFGS updates before
+c the current iteration;
+c isave(30) = the number of current iteration;
+c isave(31) = the total number of BFGS updates prior the current
+c iteration;
+c isave(33) = the number of intervals explored in the search of
+c Cauchy point in the current iteration;
+c isave(34) = the total number of function and gradient
+c evaluations;
+c isave(36) = the number of function value or gradient
+c evaluations in the current iteration;
+c if isave(37) = 0 then the subspace argmin is within the box;
+c if isave(37) = 1 then the subspace argmin is beyond the box;
+c isave(38) = the number of free variables in the current
+c iteration;
+c isave(39) = the number of active constraints in the current
+c iteration;
+c n + 1 - isave(40) = the number of variables leaving the set of
+c active constraints in the current iteration;
+c isave(41) = the number of variables entering the set of active
+c constraints in the current iteration.
+c
+c dsave is a double precision working array of dimension 29.
+c On exit with 'task' = NEW_X, the following information is
+c available:
+c dsave(1) = current 'theta' in the BFGS matrix;
+c dsave(2) = f(x) in the previous iteration;
+c dsave(3) = factr*epsmch;
+c dsave(4) = 2-norm of the line search direction vector;
+c dsave(5) = the machine precision epsmch generated by the code;
+c dsave(7) = the accumulated time spent on searching for
+c Cauchy points;
+c dsave(8) = the accumulated time spent on
+c subspace minimization;
+c dsave(9) = the accumulated time spent on line search;
+c dsave(11) = the slope of the line search function at
+c the current point of line search;
+c dsave(12) = the maximum relative step length imposed in
+c line search;
+c dsave(13) = the infinity norm of the projected gradient;
+c dsave(14) = the relative step length in the line search;
+c dsave(15) = the slope of the line search function at
+c the starting point of the line search;
+c dsave(16) = the square of the 2-norm of the line search
+c direction vector.
+c
+c Subprograms called:
+c
+c L-BFGS-B Library ... mainlb.
+c
+c
+c References:
+c
+c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+c memory algorithm for bound constrained optimization'',
+c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+c
+c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
+c limited memory FORTRAN code for solving bound constrained
+c optimization problems'', Tech. Report, NAM-11, EECS Department,
+c Northwestern University, 1994.
+c
+c (Postscript files of these papers are available via anonymous
+c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+c-jlm-jn
+ integer lws,lr,lz,lt,ld,lxp,lwa,
+ + lwy,lsy,lss,lwt,lwn,lsnd
+
+ if (task .eq. 'START') then
+ isave(1) = m*n
+ isave(2) = m**2
+ isave(3) = 4*m**2
+ isave(4) = 1 ! ws m*n
+ isave(5) = isave(4) + isave(1) ! wy m*n
+ isave(6) = isave(5) + isave(1) ! wsy m**2
+ isave(7) = isave(6) + isave(2) ! wss m**2
+ isave(8) = isave(7) + isave(2) ! wt m**2
+ isave(9) = isave(8) + isave(2) ! wn 4*m**2
+ isave(10) = isave(9) + isave(3) ! wsnd 4*m**2
+ isave(11) = isave(10) + isave(3) ! wz n
+ isave(12) = isave(11) + n ! wr n
+ isave(13) = isave(12) + n ! wd n
+ isave(14) = isave(13) + n ! wt n
+ isave(15) = isave(14) + n ! wxp n
+ isave(16) = isave(15) + n ! wa 8*m
+ endif
+ lws = isave(4)
+ lwy = isave(5)
+ lsy = isave(6)
+ lss = isave(7)
+ lwt = isave(8)
+ lwn = isave(9)
+ lsnd = isave(10)
+ lz = isave(11)
+ lr = isave(12)
+ ld = isave(13)
+ lt = isave(14)
+ lxp = isave(15)
+ lwa = isave(16)
+
+ call mainlb(n,m,x,l,u,nbd,f,g,factr,pgtol,
+ + wa(lws),wa(lwy),wa(lsy),wa(lss), wa(lwt),
+ + wa(lwn),wa(lsnd),wa(lz),wa(lr),wa(ld),wa(lt),wa(lxp),
+ + wa(lwa),
+ + iwa(1),iwa(n+1),iwa(2*n+1),task,iprint,
+ + csave,lsave,isave(22),dsave)
+
+ return
+
+ end
+
+c======================= The end of setulb =============================
+
+ subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy,
+ + sy, ss, wt, wn, snd, z, r, d, t, xp, wa,
+ + index, iwhere, indx2, task,
+ + iprint, csave, lsave, isave, dsave)
+ implicit none
+ character*60 task, csave
+ logical lsave(4)
+ integer n, m, iprint, nbd(n), index(n),
+ + iwhere(n), indx2(n), isave(23)
+ double precision f, factr, pgtol,
+ + x(n), l(n), u(n), g(n), z(n), r(n), d(n), t(n),
+c-jlm-jn
+ + xp(n),
+ + wa(8*m),
+ + ws(n, m), wy(n, m), sy(m, m), ss(m, m),
+ + wt(m, m), wn(2*m, 2*m), snd(2*m, 2*m), dsave(29)
+
+c ************
+c
+c Subroutine mainlb
+c
+c This subroutine solves bound constrained optimization problems by
+c using the compact formula of the limited memory BFGS updates.
+c
+c n is an integer variable.
+c On entry n is the number of variables.
+c On exit n is unchanged.
+c
+c m is an integer variable.
+c On entry m is the maximum number of variable metric
+c corrections allowed in the limited memory matrix.
+c On exit m is unchanged.
+c
+c x is a double precision array of dimension n.
+c On entry x is an approximation to the solution.
+c On exit x is the current approximation.
+c
+c l is a double precision array of dimension n.
+c On entry l is the lower bound of x.
+c On exit l is unchanged.
+c
+c u is a double precision array of dimension n.
+c On entry u is the upper bound of x.
+c On exit u is unchanged.
+c
+c nbd is an integer array of dimension n.
+c On entry nbd represents the type of bounds imposed on the
+c variables, and must be specified as follows:
+c nbd(i)=0 if x(i) is unbounded,
+c 1 if x(i) has only a lower bound,
+c 2 if x(i) has both lower and upper bounds,
+c 3 if x(i) has only an upper bound.
+c On exit nbd is unchanged.
+c
+c f is a double precision variable.
+c On first entry f is unspecified.
+c On final exit f is the value of the function at x.
+c
+c g is a double precision array of dimension n.
+c On first entry g is unspecified.
+c On final exit g is the value of the gradient at x.
+c
+c factr is a double precision variable.
+c On entry factr >= 0 is specified by the user. The iteration
+c will stop when
+c
+c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
+c
+c where epsmch is the machine precision, which is automatically
+c generated by the code.
+c On exit factr is unchanged.
+c
+c pgtol is a double precision variable.
+c On entry pgtol >= 0 is specified by the user. The iteration
+c will stop when
+c
+c max{|proj g_i | i = 1, ..., n} <= pgtol
+c
+c where pg_i is the ith component of the projected gradient.
+c On exit pgtol is unchanged.
+c
+c ws, wy, sy, and wt are double precision working arrays used to
+c store the following information defining the limited memory
+c BFGS matrix:
+c ws, of dimension n x m, stores S, the matrix of s-vectors;
+c wy, of dimension n x m, stores Y, the matrix of y-vectors;
+c sy, of dimension m x m, stores S'Y;
+c ss, of dimension m x m, stores S'S;
+c yy, of dimension m x m, stores Y'Y;
+c wt, of dimension m x m, stores the Cholesky factorization
+c of (theta*S'S+LD^(-1)L'); see eq.
+c (2.26) in [3].
+c
+c wn is a double precision working array of dimension 2m x 2m
+c used to store the LEL^T factorization of the indefinite matrix
+c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
+c [L_a -R_z theta*S'AA'S ]
+c
+c where E = [-I 0]
+c [ 0 I]
+c
+c snd is a double precision working array of dimension 2m x 2m
+c used to store the lower triangular part of
+c N = [Y' ZZ'Y L_a'+R_z']
+c [L_a +R_z S'AA'S ]
+c
+c z(n),r(n),d(n),t(n), xp(n),wa(8*m) are double precision working arrays.
+c z is used at different times to store the Cauchy point and
+c the Newton point.
+c xp is used to safeguard the projected Newton direction
+c
+c sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays.
+c
+c index is an integer working array of dimension n.
+c In subroutine freev, index is used to store the free and fixed
+c variables at the Generalized Cauchy Point (GCP).
+c
+c iwhere is an integer working array of dimension n used to record
+c the status of the vector x for GCP computation.
+c iwhere(i)=0 or -3 if x(i) is free and has bounds,
+c 1 if x(i) is fixed at l(i), and l(i) .ne. u(i)
+c 2 if x(i) is fixed at u(i), and u(i) .ne. l(i)
+c 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i)
+c -1 if x(i) is always free, i.e., no bounds on it.
+c
+c indx2 is an integer working array of dimension n.
+c Within subroutine cauchy, indx2 corresponds to the array iorder.
+c In subroutine freev, a list of variables entering and leaving
+c the free set is stored in indx2, and it is passed on to
+c subroutine formk with this information.
+c
+c task is a working string of characters of length 60 indicating
+c the current job when entering and leaving this subroutine.
+c
+c iprint is an INTEGER variable that must be set by the user.
+c It controls the frequency and type of output generated:
+c iprint<0 no output is generated;
+c iprint=0 print only one line at the last iteration;
+c 0<iprint<99 print also f and |proj g| every iprint iterations;
+c iprint=99 print details of every iteration except n-vectors;
+c iprint=100 print also the changes of active set and final x;
+c iprint>100 print details of every iteration including x and g;
+c When iprint > 0, the file iterate.dat will be created to
+c summarize the iteration.
+c
+c csave is a working string of characters of length 60.
+c
+c lsave is a logical working array of dimension 4.
+c
+c isave is an integer working array of dimension 23.
+c
+c dsave is a double precision working array of dimension 29.
+c
+c
+c Subprograms called
+c
+c L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk,
+c
+c errclb, prn1lb, prn2lb, prn3lb, active, projgr,
+c
+c freev, cmprlb, matupd, formt.
+c
+c Minpack2 Library ... timer
+c
+c Linpack Library ... dcopy, ddot.
+c
+c
+c References:
+c
+c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+c memory algorithm for bound constrained optimization'',
+c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+c
+c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
+c Subroutines for Large Scale Bound Constrained Optimization''
+c Tech. Report, NAM-11, EECS Department, Northwestern University,
+c 1994.
+c
+c [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of
+c Quasi-Newton Matrices and their use in Limited Memory Methods'',
+c Mathematical Programming 63 (1994), no. 4, pp. 129-156.
+c
+c (Postscript files of these papers are available via anonymous
+c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ logical prjctd,cnstnd,boxed,updatd,wrk
+ character*3 word
+ integer i,k,nintol,itfile,iback,nskip,
+ + head,col,iter,itail,iupdat,
+ + nseg,nfgv,info,ifun,
+ + iword,nfree,nact,ileave,nenter
+ double precision theta,fold,ddot,dr,rr,tol,
+ + xstep,sbgnrm,ddum,dnorm,dtd,epsmch,
+ + cpu1,cpu2,cachyt,sbtime,lnscht,time1,time2,
+ + gd,gdold,stp,stpmx,time
+ double precision one,zero
+ parameter (one=1.0d0,zero=0.0d0)
+
+ if (task .eq. 'START') then
+
+ epsmch = epsilon(one)
+
+ call timer(time1)
+
+c Initialize counters and scalars when task='START'.
+
+c for the limited memory BFGS matrices:
+ col = 0
+ head = 1
+ theta = one
+ iupdat = 0
+ updatd = .false.
+ iback = 0
+ itail = 0
+ iword = 0
+ nact = 0
+ ileave = 0
+ nenter = 0
+ fold = zero
+ dnorm = zero
+ cpu1 = zero
+ gd = zero
+ stpmx = zero
+ sbgnrm = zero
+ stp = zero
+ gdold = zero
+ dtd = zero
+
+c for operation counts:
+ iter = 0
+ nfgv = 0
+ nseg = 0
+ nintol = 0
+ nskip = 0
+ nfree = n
+ ifun = 0
+c for stopping tolerance:
+ tol = factr*epsmch
+
+c for measuring running time:
+ cachyt = 0
+ sbtime = 0
+ lnscht = 0
+
+c 'word' records the status of subspace solutions.
+ word = '---'
+
+c 'info' records the termination information.
+ info = 0
+
+ itfile = 8
+ if (iprint .ge. 1) then
+c open a summary file 'iterate.dat'
+ open (8, file = 'iterate.dat', status = 'unknown')
+ endif
+
+c Check the input arguments for errors.
+
+ call errclb(n,m,factr,l,u,nbd,task,info,k)
+ if (task(1:5) .eq. 'ERROR') then
+ call prn3lb(n,x,f,task,iprint,info,itfile,
+ + iter,nfgv,nintol,nskip,nact,sbgnrm,
+ + zero,nseg,word,iback,stp,xstep,k,
+ + cachyt,sbtime,lnscht)
+ return
+ endif
+
+ call prn1lb(n,m,l,u,x,iprint,itfile,epsmch)
+
+c Initialize iwhere & project x onto the feasible set.
+
+ call active(n,l,u,nbd,x,iwhere,iprint,prjctd,cnstnd,boxed)
+
+c The end of the initialization.
+
+ else
+c restore local variables.
+
+ prjctd = lsave(1)
+ cnstnd = lsave(2)
+ boxed = lsave(3)
+ updatd = lsave(4)
+
+ nintol = isave(1)
+ itfile = isave(3)
+ iback = isave(4)
+ nskip = isave(5)
+ head = isave(6)
+ col = isave(7)
+ itail = isave(8)
+ iter = isave(9)
+ iupdat = isave(10)
+ nseg = isave(12)
+ nfgv = isave(13)
+ info = isave(14)
+ ifun = isave(15)
+ iword = isave(16)
+ nfree = isave(17)
+ nact = isave(18)
+ ileave = isave(19)
+ nenter = isave(20)
+
+ theta = dsave(1)
+ fold = dsave(2)
+ tol = dsave(3)
+ dnorm = dsave(4)
+ epsmch = dsave(5)
+ cpu1 = dsave(6)
+ cachyt = dsave(7)
+ sbtime = dsave(8)
+ lnscht = dsave(9)
+ time1 = dsave(10)
+ gd = dsave(11)
+ stpmx = dsave(12)
+ sbgnrm = dsave(13)
+ stp = dsave(14)
+ gdold = dsave(15)
+ dtd = dsave(16)
+
+c After returning from the driver go to the point where execution
+c is to resume.
+
+ if (task(1:5) .eq. 'FG_LN') goto 666
+ if (task(1:5) .eq. 'NEW_X') goto 777
+ if (task(1:5) .eq. 'FG_ST') goto 111
+ if (task(1:4) .eq. 'STOP') then
+ if (task(7:9) .eq. 'CPU') then
+c restore the previous iterate.
+ call dcopy(n,t,1,x,1)
+ call dcopy(n,r,1,g,1)
+ f = fold
+ endif
+ goto 999
+ endif
+ endif
+
+c Compute f0 and g0.
+
+ task = 'FG_START'
+c return to the driver to calculate f and g; reenter at 111.
+ goto 1000
+ 111 continue
+ nfgv = 1
+
+c Compute the infinity norm of the (-) projected gradient.
+
+ call projgr(n,l,u,nbd,x,g,sbgnrm)
+
+ if (iprint .ge. 1) then
+ write (6,1002) iter,f,sbgnrm
+ write (itfile,1003) iter,nfgv,sbgnrm,f
+ endif
+ if (sbgnrm .le. pgtol) then
+c terminate the algorithm.
+ task = 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
+ goto 999
+ endif
+
+c ----------------- the beginning of the loop --------------------------
+
+ 222 continue
+ if (iprint .ge. 99) write (6,1001) iter + 1
+ iword = -1
+c
+ if (.not. cnstnd .and. col .gt. 0) then
+c skip the search for GCP.
+ call dcopy(n,x,1,z,1)
+ wrk = updatd
+ nseg = 0
+ goto 333
+ endif
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c Compute the Generalized Cauchy Point (GCP).
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+ call timer(cpu1)
+ call cauchy(n,x,l,u,nbd,g,indx2,iwhere,t,d,z,
+ + m,wy,ws,sy,wt,theta,col,head,
+ + wa(1),wa(2*m+1),wa(4*m+1),wa(6*m+1),nseg,
+ + iprint, sbgnrm, info, epsmch)
+ if (info .ne. 0) then
+c singular triangular system detected; refresh the lbfgs memory.
+ if(iprint .ge. 1) write (6, 1005)
+ info = 0
+ col = 0
+ head = 1
+ theta = one
+ iupdat = 0
+ updatd = .false.
+ call timer(cpu2)
+ cachyt = cachyt + cpu2 - cpu1
+ goto 222
+ endif
+ call timer(cpu2)
+ cachyt = cachyt + cpu2 - cpu1
+ nintol = nintol + nseg
+
+c Count the entering and leaving variables for iter > 0;
+c find the index set of free and active variables at the GCP.
+
+ call freev(n,nfree,index,nenter,ileave,indx2,
+ + iwhere,wrk,updatd,cnstnd,iprint,iter)
+ nact = n - nfree
+
+ 333 continue
+
+c If there are no free variables or B=theta*I, then
+c skip the subspace minimization.
+
+ if (nfree .eq. 0 .or. col .eq. 0) goto 555
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c Subspace minimization.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+ call timer(cpu1)
+
+c Form the LEL^T factorization of the indefinite
+c matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
+c [L_a -R_z theta*S'AA'S ]
+c where E = [-I 0]
+c [ 0 I]
+
+ if (wrk) call formk(n,nfree,index,nenter,ileave,indx2,iupdat,
+ + updatd,wn,snd,m,ws,wy,sy,theta,col,head,info)
+ if (info .ne. 0) then
+c nonpositive definiteness in Cholesky factorization;
+c refresh the lbfgs memory and restart the iteration.
+ if(iprint .ge. 1) write (6, 1006)
+ info = 0
+ col = 0
+ head = 1
+ theta = one
+ iupdat = 0
+ updatd = .false.
+ call timer(cpu2)
+ sbtime = sbtime + cpu2 - cpu1
+ goto 222
+ endif
+
+c compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x)
+c from 'cauchy').
+ call cmprlb(n,m,x,g,ws,wy,sy,wt,z,r,wa,index,
+ + theta,col,head,nfree,cnstnd,info)
+ if (info .ne. 0) goto 444
+
+c-jlm-jn call the direct method.
+
+ call subsm( n, m, nfree, index, l, u, nbd, z, r, xp, ws, wy,
+ + theta, x, g, col, head, iword, wa, wn, iprint, info)
+ 444 continue
+ if (info .ne. 0) then
+c singular triangular system detected;
+c refresh the lbfgs memory and restart the iteration.
+ if(iprint .ge. 1) write (6, 1005)
+ info = 0
+ col = 0
+ head = 1
+ theta = one
+ iupdat = 0
+ updatd = .false.
+ call timer(cpu2)
+ sbtime = sbtime + cpu2 - cpu1
+ goto 222
+ endif
+
+ call timer(cpu2)
+ sbtime = sbtime + cpu2 - cpu1
+ 555 continue
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c Line search and optimality tests.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c Generate the search direction d:=z-x.
+
+ do 40 i = 1, n
+ d(i) = z(i) - x(i)
+ 40 continue
+ call timer(cpu1)
+ 666 continue
+ call lnsrlb(n,l,u,nbd,x,f,fold,gd,gdold,g,d,r,t,z,stp,dnorm,
+ + dtd,xstep,stpmx,iter,ifun,iback,nfgv,info,task,
+ + boxed,cnstnd,csave,isave(22),dsave(17))
+ if (info .ne. 0 .or. iback .ge. 20) then
+c restore the previous iterate.
+ call dcopy(n,t,1,x,1)
+ call dcopy(n,r,1,g,1)
+ f = fold
+ if (col .eq. 0) then
+c abnormal termination.
+ if (info .eq. 0) then
+ info = -9
+c restore the actual number of f and g evaluations etc.
+ nfgv = nfgv - 1
+ ifun = ifun - 1
+ iback = iback - 1
+ endif
+ task = 'ABNORMAL_TERMINATION_IN_LNSRCH'
+ iter = iter + 1
+ goto 999
+ else
+c refresh the lbfgs memory and restart the iteration.
+ if(iprint .ge. 1) write (6, 1008)
+ if (info .eq. 0) nfgv = nfgv - 1
+ info = 0
+ col = 0
+ head = 1
+ theta = one
+ iupdat = 0
+ updatd = .false.
+ task = 'RESTART_FROM_LNSRCH'
+ call timer(cpu2)
+ lnscht = lnscht + cpu2 - cpu1
+ goto 222
+ endif
+ else if (task(1:5) .eq. 'FG_LN') then
+c return to the driver for calculating f and g; reenter at 666.
+ goto 1000
+ else
+c calculate and print out the quantities related to the new X.
+ call timer(cpu2)
+ lnscht = lnscht + cpu2 - cpu1
+ iter = iter + 1
+
+c Compute the infinity norm of the projected (-)gradient.
+
+ call projgr(n,l,u,nbd,x,g,sbgnrm)
+
+c Print iteration information.
+
+ call prn2lb(n,x,f,g,iprint,itfile,iter,nfgv,nact,
+ + sbgnrm,nseg,word,iword,iback,stp,xstep)
+ goto 1000
+ endif
+ 777 continue
+
+c Test for termination.
+
+ if (sbgnrm .le. pgtol) then
+c terminate the algorithm.
+ task = 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
+ goto 999
+ endif
+
+ ddum = max(abs(fold), abs(f), one)
+ if ((fold - f) .le. tol*ddum) then
+c terminate the algorithm.
+ task = 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
+ if (iback .ge. 10) info = -5
+c i.e., to issue a warning if iback>10 in the line search.
+ goto 999
+ endif
+
+c Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's.
+
+ do 42 i = 1, n
+ r(i) = g(i) - r(i)
+ 42 continue
+ rr = ddot(n,r,1,r,1)
+ if (stp .eq. one) then
+ dr = gd - gdold
+ ddum = -gdold
+ else
+ dr = (gd - gdold)*stp
+ call dscal(n,stp,d,1)
+ ddum = -gdold*stp
+ endif
+
+ if (dr .le. epsmch*ddum) then
+c skip the L-BFGS update.
+ nskip = nskip + 1
+ updatd = .false.
+ if (iprint .ge. 1) write (6,1004) dr, ddum
+ goto 888
+ endif
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c Update the L-BFGS matrix.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+ updatd = .true.
+ iupdat = iupdat + 1
+
+c Update matrices WS and WY and form the middle matrix in B.
+
+ call matupd(n,m,ws,wy,sy,ss,d,r,itail,
+ + iupdat,col,head,theta,rr,dr,stp,dtd)
+
+c Form the upper half of the pds T = theta*SS + L*D^(-1)*L';
+c Store T in the upper triangular of the array wt;
+c Cholesky factorize T to J*J' with
+c J' stored in the upper triangular of wt.
+
+ call formt(m,wt,sy,ss,col,theta,info)
+
+ if (info .ne. 0) then
+c nonpositive definiteness in Cholesky factorization;
+c refresh the lbfgs memory and restart the iteration.
+ if(iprint .ge. 1) write (6, 1007)
+ info = 0
+ col = 0
+ head = 1
+ theta = one
+ iupdat = 0
+ updatd = .false.
+ goto 222
+ endif
+
+c Now the inverse of the middle matrix in B is
+
+c [ D^(1/2) O ] [ -D^(1/2) D^(-1/2)*L' ]
+c [ -L*D^(-1/2) J ] [ 0 J' ]
+
+ 888 continue
+
+c -------------------- the end of the loop -----------------------------
+
+ goto 222
+ 999 continue
+ call timer(time2)
+ time = time2 - time1
+ call prn3lb(n,x,f,task,iprint,info,itfile,
+ + iter,nfgv,nintol,nskip,nact,sbgnrm,
+ + time,nseg,word,iback,stp,xstep,k,
+ + cachyt,sbtime,lnscht)
+ 1000 continue
+
+c Save local variables.
+
+ lsave(1) = prjctd
+ lsave(2) = cnstnd
+ lsave(3) = boxed
+ lsave(4) = updatd
+
+ isave(1) = nintol
+ isave(3) = itfile
+ isave(4) = iback
+ isave(5) = nskip
+ isave(6) = head
+ isave(7) = col
+ isave(8) = itail
+ isave(9) = iter
+ isave(10) = iupdat
+ isave(12) = nseg
+ isave(13) = nfgv
+ isave(14) = info
+ isave(15) = ifun
+ isave(16) = iword
+ isave(17) = nfree
+ isave(18) = nact
+ isave(19) = ileave
+ isave(20) = nenter
+
+ dsave(1) = theta
+ dsave(2) = fold
+ dsave(3) = tol
+ dsave(4) = dnorm
+ dsave(5) = epsmch
+ dsave(6) = cpu1
+ dsave(7) = cachyt
+ dsave(8) = sbtime
+ dsave(9) = lnscht
+ dsave(10) = time1
+ dsave(11) = gd
+ dsave(12) = stpmx
+ dsave(13) = sbgnrm
+ dsave(14) = stp
+ dsave(15) = gdold
+ dsave(16) = dtd
+
+ 1001 format (//,'ITERATION ',i5)
+ 1002 format
+ + (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5)
+ 1003 format (2(1x,i4),5x,'-',5x,'-',3x,'-',5x,'-',5x,'-',8x,'-',3x,
+ + 1p,2(1x,d10.3))
+ 1004 format (' ys=',1p,e10.3,' -gs=',1p,e10.3,' BFGS update SKIPPED')
+ 1005 format (/,
+ +' Singular triangular system detected;',/,
+ +' refresh the lbfgs memory and restart the iteration.')
+ 1006 format (/,
+ +' Nonpositive definiteness in Cholesky factorization in formk;',/,
+ +' refresh the lbfgs memory and restart the iteration.')
+ 1007 format (/,
+ +' Nonpositive definiteness in Cholesky factorization in formt;',/,
+ +' refresh the lbfgs memory and restart the iteration.')
+ 1008 format (/,
+ +' Bad direction in the line search;',/,
+ +' refresh the lbfgs memory and restart the iteration.')
+
+ return
+
+ end
+
+c======================= The end of mainlb =============================
+
+ subroutine active(n, l, u, nbd, x, iwhere, iprint,
+ + prjctd, cnstnd, boxed)
+
+ logical prjctd, cnstnd, boxed
+ integer n, iprint, nbd(n), iwhere(n)
+ double precision x(n), l(n), u(n)
+
+c ************
+c
+c Subroutine active
+c
+c This subroutine initializes iwhere and projects the initial x to
+c the feasible set if necessary.
+c
+c iwhere is an integer array of dimension n.
+c On entry iwhere is unspecified.
+c On exit iwhere(i)=-1 if x(i) has no bounds
+c 3 if l(i)=u(i)
+c 0 otherwise.
+c In cauchy, iwhere is given finer gradations.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer nbdd,i
+ double precision zero
+ parameter (zero=0.0d0)
+
+c Initialize nbdd, prjctd, cnstnd and boxed.
+
+ nbdd = 0
+ prjctd = .false.
+ cnstnd = .false.
+ boxed = .true.
+
+c Project the initial x to the easible set if necessary.
+
+ do 10 i = 1, n
+ if (nbd(i) .gt. 0) then
+ if (nbd(i) .le. 2 .and. x(i) .le. l(i)) then
+ if (x(i) .lt. l(i)) then
+ prjctd = .true.
+ x(i) = l(i)
+ endif
+ nbdd = nbdd + 1
+ else if (nbd(i) .ge. 2 .and. x(i) .ge. u(i)) then
+ if (x(i) .gt. u(i)) then
+ prjctd = .true.
+ x(i) = u(i)
+ endif
+ nbdd = nbdd + 1
+ endif
+ endif
+ 10 continue
+
+c Initialize iwhere and assign values to cnstnd and boxed.
+
+ do 20 i = 1, n
+ if (nbd(i) .ne. 2) boxed = .false.
+ if (nbd(i) .eq. 0) then
+c this variable is always free
+ iwhere(i) = -1
+
+c otherwise set x(i)=mid(x(i), u(i), l(i)).
+ else
+ cnstnd = .true.
+ if (nbd(i) .eq. 2 .and. u(i) - l(i) .le. zero) then
+c this variable is always fixed
+ iwhere(i) = 3
+ else
+ iwhere(i) = 0
+ endif
+ endif
+ 20 continue
+
+ if (iprint .ge. 0) then
+ if (prjctd) write (6,*)
+ + 'The initial X is infeasible. Restart with its projection.'
+ if (.not. cnstnd)
+ + write (6,*) 'This problem is unconstrained.'
+ endif
+
+ if (iprint .gt. 0) write (6,1001) nbdd
+
+ 1001 format (/,'At X0 ',i9,' variables are exactly at the bounds')
+
+ return
+
+ end
+
+c======================= The end of active =============================
+
+ subroutine bmv(m, sy, wt, col, v, p, info)
+
+ integer m, col, info
+ double precision sy(m, m), wt(m, m), v(2*col), p(2*col)
+
+c ************
+c
+c Subroutine bmv
+c
+c This subroutine computes the product of the 2m x 2m middle matrix
+c in the compact L-BFGS formula of B and a 2m vector v;
+c it returns the product in p.
+c
+c m is an integer variable.
+c On entry m is the maximum number of variable metric corrections
+c used to define the limited memory matrix.
+c On exit m is unchanged.
+c
+c sy is a double precision array of dimension m x m.
+c On entry sy specifies the matrix S'Y.
+c On exit sy is unchanged.
+c
+c wt is a double precision array of dimension m x m.
+c On entry wt specifies the upper triangular matrix J' which is
+c the Cholesky factor of (thetaS'S+LD^(-1)L').
+c On exit wt is unchanged.
+c
+c col is an integer variable.
+c On entry col specifies the number of s-vectors (or y-vectors)
+c stored in the compact L-BFGS formula.
+c On exit col is unchanged.
+c
+c v is a double precision array of dimension 2col.
+c On entry v specifies vector v.
+c On exit v is unchanged.
+c
+c p is a double precision array of dimension 2col.
+c On entry p is unspecified.
+c On exit p is the product Mv.
+c
+c info is an integer variable.
+c On entry info is unspecified.
+c On exit info = 0 for normal return,
+c = nonzero for abnormal return when the system
+c to be solved by dtrsl is singular.
+c
+c Subprograms called:
+c
+c Linpack ... dtrsl.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer i,k,i2
+ double precision sum
+
+ if (col .eq. 0) return
+
+c PART I: solve [ D^(1/2) O ] [ p1 ] = [ v1 ]
+c [ -L*D^(-1/2) J ] [ p2 ] [ v2 ].
+
+c solve Jp2=v2+LD^(-1)v1.
+ p(col + 1) = v(col + 1)
+ do 20 i = 2, col
+ i2 = col + i
+ sum = 0.0d0
+ do 10 k = 1, i - 1
+ sum = sum + sy(i,k)*v(k)/sy(k,k)
+ 10 continue
+ p(i2) = v(i2) + sum
+ 20 continue
+c Solve the triangular system
+ call dtrsl(wt,m,col,p(col+1),11,info)
+ if (info .ne. 0) return
+
+c solve D^(1/2)p1=v1.
+ do 30 i = 1, col
+ p(i) = v(i)/sqrt(sy(i,i))
+ 30 continue
+
+c PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ]
+c [ 0 J' ] [ p2 ] [ p2 ].
+
+c solve J^Tp2=p2.
+ call dtrsl(wt,m,col,p(col+1),01,info)
+ if (info .ne. 0) return
+
+c compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2)
+c =-D^(-1/2)p1+D^(-1)L'p2.
+ do 40 i = 1, col
+ p(i) = -p(i)/sqrt(sy(i,i))
+ 40 continue
+ do 60 i = 1, col
+ sum = 0.d0
+ do 50 k = i + 1, col
+ sum = sum + sy(k,i)*p(col+k)/sy(i,i)
+ 50 continue
+ p(i) = p(i) + sum
+ 60 continue
+
+ return
+
+ end
+
+c======================== The end of bmv ===============================
+
+ subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
+ + m, wy, ws, sy, wt, theta, col, head, p, c, wbp,
+ + v, nseg, iprint, sbgnrm, info, epsmch)
+ implicit none
+ integer n, m, head, col, nseg, iprint, info,
+ + nbd(n), iorder(n), iwhere(n)
+ double precision theta, epsmch,
+ + x(n), l(n), u(n), g(n), t(n), d(n), xcp(n),
+ + wy(n, col), ws(n, col), sy(m, m),
+ + wt(m, m), p(2*m), c(2*m), wbp(2*m), v(2*m)
+
+c ************
+c
+c Subroutine cauchy
+c
+c For given x, l, u, g (with sbgnrm > 0), and a limited memory
+c BFGS matrix B defined in terms of matrices WY, WS, WT, and
+c scalars head, col, and theta, this subroutine computes the
+c generalized Cauchy point (GCP), defined as the first local
+c minimizer of the quadratic
+c
+c Q(x + s) = g's + 1/2 s'Bs
+c
+c along the projected gradient direction P(x-tg,l,u).
+c The routine returns the GCP in xcp.
+c
+c n is an integer variable.
+c On entry n is the dimension of the problem.
+c On exit n is unchanged.
+c
+c x is a double precision array of dimension n.
+c On entry x is the starting point for the GCP computation.
+c On exit x is unchanged.
+c
+c l is a double precision array of dimension n.
+c On entry l is the lower bound of x.
+c On exit l is unchanged.
+c
+c u is a double precision array of dimension n.
+c On entry u is the upper bound of x.
+c On exit u is unchanged.
+c
+c nbd is an integer array of dimension n.
+c On entry nbd represents the type of bounds imposed on the
+c variables, and must be specified as follows:
+c nbd(i)=0 if x(i) is unbounded,
+c 1 if x(i) has only a lower bound,
+c 2 if x(i) has both lower and upper bounds, and
+c 3 if x(i) has only an upper bound.
+c On exit nbd is unchanged.
+c
+c g is a double precision array of dimension n.
+c On entry g is the gradient of f(x). g must be a nonzero vector.
+c On exit g is unchanged.
+c
+c iorder is an integer working array of dimension n.
+c iorder will be used to store the breakpoints in the piecewise
+c linear path and free variables encountered. On exit,
+c iorder(1),...,iorder(nleft) are indices of breakpoints
+c which have not been encountered;
+c iorder(nleft+1),...,iorder(nbreak) are indices of
+c encountered breakpoints; and
+c iorder(nfree),...,iorder(n) are indices of variables which
+c have no bound constraits along the search direction.
+c
+c iwhere is an integer array of dimension n.
+c On entry iwhere indicates only the permanently fixed (iwhere=3)
+c or free (iwhere= -1) components of x.
+c On exit iwhere records the status of the current x variables.
+c iwhere(i)=-3 if x(i) is free and has bounds, but is not moved
+c 0 if x(i) is free and has bounds, and is moved
+c 1 if x(i) is fixed at l(i), and l(i) .ne. u(i)
+c 2 if x(i) is fixed at u(i), and u(i) .ne. l(i)
+c 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i)
+c -1 if x(i) is always free, i.e., it has no bounds.
+c
+c t is a double precision working array of dimension n.
+c t will be used to store the break points.
+c
+c d is a double precision array of dimension n used to store
+c the Cauchy direction P(x-tg)-x.
+c
+c xcp is a double precision array of dimension n used to return the
+c GCP on exit.
+c
+c m is an integer variable.
+c On entry m is the maximum number of variable metric corrections
+c used to define the limited memory matrix.
+c On exit m is unchanged.
+c
+c ws, wy, sy, and wt are double precision arrays.
+c On entry they store information that defines the
+c limited memory BFGS matrix:
+c ws(n,m) stores S, a set of s-vectors;
+c wy(n,m) stores Y, a set of y-vectors;
+c sy(m,m) stores S'Y;
+c wt(m,m) stores the
+c Cholesky factorization of (theta*S'S+LD^(-1)L').
+c On exit these arrays are unchanged.
+c
+c theta is a double precision variable.
+c On entry theta is the scaling factor specifying B_0 = theta I.
+c On exit theta is unchanged.
+c
+c col is an integer variable.
+c On entry col is the actual number of variable metric
+c corrections stored so far.
+c On exit col is unchanged.
+c
+c head is an integer variable.
+c On entry head is the location of the first s-vector (or y-vector)
+c in S (or Y).
+c On exit col is unchanged.
+c
+c p is a double precision working array of dimension 2m.
+c p will be used to store the vector p = W^(T)d.
+c
+c c is a double precision working array of dimension 2m.
+c c will be used to store the vector c = W^(T)(xcp-x).
+c
+c wbp is a double precision working array of dimension 2m.
+c wbp will be used to store the row of W corresponding
+c to a breakpoint.
+c
+c v is a double precision working array of dimension 2m.
+c
+c nseg is an integer variable.
+c On exit nseg records the number of quadratic segments explored
+c in searching for the GCP.
+c
+c sg and yg are double precision arrays of dimension m.
+c On entry sg and yg store S'g and Y'g correspondingly.
+c On exit they are unchanged.
+c
+c iprint is an INTEGER variable that must be set by the user.
+c It controls the frequency and type of output generated:
+c iprint<0 no output is generated;
+c iprint=0 print only one line at the last iteration;
+c 0<iprint<99 print also f and |proj g| every iprint iterations;
+c iprint=99 print details of every iteration except n-vectors;
+c iprint=100 print also the changes of active set and final x;
+c iprint>100 print details of every iteration including x and g;
+c When iprint > 0, the file iterate.dat will be created to
+c summarize the iteration.
+c
+c sbgnrm is a double precision variable.
+c On entry sbgnrm is the norm of the projected gradient at x.
+c On exit sbgnrm is unchanged.
+c
+c info is an integer variable.
+c On entry info is 0.
+c On exit info = 0 for normal return,
+c = nonzero for abnormal return when the the system
+c used in routine bmv is singular.
+c
+c Subprograms called:
+c
+c L-BFGS-B Library ... hpsolb, bmv.
+c
+c Linpack ... dscal dcopy, daxpy.
+c
+c
+c References:
+c
+c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+c memory algorithm for bound constrained optimization'',
+c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+c
+c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
+c Subroutines for Large Scale Bound Constrained Optimization''
+c Tech. Report, NAM-11, EECS Department, Northwestern University,
+c 1994.
+c
+c (Postscript files of these papers are available via anonymous
+c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ logical xlower,xupper,bnded
+ integer i,j,col2,nfree,nbreak,pointr,
+ + ibp,nleft,ibkmin,iter
+ double precision f1,f2,dt,dtm,tsum,dibp,zibp,dibp2,bkmin,
+ + tu,tl,wmc,wmp,wmw,ddot,tj,tj0,neggi,sbgnrm,
+ + f2_org
+ double precision one,zero
+ parameter (one=1.0d0,zero=0.0d0)
+
+c Check the status of the variables, reset iwhere(i) if necessary;
+c compute the Cauchy direction d and the breakpoints t; initialize
+c the derivative f1 and the vector p = W'd (for theta = 1).
+
+ if (sbgnrm .le. zero) then
+ if (iprint .ge. 0) write (6,*) 'Subgnorm = 0. GCP = X.'
+ call dcopy(n,x,1,xcp,1)
+ return
+ endif
+ bnded = .true.
+ nfree = n + 1
+ nbreak = 0
+ ibkmin = 0
+ bkmin = zero
+ col2 = 2*col
+ f1 = zero
+ if (iprint .ge. 99) write (6,3010)
+
+c We set p to zero and build it up as we determine d.
+
+ do 20 i = 1, col2
+ p(i) = zero
+ 20 continue
+
+c In the following loop we determine for each variable its bound
+c status and its breakpoint, and update p accordingly.
+c Smallest breakpoint is identified.
+
+ do 50 i = 1, n
+ neggi = -g(i)
+ if (iwhere(i) .ne. 3 .and. iwhere(i) .ne. -1) then
+c if x(i) is not a constant and has bounds,
+c compute the difference between x(i) and its bounds.
+ if (nbd(i) .le. 2) tl = x(i) - l(i)
+ if (nbd(i) .ge. 2) tu = u(i) - x(i)
+
+c If a variable is close enough to a bound
+c we treat it as at bound.
+ xlower = nbd(i) .le. 2 .and. tl .le. zero
+ xupper = nbd(i) .ge. 2 .and. tu .le. zero
+
+c reset iwhere(i).
+ iwhere(i) = 0
+ if (xlower) then
+ if (neggi .le. zero) iwhere(i) = 1
+ else if (xupper) then
+ if (neggi .ge. zero) iwhere(i) = 2
+ else
+ if (abs(neggi) .le. zero) iwhere(i) = -3
+ endif
+ endif
+ pointr = head
+ if (iwhere(i) .ne. 0 .and. iwhere(i) .ne. -1) then
+ d(i) = zero
+ else
+ d(i) = neggi
+ f1 = f1 - neggi*neggi
+c calculate p := p - W'e_i* (g_i).
+ do 40 j = 1, col
+ p(j) = p(j) + wy(i,pointr)* neggi
+ p(col + j) = p(col + j) + ws(i,pointr)*neggi
+ pointr = mod(pointr,m) + 1
+ 40 continue
+ if (nbd(i) .le. 2 .and. nbd(i) .ne. 0
+ + .and. neggi .lt. zero) then
+c x(i) + d(i) is bounded; compute t(i).
+ nbreak = nbreak + 1
+ iorder(nbreak) = i
+ t(nbreak) = tl/(-neggi)
+ if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
+ bkmin = t(nbreak)
+ ibkmin = nbreak
+ endif
+ else if (nbd(i) .ge. 2 .and. neggi .gt. zero) then
+c x(i) + d(i) is bounded; compute t(i).
+ nbreak = nbreak + 1
+ iorder(nbreak) = i
+ t(nbreak) = tu/neggi
+ if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
+ bkmin = t(nbreak)
+ ibkmin = nbreak
+ endif
+ else
+c x(i) + d(i) is not bounded.
+ nfree = nfree - 1
+ iorder(nfree) = i
+ if (abs(neggi) .gt. zero) bnded = .false.
+ endif
+ endif
+ 50 continue
+
+c The indices of the nonzero components of d are now stored
+c in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n).
+c The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin.
+
+ if (theta .ne. one) then
+c complete the initialization of p for theta not= one.
+ call dscal(col,theta,p(col+1),1)
+ endif
+
+c Initialize GCP xcp = x.
+
+ call dcopy(n,x,1,xcp,1)
+
+ if (nbreak .eq. 0 .and. nfree .eq. n + 1) then
+c is a zero vector, return with the initial xcp as GCP.
+ if (iprint .gt. 100) write (6,1010) (xcp(i), i = 1, n)
+ return
+ endif
+
+c Initialize c = W'(xcp - x) = 0.
+
+ do 60 j = 1, col2
+ c(j) = zero
+ 60 continue
+
+c Initialize derivative f2.
+
+ f2 = -theta*f1
+ f2_org = f2
+ if (col .gt. 0) then
+ call bmv(m,sy,wt,col,p,v,info)
+ if (info .ne. 0) return
+ f2 = f2 - ddot(col2,v,1,p,1)
+ endif
+ dtm = -f1/f2
+ tsum = zero
+ nseg = 1
+ if (iprint .ge. 99)
+ + write (6,*) 'There are ',nbreak,' breakpoints '
+
+c If there are no breakpoints, locate the GCP and return.
+
+ if (nbreak .eq. 0) goto 888
+
+ nleft = nbreak
+ iter = 1
+
+
+ tj = zero
+
+c------------------- the beginning of the loop -------------------------
+
+ 777 continue
+
+c Find the next smallest breakpoint;
+c compute dt = t(nleft) - t(nleft + 1).
+
+ tj0 = tj
+ if (iter .eq. 1) then
+c Since we already have the smallest breakpoint we need not do
+c heapsort yet. Often only one breakpoint is used and the
+c cost of heapsort is avoided.
+ tj = bkmin
+ ibp = iorder(ibkmin)
+ else
+ if (iter .eq. 2) then
+c Replace the already used smallest breakpoint with the
+c breakpoint numbered nbreak > nlast, before heapsort call.
+ if (ibkmin .ne. nbreak) then
+ t(ibkmin) = t(nbreak)
+ iorder(ibkmin) = iorder(nbreak)
+ endif
+c Update heap structure of breakpoints
+c (if iter=2, initialize heap).
+ endif
+ call hpsolb(nleft,t,iorder,iter-2)
+ tj = t(nleft)
+ ibp = iorder(nleft)
+ endif
+
+ dt = tj - tj0
+
+ if (dt .ne. zero .and. iprint .ge. 100) then
+ write (6,4011) nseg,f1,f2
+ write (6,5010) dt
+ write (6,6010) dtm
+ endif
+
+c If a minimizer is within this interval, locate the GCP and return.
+
+ if (dtm .lt. dt) goto 888
+
+c Otherwise fix one variable and
+c reset the corresponding component of d to zero.
+
+ tsum = tsum + dt
+ nleft = nleft - 1
+ iter = iter + 1
+ dibp = d(ibp)
+ d(ibp) = zero
+ if (dibp .gt. zero) then
+ zibp = u(ibp) - x(ibp)
+ xcp(ibp) = u(ibp)
+ iwhere(ibp) = 2
+ else
+ zibp = l(ibp) - x(ibp)
+ xcp(ibp) = l(ibp)
+ iwhere(ibp) = 1
+ endif
+ if (iprint .ge. 100) write (6,*) 'Variable ',ibp,' is fixed.'
+ if (nleft .eq. 0 .and. nbreak .eq. n) then
+c all n variables are fixed,
+c return with xcp as GCP.
+ dtm = dt
+ goto 999
+ endif
+
+c Update the derivative information.
+
+ nseg = nseg + 1
+ dibp2 = dibp**2
+
+c Update f1 and f2.
+
+c temporarily set f1 and f2 for col=0.
+ f1 = f1 + dt*f2 + dibp2 - theta*dibp*zibp
+ f2 = f2 - theta*dibp2
+
+ if (col .gt. 0) then
+c update c = c + dt*p.
+ call daxpy(col2,dt,p,1,c,1)
+
+c choose wbp,
+c the row of W corresponding to the breakpoint encountered.
+ pointr = head
+ do 70 j = 1,col
+ wbp(j) = wy(ibp,pointr)
+ wbp(col + j) = theta*ws(ibp,pointr)
+ pointr = mod(pointr,m) + 1
+ 70 continue
+
+c compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'.
+ call bmv(m,sy,wt,col,wbp,v,info)
+ if (info .ne. 0) return
+ wmc = ddot(col2,c,1,v,1)
+ wmp = ddot(col2,p,1,v,1)
+ wmw = ddot(col2,wbp,1,v,1)
+
+c update p = p - dibp*wbp.
+ call daxpy(col2,-dibp,wbp,1,p,1)
+
+c complete updating f1 and f2 while col > 0.
+ f1 = f1 + dibp*wmc
+ f2 = f2 + 2.0d0*dibp*wmp - dibp2*wmw
+ endif
+
+ f2 = max(epsmch*f2_org,f2)
+ if (nleft .gt. 0) then
+ dtm = -f1/f2
+ goto 777
+c to repeat the loop for unsearched intervals.
+ else if(bnded) then
+ f1 = zero
+ f2 = zero
+ dtm = zero
+ else
+ dtm = -f1/f2
+ endif
+
+c------------------- the end of the loop -------------------------------
+
+ 888 continue
+ if (iprint .ge. 99) then
+ write (6,*)
+ write (6,*) 'GCP found in this segment'
+ write (6,4010) nseg,f1,f2
+ write (6,6010) dtm
+ endif
+ if (dtm .le. zero) dtm = zero
+ tsum = tsum + dtm
+
+c Move free variables (i.e., the ones w/o breakpoints) and
+c the variables whose breakpoints haven't been reached.
+
+ call daxpy(n,tsum,d,1,xcp,1)
+
+ 999 continue
+
+c Update c = c + dtm*p = W'(x^c - x)
+c which will be used in computing r = Z'(B(x^c - x) + g).
+
+ if (col .gt. 0) call daxpy(col2,dtm,p,1,c,1)
+ if (iprint .gt. 100) write (6,1010) (xcp(i),i = 1,n)
+ if (iprint .ge. 99) write (6,2010)
+
+ 1010 format ('Cauchy X = ',/,(4x,1p,6(1x,d11.4)))
+ 2010 format (/,'---------------- exit CAUCHY----------------------',/)
+ 3010 format (/,'---------------- CAUCHY entered-------------------')
+ 4010 format ('Piece ',i3,' --f1, f2 at start point ',1p,2(1x,d11.4))
+ 4011 format (/,'Piece ',i3,' --f1, f2 at start point ',
+ + 1p,2(1x,d11.4))
+ 5010 format ('Distance to the next break point = ',1p,d11.4)
+ 6010 format ('Distance to the stationary point = ',1p,d11.4)
+
+ return
+
+ end
+
+c====================== The end of cauchy ==============================
+
+ subroutine cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index,
+ + theta, col, head, nfree, cnstnd, info)
+
+ logical cnstnd
+ integer n, m, col, head, nfree, info, index(n)
+ double precision theta,
+ + x(n), g(n), z(n), r(n), wa(4*m),
+ + ws(n, m), wy(n, m), sy(m, m), wt(m, m)
+
+c ************
+c
+c Subroutine cmprlb
+c
+c This subroutine computes r=-Z'B(xcp-xk)-Z'g by using
+c wa(2m+1)=W'(xcp-x) from subroutine cauchy.
+c
+c Subprograms called:
+c
+c L-BFGS-B Library ... bmv.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer i,j,k,pointr
+ double precision a1,a2
+
+ if (.not. cnstnd .and. col .gt. 0) then
+ do 26 i = 1, n
+ r(i) = -g(i)
+ 26 continue
+ else
+ do 30 i = 1, nfree
+ k = index(i)
+ r(i) = -theta*(z(k) - x(k)) - g(k)
+ 30 continue
+ call bmv(m,sy,wt,col,wa(2*m+1),wa(1),info)
+ if (info .ne. 0) then
+ info = -8
+ return
+ endif
+ pointr = head
+ do 34 j = 1, col
+ a1 = wa(j)
+ a2 = theta*wa(col + j)
+ do 32 i = 1, nfree
+ k = index(i)
+ r(i) = r(i) + wy(k,pointr)*a1 + ws(k,pointr)*a2
+ 32 continue
+ pointr = mod(pointr,m) + 1
+ 34 continue
+ endif
+
+ return
+
+ end
+
+c======================= The end of cmprlb =============================
+
+ subroutine errclb(n, m, factr, l, u, nbd, task, info, k)
+
+ character*60 task
+ integer n, m, info, k, nbd(n)
+ double precision factr, l(n), u(n)
+
+c ************
+c
+c Subroutine errclb
+c
+c This subroutine checks the validity of the input data.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer i
+ double precision one,zero
+ parameter (one=1.0d0,zero=0.0d0)
+
+c Check the input arguments for errors.
+
+ if (n .le. 0) task = 'ERROR: N .LE. 0'
+ if (m .le. 0) task = 'ERROR: M .LE. 0'
+ if (factr .lt. zero) task = 'ERROR: FACTR .LT. 0'
+
+c Check the validity of the arrays nbd(i), u(i), and l(i).
+
+ do 10 i = 1, n
+ if (nbd(i) .lt. 0 .or. nbd(i) .gt. 3) then
+c return
+ task = 'ERROR: INVALID NBD'
+ info = -6
+ k = i
+ endif
+ if (nbd(i) .eq. 2) then
+ if (l(i) .gt. u(i)) then
+c return
+ task = 'ERROR: NO FEASIBLE SOLUTION'
+ info = -7
+ k = i
+ endif
+ endif
+ 10 continue
+
+ return
+
+ end
+
+c======================= The end of errclb =============================
+
+ subroutine formk(n, nsub, ind, nenter, ileave, indx2, iupdat,
+ + updatd, wn, wn1, m, ws, wy, sy, theta, col,
+ + head, info)
+
+ integer n, nsub, m, col, head, nenter, ileave, iupdat,
+ + info, ind(n), indx2(n)
+ double precision theta, wn(2*m, 2*m), wn1(2*m, 2*m),
+ + ws(n, m), wy(n, m), sy(m, m)
+ logical updatd
+
+c ************
+c
+c Subroutine formk
+c
+c This subroutine forms the LEL^T factorization of the indefinite
+c
+c matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
+c [L_a -R_z theta*S'AA'S ]
+c where E = [-I 0]
+c [ 0 I]
+c The matrix K can be shown to be equal to the matrix M^[-1]N
+c occurring in section 5.1 of [1], as well as to the matrix
+c Mbar^[-1] Nbar in section 5.3.
+c
+c n is an integer variable.
+c On entry n is the dimension of the problem.
+c On exit n is unchanged.
+c
+c nsub is an integer variable
+c On entry nsub is the number of subspace variables in free set.
+c On exit nsub is not changed.
+c
+c ind is an integer array of dimension nsub.
+c On entry ind specifies the indices of subspace variables.
+c On exit ind is unchanged.
+c
+c nenter is an integer variable.
+c On entry nenter is the number of variables entering the
+c free set.
+c On exit nenter is unchanged.
+c
+c ileave is an integer variable.
+c On entry indx2(ileave),...,indx2(n) are the variables leaving
+c the free set.
+c On exit ileave is unchanged.
+c
+c indx2 is an integer array of dimension n.
+c On entry indx2(1),...,indx2(nenter) are the variables entering
+c the free set, while indx2(ileave),...,indx2(n) are the
+c variables leaving the free set.
+c On exit indx2 is unchanged.
+c
+c iupdat is an integer variable.
+c On entry iupdat is the total number of BFGS updates made so far.
+c On exit iupdat is unchanged.
+c
+c updatd is a logical variable.
+c On entry 'updatd' is true if the L-BFGS matrix is updatd.
+c On exit 'updatd' is unchanged.
+c
+c wn is a double precision array of dimension 2m x 2m.
+c On entry wn is unspecified.
+c On exit the upper triangle of wn stores the LEL^T factorization
+c of the 2*col x 2*col indefinite matrix
+c [-D -Y'ZZ'Y/theta L_a'-R_z' ]
+c [L_a -R_z theta*S'AA'S ]
+c
+c wn1 is a double precision array of dimension 2m x 2m.
+c On entry wn1 stores the lower triangular part of
+c [Y' ZZ'Y L_a'+R_z']
+c [L_a+R_z S'AA'S ]
+c in the previous iteration.
+c On exit wn1 stores the corresponding updated matrices.
+c The purpose of wn1 is just to store these inner products
+c so they can be easily updated and inserted into wn.
+c
+c m is an integer variable.
+c On entry m is the maximum number of variable metric corrections
+c used to define the limited memory matrix.
+c On exit m is unchanged.
+c
+c ws, wy, sy, and wtyy are double precision arrays;
+c theta is a double precision variable;
+c col is an integer variable;
+c head is an integer variable.
+c On entry they store the information defining the
+c limited memory BFGS matrix:
+c ws(n,m) stores S, a set of s-vectors;
+c wy(n,m) stores Y, a set of y-vectors;
+c sy(m,m) stores S'Y;
+c wtyy(m,m) stores the Cholesky factorization
+c of (theta*S'S+LD^(-1)L')
+c theta is the scaling factor specifying B_0 = theta I;
+c col is the number of variable metric corrections stored;
+c head is the location of the 1st s- (or y-) vector in S (or Y).
+c On exit they are unchanged.
+c
+c info is an integer variable.
+c On entry info is unspecified.
+c On exit info = 0 for normal return;
+c = -1 when the 1st Cholesky factorization failed;
+c = -2 when the 2st Cholesky factorization failed.
+c
+c Subprograms called:
+c
+c Linpack ... dcopy, dpofa, dtrsl.
+c
+c
+c References:
+c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+c memory algorithm for bound constrained optimization'',
+c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+c
+c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
+c limited memory FORTRAN code for solving bound constrained
+c optimization problems'', Tech. Report, NAM-11, EECS Department,
+c Northwestern University, 1994.
+c
+c (Postscript files of these papers are available via anonymous
+c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer m2,ipntr,jpntr,iy,is,jy,js,is1,js1,k1,i,k,
+ + col2,pbegin,pend,dbegin,dend,upcl
+ double precision ddot,temp1,temp2,temp3,temp4
+ double precision one,zero
+ parameter (one=1.0d0,zero=0.0d0)
+
+c Form the lower triangular part of
+c WN1 = [Y' ZZ'Y L_a'+R_z']
+c [L_a+R_z S'AA'S ]
+c where L_a is the strictly lower triangular part of S'AA'Y
+c R_z is the upper triangular part of S'ZZ'Y.
+
+ if (updatd) then
+ if (iupdat .gt. m) then
+c shift old part of WN1.
+ do 10 jy = 1, m - 1
+ js = m + jy
+ call dcopy(m-jy,wn1(jy+1,jy+1),1,wn1(jy,jy),1)
+ call dcopy(m-jy,wn1(js+1,js+1),1,wn1(js,js),1)
+ call dcopy(m-1,wn1(m+2,jy+1),1,wn1(m+1,jy),1)
+ 10 continue
+ endif
+
+c put new rows in blocks (1,1), (2,1) and (2,2).
+ pbegin = 1
+ pend = nsub
+ dbegin = nsub + 1
+ dend = n
+ iy = col
+ is = m + col
+ ipntr = head + col - 1
+ if (ipntr .gt. m) ipntr = ipntr - m
+ jpntr = head
+ do 20 jy = 1, col
+ js = m + jy
+ temp1 = zero
+ temp2 = zero
+ temp3 = zero
+c compute element jy of row 'col' of Y'ZZ'Y
+ do 15 k = pbegin, pend
+ k1 = ind(k)
+ temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
+ 15 continue
+c compute elements jy of row 'col' of L_a and S'AA'S
+ do 16 k = dbegin, dend
+ k1 = ind(k)
+ temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
+ temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
+ 16 continue
+ wn1(iy,jy) = temp1
+ wn1(is,js) = temp2
+ wn1(is,jy) = temp3
+ jpntr = mod(jpntr,m) + 1
+ 20 continue
+
+c put new column in block (2,1).
+ jy = col
+ jpntr = head + col - 1
+ if (jpntr .gt. m) jpntr = jpntr - m
+ ipntr = head
+ do 30 i = 1, col
+ is = m + i
+ temp3 = zero
+c compute element i of column 'col' of R_z
+ do 25 k = pbegin, pend
+ k1 = ind(k)
+ temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
+ 25 continue
+ ipntr = mod(ipntr,m) + 1
+ wn1(is,jy) = temp3
+ 30 continue
+ upcl = col - 1
+ else
+ upcl = col
+ endif
+
+c modify the old parts in blocks (1,1) and (2,2) due to changes
+c in the set of free variables.
+ ipntr = head
+ do 45 iy = 1, upcl
+ is = m + iy
+ jpntr = head
+ do 40 jy = 1, iy
+ js = m + jy
+ temp1 = zero
+ temp2 = zero
+ temp3 = zero
+ temp4 = zero
+ do 35 k = 1, nenter
+ k1 = indx2(k)
+ temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
+ temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
+ 35 continue
+ do 36 k = ileave, n
+ k1 = indx2(k)
+ temp3 = temp3 + wy(k1,ipntr)*wy(k1,jpntr)
+ temp4 = temp4 + ws(k1,ipntr)*ws(k1,jpntr)
+ 36 continue
+ wn1(iy,jy) = wn1(iy,jy) + temp1 - temp3
+ wn1(is,js) = wn1(is,js) - temp2 + temp4
+ jpntr = mod(jpntr,m) + 1
+ 40 continue
+ ipntr = mod(ipntr,m) + 1
+ 45 continue
+
+c modify the old parts in block (2,1).
+ ipntr = head
+ do 60 is = m + 1, m + upcl
+ jpntr = head
+ do 55 jy = 1, upcl
+ temp1 = zero
+ temp3 = zero
+ do 50 k = 1, nenter
+ k1 = indx2(k)
+ temp1 = temp1 + ws(k1,ipntr)*wy(k1,jpntr)
+ 50 continue
+ do 51 k = ileave, n
+ k1 = indx2(k)
+ temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
+ 51 continue
+ if (is .le. jy + m) then
+ wn1(is,jy) = wn1(is,jy) + temp1 - temp3
+ else
+ wn1(is,jy) = wn1(is,jy) - temp1 + temp3
+ endif
+ jpntr = mod(jpntr,m) + 1
+ 55 continue
+ ipntr = mod(ipntr,m) + 1
+ 60 continue
+
+c Form the upper triangle of WN = [D+Y' ZZ'Y/theta -L_a'+R_z' ]
+c [-L_a +R_z S'AA'S*theta]
+
+ m2 = 2*m
+ do 70 iy = 1, col
+ is = col + iy
+ is1 = m + iy
+ do 65 jy = 1, iy
+ js = col + jy
+ js1 = m + jy
+ wn(jy,iy) = wn1(iy,jy)/theta
+ wn(js,is) = wn1(is1,js1)*theta
+ 65 continue
+ do 66 jy = 1, iy - 1
+ wn(jy,is) = -wn1(is1,jy)
+ 66 continue
+ do 67 jy = iy, col
+ wn(jy,is) = wn1(is1,jy)
+ 67 continue
+ wn(iy,iy) = wn(iy,iy) + sy(iy,iy)
+ 70 continue
+
+c Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')]
+c [(-L_a +R_z)L'^-1 S'AA'S*theta ]
+
+c first Cholesky factor (1,1) block of wn to get LL'
+c with L' stored in the upper triangle of wn.
+ call dpofa(wn,m2,col,info)
+ if (info .ne. 0) then
+ info = -1
+ return
+ endif
+c then form L^-1(-L_a'+R_z') in the (1,2) block.
+ col2 = 2*col
+ do 71 js = col+1 ,col2
+ call dtrsl(wn,m2,col,wn(1,js),11,info)
+ 71 continue
+
+c Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the
+c upper triangle of (2,2) block of wn.
+
+
+ do 72 is = col+1, col2
+ do 74 js = is, col2
+ wn(is,js) = wn(is,js) + ddot(col,wn(1,is),1,wn(1,js),1)
+ 74 continue
+ 72 continue
+
+c Cholesky factorization of (2,2) block of wn.
+
+ call dpofa(wn(col+1,col+1),m2,col,info)
+ if (info .ne. 0) then
+ info = -2
+ return
+ endif
+
+ return
+
+ end
+
+c======================= The end of formk ==============================
+
+ subroutine formt(m, wt, sy, ss, col, theta, info)
+
+ integer m, col, info
+ double precision theta, wt(m, m), sy(m, m), ss(m, m)
+
+c ************
+c
+c Subroutine formt
+c
+c This subroutine forms the upper half of the pos. def. and symm.
+c T = theta*SS + L*D^(-1)*L', stores T in the upper triangle
+c of the array wt, and performs the Cholesky factorization of T
+c to produce J*J', with J' stored in the upper triangle of wt.
+c
+c Subprograms called:
+c
+c Linpack ... dpofa.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer i,j,k,k1
+ double precision ddum
+ double precision zero
+ parameter (zero=0.0d0)
+
+
+c Form the upper half of T = theta*SS + L*D^(-1)*L',
+c store T in the upper triangle of the array wt.
+
+ do 52 j = 1, col
+ wt(1,j) = theta*ss(1,j)
+ 52 continue
+ do 55 i = 2, col
+ do 54 j = i, col
+ k1 = min(i,j) - 1
+ ddum = zero
+ do 53 k = 1, k1
+ ddum = ddum + sy(i,k)*sy(j,k)/sy(k,k)
+ 53 continue
+ wt(i,j) = ddum + theta*ss(i,j)
+ 54 continue
+ 55 continue
+
+c Cholesky factorize T to J*J' with
+c J' stored in the upper triangle of wt.
+
+ call dpofa(wt,m,col,info)
+ if (info .ne. 0) then
+ info = -3
+ endif
+
+ return
+
+ end
+
+c======================= The end of formt ==============================
+
+ subroutine freev(n, nfree, index, nenter, ileave, indx2,
+ + iwhere, wrk, updatd, cnstnd, iprint, iter)
+
+ integer n, nfree, nenter, ileave, iprint, iter,
+ + index(n), indx2(n), iwhere(n)
+ logical wrk, updatd, cnstnd
+
+c ************
+c
+c Subroutine freev
+c
+c This subroutine counts the entering and leaving variables when
+c iter > 0, and finds the index set of free and active variables
+c at the GCP.
+c
+c cnstnd is a logical variable indicating whether bounds are present
+c
+c index is an integer array of dimension n
+c for i=1,...,nfree, index(i) are the indices of free variables
+c for i=nfree+1,...,n, index(i) are the indices of bound variables
+c On entry after the first iteration, index gives
+c the free variables at the previous iteration.
+c On exit it gives the free variables based on the determination
+c in cauchy using the array iwhere.
+c
+c indx2 is an integer array of dimension n
+c On entry indx2 is unspecified.
+c On exit with iter>0, indx2 indicates which variables
+c have changed status since the previous iteration.
+c For i= 1,...,nenter, indx2(i) have changed from bound to free.
+c For i= ileave+1,...,n, indx2(i) have changed from free to bound.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer iact,i,k
+
+ nenter = 0
+ ileave = n + 1
+ if (iter .gt. 0 .and. cnstnd) then
+c count the entering and leaving variables.
+ do 20 i = 1, nfree
+ k = index(i)
+
+c write(6,*) ' k = index(i) ', k
+c write(6,*) ' index = ', i
+
+ if (iwhere(k) .gt. 0) then
+ ileave = ileave - 1
+ indx2(ileave) = k
+ if (iprint .ge. 100) write (6,*)
+ + 'Variable ',k,' leaves the set of free variables'
+ endif
+ 20 continue
+ do 22 i = 1 + nfree, n
+ k = index(i)
+ if (iwhere(k) .le. 0) then
+ nenter = nenter + 1
+ indx2(nenter) = k
+ if (iprint .ge. 100) write (6,*)
+ + 'Variable ',k,' enters the set of free variables'
+ endif
+ 22 continue
+ if (iprint .ge. 99) write (6,*)
+ + n+1-ileave,' variables leave; ',nenter,' variables enter'
+ endif
+ wrk = (ileave .lt. n+1) .or. (nenter .gt. 0) .or. updatd
+
+c Find the index set of free and active variables at the GCP.
+
+ nfree = 0
+ iact = n + 1
+ do 24 i = 1, n
+ if (iwhere(i) .le. 0) then
+ nfree = nfree + 1
+ index(nfree) = i
+ else
+ iact = iact - 1
+ index(iact) = i
+ endif
+ 24 continue
+ if (iprint .ge. 99) write (6,*)
+ + nfree,' variables are free at GCP ',iter + 1
+
+ return
+
+ end
+
+c======================= The end of freev ==============================
+
+ subroutine hpsolb(n, t, iorder, iheap)
+ integer iheap, n, iorder(n)
+ double precision t(n)
+
+c ************
+c
+c Subroutine hpsolb
+c
+c This subroutine sorts out the least element of t, and puts the
+c remaining elements of t in a heap.
+c
+c n is an integer variable.
+c On entry n is the dimension of the arrays t and iorder.
+c On exit n is unchanged.
+c
+c t is a double precision array of dimension n.
+c On entry t stores the elements to be sorted,
+c On exit t(n) stores the least elements of t, and t(1) to t(n-1)
+c stores the remaining elements in the form of a heap.
+c
+c iorder is an integer array of dimension n.
+c On entry iorder(i) is the index of t(i).
+c On exit iorder(i) is still the index of t(i), but iorder may be
+c permuted in accordance with t.
+c
+c iheap is an integer variable specifying the task.
+c On entry iheap should be set as follows:
+c iheap .eq. 0 if t(1) to t(n) is not in the form of a heap,
+c iheap .ne. 0 if otherwise.
+c On exit iheap is unchanged.
+c
+c
+c References:
+c Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT.
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c ************
+
+ integer i,j,k,indxin,indxou
+ double precision ddum,out
+
+ if (iheap .eq. 0) then
+
+c Rearrange the elements t(1) to t(n) to form a heap.
+
+ do 20 k = 2, n
+ ddum = t(k)
+ indxin = iorder(k)
+
+c Add ddum to the heap.
+ i = k
+ 10 continue
+ if (i.gt.1) then
+ j = i/2
+ if (ddum .lt. t(j)) then
+ t(i) = t(j)
+ iorder(i) = iorder(j)
+ i = j
+ goto 10
+ endif
+ endif
+ t(i) = ddum
+ iorder(i) = indxin
+ 20 continue
+ endif
+
+c Assign to 'out' the value of t(1), the least member of the heap,
+c and rearrange the remaining members to form a heap as
+c elements 1 to n-1 of t.
+
+ if (n .gt. 1) then
+ i = 1
+ out = t(1)
+ indxou = iorder(1)
+ ddum = t(n)
+ indxin = iorder(n)
+
+c Restore the heap
+ 30 continue
+ j = i+i
+ if (j .le. n-1) then
+ if (t(j+1) .lt. t(j)) j = j+1
+ if (t(j) .lt. ddum ) then
+ t(i) = t(j)
+ iorder(i) = iorder(j)
+ i = j
+ goto 30
+ endif
+ endif
+ t(i) = ddum
+ iorder(i) = indxin
+
+c Put the least member in t(n).
+
+ t(n) = out
+ iorder(n) = indxou
+ endif
+
+ return
+
+ end
+
+c====================== The end of hpsolb ==============================
+
+ subroutine lnsrlb(n, l, u, nbd, x, f, fold, gd, gdold, g, d, r, t,
+ + z, stp, dnorm, dtd, xstep, stpmx, iter, ifun,
+ + iback, nfgv, info, task, boxed, cnstnd, csave,
+ + isave, dsave)
+
+ character*60 task, csave
+ logical boxed, cnstnd
+ integer n, iter, ifun, iback, nfgv, info,
+ + nbd(n), isave(2)
+ double precision f, fold, gd, gdold, stp, dnorm, dtd, xstep,
+ + stpmx, x(n), l(n), u(n), g(n), d(n), r(n), t(n),
+ + z(n), dsave(13)
+c **********
+c
+c Subroutine lnsrlb
+c
+c This subroutine calls subroutine dcsrch from the Minpack2 library
+c to perform the line search. Subroutine dscrch is safeguarded so
+c that all trial points lie within the feasible region.
+c
+c Subprograms called:
+c
+c Minpack2 Library ... dcsrch.
+c
+c Linpack ... dtrsl, ddot.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c **********
+
+ integer i
+ double precision ddot,a1,a2
+ double precision one,zero,big
+ parameter (one=1.0d0,zero=0.0d0,big=1.0d+10)
+ double precision ftol,gtol,xtol
+ parameter (ftol=1.0d-3,gtol=0.9d0,xtol=0.1d0)
+
+ if (task(1:5) .eq. 'FG_LN') goto 556
+
+ dtd = ddot(n,d,1,d,1)
+ dnorm = sqrt(dtd)
+
+c Determine the maximum step length.
+
+ stpmx = big
+ if (cnstnd) then
+ if (iter .eq. 0) then
+ stpmx = one
+ else
+ do 43 i = 1, n
+ a1 = d(i)
+ if (nbd(i) .ne. 0) then
+ if (a1 .lt. zero .and. nbd(i) .le. 2) then
+ a2 = l(i) - x(i)
+ if (a2 .ge. zero) then
+ stpmx = zero
+ else if (a1*stpmx .lt. a2) then
+ stpmx = a2/a1
+ endif
+ else if (a1 .gt. zero .and. nbd(i) .ge. 2) then
+ a2 = u(i) - x(i)
+ if (a2 .le. zero) then
+ stpmx = zero
+ else if (a1*stpmx .gt. a2) then
+ stpmx = a2/a1
+ endif
+ endif
+ endif
+ 43 continue
+ endif
+ endif
+
+ if (iter .eq. 0 .and. .not. boxed) then
+ stp = min(one/dnorm, stpmx)
+ else
+ stp = one
+ endif
+
+ call dcopy(n,x,1,t,1)
+ call dcopy(n,g,1,r,1)
+ fold = f
+ ifun = 0
+ iback = 0
+ csave = 'START'
+ 556 continue
+ gd = ddot(n,g,1,d,1)
+ if (ifun .eq. 0) then
+ gdold=gd
+ if (gd .ge. zero) then
+c the directional derivative >=0.
+c Line search is impossible.
+ write(6,*)' ascent direction in projection gd = ', gd
+ info = -4
+ return
+ endif
+ endif
+
+ call dcsrch(f,gd,stp,ftol,gtol,xtol,zero,stpmx,csave,isave,dsave)
+
+ xstep = stp*dnorm
+ if (csave(1:4) .ne. 'CONV' .and. csave(1:4) .ne. 'WARN') then
+ task = 'FG_LNSRCH'
+ ifun = ifun + 1
+ nfgv = nfgv + 1
+ iback = ifun - 1
+ if (stp .eq. one) then
+ call dcopy(n,z,1,x,1)
+ else
+ do 41 i = 1, n
+ x(i) = stp*d(i) + t(i)
+ 41 continue
+ endif
+ else
+ task = 'NEW_X'
+ endif
+
+ return
+
+ end
+
+c======================= The end of lnsrlb =============================
+
+ subroutine matupd(n, m, ws, wy, sy, ss, d, r, itail,
+ + iupdat, col, head, theta, rr, dr, stp, dtd)
+
+ integer n, m, itail, iupdat, col, head
+ double precision theta, rr, dr, stp, dtd, d(n), r(n),
+ + ws(n, m), wy(n, m), sy(m, m), ss(m, m)
+
+c ************
+c
+c Subroutine matupd
+c
+c This subroutine updates matrices WS and WY, and forms the
+c middle matrix in B.
+c
+c Subprograms called:
+c
+c Linpack ... dcopy, ddot.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer j,pointr
+ double precision ddot
+ double precision one
+ parameter (one=1.0d0)
+
+c Set pointers for matrices WS and WY.
+
+ if (iupdat .le. m) then
+ col = iupdat
+ itail = mod(head+iupdat-2,m) + 1
+ else
+ itail = mod(itail,m) + 1
+ head = mod(head,m) + 1
+ endif
+
+c Update matrices WS and WY.
+
+ call dcopy(n,d,1,ws(1,itail),1)
+ call dcopy(n,r,1,wy(1,itail),1)
+
+c Set theta=yy/ys.
+
+ theta = rr/dr
+
+c Form the middle matrix in B.
+
+c update the upper triangle of SS,
+c and the lower triangle of SY:
+ if (iupdat .gt. m) then
+c move old information
+ do 50 j = 1, col - 1
+ call dcopy(j,ss(2,j+1),1,ss(1,j),1)
+ call dcopy(col-j,sy(j+1,j+1),1,sy(j,j),1)
+ 50 continue
+ endif
+c add new information: the last row of SY
+c and the last column of SS:
+ pointr = head
+ do 51 j = 1, col - 1
+ sy(col,j) = ddot(n,d,1,wy(1,pointr),1)
+ ss(j,col) = ddot(n,ws(1,pointr),1,d,1)
+ pointr = mod(pointr,m) + 1
+ 51 continue
+ if (stp .eq. one) then
+ ss(col,col) = dtd
+ else
+ ss(col,col) = stp*stp*dtd
+ endif
+ sy(col,col) = dr
+
+ return
+
+ end
+
+c======================= The end of matupd =============================
+
+ subroutine prn1lb(n, m, l, u, x, iprint, itfile, epsmch)
+
+ integer n, m, iprint, itfile
+ double precision epsmch, x(n), l(n), u(n)
+
+c ************
+c
+c Subroutine prn1lb
+c
+c This subroutine prints the input data, initial point, upper and
+c lower bounds of each variable, machine precision, as well as
+c the headings of the output.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer i
+
+ if (iprint .ge. 0) then
+ write (6,7001) epsmch
+ write (6,*) 'N = ',n,' M = ',m
+ if (iprint .ge. 1) then
+ write (itfile,2001) epsmch
+ write (itfile,*)'N = ',n,' M = ',m
+ write (itfile,9001)
+ if (iprint .gt. 100) then
+ write (6,1004) 'L =',(l(i),i = 1,n)
+ write (6,1004) 'X0 =',(x(i),i = 1,n)
+ write (6,1004) 'U =',(u(i),i = 1,n)
+ endif
+ endif
+ endif
+
+ 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))
+ 2001 format ('RUNNING THE L-BFGS-B CODE',/,/,
+ + 'it = iteration number',/,
+ + 'nf = number of function evaluations',/,
+ + 'nseg = number of segments explored during the Cauchy search',/,
+ + 'nact = number of active bounds at the generalized Cauchy point'
+ + ,/,
+ + 'sub = manner in which the subspace minimization terminated:'
+ + ,/,' con = converged, bnd = a bound was reached',/,
+ + 'itls = number of iterations performed in the line search',/,
+ + 'stepl = step length used',/,
+ + 'tstep = norm of the displacement (total step)',/,
+ + 'projg = norm of the projected gradient',/,
+ + 'f = function value',/,/,
+ + ' * * *',/,/,
+ + 'Machine precision =',1p,d10.3)
+ 7001 format ('RUNNING THE L-BFGS-B CODE',/,/,
+ + ' * * *',/,/,
+ + 'Machine precision =',1p,d10.3)
+ 9001 format (/,3x,'it',3x,'nf',2x,'nseg',2x,'nact',2x,'sub',2x,'itls',
+ + 2x,'stepl',4x,'tstep',5x,'projg',8x,'f')
+
+ return
+
+ end
+
+c======================= The end of prn1lb =============================
+
+ subroutine prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact,
+ + sbgnrm, nseg, word, iword, iback, stp, xstep)
+
+ character*3 word
+ integer n, iprint, itfile, iter, nfgv, nact, nseg,
+ + iword, iback
+ double precision f, sbgnrm, stp, xstep, x(n), g(n)
+
+c ************
+c
+c Subroutine prn2lb
+c
+c This subroutine prints out new information after a successful
+c line search.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer i,imod
+
+c 'word' records the status of subspace solutions.
+ if (iword .eq. 0) then
+c the subspace minimization converged.
+ word = 'con'
+ else if (iword .eq. 1) then
+c the subspace minimization stopped at a bound.
+ word = 'bnd'
+ else if (iword .eq. 5) then
+c the truncated Newton step has been used.
+ word = 'TNT'
+ else
+ word = '---'
+ endif
+ if (iprint .ge. 99) then
+ write (6,*) 'LINE SEARCH',iback,' times; norm of step = ',xstep
+ write (6,2001) iter,f,sbgnrm
+ if (iprint .gt. 100) then
+ write (6,1004) 'X =',(x(i), i = 1, n)
+ write (6,1004) 'G =',(g(i), i = 1, n)
+ endif
+ else if (iprint .gt. 0) then
+ imod = mod(iter,iprint)
+ if (imod .eq. 0) write (6,2001) iter,f,sbgnrm
+ endif
+ if (iprint .ge. 1) write (itfile,3001)
+ + iter,nfgv,nseg,nact,word,iback,stp,xstep,sbgnrm,f
+
+ 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))
+ 2001 format
+ + (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5)
+ 3001 format(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),1p,2(1x,d10.3))
+
+ return
+
+ end
+
+c======================= The end of prn2lb =============================
+
+ subroutine prn3lb(n, x, f, task, iprint, info, itfile,
+ + iter, nfgv, nintol, nskip, nact, sbgnrm,
+ + time, nseg, word, iback, stp, xstep, k,
+ + cachyt, sbtime, lnscht)
+
+ character*60 task
+ character*3 word
+ integer n, iprint, info, itfile, iter, nfgv, nintol,
+ + nskip, nact, nseg, iback, k
+ double precision f, sbgnrm, time, stp, xstep, cachyt, sbtime,
+ + lnscht, x(n)
+
+c ************
+c
+c Subroutine prn3lb
+c
+c This subroutine prints out information when either a built-in
+c convergence test is satisfied or when an error message is
+c generated.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer i
+
+ if (task(1:5) .eq. 'ERROR') goto 999
+
+ if (iprint .ge. 0) then
+ write (6,3003)
+ write (6,3004)
+ write(6,3005) n,iter,nfgv,nintol,nskip,nact,sbgnrm,f
+ if (iprint .ge. 100) then
+ write (6,1004) 'X =',(x(i),i = 1,n)
+ endif
+ if (iprint .ge. 1) write (6,*) ' F =',f
+ endif
+ 999 continue
+ if (iprint .ge. 0) then
+ write (6,3009) task
+ if (info .ne. 0) then
+ if (info .eq. -1) write (6,9011)
+ if (info .eq. -2) write (6,9012)
+ if (info .eq. -3) write (6,9013)
+ if (info .eq. -4) write (6,9014)
+ if (info .eq. -5) write (6,9015)
+ if (info .eq. -6) write (6,*)' Input nbd(',k,') is invalid.'
+ if (info .eq. -7)
+ + write (6,*)' l(',k,') > u(',k,'). No feasible solution.'
+ if (info .eq. -8) write (6,9018)
+ if (info .eq. -9) write (6,9019)
+ endif
+ if (iprint .ge. 1) write (6,3007) cachyt,sbtime,lnscht
+ write (6,3008) time
+ if (iprint .ge. 1) then
+ if (info .eq. -4 .or. info .eq. -9) then
+ write (itfile,3002)
+ + iter,nfgv,nseg,nact,word,iback,stp,xstep
+ endif
+ write (itfile,3009) task
+ if (info .ne. 0) then
+ if (info .eq. -1) write (itfile,9011)
+ if (info .eq. -2) write (itfile,9012)
+ if (info .eq. -3) write (itfile,9013)
+ if (info .eq. -4) write (itfile,9014)
+ if (info .eq. -5) write (itfile,9015)
+ if (info .eq. -8) write (itfile,9018)
+ if (info .eq. -9) write (itfile,9019)
+ endif
+ write (itfile,3008) time
+ endif
+ endif
+
+ 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))
+ 3002 format(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),6x,'-',10x,'-')
+ 3003 format (/,
+ + ' * * *',/,/,
+ + 'Tit = total number of iterations',/,
+ + 'Tnf = total number of function evaluations',/,
+ + 'Tnint = total number of segments explored during',
+ + ' Cauchy searches',/,
+ + 'Skip = number of BFGS updates skipped',/,
+ + 'Nact = number of active bounds at final generalized',
+ + ' Cauchy point',/,
+ + 'Projg = norm of the final projected gradient',/,
+ + 'F = final function value',/,/,
+ + ' * * *')
+ 3004 format (/,3x,'N',4x,'Tit',5x,'Tnf',2x,'Tnint',2x,
+ + 'Skip',2x,'Nact',5x,'Projg',8x,'F')
+ 3005 format (i5,2(1x,i6),(1x,i6),(2x,i4),(1x,i5),1p,2(2x,d10.3))
+ 3007 format (/,' Cauchy time',1p,e10.3,' seconds.',/
+ + ' Subspace minimization time',1p,e10.3,' seconds.',/
+ + ' Line search time',1p,e10.3,' seconds.')
+ 3008 format (/,' Total User time',1p,e10.3,' seconds.',/)
+ 3009 format (/,a60)
+ 9011 format (/,
+ +' Matrix in 1st Cholesky factorization in formk is not Pos. Def.')
+ 9012 format (/,
+ +' Matrix in 2st Cholesky factorization in formk is not Pos. Def.')
+ 9013 format (/,
+ +' Matrix in the Cholesky factorization in formt is not Pos. Def.')
+ 9014 format (/,
+ +' Derivative >= 0, backtracking line search impossible.',/,
+ +' Previous x, f and g restored.',/,
+ +' Possible causes: 1 error in function or gradient evaluation;',/,
+ +' 2 rounding errors dominate computation.')
+ 9015 format (/,
+ +' Warning: more than 10 function and gradient',/,
+ +' evaluations in the last line search. Termination',/,
+ +' may possibly be caused by a bad search direction.')
+ 9018 format (/,' The triangular system is singular.')
+ 9019 format (/,
+ +' Line search cannot locate an adequate point after 20 function',/
+ +,' and gradient evaluations. Previous x, f and g restored.',/,
+ +' Possible causes: 1 error in function or gradient evaluation;',/,
+ +' 2 rounding error dominate computation.')
+
+ return
+
+ end
+
+c======================= The end of prn3lb =============================
+
+ subroutine projgr(n, l, u, nbd, x, g, sbgnrm)
+
+ integer n, nbd(n)
+ double precision sbgnrm, x(n), l(n), u(n), g(n)
+
+c ************
+c
+c Subroutine projgr
+c
+c This subroutine computes the infinity norm of the projected
+c gradient.
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer i
+ double precision gi
+ double precision one,zero
+ parameter (one=1.0d0,zero=0.0d0)
+
+ sbgnrm = zero
+ do 15 i = 1, n
+ gi = g(i)
+ if (nbd(i) .ne. 0) then
+ if (gi .lt. zero) then
+ if (nbd(i) .ge. 2) gi = max((x(i)-u(i)),gi)
+ else
+ if (nbd(i) .le. 2) gi = min((x(i)-l(i)),gi)
+ endif
+ endif
+ sbgnrm = max(sbgnrm,abs(gi))
+ 15 continue
+
+ return
+
+ end
+
+c======================= The end of projgr =============================
+
+ subroutine subsm ( n, m, nsub, ind, l, u, nbd, x, d, xp, ws, wy,
+ + theta, xx, gg,
+ + col, head, iword, wv, wn, iprint, info )
+ implicit none
+ integer n, m, nsub, col, head, iword, iprint, info,
+ + ind(nsub), nbd(n)
+ double precision theta,
+ + l(n), u(n), x(n), d(n), xp(n), xx(n), gg(n),
+ + ws(n, m), wy(n, m),
+ + wv(2*m), wn(2*m, 2*m)
+
+c **********************************************************************
+c
+c This routine contains the major changes in the updated version.
+c The changes are described in the accompanying paper
+c
+c Jose Luis Morales, Jorge Nocedal
+c "Remark On Algorithm 788: L-BFGS-B: Fortran Subroutines for Large-Scale
+c Bound Constrained Optimization". Decemmber 27, 2010.
+c
+c J.L. Morales Departamento de Matematicas,
+c Instituto Tecnologico Autonomo de Mexico
+c Mexico D.F.
+c
+c J, Nocedal Department of Electrical Engineering and
+c Computer Science.
+c Northwestern University. Evanston, IL. USA
+c
+c January 17, 2011
+c
+c **********************************************************************
+c
+c
+c Subroutine subsm
+c
+c Given xcp, l, u, r, an index set that specifies
+c the active set at xcp, and an l-BFGS matrix B
+c (in terms of WY, WS, SY, WT, head, col, and theta),
+c this subroutine computes an approximate solution
+c of the subspace problem
+c
+c (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)
+c
+c subject to l<=x<=u
+c x_i=xcp_i for all i in A(xcp)
+c
+c along the subspace unconstrained Newton direction
+c
+c d = -(Z'BZ)^(-1) r.
+c
+c The formula for the Newton direction, given the L-BFGS matrix
+c and the Sherman-Morrison formula, is
+c
+c d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.
+c
+c where
+c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
+c [L_a -R_z theta*S'AA'S ]
+c
+c Note that this procedure for computing d differs
+c from that described in [1]. One can show that the matrix K is
+c equal to the matrix M^[-1]N in that paper.
+c
+c n is an integer variable.
+c On entry n is the dimension of the problem.
+c On exit n is unchanged.
+c
+c m is an integer variable.
+c On entry m is the maximum number of variable metric corrections
+c used to define the limited memory matrix.
+c On exit m is unchanged.
+c
+c nsub is an integer variable.
+c On entry nsub is the number of free variables.
+c On exit nsub is unchanged.
+c
+c ind is an integer array of dimension nsub.
+c On entry ind specifies the coordinate indices of free variables.
+c On exit ind is unchanged.
+c
+c l is a double precision array of dimension n.
+c On entry l is the lower bound of x.
+c On exit l is unchanged.
+c
+c u is a double precision array of dimension n.
+c On entry u is the upper bound of x.
+c On exit u is unchanged.
+c
+c nbd is a integer array of dimension n.
+c On entry nbd represents the type of bounds imposed on the
+c variables, and must be specified as follows:
+c nbd(i)=0 if x(i) is unbounded,
+c 1 if x(i) has only a lower bound,
+c 2 if x(i) has both lower and upper bounds, and
+c 3 if x(i) has only an upper bound.
+c On exit nbd is unchanged.
+c
+c x is a double precision array of dimension n.
+c On entry x specifies the Cauchy point xcp.
+c On exit x(i) is the minimizer of Q over the subspace of
+c free variables.
+c
+c d is a double precision array of dimension n.
+c On entry d is the reduced gradient of Q at xcp.
+c On exit d is the Newton direction of Q.
+c
+c xp is a double precision array of dimension n.
+c used to safeguard the projected Newton direction
+c
+c xx is a double precision array of dimension n
+c On entry it holds the current iterate
+c On output it is unchanged
+
+c gg is a double precision array of dimension n
+c On entry it holds the gradient at the current iterate
+c On output it is unchanged
+c
+c ws and wy are double precision arrays;
+c theta is a double precision variable;
+c col is an integer variable;
+c head is an integer variable.
+c On entry they store the information defining the
+c limited memory BFGS matrix:
+c ws(n,m) stores S, a set of s-vectors;
+c wy(n,m) stores Y, a set of y-vectors;
+c theta is the scaling factor specifying B_0 = theta I;
+c col is the number of variable metric corrections stored;
+c head is the location of the 1st s- (or y-) vector in S (or Y).
+c On exit they are unchanged.
+c
+c iword is an integer variable.
+c On entry iword is unspecified.
+c On exit iword specifies the status of the subspace solution.
+c iword = 0 if the solution is in the box,
+c 1 if some bound is encountered.
+c
+c wv is a double precision working array of dimension 2m.
+c
+c wn is a double precision array of dimension 2m x 2m.
+c On entry the upper triangle of wn stores the LEL^T factorization
+c of the indefinite matrix
+c
+c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
+c [L_a -R_z theta*S'AA'S ]
+c where E = [-I 0]
+c [ 0 I]
+c On exit wn is unchanged.
+c
+c iprint is an INTEGER variable that must be set by the user.
+c It controls the frequency and type of output generated:
+c iprint<0 no output is generated;
+c iprint=0 print only one line at the last iteration;
+c 0<iprint<99 print also f and |proj g| every iprint iterations;
+c iprint=99 print details of every iteration except n-vectors;
+c iprint=100 print also the changes of active set and final x;
+c iprint>100 print details of every iteration including x and g;
+c When iprint > 0, the file iterate.dat will be created to
+c summarize the iteration.
+c
+c info is an integer variable.
+c On entry info is unspecified.
+c On exit info = 0 for normal return,
+c = nonzero for abnormal return
+c when the matrix K is ill-conditioned.
+c
+c Subprograms called:
+c
+c Linpack dtrsl.
+c
+c
+c References:
+c
+c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
+c memory algorithm for bound constrained optimization'',
+c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
+c
+c
+c
+c * * *
+c
+c NEOS, November 1994. (Latest revision June 1996.)
+c Optimization Technology Center.
+c Argonne National Laboratory and Northwestern University.
+c Written by
+c Ciyou Zhu
+c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
+c
+c
+c ************
+
+ integer pointr,m2,col2,ibd,jy,js,i,j,k
+ double precision alpha, xk, dk, temp1, temp2
+ double precision one,zero
+ parameter (one=1.0d0,zero=0.0d0)
+c
+ double precision dd_p
+
+ if (nsub .le. 0) return
+ if (iprint .ge. 99) write (6,1001)
+
+c Compute wv = W'Zd.
+
+ pointr = head
+ do 20 i = 1, col
+ temp1 = zero
+ temp2 = zero
+ do 10 j = 1, nsub
+ k = ind(j)
+ temp1 = temp1 + wy(k,pointr)*d(j)
+ temp2 = temp2 + ws(k,pointr)*d(j)
+ 10 continue
+ wv(i) = temp1
+ wv(col + i) = theta*temp2
+ pointr = mod(pointr,m) + 1
+ 20 continue
+
+c Compute wv:=K^(-1)wv.
+
+ m2 = 2*m
+ col2 = 2*col
+ call dtrsl(wn,m2,col2,wv,11,info)
+ if (info .ne. 0) return
+ do 25 i = 1, col
+ wv(i) = -wv(i)
+ 25 continue
+ call dtrsl(wn,m2,col2,wv,01,info)
+ if (info .ne. 0) return
+
+c Compute d = (1/theta)d + (1/theta**2)Z'W wv.
+
+ pointr = head
+ do 40 jy = 1, col
+ js = col + jy
+ do 30 i = 1, nsub
+ k = ind(i)
+ d(i) = d(i) + wy(k,pointr)*wv(jy)/theta
+ + + ws(k,pointr)*wv(js)
+ 30 continue
+ pointr = mod(pointr,m) + 1
+ 40 continue
+
+ call dscal( nsub, one/theta, d, 1 )
+c
+c-----------------------------------------------------------------
+c Let us try the projection, d is the Newton direction
+
+ iword = 0
+
+ call dcopy ( n, x, 1, xp, 1 )
+c
+ do 50 i=1, nsub
+ k = ind(i)
+ dk = d(i)
+ xk = x(k)
+ if ( nbd(k) .ne. 0 ) then
+c
+ if ( nbd(k).eq.1 ) then ! lower bounds only
+ x(k) = max( l(k), xk + dk )
+ if ( x(k).eq.l(k) ) iword = 1
+ else
+c
+ if ( nbd(k).eq.2 ) then ! upper and lower bounds
+ xk = max( l(k), xk + dk )
+ x(k) = min( u(k), xk )
+ if ( x(k).eq.l(k) .or. x(k).eq.u(k) ) iword = 1
+ else
+c
+ if ( nbd(k).eq.3 ) then ! upper bounds only
+ x(k) = min( u(k), xk + dk )
+ if ( x(k).eq.u(k) ) iword = 1
+ end if
+ end if
+ end if
+c
+ else ! free variables
+ x(k) = xk + dk
+ end if
+ 50 continue
+c
+ if ( iword.eq.0 ) then
+ go to 911
+ end if
+c
+c check sign of the directional derivative
+c
+ dd_p = zero
+ do 55 i=1, n
+ dd_p = dd_p + (x(i) - xx(i))*gg(i)
+ 55 continue
+ if ( dd_p .gt.zero ) then
+ call dcopy( n, xp, 1, x, 1 )
+ write(6,*) ' Positive dir derivative in projection '
+ write(6,*) ' Using the backtracking step '
+ else
+ go to 911
+ endif
+c
+c-----------------------------------------------------------------
+c
+ alpha = one
+ temp1 = alpha
+ ibd = 0
+ do 60 i = 1, nsub
+ k = ind(i)
+ dk = d(i)
+ if (nbd(k) .ne. 0) then
+ if (dk .lt. zero .and. nbd(k) .le. 2) then
+ temp2 = l(k) - x(k)
+ if (temp2 .ge. zero) then
+ temp1 = zero
+ else if (dk*alpha .lt. temp2) then
+ temp1 = temp2/dk
+ endif
+ else if (dk .gt. zero .and. nbd(k) .ge. 2) then
+ temp2 = u(k) - x(k)
+ if (temp2 .le. zero) then
+ temp1 = zero
+ else if (dk*alpha .gt. temp2) then
+ temp1 = temp2/dk
+ endif
+ endif
+ if (temp1 .lt. alpha) then
+ alpha = temp1
+ ibd = i
+ endif
+ endif
+ 60 continue
+
+ if (alpha .lt. one) then
+ dk = d(ibd)
+ k = ind(ibd)
+ if (dk .gt. zero) then
+ x(k) = u(k)
+ d(ibd) = zero
+ else if (dk .lt. zero) then
+ x(k) = l(k)
+ d(ibd) = zero
+ endif
+ endif
+ do 70 i = 1, nsub
+ k = ind(i)
+ x(k) = x(k) + alpha*d(i)
+ 70 continue
+cccccc
+ 911 continue
+
+ if (iprint .ge. 99) write (6,1004)
+
+ 1001 format (/,'----------------SUBSM entered-----------------',/)
+ 1004 format (/,'----------------exit SUBSM --------------------',/)
+
+ return
+
+ end
+c====================== The end of subsm ===============================
+
+ subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
+ + task,isave,dsave)
+ character*(*) task
+ integer isave(2)
+ double precision f,g,stp,ftol,gtol,xtol,stpmin,stpmax
+ double precision dsave(13)
+c **********
+c
+c Subroutine dcsrch
+c
+c This subroutine finds a step that satisfies a sufficient
+c decrease condition and a curvature condition.
+c
+c Each call of the subroutine updates an interval with
+c endpoints stx and sty. The interval is initially chosen
+c so that it contains a minimizer of the modified function
+c
+c psi(stp) = f(stp) - f(0) - ftol*stp*f'(0).
+c
+c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
+c interval is chosen so that it contains a minimizer of f.
+c
+c The algorithm is designed to find a step that satisfies
+c the sufficient decrease condition
+c
+c f(stp) <= f(0) + ftol*stp*f'(0),
+c
+c and the curvature condition
+c
+c abs(f'(stp)) <= gtol*abs(f'(0)).
+c
+c If ftol is less than gtol and if, for example, the function
+c is bounded below, then there is always a step which satisfies
+c both conditions.
+c
+c If no step can be found that satisfies both conditions, then
+c the algorithm stops with a warning. In this case stp only
+c satisfies the sufficient decrease condition.
+c
+c A typical invocation of dcsrch has the following outline:
+c
+c task = 'START'
+c 10 continue
+c call dcsrch( ... )
+c if (task .eq. 'FG') then
+c Evaluate the function and the gradient at stp
+c goto 10
+c end if
+c
+c NOTE: The user must no alter work arrays between calls.
+c
+c The subroutine statement is
+c
+c subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
+c task,isave,dsave)
+c where
+c
+c f is a double precision variable.
+c On initial entry f is the value of the function at 0.
+c On subsequent entries f is the value of the
+c function at stp.
+c On exit f is the value of the function at stp.
+c
+c g is a double precision variable.
+c On initial entry g is the derivative of the function at 0.
+c On subsequent entries g is the derivative of the
+c function at stp.
+c On exit g is the derivative of the function at stp.
+c
+c stp is a double precision variable.
+c On entry stp is the current estimate of a satisfactory
+c step. On initial entry, a positive initial estimate
+c must be provided.
+c On exit stp is the current estimate of a satisfactory step
+c if task = 'FG'. If task = 'CONV' then stp satisfies
+c the sufficient decrease and curvature condition.
+c
+c ftol is a double precision variable.
+c On entry ftol specifies a nonnegative tolerance for the
+c sufficient decrease condition.
+c On exit ftol is unchanged.
+c
+c gtol is a double precision variable.
+c On entry gtol specifies a nonnegative tolerance for the
+c curvature condition.
+c On exit gtol is unchanged.
+c
+c xtol is a double precision variable.
+c On entry xtol specifies a nonnegative relative tolerance
+c for an acceptable step. The subroutine exits with a
+c warning if the relative difference between sty and stx
+c is less than xtol.
+c On exit xtol is unchanged.
+c
+c stpmin is a double precision variable.
+c On entry stpmin is a nonnegative lower bound for the step.
+c On exit stpmin is unchanged.
+c
+c stpmax is a double precision variable.
+c On entry stpmax is a nonnegative upper bound for the step.
+c On exit stpmax is unchanged.
+c
+c task is a character variable of length at least 60.
+c On initial entry task must be set to 'START'.
+c On exit task indicates the required action:
+c
+c If task(1:2) = 'FG' then evaluate the function and
+c derivative at stp and call dcsrch again.
+c
+c If task(1:4) = 'CONV' then the search is successful.
+c
+c If task(1:4) = 'WARN' then the subroutine is not able
+c to satisfy the convergence conditions. The exit value of
+c stp contains the best point found during the search.
+c
+c If task(1:5) = 'ERROR' then there is an error in the
+c input arguments.
+c
+c On exit with convergence, a warning or an error, the
+c variable task contains additional information.
+c
+c isave is an integer work array of dimension 2.
+c
+c dsave is a double precision work array of dimension 13.
+c
+c Subprograms called
+c
+c MINPACK-2 ... dcstep
+c
+c MINPACK-1 Project. June 1983.
+c Argonne National Laboratory.
+c Jorge J. More' and David J. Thuente.
+c
+c MINPACK-2 Project. October 1993.
+c Argonne National Laboratory and University of Minnesota.
+c Brett M. Averick, Richard G. Carter, and Jorge J. More'.
+c
+c **********
+ double precision zero,p5,p66
+ parameter(zero=0.0d0,p5=0.5d0,p66=0.66d0)
+ double precision xtrapl,xtrapu
+ parameter(xtrapl=1.1d0,xtrapu=4.0d0)
+
+ logical brackt
+ integer stage
+ double precision finit,ftest,fm,fx,fxm,fy,fym,ginit,gtest,
+ + gm,gx,gxm,gy,gym,stx,sty,stmin,stmax,width,width1
+
+c Initialization block.
+
+ if (task(1:5) .eq. 'START') then
+
+c Check the input arguments for errors.
+
+ if (stp .lt. stpmin) task = 'ERROR: STP .LT. STPMIN'
+ if (stp .gt. stpmax) task = 'ERROR: STP .GT. STPMAX'
+ if (g .ge. zero) task = 'ERROR: INITIAL G .GE. ZERO'
+ if (ftol .lt. zero) task = 'ERROR: FTOL .LT. ZERO'
+ if (gtol .lt. zero) task = 'ERROR: GTOL .LT. ZERO'
+ if (xtol .lt. zero) task = 'ERROR: XTOL .LT. ZERO'
+ if (stpmin .lt. zero) task = 'ERROR: STPMIN .LT. ZERO'
+ if (stpmax .lt. stpmin) task = 'ERROR: STPMAX .LT. STPMIN'
+
+c Exit if there are errors on input.
+
+ if (task(1:5) .eq. 'ERROR') return
+
+c Initialize local variables.
+
+ brackt = .false.
+ stage = 1
+ finit = f
+ ginit = g
+ gtest = ftol*ginit
+ width = stpmax - stpmin
+ width1 = width/p5
+
+c The variables stx, fx, gx contain the values of the step,
+c function, and derivative at the best step.
+c The variables sty, fy, gy contain the value of the step,
+c function, and derivative at sty.
+c The variables stp, f, g contain the values of the step,
+c function, and derivative at stp.
+
+ stx = zero
+ fx = finit
+ gx = ginit
+ sty = zero
+ fy = finit
+ gy = ginit
+ stmin = zero
+ stmax = stp + xtrapu*stp
+ task = 'FG'
+
+ goto 1000
+
+ else
+
+c Restore local variables.
+
+ if (isave(1) .eq. 1) then
+ brackt = .true.
+ else
+ brackt = .false.
+ endif
+ stage = isave(2)
+ ginit = dsave(1)
+ gtest = dsave(2)
+ gx = dsave(3)
+ gy = dsave(4)
+ finit = dsave(5)
+ fx = dsave(6)
+ fy = dsave(7)
+ stx = dsave(8)
+ sty = dsave(9)
+ stmin = dsave(10)
+ stmax = dsave(11)
+ width = dsave(12)
+ width1 = dsave(13)
+
+ endif
+
+c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
+c algorithm enters the second stage.
+
+ ftest = finit + stp*gtest
+ if (stage .eq. 1 .and. f .le. ftest .and. g .ge. zero)
+ + stage = 2
+
+c Test for warnings.
+
+ if (brackt .and. (stp .le. stmin .or. stp .ge. stmax))
+ + task = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS'
+ if (brackt .and. stmax - stmin .le. xtol*stmax)
+ + task = 'WARNING: XTOL TEST SATISFIED'
+ if (stp .eq. stpmax .and. f .le. ftest .and. g .le. gtest)
+ + task = 'WARNING: STP = STPMAX'
+ if (stp .eq. stpmin .and. (f .gt. ftest .or. g .ge. gtest))
+ + task = 'WARNING: STP = STPMIN'
+
+c Test for convergence.
+
+ if (f .le. ftest .and. abs(g) .le. gtol*(-ginit))
+ + task = 'CONVERGENCE'
+
+c Test for termination.
+
+ if (task(1:4) .eq. 'WARN' .or. task(1:4) .eq. 'CONV') goto 1000
+
+c A modified function is used to predict the step during the
+c first stage if a lower function value has been obtained but
+c the decrease is not sufficient.
+
+ if (stage .eq. 1 .and. f .le. fx .and. f .gt. ftest) then
+
+c Define the modified function and derivative values.
+
+ fm = f - stp*gtest
+ fxm = fx - stx*gtest
+ fym = fy - sty*gtest
+ gm = g - gtest
+ gxm = gx - gtest
+ gym = gy - gtest
+
+c Call dcstep to update stx, sty, and to compute the new step.
+
+ call dcstep(stx,fxm,gxm,sty,fym,gym,stp,fm,gm,
+ + brackt,stmin,stmax)
+
+c Reset the function and derivative values for f.
+
+ fx = fxm + stx*gtest
+ fy = fym + sty*gtest
+ gx = gxm + gtest
+ gy = gym + gtest
+
+ else
+
+c Call dcstep to update stx, sty, and to compute the new step.
+
+ call dcstep(stx,fx,gx,sty,fy,gy,stp,f,g,
+ + brackt,stmin,stmax)
+
+ endif
+
+c Decide if a bisection step is needed.
+
+ if (brackt) then
+ if (abs(sty-stx) .ge. p66*width1) stp = stx + p5*(sty - stx)
+ width1 = width
+ width = abs(sty-stx)
+ endif
+
+c Set the minimum and maximum steps allowed for stp.
+
+ if (brackt) then
+ stmin = min(stx,sty)
+ stmax = max(stx,sty)
+ else
+ stmin = stp + xtrapl*(stp - stx)
+ stmax = stp + xtrapu*(stp - stx)
+ endif
+
+c Force the step to be within the bounds stpmax and stpmin.
+
+ stp = max(stp,stpmin)
+ stp = min(stp,stpmax)
+
+c If further progress is not possible, let stp be the best
+c point obtained during the search.
+
+ if (brackt .and. (stp .le. stmin .or. stp .ge. stmax)
+ + .or. (brackt .and. stmax-stmin .le. xtol*stmax)) stp = stx
+
+c Obtain another function and derivative.
+
+ task = 'FG'
+
+ 1000 continue
+
+c Save local variables.
+
+ if (brackt) then
+ isave(1) = 1
+ else
+ isave(1) = 0
+ endif
+ isave(2) = stage
+ dsave(1) = ginit
+ dsave(2) = gtest
+ dsave(3) = gx
+ dsave(4) = gy
+ dsave(5) = finit
+ dsave(6) = fx
+ dsave(7) = fy
+ dsave(8) = stx
+ dsave(9) = sty
+ dsave(10) = stmin
+ dsave(11) = stmax
+ dsave(12) = width
+ dsave(13) = width1
+
+ return
+ end
+
+c====================== The end of dcsrch ==============================
+
+ subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,
+ + stpmin,stpmax)
+ logical brackt
+ double precision stx,fx,dx,sty,fy,dy,stp,fp,dp,stpmin,stpmax
+c **********
+c
+c Subroutine dcstep
+c
+c This subroutine computes a safeguarded step for a search
+c procedure and updates an interval that contains a step that
+c satisfies a sufficient decrease and a curvature condition.
+c
+c The parameter stx contains the step with the least function
+c value. If brackt is set to .true. then a minimizer has
+c been bracketed in an interval with endpoints stx and sty.
+c The parameter stp contains the current step.
+c The subroutine assumes that if brackt is set to .true. then
+c
+c min(stx,sty) < stp < max(stx,sty),
+c
+c and that the derivative at stx is negative in the direction
+c of the step.
+c
+c The subroutine statement is
+c
+c subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,
+c stpmin,stpmax)
+c
+c where
+c
+c stx is a double precision variable.
+c On entry stx is the best step obtained so far and is an
+c endpoint of the interval that contains the minimizer.
+c On exit stx is the updated best step.
+c
+c fx is a double precision variable.
+c On entry fx is the function at stx.
+c On exit fx is the function at stx.
+c
+c dx is a double precision variable.
+c On entry dx is the derivative of the function at
+c stx. The derivative must be negative in the direction of
+c the step, that is, dx and stp - stx must have opposite
+c signs.
+c On exit dx is the derivative of the function at stx.
+c
+c sty is a double precision variable.
+c On entry sty is the second endpoint of the interval that
+c contains the minimizer.
+c On exit sty is the updated endpoint of the interval that
+c contains the minimizer.
+c
+c fy is a double precision variable.
+c On entry fy is the function at sty.
+c On exit fy is the function at sty.
+c
+c dy is a double precision variable.
+c On entry dy is the derivative of the function at sty.
+c On exit dy is the derivative of the function at the exit sty.
+c
+c stp is a double precision variable.
+c On entry stp is the current step. If brackt is set to .true.
+c then on input stp must be between stx and sty.
+c On exit stp is a new trial step.
+c
+c fp is a double precision variable.
+c On entry fp is the function at stp
+c On exit fp is unchanged.
+c
+c dp is a double precision variable.
+c On entry dp is the the derivative of the function at stp.
+c On exit dp is unchanged.
+c
+c brackt is an logical variable.
+c On entry brackt specifies if a minimizer has been bracketed.
+c Initially brackt must be set to .false.
+c On exit brackt specifies if a minimizer has been bracketed.
+c When a minimizer is bracketed brackt is set to .true.
+c
+c stpmin is a double precision variable.
+c On entry stpmin is a lower bound for the step.
+c On exit stpmin is unchanged.
+c
+c stpmax is a double precision variable.
+c On entry stpmax is an upper bound for the step.
+c On exit stpmax is unchanged.
+c
+c MINPACK-1 Project. June 1983
+c Argonne National Laboratory.
+c Jorge J. More' and David J. Thuente.
+c
+c MINPACK-2 Project. October 1993.
+c Argonne National Laboratory and University of Minnesota.
+c Brett M. Averick and Jorge J. More'.
+c
+c **********
+ double precision zero,p66,two,three
+ parameter(zero=0.0d0,p66=0.66d0,two=2.0d0,three=3.0d0)
+
+ double precision gamma,p,q,r,s,sgnd,stpc,stpf,stpq,theta
+
+ sgnd = dp*(dx/abs(dx))
+
+c First case: A higher function value. The minimum is bracketed.
+c If the cubic step is closer to stx than the quadratic step, the
+c cubic step is taken, otherwise the average of the cubic and
+c quadratic steps is taken.
+
+ if (fp .gt. fx) then
+ theta = three*(fx - fp)/(stp - stx) + dx + dp
+ s = max(abs(theta),abs(dx),abs(dp))
+ gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp/s))
+ if (stp .lt. stx) gamma = -gamma
+ p = (gamma - dx) + theta
+ q = ((gamma - dx) + gamma) + dp
+ r = p/q
+ stpc = stx + r*(stp - stx)
+ stpq = stx + ((dx/((fx - fp)/(stp - stx) + dx))/two)*
+ + (stp - stx)
+ if (abs(stpc-stx) .lt. abs(stpq-stx)) then
+ stpf = stpc
+ else
+ stpf = stpc + (stpq - stpc)/two
+ endif
+ brackt = .true.
+
+c Second case: A lower function value and derivatives of opposite
+c sign. The minimum is bracketed. If the cubic step is farther from
+c stp than the secant step, the cubic step is taken, otherwise the
+c secant step is taken.
+
+ else if (sgnd .lt. zero) then
+ theta = three*(fx - fp)/(stp - stx) + dx + dp
+ s = max(abs(theta),abs(dx),abs(dp))
+ gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp/s))
+ if (stp .gt. stx) gamma = -gamma
+ p = (gamma - dp) + theta
+ q = ((gamma - dp) + gamma) + dx
+ r = p/q
+ stpc = stp + r*(stx - stp)
+ stpq = stp + (dp/(dp - dx))*(stx - stp)
+ if (abs(stpc-stp) .gt. abs(stpq-stp)) then
+ stpf = stpc
+ else
+ stpf = stpq
+ endif
+ brackt = .true.
+
+c Third case: A lower function value, derivatives of the same sign,
+c and the magnitude of the derivative decreases.
+
+ else if (abs(dp) .lt. abs(dx)) then
+
+c The cubic step is computed only if the cubic tends to infinity
+c in the direction of the step or if the minimum of the cubic
+c is beyond stp. Otherwise the cubic step is defined to be the
+c secant step.
+
+ theta = three*(fx - fp)/(stp - stx) + dx + dp
+ s = max(abs(theta),abs(dx),abs(dp))
+
+c The case gamma = 0 only arises if the cubic does not tend
+c to infinity in the direction of the step.
+
+ gamma = s*sqrt(max(zero,(theta/s)**2-(dx/s)*(dp/s)))
+ if (stp .gt. stx) gamma = -gamma
+ p = (gamma - dp) + theta
+ q = (gamma + (dx - dp)) + gamma
+ r = p/q
+ if (r .lt. zero .and. gamma .ne. zero) then
+ stpc = stp + r*(stx - stp)
+ else if (stp .gt. stx) then
+ stpc = stpmax
+ else
+ stpc = stpmin
+ endif
+ stpq = stp + (dp/(dp - dx))*(stx - stp)
+
+ if (brackt) then
+
+c A minimizer has been bracketed. If the cubic step is
+c closer to stp than the secant step, the cubic step is
+c taken, otherwise the secant step is taken.
+
+ if (abs(stpc-stp) .lt. abs(stpq-stp)) then
+ stpf = stpc
+ else
+ stpf = stpq
+ endif
+ if (stp .gt. stx) then
+ stpf = min(stp+p66*(sty-stp),stpf)
+ else
+ stpf = max(stp+p66*(sty-stp),stpf)
+ endif
+ else
+
+c A minimizer has not been bracketed. If the cubic step is
+c farther from stp than the secant step, the cubic step is
+c taken, otherwise the secant step is taken.
+
+ if (abs(stpc-stp) .gt. abs(stpq-stp)) then
+ stpf = stpc
+ else
+ stpf = stpq
+ endif
+ stpf = min(stpmax,stpf)
+ stpf = max(stpmin,stpf)
+ endif
+
+c Fourth case: A lower function value, derivatives of the same sign,
+c and the magnitude of the derivative does not decrease. If the
+c minimum is not bracketed, the step is either stpmin or stpmax,
+c otherwise the cubic step is taken.
+
+ else
+ if (brackt) then
+ theta = three*(fp - fy)/(sty - stp) + dy + dp
+ s = max(abs(theta),abs(dy),abs(dp))
+ gamma = s*sqrt((theta/s)**2 - (dy/s)*(dp/s))
+ if (stp .gt. sty) gamma = -gamma
+ p = (gamma - dp) + theta
+ q = ((gamma - dp) + gamma) + dy
+ r = p/q
+ stpc = stp + r*(sty - stp)
+ stpf = stpc
+ else if (stp .gt. stx) then
+ stpf = stpmax
+ else
+ stpf = stpmin
+ endif
+ endif
+
+c Update the interval which contains a minimizer.
+
+ if (fp .gt. fx) then
+ sty = stp
+ fy = fp
+ dy = dp
+ else
+ if (sgnd .lt. zero) then
+ sty = stx
+ fy = fx
+ dy = dx
+ endif
+ stx = stp
+ fx = fp
+ dx = dp
+ endif
+
+c Compute the new step.
+
+ stp = stpf
+
+ return
+ end
+
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
+