ACLOCAL_AMFLAGS=-I ./m4
-SUBDIRS= util lbfgs subplex direct cdirect stogo praxis api octave . test
+SUBDIRS= util subplex direct cdirect stogo praxis luksan api octave . test
EXTRA_DIST=COPYRIGHT autogen.sh nlopt.pc.in m4
libnlopt_la_SOURCES =
-libnlopt_la_LIBADD = lbfgs/liblbfgs.la subplex/libsubplex.la \
-direct/libdirect.la cdirect/libcdirect.la stogo/libstogo.la \
-praxis/libpraxis.la api/libapi.la util/libutil.la
+libnlopt_la_LIBADD = subplex/libsubplex.la direct/libdirect.la \
+cdirect/libcdirect.la stogo/libstogo.la praxis/libpraxis.la \
+luksan/libluksan.la api/libapi.la util/libutil.la
libnlopt_la_LDFLAGS = -no-undefined -version-info @SHARED_VERSION_INFO@
-AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/lbfgs -I$(top_srcdir)/subplex -I$(top_srcdir)/praxis -I$(top_srcdir)/util
+AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/subplex -I$(top_srcdir)/praxis -I$(top_srcdir)/luksan -I$(top_srcdir)/util
include_HEADERS = nlopt.h
noinst_LTLIBRARIES = libapi.la
#include "stogo.h"
-#include "l-bfgs-b.h"
-
#include "cdirect.h"
+#include "luksan.h"
+
/*************************************************************************/
/* for "hybrid" algorithms that combine global search with some
&stop, minf);
}
- case NLOPT_LD_LBFGS: {
- int iret, *nbd = (int *) malloc(sizeof(int) * n);
- if (!nbd) return NLOPT_OUT_OF_MEMORY;
- for (i = 0; i < n; ++i) {
- int linf = my_isinf(lb[i]) && lb[i] < 0;
- int uinf = my_isinf(ub[i]) && ub[i] > 0;
- nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
- }
- iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,
- MIN(n, 5), 0.0, ftol_rel,
- xtol_abs ? *xtol_abs : xtol_rel,
- maxeval);
- free(nbd);
- if (iret <= 0) {
- switch (iret) {
- case -1: return NLOPT_INVALID_ARGS;
- case -2: default: return NLOPT_FAILURE;
- }
- }
- else {
- *minf = f(n, x, NULL, f_data);
- switch (iret) {
- case 5: return NLOPT_MAXEVAL_REACHED;
- case 2: return NLOPT_XTOL_REACHED;
- case 1: return NLOPT_FTOL_REACHED;
- default: return NLOPT_SUCCESS;
- }
- }
- break;
- }
+ case NLOPT_LD_LBFGS:
+ return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop);
default:
return NLOPT_INVALID_ARGS;
direct/Makefile
cdirect/Makefile
stogo/Makefile
- lbfgs/Makefile
subplex/Makefile
praxis/Makefile
+ luksan/Makefile
test/Makefile
])
noinst_LTLIBRARIES = liblbfgs.la
liblbfgs_la_SOURCES = ap.cpp ap.h l-bfgs-b.cpp l-bfgs-b.h
-EXTRA_DIST =
+EXTRA_DIST = README
--- /dev/null
+Limited-memory BFGS code by Ciyou Zhu, Richard Byrd, Peihuang Lu, and
+Jorge Nocedal, which can be found (in various translations) at various
+online sites:
+
+ http://www.alglib.net/optimization/lbfgsb.php
+ http://plato.asu.edu/ftp/other_software/toms778_f90.tar.gz
+ http://www-fp.mcs.anl.gov/otc/Tools/LBFGS-B/
+
+Unfortunately, since it was published in ACM TOMS, it seems to be
+under the restrictive (semi-free) ACM copyright conditions. This was
+confirmed in a private email from the author, Prof. Nocedal.
+
+Therefore, we no longer use this code in NLopt.
+
+Fortunately, Prof. Luksan's PLIN.FOR code (see the luksan/ directory)
+is an independent implementation of the same algorithm, and is licensed
+under the LGPL.
--- /dev/null
+---------------------------------------------------------------------------
+
+Copyright:
+
+ Subroutines PBUN, PNEW, PVAR, PSEN, Copyright ACM, 2001. The original
+ versions were published in Transactions on Mathematical Software,
+ Vol.27, 2001, pp.193-213. Here are the author's modifications. They
+ are posted here by permission of ACM for your personal use. Not for
+ redistribution. Subroutines PLIP, PSEN, Copyright Jan Vlcek, 2007.
+ The remaining subroutines, Copyright Ladislav Luksan, 2007. Many of
+ sparse matrix modules were prepared by Miroslav Tuma.
+
+License:
+
+ This library (with exception of PBUN, PNEW, PVAR, PSEN) is a free
+ software; you can redistribute it and/or modify it under the terms
+ of the GNU Lesser General Public License as published by the Free
+ Software Foundation; either version 2.1 of the License, or (at your
+ option) any later version (see http://www.gnu.org/copyleft/gpl.html).
+
+ This library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ Permission is hereby granted to use or copy this program under the
+ terms of the GNU LGPL, provided that the Copyright, this License,
+ and the Availability of the original version is retained on all copies.
+ User documentation of any code that uses this code or any modified
+ version of this code must cite the Copyright, this License, the
+ Availability note, and "Used by permission." Permission to modify
+ the code and to distribute modified code is granted, provided the
+ Copyright, this License, and the Availability note are retained,
+ and a notice that the code was modified is included.
+
+Availability:
+
+ http://www.cs.cas.cz/~luksan/subroutines.html
+
+Acknowledgements:
+
+ This work was supported by the Grant Agency of the Czech Academy of
+ Sciences, under grant IAA1030405.
+
+
--- /dev/null
+AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api
+
+noinst_LTLIBRARIES = libluksan.la
+libluksan_la_SOURCES = plis.c mssubs.c pssubs.c luksan.h
+
+EXTRA_DIST = README COPYRIGHT plis.txt v999-07.pdf
--- /dev/null
+These routines for nonlinear optimization are based on the Fortran 77
+code placed online by Professor Ladislav Luksan at his web site:
+
+ http://www.uivt.cas.cz/~luksan/subroutines.html
+
+and graciously licensed under the GNU Lesser General Public License
+(LGPL). See also the COPYRIGHT file for details.
+
+The C conversions were done via f2c by S. G. Johnson, with the f2c
+output manually cleaned up somewhat, and then converted to use the
+NLopt termination conditions, C dynamic allocation, etc.
+
+In particular, we converted the subroutines PLIS, PLIP, and PNET (as
+well as various auxiliary routines required by this code), all of
+which solve general nonlinear unconstrained or box-constrained
+optimization problems. Prof. Luksan also provides a large number of
+more specialized routines for cases where one has additional knowledge
+of the objective function, but these cases seem outside the scope of
+NLopt.
+
+[We also did not convert the PMIN, PBUN, PNEW, or PVAR subroutines
+since these were published in ACM Trans. on Math. Software (TOMS) and
+are subject to the non-free ACM licensing conditions.]
--- /dev/null
+#ifndef LUKSAN_H
+#define LUKSAN_H
+
+#include <nlopt.h>
+#include <nlopt-util.h>
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif /* __cplusplus */
+
+nlopt_result luksan_plis(int n, nlopt_func f, void *f_data,
+ const double *lb, const double *ub, /* bounds */
+ double *x, /* in: initial guess, out: minimizer */
+ double *minf,
+ nlopt_stopping *stop);
+
+/***************************** internal routines *************************/
+
+/* mssubs.c: */
+void luksan_mxdcmd__(int *n, int *m, double *a,
+ double *x, double *alf, double *y, double *z__);
+void luksan_mxdrcb__(int *n, int *m, double *a,
+ double *b, double *u, double *v, double *x, int *
+ ix, int *job);
+void luksan_mxdrcf__(int *n, int *m, double *a,
+ double *b, double *u, double *v, double *x, int *
+ ix, int *job);
+void luksan_mxdrmm__(int *n, int *m, double *a,
+ double *x, double *y);
+void luksan_mxdrsu__(int *n, int *m, double *a,
+ double *b, double *u);
+void luksan_mxucop__(int *n, double *x, double *y,
+ int *ix, int *job);
+void luksan_mxudir__(int *n, double *a, double *x,
+ double *y, double *z__, int *ix, int *job);
+void luksan_mxuneg__(int *n, double *x, double *y,
+ int *ix, int *job);
+void luksan_mxuzer__(int *n, double *x, int *ix,
+ int *job);
+void luksan_mxvcop__(int *n, double *x, double *y);
+void luksan_mxvdif__(int *n, double *x, double *y,
+ double *z__);
+void luksan_mxvneg__(int *n, double *x, double *y);
+void luksan_mxvscl__(int *n, double *a, double *x,
+ double *y);
+void luksan_mxvset__(int *n, double *a, double *x);
+double luksan_mxudot__(int *n, double *x, double *y, int *ix,
+ int *job);
+double luksan_mxvdot__(int *n, double *x, double *y);
+void luksan_mxvdir__(int *n, double *a, double *x,
+ double *y, double *z__);
+void luksan_mxdcmu__(int *n, int *m, double *a,
+ double *alf, double *x, double *y);
+void luksan_mxvlin__(int *n, double *a, double *x,
+ double *b, double *y, double *z__);
+void luksan_mxdcmv__(int *n, int *m, double *a,
+ double *alf, double *x, double *u, double *bet,
+ double *y, double *v);
+void luksan_mxvsav__(int *n, double *x, double *y);
+void luksan_mxvine__(int *n, int *ix);
+double luksan_mxvmax__(int *n, double *x);
+
+/* pssubs.c: */
+void luksan_pcbs04__(int *nf, double *x, int *ix,
+ double *xl, double *xu, double *eps9, int *kbf);
+void luksan_ps1l01__(double *r__, double *rp,
+ double *f, double *fo, double *fp, double *p,
+ double *po, double *pp, double *minf, double *fmax,
+ double *rmin, double *rmax, double *tols, double *
+ tolp, double *par1, double *par2, int *kd, int *ld,
+ int *nit, int *kit, int *nred, int *mred, int *
+ maxst, int *iest, int *inits, int *iters, int *kters,
+ int *mes, int *isys);
+void luksan_pulsp3__(int *n, int *m, int *mf,
+ double *xm, double *gr, double *xo, double *go,
+ double *r__, double *po, double *sig, int *iterh,
+ int *met3);
+void luksan_pulvp3__(int *n, int *m, double *xm,
+ double *xr, double *gr, double *s, double *so,
+ double *xo, double *go, double *r__, double *po,
+ double *sig, int *iterh, int *met2, int *met3,
+ int *met5);
+void luksan_pyadc0__(int *nf, int *n, double *x,
+ int *ix, double *xl, double *xu, int *inew);
+void luksan_pyfut1__(int *n, double *f, double *
+ fo, double *umax, double *gmax, double *dmax__,
+ double *tolx, double *tolf, double *tolb, double *
+ tolg, int *kd, int *nit, int *kit, int *mit, int *
+ nfv, int *mfv, int *nfg, int *mfg, int *ntesx,
+ int *mtesx, int *ntesf, int *mtesf, int *ites,
+ int *ires1, int *ires2, int *irest, int *iters,
+ int *iterm);
+void luksan_pyrmc0__(int *nf, int *n, int *ix,
+ double *g, double *eps8, double *umax, double *gmax,
+ double *rmax, int *iold, int *irest);
+void luksan_pytrcd__(int *nf, double *x, int *ix,
+ double *xo, double *g, double *go, double *r__,
+ double *f, double *fo, double *p, double *po,
+ double *dmax__, int *kbf, int *kd, int *ld, int *
+ iters);
+void luksan_pytrcg__(int *nf, int *n, int *ix,
+ double *g, double *umax, double *gmax, int *kbf,
+ int *iold);
+void luksan_pytrcs__(int *nf, double *x, int *ix,
+ double *xo, double *xl, double *xu, double *g,
+ double *go, double *s, double *ro, double *fp,
+ double *fo, double *f, double *po, double *p,
+ double *rmax, double *eta9, int *kbf);
+void luksan_pnint1__(double *rl, double *ru, double *fl,
+ double *fu, double *pl, double *pu, double *r__,
+ int *mode, int *mtyp, int *merr);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif /* __cplusplus */
+
+#endif /* LUKSAN_H */
--- /dev/null
+#include <math.h>
+#include "luksan.h"
+
+#define max(a,b) ((a) > (b) ? (a) : (b))
+#define iabs(a) ((a) < 0 ? -(a) : (a))
+
+/* subroutines extracted from mssubs.for */
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* FUNCTION MXVMAX ALL SYSTEMS 91/12/01 */
+/* PURPOSE : */
+/* L-INFINITY NORM OF A VECTOR. */
+
+/* PARAMETERS : */
+/* II N VECTOR DIMENSION. */
+/* RI X(N) INPUT VECTOR. */
+/* RR MXVMAX L-INFINITY NORM OF THE VECTOR X. */
+
+double luksan_mxvmax__(int *n, double *x)
+{
+ /* System generated locals */
+ int i__1;
+ double d__1, d__2, d__3;
+
+ /* Local variables */
+ int i__;
+ double mxvmax;
+
+ /* Parameter adjustments */
+ --x;
+
+ /* Function Body */
+ mxvmax = 0.;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+/* Computing MAX */
+ d__2 = mxvmax, d__3 = (d__1 = x[i__], fabs(d__1));
+ mxvmax = max(d__2,d__3);
+/* L1: */
+ }
+ return mxvmax;
+} /* luksan_mxvmax__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXVINE ALL SYSTEMS 94/12/01 */
+/* PURPOSE : */
+/* ELEMENTS OF THE INTEGER VECTOR ARE REPLACED BY THEIR ABSOLUTE VALUES. */
+
+/* PARAMETERS : */
+/* II N DIMENSION OF THE INTEGER VECTOR. */
+/* IU IX(N) INTEGER VECTOR WHICH IS UPDATED SO THAT IX(I):=ABS(IX(I)) */
+/* FOR ALL I. */
+
+void luksan_mxvine__(int *n, int *ix)
+{
+ /* System generated locals */
+ int i__1, i__2;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --ix;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ ix[i__] = (i__2 = ix[i__], iabs(i__2));
+/* L1: */
+ }
+ return;
+} /* luksan_mxvine__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXDCMV ALL SYSTEMS 91/12/01 */
+/* PURPOSE : */
+/* RANK-TWO UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A. */
+/* THIS MATRIX IS UPDATED BY THE RULE A:=A+ALF*X*TRANS(U)+BET*Y*TRANS(V). */
+
+/* PARAMETERS : */
+/* II N NUMBER OF ROWS OF THE MATRIX A. */
+/* II M NUMBER OF COLUMNS OF THE MATRIX A. */
+/* RU A(N*M) RECTANGULAR MATRIX STORED COLUMNWISE IN THE */
+/* ONE-DIMENSIONAL ARRAY. */
+/* RI ALF SCALAR PARAMETER. */
+/* RI X(N) INPUT VECTOR. */
+/* RI U(M) INPUT VECTOR. */
+/* RI BET SCALAR PARAMETER. */
+/* RI Y(N) INPUT VECTOR. */
+/* RI V(M) INPUT VECTOR. */
+
+void luksan_mxdcmv__(int *n, int *m, double *a,
+ double *alf, double *x, double *u, double *bet,
+ double *y, double *v)
+{
+ /* System generated locals */
+ int i__1, i__2;
+
+ /* Local variables */
+ int i__, j, k;
+ double tempa, tempb;
+
+ /* Parameter adjustments */
+ --v;
+ --y;
+ --u;
+ --x;
+ --a;
+
+ /* Function Body */
+ k = 0;
+ i__1 = *m;
+ for (j = 1; j <= i__1; ++j) {
+ tempa = *alf * u[j];
+ tempb = *bet * v[j];
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ a[k + i__] = a[k + i__] + tempa * x[i__] + tempb * y[i__];
+/* L1: */
+ }
+ k += *n;
+/* L2: */
+ }
+ return;
+} /* luksan_mxdcmv__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXVSAV ALL SYSTEMS 91/12/01 */
+/* PORTABILITY : ALL SYSTEMS */
+/* 91/12/01 LU : ORIGINAL VERSION */
+
+/* PURPOSE : */
+/* DIFFERENCE OF TWO VECTORS RETURNED IN THE SUBSTRACTED ONE. */
+
+/* PARAMETERS : */
+/* II N VECTOR DIMENSION. */
+/* RI X(N) INPUT VECTOR. */
+/* RU Y(N) UPDATE VECTOR WHERE Y:= X - Y. */
+
+void luksan_mxvsav__(int *n, double *x, double *y)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+ double temp;
+
+ /* Parameter adjustments */
+ --y;
+ --x;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = y[i__];
+ y[i__] = x[i__] - y[i__];
+ x[i__] = temp;
+/* L10: */
+ }
+ return;
+} /* luksan_mxvsav__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXVLIN ALL SYSTEMS 92/12/01 */
+/* PURPOSE : */
+/* LINEAR COMBINATION OF TWO VECTORS. */
+
+/* PARAMETERS : */
+/* II N VECTOR DIMENSION. */
+/* RI A SCALING FACTOR. */
+/* RI X(N) INPUT VECTOR. */
+/* RI B SCALING FACTOR. */
+/* RI Y(N) INPUT VECTOR. */
+/* RO Z(N) OUTPUT VECTOR WHERE Z:= A*X + B*Y. */
+
+void luksan_mxvlin__(int *n, double *a, double *x,
+ double *b, double *y, double *z__)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --z__;
+ --y;
+ --x;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ z__[i__] = *a * x[i__] + *b * y[i__];
+/* L1: */
+ }
+ return;
+} /* luksan_mxvlin__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXDCMU ALL SYSTEMS 91/12/01 */
+/* PURPOSE : */
+/* UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A. THIS MATRIX */
+/* IS UPDATED BY THE RULE A:=A+ALF*X*TRANS(Y). */
+
+/* PARAMETERS : */
+/* II N NUMBER OF ROWS OF THE MATRIX A. */
+/* II M NUMBER OF COLUMNS OF THE MATRIX A. */
+/* RU A(N*M) RECTANGULAR MATRIX STORED COLUMNWISE IN THE */
+/* ONE-DIMENSIONAL ARRAY. */
+/* RI ALF SCALAR PARAMETER. */
+/* RI X(N) INPUT VECTOR. */
+/* RI Y(M) INPUT VECTOR. */
+
+void luksan_mxdcmu__(int *n, int *m, double *a,
+ double *alf, double *x, double *y)
+{
+ /* System generated locals */
+ int i__1, i__2;
+
+ /* Local variables */
+ int i__, j, k;
+ double temp;
+
+ /* Parameter adjustments */
+ --y;
+ --x;
+ --a;
+
+ /* Function Body */
+ k = 0;
+ i__1 = *m;
+ for (j = 1; j <= i__1; ++j) {
+ temp = *alf * y[j];
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ a[k + i__] += temp * x[i__];
+/* L1: */
+ }
+ k += *n;
+/* L2: */
+ }
+ return;
+} /* luksan_mxdcmu__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXVDIR ALL SYSTEMS 91/12/01 */
+/* PURPOSE : */
+/* VECTOR AUGMENTED BY THE SCALED VECTOR. */
+
+/* PARAMETERS : */
+/* II N VECTOR DIMENSION. */
+/* RI A SCALING FACTOR. */
+/* RI X(N) INPUT VECTOR. */
+/* RI Y(N) INPUT VECTOR. */
+/* RO Z(N) OUTPUT VECTOR WHERE Z:= Y + A*X. */
+
+void luksan_mxvdir__(int *n, double *a, double *x,
+ double *y, double *z__)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --z__;
+ --y;
+ --x;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ z__[i__] = y[i__] + *a * x[i__];
+/* L10: */
+ }
+} /* luksan_mxvdir__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXDCMD ALL SYSTEMS 91/12/01
+* PURPOSE :
+* MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A
+* BY A VECTOR X AND ADDITION OF THE SCALED VECTOR ALF*Y.
+*
+* PARAMETERS :
+* II N NUMBER OF ROWS OF THE MATRIX A.
+* II M NUMBER OF COLUMNS OF THE MATRIX A.
+* RI A(N*M) RECTANGULAR MATRIX STORED COLUMNWISE IN THE
+* ONE-DIMENSIONAL ARRAY.
+* RI X(M) INPUT VECTOR.
+* RI ALF SCALING FACTOR.
+* RI Y(N) INPUT VECTOR.
+* RO Z(N) OUTPUT VECTOR EQUAL TO A*X+ALF*Y.
+*
+* SUBPROGRAMS USED :
+* S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR.
+* S MXVSCL SCALING OF A VECTOR.
+*/
+void luksan_mxdcmd__(int *n, int *m, double *a,
+ double *x, double *alf, double *y, double *z__)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int j, k;
+
+ /* Parameter adjustments */
+ --z__;
+ --y;
+ --x;
+ --a;
+
+ /* Function Body */
+ luksan_mxvscl__(n, alf, &y[1], &z__[1]);
+ k = 0;
+ i__1 = *m;
+ for (j = 1; j <= i__1; ++j) {
+ luksan_mxvdir__(n, &x[j], &a[k + 1], &z__[1], &z__[1]);
+ k += *n;
+/* L1: */
+ }
+ return;
+} /* luksan_mxdcmd__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXDRCB ALL SYSTEMS 91/12/01
+* PURPOSE :
+* BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION OF
+* THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
+*
+* PARAMETERS :
+* II N NUMBER OF ROWS OF THE MATRICES A AND B.
+* II M NUMBER OF COLUMNS OF THE MATRICES A AND B.
+* RI A(N*M) RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
+* RI B(N*M) RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
+* RI U(M) VECTOR OF SCALAR COEFFICIENTS.
+* RO V(M) VECTOR OF SCALAR COEFFICIENTS.
+* RU X(N) PREMULTIPLIED VECTOR.
+* II IX(N) VECTOR CONTAINING TYPES OF BOUNDS.
+* II JOB OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).EQ.-5.
+*
+* SUBPROGRAM USED :
+* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR.
+* RF MXUDOT DOT PRODUCT OF VECTORS.
+*
+* METHOD :
+* H.MATTHIES, G.STRANG: THE SOLUTION OF NONLINEAR FINITE ELEMENT
+* EQUATIONS. INT.J.NUMER. METHODS ENGN. 14 (1979) 1613-1626.
+*/
+void luksan_mxdrcb__(int *n, int *m, double *a,
+ double *b, double *u, double *v, double *x, int *
+ ix, int *job)
+{
+ /* System generated locals */
+ int i__1;
+ double d__1;
+
+ /* Local variables */
+ int i__, k;
+
+ /* Parameter adjustments */
+ --ix;
+ --x;
+ --v;
+ --u;
+ --b;
+ --a;
+
+ /* Function Body */
+ k = 1;
+ i__1 = *m;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ v[i__] = u[i__] * luksan_mxudot__(n, &x[1], &a[k], &ix[1], job);
+ d__1 = -v[i__];
+ luksan_mxudir__(n, &d__1, &b[k], &x[1], &x[1], &ix[1], job);
+ k += *n;
+/* L1: */
+ }
+ return;
+} /* luksan_mxdrcb__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXDRCF ALL SYSTEMS 91/12/01
+* PURPOSE :
+* FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION OF
+* THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
+*
+* PARAMETERS :
+* II N NUMBER OF ROWS OF THE MATRICES A AND B.
+* II M NUMBER OF COLUMNS OF THE MATRICES A AND B.
+* RI A(N*M) RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
+* RI B(N*M) RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
+* RI U(M) VECTOR OF SCALAR COEFFICIENTS.
+* RI V(M) VECTOR OF SCALAR COEFFICIENTS.
+* RU X(N) PREMULTIPLIED VECTOR.
+* II IX(N) VECTOR CONTAINING TYPES OF BOUNDS.
+* II JOB OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).EQ.-5.
+*
+* SUBPROGRAM USED :
+* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR.
+* RF MXUDOT DOT PRODUCT OF VECTORS.
+*
+* METHOD :
+* H.MATTHIES, G.STRANG: THE SOLUTION OF NONLINEAR FINITE ELEMENT
+* EQUATIONS. INT.J.NUMER. METHODS ENGN. 14 (1979) 1613-1626.
+*/
+void luksan_mxdrcf__(int *n, int *m, double *a,
+ double *b, double *u, double *v, double *x, int *
+ ix, int *job)
+{
+ /* System generated locals */
+ double d__1;
+
+ /* Local variables */
+ int i__, k;
+ double temp;
+
+ /* Parameter adjustments */
+ --ix;
+ --x;
+ --v;
+ --u;
+ --b;
+ --a;
+
+ /* Function Body */
+ k = (*m - 1) * *n + 1;
+ for (i__ = *m; i__ >= 1; --i__) {
+ temp = u[i__] * luksan_mxudot__(n, &x[1], &b[k], &ix[1], job);
+ d__1 = v[i__] - temp;
+ luksan_mxudir__(n, &d__1, &a[k], &x[1], &x[1], &ix[1], job);
+ k -= *n;
+/* L1: */
+ }
+ return;
+} /* luksan_mxdrcf__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXDRMM ALL SYSTEMS 91/12/01
+* PURPOSE :
+* MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY
+* A VECTOR X.
+*
+* PARAMETERS :
+* II N NUMBER OF COLUMNS OF THE MATRIX A.
+* II M NUMBER OF ROWS OF THE MATRIX A.
+* RI A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE
+* ONE-DIMENSIONAL ARRAY.
+* RI X(N) INPUT VECTOR.
+* RO Y(M) OUTPUT VECTOR EQUAL TO A*X.
+*/
+void luksan_mxdrmm__(int *n, int *m, double *a,
+ double *x, double *y)
+{
+ /* System generated locals */
+ int i__1, i__2;
+
+ /* Local variables */
+ int i__, j, k;
+ double temp;
+
+ /* Parameter adjustments */
+ --y;
+ --x;
+ --a;
+
+ /* Function Body */
+ k = 0;
+ i__1 = *m;
+ for (j = 1; j <= i__1; ++j) {
+ temp = 0.;
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp += a[k + i__] * x[i__];
+/* L1: */
+ }
+ y[j] = temp;
+ k += *n;
+/* L2: */
+ }
+ return;
+} /* luksan_mxdrmm__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXDRSU ALL SYSTEMS 91/12/01
+* PURPOSE :
+* SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B. SHIFT OF
+* ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN THE LIMITED
+* MEMORY BFGS METHOD.
+*
+* PARAMETERS :
+* II N NUMBER OF ROWS OF THE MATRIX A AND B.
+* II M NUMBER OF COLUMNS OF THE MATRIX A AND B.
+* RU A(N*M) RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
+* RU B(N*M) RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
+* RU U(M) VECTOR.
+*/
+void luksan_mxdrsu__(int *n, int *m, double *a,
+ double *b, double *u)
+{
+ int i__, k, l;
+
+ /* Parameter adjustments */
+ --u;
+ --b;
+ --a;
+
+ /* Function Body */
+ k = (*m - 1) * *n + 1;
+ for (i__ = *m - 1; i__ >= 1; --i__) {
+ l = k - *n;
+ luksan_mxvcop__(n, &a[l], &a[k]);
+ luksan_mxvcop__(n, &b[l], &b[k]);
+ u[i__ + 1] = u[i__];
+ k = l;
+/* L1: */
+ }
+ return;
+} /* luksan_mxdrsu__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXUCOP ALL SYSTEMS 99/12/01
+* PURPOSE :
+* COPY OF THE VECTOR WITH INITIATION OF THE ACTIVE PART.
+*
+* PARAMETERS :
+* II N VECTOR DIMENSION.
+* RI X(N) INPUT VECTOR.
+* RO Y(N) OUTPUT VECTOR WHERE Y:= X.
+* II IX(N) VECTOR CONTAINING TYPES OF BOUNDS.
+* II JOB OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).EQ.-5.
+*/
+void luksan_mxucop__(int *n, double *x, double *y,
+ int *ix, int *job)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --ix;
+ --y;
+ --x;
+
+ /* Function Body */
+ if (*job == 0) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ y[i__] = x[i__];
+/* L1: */
+ }
+ } else if (*job > 0) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] >= 0) {
+ y[i__] = x[i__];
+ } else {
+ y[i__] = 0.;
+ }
+/* L2: */
+ }
+ } else {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] != -5) {
+ y[i__] = x[i__];
+ } else {
+ y[i__] = 0.;
+ }
+/* L3: */
+ }
+ }
+ return;
+} /* luksan_mxucop__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXUDIR ALL SYSTEMS 99/12/01
+* PURPOSE :
+* VECTOR AUGMENTED BY THE SCALED VECTOR IN A BOUND CONSTRAINED CASE.
+*
+* PARAMETERS :
+* II N VECTOR DIMENSION.
+* RI A SCALING FACTOR.
+* RI X(N) INPUT VECTOR.
+* RI Y(N) INPUT VECTOR.
+* RO Z(N) OUTPUT VECTOR WHERE Z:= Y + A*X.
+* II IX(N) VECTOR CONTAINING TYPES OF BOUNDS.
+* II JOB OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).EQ.-5.
+*/
+void luksan_mxudir__(int *n, double *a, double *x,
+ double *y, double *z__, int *ix, int *job)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --ix;
+ --z__;
+ --y;
+ --x;
+
+ /* Function Body */
+ if (*job == 0) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ z__[i__] = y[i__] + *a * x[i__];
+/* L1: */
+ }
+ } else if (*job > 0) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] >= 0) {
+ z__[i__] = y[i__] + *a * x[i__];
+ }
+/* L2: */
+ }
+ } else {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] != -5) {
+ z__[i__] = y[i__] + *a * x[i__];
+ }
+/* L3: */
+ }
+ }
+ return;
+} /* luksan_mxudir__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* FUNCTION MXVDOT ALL SYSTEMS 91/12/01 */
+/* PURPOSE : */
+/* DOT PRODUCT OF TWO VECTORS. */
+
+/* PARAMETERS : */
+/* II N VECTOR DIMENSION. */
+/* RI X(N) INPUT VECTOR. */
+/* RI Y(N) INPUT VECTOR. */
+/* RR MXVDOT VALUE OF DOT PRODUCT MXVDOT=TRANS(X)*Y. */
+
+double luksan_mxvdot__(int *n, double *x, double *y)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+ double temp;
+
+ /* Parameter adjustments */
+ --y;
+ --x;
+
+ /* Function Body */
+ temp = 0.;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp += x[i__] * y[i__];
+/* L10: */
+ }
+ return temp;
+} /* luksan_mxvdot__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* FUNCTION MXUDOT ALL SYSTEMS 99/12/01
+* PURPOSE :
+* DOT PRODUCT OF VECTORS IN A BOUND CONSTRAINED CASE.
+*
+* PARAMETERS :
+* II N VECTOR DIMENSION.
+* RI X(N) INPUT VECTOR.
+* RI Y(N) INPUT VECTOR.
+* II IX(N) VECTOR CONTAINING TYPES OF BOUNDS.
+* II JOB OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).EQ.-5.
+* RR MXUDOT VALUE OF DOT PRODUCT MXUDOT=TRANS(X)*Y.
+*/
+double luksan_mxudot__(int *n, double *x, double *y, int *ix,
+ int *job)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+ double temp;
+
+ /* Parameter adjustments */
+ --ix;
+ --y;
+ --x;
+
+ /* Function Body */
+ temp = 0.;
+ if (*job == 0) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp += x[i__] * y[i__];
+/* L1: */
+ }
+ } else if (*job > 0) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] >= 0) {
+ temp += x[i__] * y[i__];
+ }
+/* L2: */
+ }
+ } else {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] != -5) {
+ temp += x[i__] * y[i__];
+ }
+/* L3: */
+ }
+ }
+ return temp;
+} /* luksan_mxudot__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXUNEG ALL SYSTEMS 00/12/01
+* PURPOSE :
+* COPY OF THE VECTOR WITH INITIATION OF THE ACTIVE PART.
+*
+* PARAMETERS :
+* II N VECTOR DIMENSION.
+* RI X(N) INPUT VECTOR.
+* RO Y(N) OUTPUT VECTOR WHERE Y:= X.
+* II IX(N) VECTOR CONTAINING TYPES OF BOUNDS.
+* II JOB OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).EQ.-5.
+*/
+void luksan_mxuneg__(int *n, double *x, double *y,
+ int *ix, int *job)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --ix;
+ --y;
+ --x;
+
+ /* Function Body */
+ if (*job == 0) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ y[i__] = -x[i__];
+/* L1: */
+ }
+ } else if (*job > 0) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] >= 0) {
+ y[i__] = -x[i__];
+ } else {
+ y[i__] = 0.;
+ }
+/* L2: */
+ }
+ } else {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] != -5) {
+ y[i__] = -x[i__];
+ } else {
+ y[i__] = 0.;
+ }
+/* L3: */
+ }
+ }
+ return;
+} /* luksan_mxuneg__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXUZER ALL SYSTEMS 99/12/01
+* PURPOSE :
+* VECTOR ELEMENTS CORRESPONDING TO ACTIVE BOUNDS ARE SET TO ZERO.
+*
+* PARAMETERS :
+* II N VECTOR DIMENSION.
+* RO X(N) OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I.
+* II IX(N) VECTOR CONTAINING TYPES OF BOUNDS.
+* II JOB OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
+* IX(I).EQ.-5.
+*/
+void luksan_mxuzer__(int *n, double *x, int *ix,
+ int *job)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --ix;
+ --x;
+
+ /* Function Body */
+ if (*job == 0) {
+ return;
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] < 0) {
+ x[i__] = 0.;
+ }
+/* L1: */
+ }
+ return;
+} /* luksan_mxuzer__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXVCOP ALL SYSTEMS 88/12/01
+* PURPOSE :
+* COPYING OF A VECTOR.
+*
+* PARAMETERS :
+* II N VECTOR DIMENSION.
+* RI X(N) INPUT VECTOR.
+* RO Y(N) OUTPUT VECTOR WHERE Y:= X.
+*/
+void luksan_mxvcop__(int *n, double *x, double *y)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --y;
+ --x;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ y[i__] = x[i__];
+/* L10: */
+ }
+ return;
+} /* luksan_mxvcop__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXVDIF ALL SYSTEMS 88/12/01
+* PURPOSE :
+* VECTOR DIFFERENCE.
+*
+* PARAMETERS :
+* RI X(N) INPUT VECTOR.
+* RI Y(N) INPUT VECTOR.
+* RO Z(N) OUTPUT VECTOR WHERE Z:= X - Y.
+*/
+void luksan_mxvdif__(int *n, double *x, double *y,
+ double *z__)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --z__;
+ --y;
+ --x;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ z__[i__] = x[i__] - y[i__];
+/* L10: */
+ }
+ return;
+} /* luksan_mxvdif__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXVNEG ALL SYSTEMS 88/12/01
+* PURPOSE :
+* CHANGE THE SIGNS OF VECTOR ELEMENTS.
+*
+* PARAMETERS :
+* II N VECTOR DIMENSION.
+* RI X(N) INPUT VECTOR.
+* RO Y(N) OUTPUT VECTOR WHERE Y:= - X.
+*/
+void luksan_mxvneg__(int *n, double *x, double *y)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --y;
+ --x;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ y[i__] = -x[i__];
+/* L10: */
+ }
+ return;
+} /* luksan_mxvneg__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXVSCL ALL SYSTEMS 88/12/01
+* PURPOSE :
+* SCALING OF A VECTOR.
+*
+* PARAMETERS :
+* II N VECTOR DIMENSION.
+* RI X(N) INPUT VECTOR.
+* RI A SCALING FACTOR.
+* RO Y(N) OUTPUT VECTOR WHERE Y:= A*X.
+*/
+void luksan_mxvscl__(int *n, double *a, double *x,
+ double *y)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --y;
+ --x;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ y[i__] = *a * x[i__];
+/* L1: */
+ }
+ return;
+} /* luksan_mxvscl__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE MXVSET ALL SYSTEMS 88/12/01
+* PURPOSE :
+* A SCALAR IS SET TO ALL THE ELEMENTS OF A VECTOR.
+*
+* PARAMETERS :
+* II N VECTOR DIMENSION.
+* RI A INITIAL VALUE.
+* RO X(N) OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I.
+*/
+void luksan_mxvset__(int *n, double *a, double *x)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --x;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ x[i__] = *a;
+/* L10: */
+ }
+ return;
+} /* luksan_mxvset__ */
+
--- /dev/null
+#include <limits.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include "luksan.h"
+
+#define max(a,b) ((a) > (b) ? (a) : (b))
+#define min(a,b) ((a) < (b) ? (a) : (b))
+
+/* Common Block Declarations */
+
+typedef struct {
+ int nres, ndec, nin, nit;
+ /* int nfv; -- now stored in stop->nevals */
+ int nfg, nfh;
+} stat_common;
+
+/* *********************************************************************** */
+/* SUBROUTINE PLIS ALL SYSTEMS 01/09/22 */
+/* PURPOSE : */
+/* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT */
+/* USE THE LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
+/* RECURRENCES. */
+
+/* PARAMETERS : */
+/* II NF NUMBER OF VARIABLES. */
+/* II NB CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. */
+/* NB>0-SIMPLE BOUNDS ACCEPTED. */
+/* RI X(NF) VECTOR OF VARIABLES. */
+/* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE */
+/* X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I). */
+/* IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND */
+/* XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. */
+/* RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. */
+/* RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. */
+/* RO GF(NF) GRADIENT OF THE OBJECTIVE FUNCTION. */
+/* RO S(NF) DIRECTION VECTOR. */
+/* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. */
+/* RI GO(NF) GRADIENTS DIFFERENCE. */
+/* RA UO(NF) AUXILIARY VECTOR. */
+/* RA VO(NF) AUXILIARY VECTOR. */
+/* RI XMAX MAXIMUM STEPSIZE. */
+/* RI TOLX TOLERANCE FOR CHANGE OF VARIABLES. */
+/* RI TOLF TOLERANCE FOR CHANGE OF FUNCTION VALUES. */
+/* RI TOLB TOLERANCE FOR THE FUNCTION VALUE. */
+/* RI TOLG TOLERANCE FOR THE GRADIENT NORM. */
+/* RI MINF_EST ESTIMATION OF THE MINIMUM FUNCTION VALUE. */
+/* RO GMAX MAXIMUM PARTIAL DERIVATIVE. */
+/* RO F VALUE OF THE OBJECTIVE FUNCTION. */
+/* II MIT MAXIMUM NUMBER OF ITERATIONS. */
+/* II MFV MAXIMUM NUMBER OF FUNCTION EVALUATIONS. */
+/* II IEST ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. */
+/* IEST=1-MINIMUM IS ESTIMATED BY THE VALUE MINF_EST. */
+/* II MF NUMBER OF LIMITED MEMORY STEPS. */
+/* IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. */
+/* ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN */
+/* MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
+/* ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN */
+/* MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
+/* ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. */
+/* ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. */
+/* ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, */
+/* BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. */
+/* ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. */
+/* ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. */
+
+/* VARIABLES IN COMMON /STAT/ (STATISTICS) : */
+/* IO NRES NUMBER OF RESTARTS. */
+/* IO NDEC NUMBER OF MATRIX DECOMPOSITION. */
+/* IO NIN NUMBER OF INNER ITERATIONS. */
+/* IO NIT NUMBER OF ITERATIONS. */
+/* IO NFV NUMBER OF FUNCTION EVALUATIONS. */
+/* IO NFG NUMBER OF GRADIENT EVALUATIONS. */
+/* IO NFH NUMBER OF HESSIAN EVALUATIONS. */
+
+/* SUBPROGRAMS USED : */
+/* S PCBS04 ELIMINATION OF BOX CONSTRAINT VIOLATIONS. */
+/* S PS1L01 STEPSIZE SELECTION USING LINE SEARCH. */
+/* S PYADC0 ADDITION OF A BOX CONSTRAINT. */
+/* S PYFUT1 TEST ON TERMINATION. */
+/* S PYRMC0 DELETION OF A BOX CONSTRAINT. */
+/* S PYTRCD COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC */
+/* UPDATE. */
+/* S PYTRCG COMPUTATION OF THE PROJECTED GRADIENT. */
+/* S PYTRCS COMPUTATION OF THE PROJECTED DIRECTION VECTOR. */
+/* S MXDRCB BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
+/* OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
+/* S MXDRCF FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
+/* OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
+/* S MXDRSU SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B. */
+/* SHIFT OF ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN */
+/* THE LIMITED MEMORY BFGS METHOD. */
+/* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. */
+/* RF MXUDOT DOT PRODUCT OF TWO VECTORS. */
+/* S MXUNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. */
+/* S MXVCOP COPYING OF A VECTOR. */
+/* S MXVSCL SCALING OF A VECTOR. */
+
+/* EXTERNAL SUBROUTINES : */
+/* SE OBJ COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. */
+/* CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER */
+/* OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS */
+/* THE VALUE OF THE OBJECTIVE FUNCTION. */
+/* SE DOBJ COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. */
+/* CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER */
+/* OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF) */
+/* IS THE GRADIENT OF THE OBJECTIVE FUNCTION. */
+/* -- OBJ and DOBJ are replaced by a single function, objgrad, in NLopt */
+
+/* METHOD : */
+/* LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
+/* RECURRENCES. */
+
+static void plis_(int *nf, int *nb, double *x, int *
+ ix, double *xl, double *xu, double *gf, double *s,
+ double *xo, double *go, double *uo, double *vo,
+ double *xmax, double *tolx, double *tolf, double *
+ tolb, double *tolg, nlopt_stopping *stop,
+ double *minf_est, double *gmax,
+ double *f, int *mit, int *mfv, int *iest, int *mf,
+ int *iterm, stat_common *stat_1,
+ nlopt_func objgrad, void *objgrad_data)
+{
+ /* System generated locals */
+ int i__1;
+ double d__1, d__2;
+
+ /* Builtin functions */
+
+ /* Local variables */
+ double a, b;
+ int i__, k, n;
+ double p, r__;
+ int kd, ld;
+ double fo, fp, po, pp, ro, rp;
+ int kbf, mfg;
+ int mes, kit;
+ double alf1, alf2, eta0, eta9, par1, par2;
+ int mes1, mes2, mes3;
+ double eps8, eps9;
+ int mred, iold, nred;
+ double maxf, dmax__;
+ int inew;
+ double told;
+ int ites;
+ double rmin, rmax, umax, tolp, tols;
+ int isys;
+ int ires1, ires2;
+ int iterd, mtesf, ntesf;
+ double gnorm;
+ int iters, irest, inits, kters, maxst;
+ double snorm;
+ int mtesx, ntesx;
+
+/* INITIATION */
+
+ /* Parameter adjustments */
+ --vo;
+ --uo;
+ --go;
+ --xo;
+ --s;
+ --gf;
+ --xu;
+ --xl;
+ --ix;
+ --x;
+
+ /* Function Body */
+ kbf = 0;
+ if (*nb > 0) {
+ kbf = 2;
+ }
+ stat_1->nres = 0;
+ stat_1->ndec = 0;
+ stat_1->nin = 0;
+ stat_1->nit = 0;
+ stat_1->nfg = 0;
+ stat_1->nfh = 0;
+ isys = 0;
+ ites = 1;
+ mtesx = 2;
+ mtesf = 2;
+ inits = 2;
+ *iterm = 0;
+ iterd = 0;
+ iters = 2;
+ kters = 3;
+ irest = 0;
+ ires1 = 999;
+ ires2 = 0;
+ mred = 10;
+ mes = 4;
+ mes1 = 2;
+ mes2 = 2;
+ mes3 = 2;
+ eta0 = 1e-15;
+ eta9 = 1e120;
+ eps8 = 1.;
+ eps9 = 1e-8;
+ alf1 = 1e-10;
+ alf2 = 1e10;
+ rmax = eta9;
+ dmax__ = eta9;
+ maxf = 1e20;
+ if (*iest <= 0) {
+ *minf_est = -HUGE_VAL;
+ }
+ if (*iest > 0) {
+ *iest = 1;
+ }
+ if (*xmax <= 0.) {
+ *xmax = 1e16;
+ }
+ if (*tolx <= 0.) {
+ *tolx = 1e-16;
+ }
+ if (*tolf <= 0.) {
+ *tolf = 1e-14;
+ }
+ if (*tolg <= 0.) {
+ *tolg = 1e-8; /* SGJ: was 1e-6, but this sometimes stops too soon */
+ }
+#if 0
+ /* removed by SGJ: this check prevented us from using minf_max <= 0,
+ which doesn't make sense. Instead, if you don't want to have a
+ lower limit, you should set minf_max = -HUGE_VAL */
+ if (*tolb <= 0.) {
+ *tolb = *minf_est + 1e-16;
+ }
+#endif
+ told = 1e-4;
+ tols = 1e-4;
+ tolp = .8;
+ /* changed by SGJ: default is no limit (INT_MAX) on # iterations/fevals */
+ if (*mit <= 0) {
+ *mit = INT_MAX;
+ }
+ if (*mfv <= 0) {
+ *mfv = INT_MAX;
+ }
+ mfg = *mfv;
+ kd = 1;
+ ld = -1;
+ kit = -(ires1 * *nf + ires2);
+ fo = *minf_est;
+
+/* INITIAL OPERATIONS WITH SIMPLE BOUNDS */
+
+ if (kbf > 0) {
+ i__1 = *nf;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
+ xu[i__] = xl[i__];
+ ix[i__] = 5;
+ } else if (ix[i__] == 5 || ix[i__] == 6) {
+ xl[i__] = x[i__];
+ xu[i__] = x[i__];
+ ix[i__] = 5;
+ }
+/* L2: */
+ }
+ luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
+ luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
+ }
+ if (*iterm != 0) {
+ goto L11190;
+ }
+ *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
+ ++stop->nevals;
+ ++stat_1->nfg;
+L11120:
+ luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
+ luksan_pyfut1__(nf, f, &fo, &umax, gmax, &dmax__, tolx, tolf, tolb, tolg,
+ &kd, &stat_1->nit, &kit, mit, &stop->nevals, mfv, &stat_1->nfg, &mfg,
+ &ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
+ iters, iterm);
+ if (*iterm != 0) {
+ goto L11190;
+ }
+ if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
+ if (kbf > 0 && rmax > 0.) {
+ luksan_pyrmc0__(nf, &n, &ix[1], &gf[1], &eps8, &umax, gmax, &rmax, &
+ iold, &irest);
+ }
+L11130:
+
+/* DIRECTION DETERMINATION */
+
+ gnorm = sqrt(luksan_mxudot__(nf, &gf[1], &gf[1], &ix[1], &kbf));
+ if (irest != 0) {
+ goto L12620;
+ }
+/* Computing MIN */
+ i__1 = stat_1->nit - kit;
+ k = min(i__1,*mf);
+ if (k <= 0) {
+ irest = max(irest,1);
+ goto L12620;
+ }
+
+/* DETERMINATION OF THE PARAMETER B */
+
+ b = luksan_mxudot__(nf, &xo[1], &go[1], &ix[1], &kbf);
+ if (b <= 0.) {
+ irest = max(irest,1);
+ goto L12620;
+ }
+ uo[1] = 1. / b;
+ luksan_mxuneg__(nf, &gf[1], &s[1], &ix[1], &kbf);
+ luksan_mxdrcb__(nf, &k, &xo[1], &go[1], &uo[1], &vo[1], &s[1], &ix[1], &
+ kbf);
+ a = luksan_mxudot__(nf, &go[1], &go[1], &ix[1], &kbf);
+ if (a > 0.) {
+ d__1 = b / a;
+ luksan_mxvscl__(nf, &d__1, &s[1], &s[1]);
+ }
+ luksan_mxdrcf__(nf, &k, &xo[1], &go[1], &uo[1], &vo[1], &s[1], &ix[1], &
+ kbf);
+ snorm = sqrt(luksan_mxudot__(nf, &s[1], &s[1], &ix[1], &kbf));
+/* Computing MIN */
+ i__1 = k + 1;
+ k = min(i__1,*mf);
+ luksan_mxdrsu__(nf, &k, &xo[1], &go[1], &uo[1]);
+L12620:
+ iterd = 0;
+ if (irest != 0) {
+
+/* STEEPEST DESCENT DIRECTION */
+
+ luksan_mxuneg__(nf, &gf[1], &s[1], &ix[1], &kbf);
+ snorm = gnorm;
+ if (kit < stat_1->nit) {
+ ++stat_1->nres;
+ kit = stat_1->nit;
+ } else {
+ *iterm = -10;
+ if (iters < 0) {
+ *iterm = iters - 5;
+ }
+ }
+ }
+
+/* TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH */
+
+ if (kd > 0) {
+ p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
+ }
+ if (iterd < 0) {
+ *iterm = iterd;
+ } else {
+
+/* TEST ON DESCENT DIRECTION */
+
+ if (snorm <= 0.) {
+ irest = max(irest,1);
+ } else if (p + told * gnorm * snorm <= 0.) {
+ irest = 0;
+ } else {
+
+/* UNIFORM DESCENT CRITERION */
+
+ irest = max(irest,1);
+ }
+ if (irest == 0) {
+
+/* PREPARATION OF LINE SEARCH */
+
+ nred = 0;
+ rmin = alf1 * gnorm / snorm;
+/* Computing MIN */
+ d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
+ rmax = min(d__1,d__2);
+ }
+ }
+ if (*iterm != 0) {
+ goto L11190;
+ }
+ if (irest != 0) {
+ goto L11130;
+ }
+ luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
+ &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
+ if (rmax == 0.) {
+ goto L11175;
+ }
+L11170:
+ luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, minf_est, &maxf, &rmin,
+ &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1->nit, &kit, &
+ nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
+ if (isys == 0) {
+ goto L11174;
+ }
+ luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
+ luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
+ *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
+ ++stop->nevals;
+ ++stat_1->nfg;
+ p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
+ goto L11170;
+L11174:
+ if (iters <= 0) {
+ r__ = 0.;
+ *f = fo;
+ p = po;
+ luksan_mxvcop__(nf, &xo[1], &x[1]);
+ luksan_mxvcop__(nf, &go[1], &gf[1]);
+ irest = max(irest,1);
+ ld = kd;
+ goto L11130;
+ }
+ luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
+ p, &po, &dmax__, &kbf, &kd, &ld, &iters);
+L11175:
+ if (kbf > 0) {
+ luksan_mxvine__(nf, &ix[1]);
+ luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
+ }
+ goto L11120;
+L11190:
+ return;
+} /* plis_ */
+
+/* number of double variables that can be stored in scratch memory
+ ... it's >= 2007, and this is in the context of scientific computation,
+ so assume that at least 10M are available, and that sizeof(double)==8 */
+#define MEMAVAIL 1310720
+
+/* NLopt wrapper around plis_, handling dynamic allocation etc. */
+nlopt_result luksan_plis(int n, nlopt_func f, void *f_data,
+ const double *lb, const double *ub, /* bounds */
+ double *x, /* in: initial guess, out: minimizer */
+ double *minf,
+ nlopt_stopping *stop)
+{
+ int i, *ix;
+ double *work, *xl, *xu, *xo, *gf, *s, *go, *uo, *vo;
+ double gmax, minf_est;
+ double xmax = 0; /* no maximum */
+ double tolg = 0; /* default gradient tolerance */
+ int iest = 0; /* we have no estimate of min function value */
+ int mit = 0; /* default no limit on #iterations */
+ int mfv = stop->maxeval;
+ stat_common stat;
+ int iterm;
+ int mf;
+
+ ix = (int*) malloc(sizeof(int) * n);
+ if (!ix) return NLOPT_OUT_OF_MEMORY;
+
+ /* FIXME: what should we set mf to? The example program tlis.for
+ sets it to zero as far as I can tell, but it seems to greatly
+ improve convergence to make it > 0. The computation time
+ per iteration, and of course the memory, seem to go as O(n * mf),
+ and we'll assume that the main limiting factor is the memory.
+ We'll assume that at least MEMAVAIL memory, or 4*n memory, whichever
+ is bigger, is available. */
+ mf = max(MEMAVAIL/n, 4);
+ if (stop->maxeval && stop->maxeval <= mf)
+ mf = max(stop->maxeval - 5, 1); /* mf > maxeval seems not good */
+
+ retry_alloc:
+ work = (double*) malloc(sizeof(double) * (n * 4 + max(n,n*mf)*2 +
+ max(n,mf)*2));
+ if (!work) {
+ if (mf > 0) {
+ mf /= 4;
+ goto retry_alloc;
+ }
+ free(ix);
+ return NLOPT_OUT_OF_MEMORY;
+ }
+
+ xl = work; xu = xl + n; gf = xu + n; s = gf + n;
+ xo = s + n; go = xo + max(n,n*mf);
+ uo = go + max(n,n*mf); vo = uo + max(n,mf);
+
+ for (i = 0; i < n; ++i) {
+ int lbu = lb[i] <= -0.99 * HUGE_VAL; /* lb unbounded */
+ int ubu = ub[i] >= 0.99 * HUGE_VAL; /* ub unbounded */
+ ix[i] = lbu ? (ubu ? 0 : 2) : (ubu ? 1 : (lb[i] == ub[i] ? 5 : 3));
+ xl[i] = lb[i];
+ xu[i] = ub[i];
+ }
+
+ /* ? xo does not seem to be initialized in the
+ original Fortran code, but it is used upon
+ input to plis if mf > 0 ... perhaps ALLOCATE initializes
+ arrays to zero by default? */
+ memset(xo, 0, sizeof(double) * max(n,n*mf));
+
+ plis_(&n, &n, x, ix, xl, xu,
+ gf, s, xo, go, uo, vo,
+ &xmax,
+
+ /* fixme: pass tol_rel and tol_abs and use NLopt check */
+ &stop->xtol_rel,
+ &stop->ftol_rel,
+ &stop->minf_max,
+ &tolg,
+ stop,
+
+ &minf_est, &gmax,
+ minf,
+ &mit, &mfv,
+ &iest,
+ &mf,
+ &iterm, &stat,
+ f, f_data);
+
+ free(work);
+ free(ix);
+
+ switch (iterm) {
+ case 1: return NLOPT_XTOL_REACHED;
+ case 2: return NLOPT_FTOL_REACHED;
+ case 3: return NLOPT_MINF_MAX_REACHED;
+ case 4: return NLOPT_SUCCESS; /* gradient tolerance reached */
+ case 6: return NLOPT_SUCCESS;
+ case 12: case 13: return NLOPT_MAXEVAL_REACHED;
+ default: return NLOPT_FAILURE;
+ }
+}
--- /dev/null
+***********************************************************************
+* *
+* PLIS - A LIMITED MEMORY VARIABLE METRIC ALGORITHM FOR *
+* LARGE-SCALE OPTIMIZATION. *
+* *
+***********************************************************************
+
+
+1. Introduction:
+----------------
+
+ The double-precision FORTRAN 77 basic subroutine PLIS is designed
+to find a close approximation to a local minimum of a nonlinear
+function F(X) with simple bounds on variables. Here X is a vector of NF
+variables and F(X) is a smooth function. We suppose that NF is large
+but the sparsity pattern of the Hessian matrix is not known (or the
+Hessian matrix is dense). Simple bounds are assumed in the form
+
+ X(I) unbounded if IX(I) = 0,
+ XL(I) <= X(I) if IX(I) = 1,
+ X(I) <= XU(I) if IX(I) = 2,
+ XL(I) <= X(I) <= XU(I) if IX(I) = 3,
+ XL(I) = X(I) = XU(I) if IX(I) = 5,
+
+where 1 <= I <= NF. To simplify user's work, two additional easy to use
+subroutines are added. They call the basic general subroutine PLIS:
+
+ PLISU - unconstrained large-scale optimization,
+ PLISS - large-scale optimization with simple bounds.
+
+All subroutines contain a description of formal parameters and
+extensive comments. Furthermore, two test programs TLISU and TLISS are
+included, which contain several test problems (see e.g. [2]). These
+test programs serve as examples for using the subroutines, verify their
+correctness and demonstrate their efficiency.
+ In this short guide, we describe all subroutines which can be
+called from the user's program. A detailed description of the method is
+given in [1]. In the description of formal parameters, we introduce a
+type of the argument that specifies whether the argument must have a
+value defined on entry to the subroutine (I), whether it is a value
+which will be returned (O), or both (U), or whether it is an auxiliary
+value (A). Note that the arguments of the type I can be changed on
+output under some circumstances, especially if improper input values
+were given. Besides formal parameters, we can use a COMMON /STAT/ block
+containing statistical information. This block, used in each subroutine
+has the following form:
+
+ COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
+
+The arguments have the following meaning:
+
+ Argument Type Significance
+ ----------------------------------------------------------------------
+ NRES O Positive INTEGER variable that indicates the number of
+ restarts.
+ NDEC O Positive INTEGER variable that indicates the number of
+ matrix decompositions.
+ NIN O Positive INTEGER variable that indicates the number of
+ inner iterations (for solving linear systems).
+ NIT O Positive INTEGER variable that indicates the number of
+ iterations.
+ NFV O Positive INTEGER variable that indicates the number of
+ function evaluations.
+ NFG O Positive INTEGER variable that specifies the number of
+ gradient evaluations.
+ NFH O Positive INTEGER variable that specifies the number of
+ Hessian evaluations.
+
+
+2. Subroutines PLISU, PLISS:
+----------------------------
+
+The calling sequences are
+
+ CALL PLISU(NF,X,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
+ CALL PLISS(NF,X,IX,XL,XU,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
+
+The arguments have the following meaning.
+
+ Argument Type Significance
+ ----------------------------------------------------------------------
+ NF I Positive INTEGER variable that specifies the number of
+ variables of the objective function.
+ X(NF) U On input, DOUBLE PRECISION vector with the initial
+ estimate to the solution. On output, the approximation
+ to the minimum.
+ IX(NF) I On input (significant only if NB>0) INTEGER vector
+ containing the simple bounds types:
+ IX(I)=0 - the variable X(I) is unbounded,
+ IX(I)=1 - the lower bound X(I) >= XL(I),
+ IX(I)=2 - the upper bound X(I) <= XU(I),
+ IX(I)=3 - the two side bound XL(I) <= X(I) <= XU(I),
+ IX(I)=5 - the variable X(I) is fixed (given by its
+ initial estimate).
+ XL(NF) I DOUBLE PRECISION vector with lower bounds for variables
+ (significant only if NB>0).
+ XU(NF) I DOUBLE PRECISION vector with upper bounds for variables
+ (significant only if NB>0).
+ IPAR(7) I INTEGER parameters:
+ IPAR(1)=MIT, IPAR(2)=MFV, IPAR(3)-NONE,
+ IPAR(4)=IEST, IPAR(5)-NONE, IPAR(6)-NONE,
+ IPAR(7)=MF.
+ Parameters MIT, MFV, IEST, MF are described in Section 3
+ together with other parameters of the subroutine PLIS.
+ RPAR(9) I DOUBLE PRECISION parameters:
+ RPAR(1)=XMAX, RPAR(2)=TOLX, RPAR(3)=TOLF,
+ RPAR(4)=TOLB, RPAR(5)=TOLG, RPAR(6)=FMIN,
+ RPAR(7)-NONE, RPAR(6)-NONE, RPAR(9)-NONE.
+ Parameters XMAX, TOLX, TOLF, TOLB, TOLG, FMIN are
+ described in Section 3 together with other parameters
+ of the subroutine PLIS.
+ F O DOUBLE PRECISION value of the objective function at the
+ solution X.
+ GMAX O DOUBLE PRECISION maximum absolute value of a partial
+ derivative of the Lagrangian function.
+ IPRNT I INTEGER variable that specifies PRINT:
+ IPRNT= 0 - print is suppressed,
+ IPRNT= 1 - basic print of final results,
+ IPRNT=-1 - extended print of final results,
+ IPRNT= 2 - basic print of intermediate and final
+ results,
+ IPRNT=-2 - extended print of intermediate and final
+ results.
+ ITERM O INTEGER variable that indicates the cause of termination:
+ ITERM= 1 - if |X - XO| was less than or equal to TOLX
+ in two subsequent iterations,
+ ITERM= 2 - if |F - FO| was less than or equal to TOLF
+ in two subsequent iterations,
+ ITERM= 3 - if F is less than or equal to TOLB,
+ ITERM= 4 - if GMAX is less than or equal to TOLG,
+ ITERM= 6 - if termination criterion was not satisfied,
+ but the solution is probably acceptable,
+ ITERM=11 - if NIT exceeded MIT,
+ ITERM=12 - if NFV exceeded MFV,
+ ITERM< 0 - if the method failed.
+
+ The subroutines PLISU, PLISS require the user supplied subroutines
+OBJ and DOBJ that define the objective function and its gradient and
+have the form
+
+ SUBROUTINE OBJ(NF,X,F)
+ SUBROUTINE DOBJ(NF,X,G)
+
+The arguments of the user supplied subroutines have the following
+meaning.
+
+ Argument Type Significance
+ ----------------------------------------------------------------------
+ NF I Positive INTEGER variable that specifies the number of
+ variables of the objective function.
+ X(NF) I DOUBLE PRECISION an estimate to the solution.
+ F O DOUBLE PRECISION value of the objective function at the
+ point X.
+ G(NF) O DOUBLE PRECISION gradient of the objective function
+ at the point X.
+
+
+3. Subroutine PLIS:
+-------------------
+
+ This general subroutine is called from all subroutines described
+in Section 2. The calling sequence is
+
+ CALL PLIS(NF,NB,X,IX,XL,XU,GF,S,XO,GO,UO,VO,XMAX,TOLX,TOLF,TOLB,
+ & TOLG,FMIN,GMAX,F,MIT,MFV,IEST,MF,IPRNT,ITERM)
+
+The arguments NF, NB, X, IX, XL, XU, GMAX, F, IPRNT, ITERM, have the
+same meaning as in Section 2. Other arguments have the following meaning:
+
+ Argument Type Significance
+ ----------------------------------------------------------------------
+ GF(NF) A DOUBLE PRECISION gradient of the objective function.
+ S(NF) A DOUBLE PRECISION direction vector.
+ XO(NF*MF) A DOUBLE PRECISION array which contains increments of
+ variables.
+ GO(NF*MF) A DOUBLE PRECISION array which contains increments of
+ gradients.
+ UO(MF) A DOUBLE PRECISION Auxiliary array.
+ VO(MF) A DOUBLE PRECISION Auxiliary array.
+ XMAX I DOUBLE PRECISION maximum stepsize; the choice XMAX=0
+ causes that the default value 1.0D+16 will be taken.
+ TOLX I DOUBLE PRECISION tolerance for the change of the
+ coordinate vector X; the choice TOLX=0 causes that the
+ default value TOLX=1.0D-16 will be taken.
+ TOLF I DOUBLE PRECISION tolerance for the change of function
+ values; the choice TOLF=0 causes that the default
+ value TOLF=1.0D-14 will be taken.
+ TOLB I DOUBLE PRECISION minimum acceptable function value;
+ the choice TOLB=0 causes that the default value
+ TOLB=FMIN+1.0D-16 will be taken.
+ TOLG I DOUBLE PRECISION tolerance for the Lagrangian function
+ gradient; the choice TOLG=0 causes that the default
+ value TOLG=1.0D-6 will be taken.
+ FMIN I DOUBLE PRECISION lower bound for the minimum function
+ value.
+ MIT I INTEGER variable that specifies the maximum number of
+ iterations; the choice MIT=0 causes that the default
+ value 9000 will be taken.
+ MFV I INTEGER variable that specifies the maximum number of
+ function evaluations; the choice MFV=0 causes that
+ the default value 9000 will be taken.
+ IEST I INTEGER estimation of the minimum functiom value for
+ the line search:
+ IEST=0 - estimation is not used,
+ IEST=1 - lower bound FMIN is used as an estimation
+ for the minimum function value.
+ MF I The number of limited-memory variable metric updates
+ in each iteration (they use 2*MF stored vectors).
+
+The choice of parameter XMAX can be sensitive in many cases. First, the
+objective function can be evaluated only in a relatively small region
+(if it contains exponentials) so that the maximum stepsize is necessary.
+Secondly, the problem can be very ill-conditioned far from the solution
+point so that large steps can be unsuitable. Finally, if the problem has
+more local solutions, a suitably chosen maximum stepsize can lead to
+obtaining a better local solution.
+ The subroutine PLIS requires the user supplied subroutines OBJ
+and DOBJ which are described in Section 2.
+
+4. Verification of the subroutines:
+-----------------------------------
+
+ Subroutine PLISU can be verified and tested using the program
+TLISU. This program calls the subroutines TIUD14 (initiation), TFFU14
+(function evaluation) and TFGU14 (gradient evaluation) containing
+22 unconstrained test problems with at most 1000 variables [2]. The
+results obtained by the program TLISU on a PC computer with Microsoft
+Power Station Fortran compiler have the following form.
+
+NIT= 4988 NFV= 5554 NFG= 5554 F= 0.963780013E-14 G= 0.891E-06 ITERM= 4
+NIT= 425 NFV= 454 NFG= 454 F= 14.9944763 G= 0.773E-05 ITERM= 2
+NIT= 74 NFV= 78 NFG= 78 F= 0.655101686E-09 G= 0.539E-06 ITERM= 4
+NIT= 103 NFV= 112 NFG= 112 F= 269.499543 G= 0.899E-06 ITERM= 4
+NIT= 24 NFV= 26 NFG= 26 F= 0.130639280E-11 G= 0.671E-06 ITERM= 4
+NIT= 30 NFV= 31 NFG= 31 F= 0.216102227E-10 G= 0.946E-06 ITERM= 4
+NIT= 38 NFV= 43 NFG= 43 F= 335.137433 G= 0.730E-06 ITERM= 4
+NIT= 29 NFV= 33 NFG= 33 F= 761774.954 G= 0.432E-03 ITERM= 2
+NIT= 13 NFV= 16 NFG= 16 F= 316.436141 G= 0.369E-06 ITERM= 4
+NIT= 1540 NFV= 1582 NFG= 1582 F= -124.630000 G= 0.124E-04 ITERM= 2
+NIT= 114 NFV= 138 NFG= 138 F= 10.7765879 G= 0.380E-06 ITERM= 4
+NIT= 248 NFV= 267 NFG= 267 F= 982.273617 G= 0.123E-04 ITERM= 2
+NIT= 7 NFV= 8 NFG= 8 F= 0.165734137E-12 G= 0.453E-06 ITERM= 4
+NIT= 10 NFV= 12 NFG= 12 F= 0.128729169E-08 G= 0.916E-06 ITERM= 4
+NIT= 2830 NFV= 2929 NFG= 2929 F= 1.92401599 G= 0.936E-06 ITERM= 4
+NIT= 196 NFV= 210 NFG= 210 F= -427.404476 G= 0.991E-05 ITERM= 2
+NIT= 1007 NFV= 1032 NFG= 1032 F=-0.379921091E-01 G= 0.876E-06 ITERM= 4
+NIT= 1449 NFV= 1474 NFG= 1474 F=-0.245741193E-01 G= 0.862E-06 ITERM= 4
+NIT= 1393 NFV= 1431 NFG= 1431 F= 59.5986241 G= 0.259E-05 ITERM= 2
+NIT= 2129 NFV= 2191 NFG= 2191 F= -1.00013520 G= 0.908E-06 ITERM= 4
+NIT= 2120 NFV= 2169 NFG= 2169 F= 2.13866377 G= 0.927E-06 ITERM= 4
+NIT= 1305 NFV= 1346 NFG= 1346 F= 1.00000000 G= 0.982E-06 ITERM= 4
+NITER =20072 NFVAL =21136 NSUCC = 22
+TIME= 0:00:10.78
+
+The rows corresponding to individual test problems contain the number of
+iterations NIT, the number of function evaluations NFV, the number of
+gradient evaluations NFG, the final value of the objective function F,
+the norm of gradient G and the cause of termination ITERM.
+ Subroutine PLISS can be verified and tested using the program
+TLISS. This program calls the subroutines TIUD14 (initiation), TFFU14
+(function evaluation), TFGU14 (gradient evaluation) containing 22 box
+constrained test problems with at most 1000 variables [2]. The results
+obtained by the program TLISS on a PC computer with Microsoft Power
+Station Fortran compiler have the following form.
+
+NIT= 5063 NFV= 5738 NFG= 5738 F= 0.00000000 G= 0.000E+00 ITERM= 3
+NIT= 3167 NFV= 4664 NFG= 4664 F= 3926.45961 G= 0.626E-04 ITERM= 2
+NIT= 113 NFV= 124 NFG= 124 F= 0.459503394E-12 G= 0.600E-06 ITERM= 4
+NIT= 59 NFV= 64 NFG= 64 F= 269.522686 G= 0.838E-06 ITERM= 4
+NIT= 24 NFV= 26 NFG= 26 F= 0.130639280E-11 G= 0.671E-06 ITERM= 4
+NIT= 30 NFV= 31 NFG= 31 F= 0.216102227E-10 G= 0.946E-06 ITERM= 4
+NIT= 33 NFV= 40 NFG= 40 F= 337.722479 G= 0.592E-06 ITERM= 4
+NIT= 50 NFV= 55 NFG= 55 F= 761925.725 G= 0.240E-03 ITERM= 2
+NIT= 505 NFV= 508 NFG= 508 F= 428.056916 G= 0.334E-07 ITERM= 4
+NIT= 1167 NFV= 1227 NFG= 1227 F= -81.0913589 G= 0.100E-04 ITERM= 2
+NIT= 20 NFV= 26 NFG= 26 F= 96517.2947 G= 0.745E-05 ITERM= 2
+NIT= 91 NFV= 109 NFG= 109 F= 4994.21410 G= 0.104E-04 ITERM= 2
+NIT= 7 NFV= 8 NFG= 8 F= 0.165734137E-12 G= 0.453E-06 ITERM= 4
+NIT= 10 NFV= 12 NFG= 12 F= 0.128729169E-08 G= 0.916E-06 ITERM= 4
+NIT= 2830 NFV= 2929 NFG= 2929 F= 1.92401599 G= 0.936E-06 ITERM= 4
+NIT= 178 NFV= 184 NFG= 184 F= -427.391653 G= 0.107E-04 ITERM= 2
+NIT= 1007 NFV= 1032 NFG= 1032 F=-0.379921091E-01 G= 0.876E-06 ITERM= 4
+NIT= 1449 NFV= 1474 NFG= 1474 F=-0.245741193E-01 G= 0.862E-06 ITERM= 4
+NIT= 1561 NFV= 1595 NFG= 1595 F= 1654.94525 G= 0.112E-04 ITERM= 2
+NIT= 2075 NFV= 2121 NFG= 2121 F= -1.00013520 G= 0.916E-06 ITERM= 4
+NIT= 1361 NFV= 1389 NFG= 1389 F= 2.41354873 G= 0.709E-06 ITERM= 4
+NIT= 1562 NFV= 1598 NFG= 1598 F= 1.00000000 G= 0.786E-06 ITERM= 4
+NITER =22362 NFVAL =24954 NSUCC = 22
+TIME= 0:00:12.39
+
+References:
+-----------
+
+[1] Luksan L., Matonoha C., Vlcek J.: LSA: Algorithms for large-scale
+ unconstrained and box constrained optimization Technical Report V-896.
+ Prague, ICS AS CR, 2004.
+
+[2] Luksan L., Vlcek J.: Sparse and partially separable test problems
+ for unconstrained and equality constrained optimization. Research
+ Report V-767, Institute of Computer Science, Academy of Sciences
+ of the Czech Republic, Prague, Czech Republic, 1998.
+\1a
\ No newline at end of file
--- /dev/null
+#include <math.h>
+#include "luksan.h"
+
+#define FALSE_ 0
+#define max(a,b) ((a) > (b) ? (a) : (b))
+#define min(a,b) ((a) < (b) ? (a) : (b))
+#define iabs(a) ((a) < 0 ? -(a) : (a))
+
+/* subroutines extracted from pssubs.for */
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PCBS04 ALL SYSTEMS 98/12/01
+* PURPOSE :
+* INITIATION OF THE VECTOR CONTAINING TYPES OF CONSTRAINTS.
+*
+* PARAMETERS :
+* II NF NUMBER OF VARIABLES.
+* RI X(NF) VECTOR OF VARIABLES.
+* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
+* RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
+* RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
+* RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS.
+* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
+* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
+*/
+void luksan_pcbs04__(int *nf, double *x, int *ix,
+ double *xl, double *xu, double *eps9, int *kbf)
+{
+ /* System generated locals */
+ int i__1, i__2;
+ double d__1, d__2;
+
+ /* Local variables */
+ int i__, ixi;
+ double temp;
+
+ /* Parameter adjustments */
+ --xu;
+ --xl;
+ --ix;
+ --x;
+
+ /* Function Body */
+ if (*kbf > 0) {
+ i__1 = *nf;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = 1.;
+ ixi = (i__2 = ix[i__], iabs(i__2));
+/* Computing MAX */
+ d__2 = (d__1 = xl[i__], fabs(d__1));
+ if ((ixi == 1 || ixi == 3 || ixi == 4) && x[i__] <= xl[i__] + *
+ eps9 * max(d__2,temp)) {
+ x[i__] = xl[i__];
+ }
+/* Computing MAX */
+ d__2 = (d__1 = xu[i__], fabs(d__1));
+ if ((ixi == 2 || ixi == 3 || ixi == 4) && x[i__] >= xu[i__] - *
+ eps9 * max(d__2,temp)) {
+ x[i__] = xu[i__];
+ }
+/* L1: */
+ }
+ }
+ return;
+} /* luksan_pcbs04__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PNINT1 ALL SYSTEMS 91/12/01 */
+/* PURPOSE : */
+/* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITH DIRECTIONAL */
+/* DERIVATIVES. */
+
+/* PARAMETERS : */
+/* RI RL LOWER VALUE OF THE STEPSIZE PARAMETER. */
+/* RI RU UPPER VALUE OF THE STEPSIZE PARAMETER. */
+/* RI FL VALUE OF THE OBJECTIVE FUNCTION FOR R=RL. */
+/* RI FU VALUE OF THE OBJECTIVE FUNCTION FOR R=RU. */
+/* RI PL DIRECTIONAL DERIVATIVE FOR R=RL. */
+/* RI PU DIRECTIONAL DERIVATIVE FOR R=RU. */
+/* RO R VALUE OF THE STEPSIZE PARAMETER OBTAINED. */
+/* II MODE MODE OF LINE SEARCH. */
+/* II MTYP METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-QUADRATIC */
+/* INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE). */
+/* MTYP=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL */
+/* DERIVATIVES). MTYP=4-CUBIC INTERPOLATION. MTYP=5-CONIC */
+/* INTERPOLATION. */
+/* IO MERR ERROR INDICATOR. MERR=0 FOR NORMAL RETURN. */
+
+/* METHOD : */
+/* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS. */
+
+void luksan_pnint1__(double *rl, double *ru, double *fl,
+ double *fu, double *pl, double *pu, double *r__,
+ int *mode, int *mtyp, int *merr)
+{
+ /* System generated locals */
+ double d__1, d__2;
+
+ /* Local variables */
+ double a, b, c__, d__, den, dis;
+ int ntyp;
+
+ *merr = 0;
+ if (*mode <= 0) {
+ return;
+ }
+ if (*pl >= 0.) {
+ *merr = 2;
+ return;
+ } else if (*ru <= *rl) {
+ *merr = 3;
+ return;
+ }
+ for (ntyp = *mtyp; ntyp >= 1; --ntyp) {
+ if (ntyp == 1) {
+
+/* BISECTION */
+
+ if (*mode == 1) {
+ *r__ = *ru * 4.;
+ return;
+ } else {
+ *r__ = (*rl + *ru) * .5;
+ return;
+ }
+ } else if (ntyp == *mtyp) {
+ a = (*fu - *fl) / (*pl * (*ru - *rl));
+ b = *pu / *pl;
+ }
+ if (ntyp == 2) {
+
+/* QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL */
+/* DERIVATIVE */
+
+ den = (1. - a) * 2.;
+ } else if (ntyp == 3) {
+
+/* QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL */
+/* DERIVATIVES */
+
+ den = 1. - b;
+ } else if (ntyp == 4) {
+
+/* CUBIC EXTRAPOLATION OR INTERPOLATION */
+
+ c__ = b - a * 2. + 1.;
+ d__ = b - a * 3. + 2.;
+ dis = d__ * d__ - c__ * 3.;
+ if (dis < 0.) {
+ goto L1;
+ }
+ den = d__ + sqrt(dis);
+ } else if (ntyp == 5) {
+
+/* CONIC EXTRAPOLATION OR INTERPOLATION */
+
+ dis = a * a - b;
+ if (dis < 0.) {
+ goto L1;
+ }
+ den = a + sqrt(dis);
+ if (den <= 0.) {
+ goto L1;
+ }
+/* Computing 3rd power */
+ d__1 = 1. / den;
+ den = 1. - b * (d__1 * (d__1 * d__1));
+ }
+ if (*mode == 1 && den > 0. && den < 1.) {
+
+/* EXTRAPOLATION ACCEPTED */
+
+ *r__ = *rl + (*ru - *rl) / den;
+/* Computing MAX */
+ d__1 = *r__, d__2 = *ru * 1.1;
+ *r__ = max(d__1,d__2);
+/* Computing MIN */
+ d__1 = *r__, d__2 = *ru * 1e3;
+ *r__ = min(d__1,d__2);
+ return;
+ } else if (*mode == 2 && den > 1.) {
+
+/* INTERPOLATION ACCEPTED */
+
+ *r__ = *rl + (*ru - *rl) / den;
+ if (*rl == 0.) {
+/* Computing MAX */
+ d__1 = *r__, d__2 = *rl + (*ru - *rl) * .01;
+ *r__ = max(d__1,d__2);
+ } else {
+/* Computing MAX */
+ d__1 = *r__, d__2 = *rl + (*ru - *rl) * .1;
+ *r__ = max(d__1,d__2);
+ }
+/* Computing MIN */
+ d__1 = *r__, d__2 = *rl + (*ru - *rl) * .9;
+ *r__ = min(d__1,d__2);
+ return;
+ }
+L1:
+ ;
+ }
+ return;
+} /* luksan_pnint1__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PS1L01 ALL SYSTEMS 97/12/01
+* PURPOSE :
+* STANDARD LINE SEARCH WITH DIRECTIONAL DERIVATIVES.
+*
+* PARAMETERS :
+* RO R VALUE OF THE STEPSIZE PARAMETER.
+* RO RP PREVIOUS VALUE OF THE STEPSIZE PARAMETER.
+* RO F VALUE OF THE OBJECTIVE FUNCTION.
+* RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION.
+* RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
+* RO P VALUE OF THE DIRECTIONAL DERIVATIVE.
+* RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
+* RO PP PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
+* RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
+* RI MAXF UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
+* RI RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER
+* RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER
+* RI TOLS TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
+* CHANGE OF THE FUNCTION VALUE).
+* RI TOLP TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
+* CHANGE OF THE DIRECTIONAL DERIVATIVE).
+* RO PAR1 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
+* UPDATES.
+* RO PAR2 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
+* UPDATES.
+* II KD DEGREE OF REQUIRED DERIVATIVES.
+* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES OF OBJECTIVE
+* II NIT ACTUAL NUMBER OF ITERATIONS.
+* II KIT NUMBER OF THE ITERATION AFTER LAST RESTART.
+* IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
+* II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
+* IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
+* STEPSIZE WAS NOT OR WAS REACHED.
+* II IEST LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND
+* IS NOT OR IS GIVEN.
+* II INITS CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE
+* IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL
+* STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED
+* INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL
+* STEPSIZE.
+* IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
+* LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
+* STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
+* ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
+* ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
+* ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
+* DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
+* II KTERS TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH.
+* KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE.
+* KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE.
+* KTERS=6-FIRST STEPSIZE.
+* II MES METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC
+* INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
+* MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
+* DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC
+* INTERPOLATION.
+* IU ISYS CONTROL PARAMETER.
+*
+* SUBPROGRAM USED :
+* S PNINT1 EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL
+* DERIVATIVES.
+*
+* METHOD :
+* SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH STANDARD TERMINATION
+* CRITERIA.
+*/
+void luksan_ps1l01__(double *r__, double *rp,
+ double *f, double *fo, double *fp, double *p,
+ double *po, double *pp, double *minf, double *maxf,
+ double *rmin, double *rmax, double *tols, double *
+ tolp, double *par1, double *par2, int *kd, int *ld,
+ int *nit, int *kit, int *nred, int *mred, int *
+ maxst, int *iest, int *inits, int *iters, int *kters,
+ int *mes, int *isys)
+{
+ /* System generated locals */
+ double d__1, d__2;
+
+ /* Local variables */
+ unsigned l1, l2, l3, m1, l5, m2, l7, m3;
+ static double fl, fu, pl, rl, pu, ru;
+ static int mes1, mes2, mes3, mode;
+ int merr;
+ static int mtyp;
+ int init1;
+ double rtemp;
+
+ if (*isys == 1) {
+ goto L3;
+ }
+ mes1 = 2;
+ mes2 = 2;
+ mes3 = 2;
+ *iters = 0;
+ if (*po >= 0.) {
+ *r__ = 0.;
+ *iters = -2;
+ goto L4;
+ }
+ if (*rmax <= 0.) {
+ *iters = 0;
+ goto L4;
+ }
+
+/* INITIAL STEPSIZE SELECTION */
+
+ if (*inits > 0) {
+ rtemp = *minf - *f;
+ } else if (*iest == 0) {
+ rtemp = *f - *fp;
+ } else {
+/* Computing MAX */
+ d__1 = *f - *fp, d__2 = *minf - *f;
+ rtemp = max(d__1,d__2);
+ }
+ init1 = iabs(*inits);
+ *rp = 0.;
+ *fp = *fo;
+ *pp = *po;
+ if (init1 == 0) {
+ } else if (init1 == 1 || (*inits >= 1 && *iest == 0)) {
+ *r__ = 1.;
+ } else if (init1 == 2) {
+/* Computing MIN */
+ d__1 = 1., d__2 = rtemp * 4. / *po;
+ *r__ = min(d__1,d__2);
+ } else if (init1 == 3) {
+/* Computing MIN */
+ d__1 = 1., d__2 = rtemp * 2. / *po;
+ *r__ = min(d__1,d__2);
+ } else if (init1 == 4) {
+ *r__ = rtemp * 2. / *po;
+ }
+ *r__ = max(*r__,*rmin);
+ *r__ = min(*r__,*rmax);
+ mode = 0;
+ ru = 0.;
+ fu = *fo;
+ pu = *po;
+
+/* NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) */
+
+L2:
+ luksan_pnint1__(&rl, &ru, &fl, &fu, &pl, &pu, r__, &mode, &mtyp, &merr);
+ if (merr > 0) {
+ *iters = -merr;
+ goto L4;
+ } else if (mode == 1) {
+ --(*nred);
+ *r__ = min(*r__,*rmax);
+ } else if (mode == 2) {
+ ++(*nred);
+ }
+
+/* COMPUTATION OF THE NEW FUNCTION VALUE AND THE NEW DIRECTIONAL */
+/* DERIVATIVE */
+
+ *kd = 1;
+ *ld = -1;
+ *isys = 1;
+ return;
+L3:
+ if (mode == 0) {
+ *par1 = *p / *po;
+ *par2 = *f - *fo;
+ }
+ if (*iters != 0) {
+ goto L4;
+ }
+ if (*f <= *minf) {
+ *iters = 7;
+ goto L4;
+ } else {
+ l1 = *r__ <= *rmin && *nit != *kit;
+ l2 = *r__ >= *rmax;
+ l3 = *f - *fo <= *tols * *r__ * *po;
+ l5 = *p >= *tolp * *po || (mes2 == 2 && mode == 2);
+ l7 = mes2 <= 2 || mode != 0;
+ m1 = FALSE_;
+ m2 = FALSE_;
+ m3 = l3;
+ if (mes3 >= 1) {
+ m1 = fabs(*p) <= fabs(*po) * .01 && *fo - *f >= fabs(*fo) *
+ 9.9999999999999994e-12;
+ l3 = l3 || m1;
+ }
+ if (mes3 >= 2) {
+ m2 = fabs(*p) <= fabs(*po) * .5 && (d__1 = *fo - *f, fabs(d__1)) <=
+ fabs(*fo) * 2.0000000000000001e-13;
+ l3 = l3 || m2;
+ }
+ *maxst = 0;
+ if (l2) {
+ *maxst = 1;
+ }
+ }
+
+/* TEST ON TERMINATION */
+
+ if (l1 && ! l3) {
+ *iters = 0;
+ goto L4;
+ } else if (l2 && l3 && ! l5) {
+ *iters = 7;
+ goto L4;
+ } else if (m3 && mes1 == 3) {
+ *iters = 5;
+ goto L4;
+ } else if (l3 && l5 && l7) {
+ *iters = 4;
+ goto L4;
+ } else if (*kters < 0 || (*kters == 6 && l7)) {
+ *iters = 6;
+ goto L4;
+ } else if (iabs(*nred) >= *mred) {
+ *iters = -1;
+ goto L4;
+ } else {
+ *rp = *r__;
+ *fp = *f;
+ *pp = *p;
+ mode = max(mode,1);
+ mtyp = iabs(*mes);
+ if (*f >= *maxf) {
+ mtyp = 1;
+ }
+ }
+ if (mode == 1) {
+
+/* INTERVAL CHANGE AFTER EXTRAPOLATION */
+
+ rl = ru;
+ fl = fu;
+ pl = pu;
+ ru = *r__;
+ fu = *f;
+ pu = *p;
+ if (! l3) {
+ *nred = 0;
+ mode = 2;
+ } else if (mes1 == 1) {
+ mtyp = 1;
+ }
+ } else {
+
+/* INTERVAL CHANGE AFTER INTERPOLATION */
+
+ if (! l3) {
+ ru = *r__;
+ fu = *f;
+ pu = *p;
+ } else {
+ rl = *r__;
+ fl = *f;
+ pl = *p;
+ }
+ }
+ goto L2;
+L4:
+ *isys = 0;
+ return;
+} /* luksan_ps1l01__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PULSP3 ALL SYSTEMS 02/12/01
+* PURPOSE :
+* LIMITED STORAGE VARIABLE METRIC UPDATE.
+*
+* PARAMETERS :
+* II N NUMBER OF VARIABLES (NUMBER OF ROWS OF XM).
+* II M NUMBER OF COLUMNS OF XM.
+* II MF MAXIMUM NUMBER OF COLUMNS OF XM.
+* RI XM(N*M) RECTANGULAR MATRIX IN THE PRODUCT FORM SHIFTED BROYDEN
+* METHOD (STORED COLUMNWISE): H-SIGMA*I=XM*TRANS(XM)
+* RO GR(M) MATRIX TRANS(XM)*GO.
+* RU XO(N) VECTORS OF VARIABLES DIFFERENCE XO AND VECTOR XO-TILDE.
+* RU GO(N) GRADIENT DIFFERENCE GO AND VECTOR XM*TRANS(XM)*GO.
+* RI R STEPSIZE PARAMETER.
+* RI PO OLD DIRECTIONAL DERIVATIVE (MULTIPLIED BY R)
+* RU SIG SCALING PARAMETER (ZETA AND SIGMA).
+* IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
+* ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
+* II MET3 CHOICE OF SIGMA (1-CONSTANT, 2-QUADRATIC EQUATION).
+*
+* SUBPROGRAMS USED :
+* S MXDRMM MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR
+* MATRIX BY A VECTOR.
+* S MXDCMU UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX.
+* WITH CONTROLLING OF POSITIVE DEFINITENESS.
+* S MXVDIR VECTOR AUGMENTED BY A SCALED VECTOR.
+* RF MXVDOT DOT PRODUCT OF VECTORS.
+* S MXVSCL SCALING OF A VECTOR.
+*
+* METHOD :
+* SHIFTED BFGS METHOD IN THE PRODUCT FORM.
+*/
+void luksan_pulsp3__(int *n, int *m, int *mf,
+ double *xm, double *gr, double *xo, double *go,
+ double *r__, double *po, double *sig, int *iterh,
+ int *met3)
+{
+ /* System generated locals */
+ double d__1, d__2, d__3, d__4;
+
+ /* Builtin functions */
+
+ /* Local variables */
+ double a, b, c__, aa, bb, ah, den, par, pom;
+
+ /* Parameter adjustments */
+ --go;
+ --xo;
+ --gr;
+ --xm;
+
+ /* Function Body */
+ if (*m >= *mf) {
+ return;
+ }
+ b = luksan_mxvdot__(n, &xo[1], &go[1]);
+ if (b <= 0.) {
+ *iterh = 2;
+ goto L22;
+ }
+ luksan_mxdrmm__(n, m, &xm[1], &go[1], &gr[1]);
+ ah = luksan_mxvdot__(n, &go[1], &go[1]);
+ aa = luksan_mxvdot__(m, &gr[1], &gr[1]);
+ a = aa + ah * *sig;
+ c__ = -(*r__) * *po;
+
+/* DETERMINATION OF THE PARAMETER SIG (SHIFT) */
+
+ par = 1.;
+ pom = b / ah;
+ if (a > 0.) {
+ den = luksan_mxvdot__(n, &xo[1], &xo[1]);
+ if (*met3 <= 4) {
+/* Computing MAX */
+ d__1 = 0., d__2 = 1. - aa / a;
+/* Computing MAX */
+ d__3 = 0., d__4 = 1. - b * b / (den * ah);
+ *sig = sqrt((max(d__1,d__2))) / (sqrt((max(d__3,d__4))) + 1.) *
+ pom;
+ } else {
+/* Computing MAX */
+ d__1 = 0., d__2 = *sig * ah / a;
+/* Computing MAX */
+ d__3 = 0., d__4 = 1. - b * b / (den * ah);
+ *sig = sqrt((max(d__1,d__2))) / (sqrt((max(d__3,d__4))) + 1.) *
+ pom;
+ }
+/* Computing MAX */
+ d__1 = *sig, d__2 = pom * .2;
+ *sig = max(d__1,d__2);
+/* Computing MIN */
+ d__1 = *sig, d__2 = pom * .8;
+ *sig = min(d__1,d__2);
+ } else {
+ *sig = pom * .25;
+ }
+
+/* COMPUTATION OF SHIFTED XO AND SHIFTED B */
+
+ bb = b - ah * *sig;
+ d__1 = -(*sig);
+ luksan_mxvdir__(n, &d__1, &go[1], &xo[1], &xo[1]);
+
+/* BFGS-BASED SHIFTED BFGS UPDATE */
+
+ pom = 1.;
+ d__1 = -1. / bb;
+ luksan_mxdcmu__(n, m, &xm[1], &d__1, &xo[1], &gr[1]);
+ d__1 = sqrt(par / bb);
+ luksan_mxvscl__(n, &d__1, &xo[1], &xm[*n * *m + 1]);
+ ++(*m);
+L22:
+ *iterh = 0;
+ return;
+} /* luksan_pulsp3__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PULVP3 ALL SYSTEMS 03/12/01
+* PURPOSE :
+* RANK-TWO LIMITED-STORAGE VARIABLE-METRIC METHODS IN THE PRODUCT FORM.
+*
+* PARAMETERS :
+* II N NUMBER OF VARIABLES (NUMBER OF ROWS OF XM).
+* II M NUMBER OF COLUMNS OF XM.
+* RI XM(N*M) RECTANGULAR MATRIX IN THE PRODUCT FORM SHIFTED BROYDEN
+* METHOD (STORED COLUMNWISE): H-SIGMA*I=XM*TRANS(XM)
+* RO XR(M) VECTOR TRANS(XM)*H**(-1)*XO.
+* RO GR(M) MATRIX TRANS(XM)*GO.
+* RA S(N) AUXILIARY VECTORS (H**(-1)*XO AND U).
+* RA SO(N) AUXILIARY VECTORS ((H-SIGMA*I)*H**(-1)*XO AND V).
+* RU XO(N) VECTORS OF VARIABLES DIFFERENCE XO AND VECTOR XO-TILDE.
+* RU GO(N) GRADIENT DIFFERENCE GO AND VECTOR XM*TRANS(XM)*GO.
+* RI R STEPSIZE PARAMETER.
+* RI PO OLD DIRECTIONAL DERIVATIVE (MULTIPLIED BY R)
+* RU SIG SCALING PARAMETER (ZETA AND SIGMA).
+* IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
+* ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
+* II MET2 CHOICE OF THE CORRECTION PARAMETER (1-THE UNIT VALUE,
+* 2-THE BALANCING VALUE, 3-THE SQUARE ROOT, 4-THE GEOMETRIC
+* MEAN).
+* II MET3 CHOICE OF THE SHIFT PARAMETER (4-THE FIRST FORMULA,
+* 5-THE SECOND FORMULA).
+* II MET5 CHOICE OF THE METHOD (1-RANK-ONE METHOD, 2-RANK-TWO
+* METHOD).
+*
+* SUBPROGRAMS USED :
+* S MXDRMM MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR
+* MATRIX BY A VECTOR.
+* S MXDCMU UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX.
+* WITH CONTROLLING OF POSITIVE DEFINITENESS. RANK-ONE FORMULA.
+* S MXDCMV UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX.
+* WITH CONTROLLING OF POSITIVE DEFINITENESS. RANK-TWO FORMULA.
+* S MXVDIR VECTOR AUGMENTED BY A SCALED VECTOR.
+* RF MXVDOT DOT PRODUCT OF VECTORS.
+* S MXVLIN LINEAR COMBINATION OF TWO VECTORS.
+* S MXVSCL SCALING OF A VECTOR.
+*
+* METHOD :
+* RANK-ONE LIMITED-STORAGE VARIABLE-METRIC METHOD IN THE PRODUCT FORM.
+*/
+void luksan_pulvp3__(int *n, int *m, double *xm,
+ double *xr, double *gr, double *s, double *so,
+ double *xo, double *go, double *r__, double *po,
+ double *sig, int *iterh, int *met2, int *met3,
+ int *met5)
+{
+ /* System generated locals */
+ double d__1, d__2, d__3, d__4;
+
+ /* Builtin functions */
+
+ /* Local variables */
+ double a, b, c__, aa, bb, cc, ah, den, par, pom, zet;
+
+ /* Parameter adjustments */
+ --go;
+ --xo;
+ --so;
+ --s;
+ --gr;
+ --xr;
+ --xm;
+
+ /* Function Body */
+ zet = *sig;
+
+/* COMPUTATION OF B */
+
+ b = luksan_mxvdot__(n, &xo[1], &go[1]);
+ if (b <= 0.) {
+ *iterh = 2;
+ goto L22;
+ }
+
+/* COMPUTATION OF GR=TRANS(XM)*GO, XR=TRANS(XM)*H**(-1)*XO */
+/* AND S=H**(-1)*XO, SO=(H-SIGMA*I)*H**(-1)*XO. COMPUTATION */
+/* OF AA=GR*GR, BB=GR*XR, CC=XR*XR. COMPUTATION OF A AND C. */
+
+ luksan_mxdrmm__(n, m, &xm[1], &go[1], &gr[1]);
+ luksan_mxvscl__(n, r__, &s[1], &s[1]);
+ luksan_mxdrmm__(n, m, &xm[1], &s[1], &xr[1]);
+ d__1 = -(*sig);
+ luksan_mxvdir__(n, &d__1, &s[1], &xo[1], &so[1]);
+ ah = luksan_mxvdot__(n, &go[1], &go[1]);
+ aa = luksan_mxvdot__(m, &gr[1], &gr[1]);
+ bb = luksan_mxvdot__(m, &gr[1], &xr[1]);
+ cc = luksan_mxvdot__(m, &xr[1], &xr[1]);
+ a = aa + ah * *sig;
+ c__ = -(*r__) * *po;
+
+/* DETERMINATION OF THE PARAMETER SIG (SHIFT) */
+
+ pom = b / ah;
+ if (a > 0.) {
+ den = luksan_mxvdot__(n, &xo[1], &xo[1]);
+ if (*met3 <= 4) {
+/* Computing MAX */
+ d__1 = 0., d__2 = 1. - aa / a;
+/* Computing MAX */
+ d__3 = 0., d__4 = 1. - b * b / (den * ah);
+ *sig = sqrt((max(d__1,d__2))) / (sqrt((max(d__3,d__4))) + 1.) *
+ pom;
+ } else {
+/* Computing MAX */
+ d__1 = 0., d__2 = *sig * ah / a;
+/* Computing MAX */
+ d__3 = 0., d__4 = 1. - b * b / (den * ah);
+ *sig = sqrt((max(d__1,d__2))) / (sqrt((max(d__3,d__4))) + 1.) *
+ pom;
+ }
+/* Computing MAX */
+ d__1 = *sig, d__2 = pom * .2;
+ *sig = max(d__1,d__2);
+/* Computing MIN */
+ d__1 = *sig, d__2 = pom * .8;
+ *sig = min(d__1,d__2);
+ } else {
+ *sig = pom * .25;
+ }
+
+/* COMPUTATION OF SHIFTED XO AND SHIFTED B */
+
+ b -= ah * *sig;
+ d__1 = -(*sig);
+ luksan_mxvdir__(n, &d__1, &go[1], &xo[1], &xo[1]);
+
+/* COMPUTATION OF THE PARAMETER RHO (CORRECTION) */
+
+ if (*met2 <= 1) {
+ par = 1.;
+ } else if (*met2 == 2) {
+ par = *sig * ah / b;
+ } else if (*met2 == 3) {
+ par = sqrt(1. - aa / a);
+ } else if (*met2 == 4) {
+ par = sqrt(sqrt(1. - aa / a) * (*sig * ah / b));
+ } else {
+ par = zet / (zet + *sig);
+ }
+
+/* COMPUTATION OF THE PARAMETER THETA (BFGS) */
+
+ d__1 = sqrt(par * b / cc);
+ pom = copysign(d__1, bb);
+
+/* COMPUTATION OF Q AND P */
+
+ if (*met5 == 1) {
+
+/* RANK ONE UPDATE OF XM */
+
+ luksan_mxvdir__(m, &pom, &xr[1], &gr[1], &xr[1]);
+ luksan_mxvlin__(n, &par, &xo[1], &pom, &so[1], &s[1]);
+ d__1 = -1. / (par * b + pom * bb);
+ luksan_mxdcmu__(n, m, &xm[1], &d__1, &s[1], &xr[1]);
+ } else {
+
+/* RANK TWO UPDATE OF XM */
+
+ d__1 = par / pom - bb / b;
+ luksan_mxvdir__(n, &d__1, &xo[1], &so[1], &s[1]);
+ d__1 = -1. / b;
+ d__2 = -1. / cc;
+ luksan_mxdcmv__(n, m, &xm[1], &d__1, &xo[1], &gr[1], &d__2, &s[1], &xr[1]);
+ }
+L22:
+ *iterh = 0;
+ return;
+} /* luksan_pulvp3__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PYADC0 ALL SYSTEMS 98/12/01
+* PURPOSE :
+* NEW SIMPLE BOUNDS ARE ADDED TO THE ACTIVE SET.
+*
+* PARAMETERS :
+* II NF DECLARED NUMBER OF VARIABLES.
+* II N REDUCED NUMBER OF VARIABLES.
+* RI X(NF) VECTOR OF VARIABLES.
+* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
+* RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
+* RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
+* IO INEW NUMBER OF ACTIVE CONSTRAINTS.
+*/
+void luksan_pyadc0__(int *nf, int *n, double *x,
+ int *ix, double *xl, double *xu, int *inew)
+{
+ /* System generated locals */
+ int i__1;
+
+ /* Local variables */
+ int i__, ii, ixi;
+
+ /* Parameter adjustments */
+ --ix;
+ --x;
+ --xl;
+ --xu;
+
+ /* Function Body */
+ *n = *nf;
+ *inew = 0;
+ i__1 = *nf;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ ii = ix[i__];
+ ixi = iabs(ii);
+ if (ixi >= 5) {
+ ix[i__] = -ixi;
+ } else if ((ixi == 1 || ixi == 3 || ixi == 4) && x[i__] <= xl[i__]) {
+ x[i__] = xl[i__];
+ if (ixi == 4) {
+ ix[i__] = -3;
+ } else {
+ ix[i__] = -ixi;
+ }
+ --(*n);
+ if (ii > 0) {
+ ++(*inew);
+ }
+ } else if ((ixi == 2 || ixi == 3 || ixi == 4) && x[i__] >= xu[i__]) {
+ x[i__] = xu[i__];
+ if (ixi == 3) {
+ ix[i__] = -4;
+ } else {
+ ix[i__] = -ixi;
+ }
+ --(*n);
+ if (ii > 0) {
+ ++(*inew);
+ }
+ }
+/* L1: */
+ }
+ return;
+} /* luksan_pyadc0__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PYFUT1 ALL SYSTEMS 98/12/01
+* PURPOSE :
+* TERMINATION CRITERIA AND TEST ON RESTART.
+*
+* PARAMETERS :
+* II N ACTUAL NUMBER OF VARIABLES.
+* RI F NEW VALUE OF THE OBJECTIVE FUNCTION.
+* RI FO OLD VALUE OF THE OBJECTIVE FUNCTION.
+* RI UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
+* RO GMAX NORM OF THE TRANSFORMED GRADIENT.
+* RI DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
+* RI TOLX LOWER BOUND FOR STEPLENGTH.
+* RI TOLF LOWER BOUND FOR FUNCTION DECREASE.
+* RI TOLB LOWER BOUND FOR FUNCTION VALUE.
+* RI TOLG LOWER BOUND FOR GRADIENT.
+* II KD DEGREE OF REQUIRED DERIVATIVES.
+* IU NIT ACTUAL NUMBER OF ITERATIONS.
+* II KIT NUMBER OF THE ITERATION AFTER RESTART.
+* II MIT MAXIMUM NUMBER OF ITERATIONS.
+* IU NFV ACTUAL NUMBER OF COMPUTED FUNCTION VALUES.
+* II MFV MAXIMUM NUMBER OF COMPUTED FUNCTION VALUES.
+* IU NFG ACTUAL NUMBER OF COMPUTED GRADIENT VALUES.
+* II MFG MAXIMUM NUMBER OF COMPUTED GRADIENT VALUES.
+* IU NTESX ACTUAL NUMBER OF TESTS ON STEPLENGTH.
+* II MTESX MAXIMUM NUMBER OF TESTS ON STEPLENGTH.
+* IU NTESF ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE.
+* II MTESF MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE.
+* II IRES1 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
+* IRES1*N+IRES2 ITERATIONS.
+* II IRES2 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
+* IRES1*N+IRES2 ITERATIONS.
+* IU IREST RESTART INDICATOR. RESTART IS PERFORMED IF IREST>0.
+* II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
+* ITERS=0 FOR ZERO STEP.
+* IO ITERM TERMINATION INDICATOR. ITERM=1-TERMINATION AFTER MTESX
+* UNSUFFICIENT STEPLENGTHS. ITERM=2-TERMINATION AFTER MTESF
+* UNSUFFICIENT FUNCTION DECREASES. ITERM=3-TERMINATION ON LOWER
+* BOUND FOR FUNCTION VALUE. ITERM=4-TERMINATION ON LOWER BOUND
+* FOR GRADIENT. ITERM=11-TERMINATION AFTER MAXIMUM NUMBER OF
+* ITERATIONS. ITERM=12-TERMINATION AFTER MAXIMUM NUMBER OF
+* COMPUTED FUNCTION VALUES.
+*/
+void luksan_pyfut1__(int *n, double *f, double *
+ fo, double *umax, double *gmax, double *dmax__,
+ double *tolx, double *tolf, double *tolb, double *
+ tolg, int *kd, int *nit, int *kit, int *mit, int *
+ nfv, int *mfv, int *nfg, int *mfg, int *ntesx,
+ int *mtesx, int *ntesf, int *mtesf, int *ites,
+ int *ires1, int *ires2, int *irest, int *iters,
+ int *iterm)
+{
+ /* System generated locals */
+ double d__1, d__2;
+
+ /* Builtin functions */
+
+ /* Local variables */
+ double temp;
+
+ if (*iterm < 0) {
+ return;
+ }
+ if (*ites <= 0) {
+ goto L1;
+ }
+ if (*iters == 0) {
+ goto L1;
+ }
+ if (*nit <= 0) {
+/* Computing MIN */
+ d__1 = sqrt((fabs(*f))), d__2 = fabs(*f) / 10.;
+ *fo = *f + min(d__1,d__2);
+ }
+ if (*f <= *tolb) {
+ *iterm = 3;
+ return;
+ }
+ if (*kd > 0) {
+ if (*gmax <= *tolg && *umax <= *tolg) {
+ *iterm = 4;
+ return;
+ }
+ }
+ if (*nit <= 0) {
+ *ntesx = 0;
+ *ntesf = 0;
+ }
+ if (*dmax__ <= *tolx) {
+ *iterm = 1;
+ ++(*ntesx);
+ if (*ntesx >= *mtesx) {
+ return;
+ }
+ } else {
+ *ntesx = 0;
+ }
+/* Computing MAX */
+ d__2 = fabs(*f);
+ temp = (d__1 = *fo - *f, fabs(d__1)) / max(d__2,1.);
+ if (temp <= *tolf) {
+ *iterm = 2;
+ ++(*ntesf);
+ if (*ntesf >= *mtesf) {
+ return;
+ }
+ } else {
+ *ntesf = 0;
+ }
+L1:
+ if (*nit >= *mit) {
+ *iterm = 11;
+ return;
+ }
+ if (*nfv >= *mfv) {
+ *iterm = 12;
+ return;
+ }
+ if (*nfg >= *mfg) {
+ *iterm = 13;
+ return;
+ }
+ *iterm = 0;
+ if (*n > 0 && *nit - *kit >= *ires1 * *n + *ires2) {
+ *irest = max(*irest,1);
+ }
+ ++(*nit);
+ return;
+} /* luksan_pyfut1__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PYRMC0 ALL SYSTEMS 98/12/01
+* PURPOSE :
+* OLD SIMPLE BOUND IS REMOVED FROM THE ACTIVE SET. TRANSFORMED
+* GRADIENT OF THE OBJECTIVE FUNCTION IS UPDATED.
+*
+* PARAMETERS :
+* II NF DECLARED NUMBER OF VARIABLES.
+* II N REDUCED NUMBER OF VARIABLES.
+* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
+* RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION.
+* RI EPS8 TOLERANCE FOR CONSTRAINT TO BE REMOVED.
+* RI UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
+* RI GMAX NORM OF THE TRANSFORMED GRADIENT.
+* RO RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
+* II IOLD NUMBER OF REMOVED CONSTRAINTS.
+* IU IREST RESTART INDICATOR.
+*/
+void luksan_pyrmc0__(int *nf, int *n, int *ix,
+ double *g, double *eps8, double *umax, double *gmax,
+ double *rmax, int *iold, int *irest)
+{
+ /* System generated locals */
+ int i__1, i__2, i__3;
+
+ /* Local variables */
+ int i__, ixi;
+
+ /* Parameter adjustments */
+ --g;
+ --ix;
+
+ /* Function Body */
+ if (*n == 0 || *rmax > 0.) {
+ if (*umax > *eps8 * *gmax) {
+ *iold = 0;
+ i__1 = *nf;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ ixi = ix[i__];
+ if (ixi >= 0) {
+ } else if (ixi <= -5) {
+ } else if ((ixi == -1 || ixi == -3) && -g[i__] <= 0.) {
+ } else if ((ixi == -2 || ixi == -4) && g[i__] <= 0.) {
+ } else {
+ ++(*iold);
+/* Computing MIN */
+ i__3 = (i__2 = ix[i__], iabs(i__2));
+ ix[i__] = min(i__3,3);
+ if (*rmax == 0.) {
+ goto L2;
+ }
+ }
+/* L1: */
+ }
+L2:
+ if (*iold > 1) {
+ *irest = max(*irest,1);
+ }
+ }
+ }
+ return;
+} /* luksan_pyrmc0__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PYTRCD ALL SYSTEMS 98/12/01
+* PURPOSE :
+* VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
+* AND SCALED AND REDUCED. TEST VALUE DMAX IS DETERMINED.
+*
+* PARAMETERS :
+* II NF DECLARED NUMBER OF VARIABLES.
+* RI X(NF) VECTOR OF VARIABLES.
+* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
+* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE.
+* RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION.
+* RU GO(NF) GRADIENTS DIFFERENCE.
+* RO R VALUE OF THE STEPSIZE PARAMETER.
+* RO F NEW VALUE OF THE OBJECTIVE FUNCTION.
+* RI FO OLD VALUE OF THE OBJECTIVE FUNCTION.
+* RO P NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
+* RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
+* RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
+* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
+* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
+* IO KD DEGREE OF REQUIRED DERIVATIVES.
+* IO LD DEGREE OF COMPUTED DERIVATIVES.
+* II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
+* ITERS=0 FOR ZERO STEP.
+*
+* SUBPROGRAMS USED :
+* S MXVDIF DIFFERENCE OF TWO VECTORS.
+* S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
+* SUBSTRACTED ONE.
+*/
+void luksan_pytrcd__(int *nf, double *x, int *ix,
+ double *xo, double *g, double *go, double *r__,
+ double *f, double *fo, double *p, double *po,
+ double *dmax__, int *kbf, int *kd, int *ld, int *
+ iters)
+{
+ /* System generated locals */
+ int i__1;
+ double d__1, d__2, d__3, d__4, d__5;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --go;
+ --g;
+ --xo;
+ --ix;
+ --x;
+
+ /* Function Body */
+ if (*iters > 0) {
+ luksan_mxvdif__(nf, &x[1], &xo[1], &xo[1]);
+ luksan_mxvdif__(nf, &g[1], &go[1], &go[1]);
+ *po = *r__ * *po;
+ *p = *r__ * *p;
+ } else {
+ *f = *fo;
+ *p = *po;
+ luksan_mxvsav__(nf, &x[1], &xo[1]);
+ luksan_mxvsav__(nf, &g[1], &go[1]);
+ *ld = *kd;
+ }
+ *dmax__ = 0.;
+ i__1 = *nf;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (*kbf > 0) {
+ if (ix[i__] < 0) {
+ xo[i__] = 0.;
+ go[i__] = 0.;
+ goto L1;
+ }
+ }
+/* Computing MAX */
+/* Computing MAX */
+ d__5 = (d__2 = x[i__], fabs(d__2));
+ d__3 = *dmax__, d__4 = (d__1 = xo[i__], fabs(d__1)) / max(d__5,1.);
+ *dmax__ = max(d__3,d__4);
+L1:
+ ;
+ }
+ return;
+} /* luksan_pytrcd__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PYTRCG ALL SYSTEMS 99/12/01
+* PURPOSE :
+* GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. TEST VALUES
+* GMAX AND UMAX ARE COMPUTED.
+*
+* PARAMETERS :
+* II NF DECLARED NUMBER OF VARIABLES.
+* II N ACTUAL NUMBER OF VARIABLES.
+* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
+* RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION.
+* RI UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
+* RI GMAX NORM OF THE TRANSFORMED GRADIENT.
+* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
+* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
+* II IOLD INDEX OF THE REMOVED CONSTRAINT.
+*
+* SUBPROGRAMS USED :
+* RF MXVMAX L-INFINITY NORM OF A VECTOR.
+*/
+void luksan_pytrcg__(int *nf, int *n, int *ix,
+ double *g, double *umax, double *gmax, int *kbf,
+ int *iold)
+{
+ /* System generated locals */
+ int i__1;
+ double d__1, d__2;
+
+ /* Local variables */
+ int i__;
+ double temp;
+
+ /* Parameter adjustments */
+ --g;
+ --ix;
+
+ /* Function Body */
+ if (*kbf > 0) {
+ *gmax = 0.;
+ *umax = 0.;
+ *iold = 0;
+ i__1 = *nf;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = g[i__];
+ if (ix[i__] >= 0) {
+/* Computing MAX */
+ d__1 = *gmax, d__2 = fabs(temp);
+ *gmax = max(d__1,d__2);
+ } else if (ix[i__] <= -5) {
+ } else if ((ix[i__] == -1 || ix[i__] == -3) && *umax + temp >= 0.)
+ {
+ } else if ((ix[i__] == -2 || ix[i__] == -4) && *umax - temp >= 0.)
+ {
+ } else {
+ *iold = i__;
+ *umax = fabs(temp);
+ }
+/* L1: */
+ }
+ } else {
+ *umax = 0.;
+ *gmax = luksan_mxvmax__(nf, &g[1]);
+ }
+ *n = *nf;
+ return;
+} /* luksan_pytrcg__ */
+
+/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
+/* SUBROUTINE PYTRCS ALL SYSTEMS 98/12/01
+* PURPOSE :
+* SCALED AND REDUCED DIRECTION VECTOR IS BACK TRANSFORMED. VECTORS
+* X,G AND VALUES F,P ARE SAVED.
+*
+* PARAMETERS :
+* II NF DECLARED NUMBER OF VARIABLES.
+* RI X(NF) VECTOR OF VARIABLES.
+* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
+* RO XO(NF) SAVED VECTOR OF VARIABLES.
+* RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
+* RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
+* RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION.
+* RO GO(NF) SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
+* RO S(NF) DIRECTION VECTOR.
+* RO RO SAVED VALUE OF THE STEPSIZE PARAMETER.
+* RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
+* RU FO SAVED VALUE OF THE OBJECTIVE FUNCTION.
+* RI F VALUE OF THE OBJECTIVE FUNCTION.
+* RO PO SAVED VALUE OF THE DIRECTIONAL DERIVATIVE.
+* RI P VALUE OF THE DIRECTIONAL DERIVATIVE.
+* RO RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
+* RI ETA9 MAXIMUM FOR REAL NUMBERS.
+* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
+* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
+*
+* SUBPROGRAMS USED :
+* S MXVCOP COPYING OF A VECTOR.
+*/
+void luksan_pytrcs__(int *nf, double *x, int *ix,
+ double *xo, double *xl, double *xu, double *g,
+ double *go, double *s, double *ro, double *fp,
+ double *fo, double *f, double *po, double *p,
+ double *rmax, double *eta9, int *kbf)
+{
+ /* System generated locals */
+ int i__1;
+ double d__1, d__2;
+
+ /* Local variables */
+ int i__;
+
+ /* Parameter adjustments */
+ --s;
+ --go;
+ --g;
+ --xu;
+ --xl;
+ --xo;
+ --ix;
+ --x;
+
+ /* Function Body */
+ *fp = *fo;
+ *ro = 0.;
+ *fo = *f;
+ *po = *p;
+ luksan_mxvcop__(nf, &x[1], &xo[1]);
+ luksan_mxvcop__(nf, &g[1], &go[1]);
+ if (*kbf > 0) {
+ i__1 = *nf;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (ix[i__] < 0) {
+ s[i__] = 0.;
+ } else {
+ if (ix[i__] == 1 || ix[i__] >= 3) {
+ if (s[i__] < -1. / *eta9) {
+/* Computing MIN */
+ d__1 = *rmax, d__2 = (xl[i__] - x[i__]) / s[i__];
+ *rmax = min(d__1,d__2);
+ }
+ }
+ if (ix[i__] == 2 || ix[i__] >= 3) {
+ if (s[i__] > 1. / *eta9) {
+/* Computing MIN */
+ d__1 = *rmax, d__2 = (xu[i__] - x[i__]) / s[i__];
+ *rmax = min(d__1,d__2);
+ }
+ }
+ }
+/* L1: */
+ }
+ }
+ return;
+} /* luksan_pytrcs__ */
+
~ianmdlvl / nlopt.git / commitdiff