From: stevenj Date: Mon, 3 Sep 2007 19:38:15 +0000 (-0400) Subject: added luksan/ directory, which will eventually contain several routines written by... X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=cab3a453b4757cdbf14c4e98076cd3cdc69e2dc9;p=nlopt.git added luksan/ directory, which will eventually contain several routines written by Prof. Luksan and made available under the LGPL. currently, I've converted the PLIN.FOR subroutine, which is a reimplementation of the LBFGS algorithm (and removed the lbfgs/ directory, since the code turns out to be non-free) darcs-hash:20070903193815-c8de0-7428452045bf7147a8424cbd56af89de2e3efa66.gz --- diff --git a/Makefile.am b/Makefile.am index af4ce1c..86cbf1d 100644 --- a/Makefile.am +++ b/Makefile.am @@ -3,13 +3,13 @@ lib_LTLIBRARIES = libnlopt.la 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@ diff --git a/api/Makefile.am b/api/Makefile.am index 08bffdc..30ba1bf 100644 --- a/api/Makefile.am +++ b/api/Makefile.am @@ -1,4 +1,4 @@ -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 diff --git a/api/nlopt.c b/api/nlopt.c index 667a5b1..f3aaca5 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -115,10 +115,10 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) #include "stogo.h" -#include "l-bfgs-b.h" - #include "cdirect.h" +#include "luksan.h" + /*************************************************************************/ /* for "hybrid" algorithms that combine global search with some @@ -301,36 +301,8 @@ static nlopt_result nlopt_minimize_( &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; diff --git a/configure.ac b/configure.ac index 534abab..0b1e504 100644 --- a/configure.ac +++ b/configure.ac @@ -137,9 +137,9 @@ AC_CONFIG_FILES([ direct/Makefile cdirect/Makefile stogo/Makefile - lbfgs/Makefile subplex/Makefile praxis/Makefile + luksan/Makefile test/Makefile ]) diff --git a/lbfgs/Makefile.am b/lbfgs/Makefile.am index fbe9986..a8db963 100644 --- a/lbfgs/Makefile.am +++ b/lbfgs/Makefile.am @@ -1,4 +1,4 @@ noinst_LTLIBRARIES = liblbfgs.la liblbfgs_la_SOURCES = ap.cpp ap.h l-bfgs-b.cpp l-bfgs-b.h -EXTRA_DIST = +EXTRA_DIST = README diff --git a/lbfgs/README b/lbfgs/README new file mode 100644 index 0000000..80c0d27 --- /dev/null +++ b/lbfgs/README @@ -0,0 +1,17 @@ +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. diff --git a/luksan/COPYRIGHT b/luksan/COPYRIGHT new file mode 100644 index 0000000..8e2f44e --- /dev/null +++ b/luksan/COPYRIGHT @@ -0,0 +1,45 @@ +--------------------------------------------------------------------------- + +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. + + diff --git a/luksan/Makefile.am b/luksan/Makefile.am new file mode 100644 index 0000000..2cfa02a --- /dev/null +++ b/luksan/Makefile.am @@ -0,0 +1,6 @@ +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 diff --git a/luksan/README b/luksan/README new file mode 100644 index 0000000..3ef4864 --- /dev/null +++ b/luksan/README @@ -0,0 +1,23 @@ +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.] diff --git a/luksan/luksan.h b/luksan/luksan.h new file mode 100644 index 0000000..6d36e68 --- /dev/null +++ b/luksan/luksan.h @@ -0,0 +1,118 @@ +#ifndef LUKSAN_H +#define LUKSAN_H + +#include +#include + +#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 */ diff --git a/luksan/mssubs.c b/luksan/mssubs.c new file mode 100644 index 0000000..0d2cd14 --- /dev/null +++ b/luksan/mssubs.c @@ -0,0 +1,991 @@ +#include +#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__ */ + diff --git a/luksan/plis.c b/luksan/plis.c new file mode 100644 index 0000000..e5c8083 --- /dev/null +++ b/luksan/plis.c @@ -0,0 +1,523 @@ +#include +#include +#include +#include +#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; + } +} diff --git a/luksan/plis.txt b/luksan/plis.txt new file mode 100644 index 0000000..d1eb798 --- /dev/null +++ b/luksan/plis.txt @@ -0,0 +1,302 @@ +*********************************************************************** +* * +* 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. + \ No newline at end of file diff --git a/luksan/pssubs.c b/luksan/pssubs.c new file mode 100644 index 0000000..c6ae5b4 --- /dev/null +++ b/luksan/pssubs.c @@ -0,0 +1,1257 @@ +#include +#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__ */ + diff --git a/luksan/v999-07.pdf b/luksan/v999-07.pdf new file mode 100644 index 0000000..a25fd90 Binary files /dev/null and b/luksan/v999-07.pdf differ