From 052a636f9c21466d03ee54037fb9348c6fafb38e Mon Sep 17 00:00:00 2001 From: relf Date: Sun, 8 Oct 2023 21:14:39 +0200 Subject: [PATCH] Initial NLopt 2.7.1 import --- .github/dependabot.yml | 10 + .github/workflows/tests.yml | 35 + .gitignore | 12 + Cargo.toml | 21 + README.md | 12 + slsqp/CMakeLists.txt | 44 + slsqp/README.md | 11 + slsqp/compile_commands.json | 7 + slsqp/nlopt-util.h | 140 ++ slsqp/nlopt.h | 350 +++++ slsqp/nlopt_config.h | 139 ++ slsqp/nlopt_config.h.in | 139 ++ slsqp/rescale.c | 96 ++ slsqp/slsqp.c | 2638 ++++++++++++++++++++++++++++++++ slsqp/slsqp.h | 62 + slsqp/slsqp.rs | 2868 +++++++++++++++++++++++++++++++++++ slsqp/stop.c | 266 ++++ slsqp/timer.c | 93 ++ src/main.rs | 3 + 19 files changed, 6946 insertions(+) create mode 100644 .github/dependabot.yml create mode 100644 .github/workflows/tests.yml create mode 100644 .gitignore create mode 100644 Cargo.toml create mode 100644 README.md create mode 100644 slsqp/CMakeLists.txt create mode 100644 slsqp/README.md create mode 100644 slsqp/compile_commands.json create mode 100644 slsqp/nlopt-util.h create mode 100644 slsqp/nlopt.h create mode 100644 slsqp/nlopt_config.h create mode 100644 slsqp/nlopt_config.h.in create mode 100644 slsqp/rescale.c create mode 100644 slsqp/slsqp.c create mode 100644 slsqp/slsqp.h create mode 100644 slsqp/slsqp.rs create mode 100644 slsqp/stop.c create mode 100644 slsqp/timer.c create mode 100644 src/main.rs diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..dfd0e30 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,10 @@ +# Set update schedule for GitHub Actions + +version: 2 +updates: + + - package-ecosystem: "github-actions" + directory: "/" + schedule: + # Check for updates to GitHub Actions every week + interval: "weekly" diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..11da360 --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,35 @@ +name: tests + +on: [push, pull_request] + +env: + CARGO_TERM_COLOR: always + +jobs: + testing: + name: testing-${{ matrix.toolchain }}-${{ matrix.os }} + runs-on: ${{ matrix.os }} + strategy: + matrix: + toolchain: + - stable + os: + - ubuntu-latest + - windows-latest + + steps: + - name: Checkout sources + uses: actions/checkout@v4 + + - name: Install toolchain + uses: actions-rs/toolchain@v1 + with: + profile: minimal + toolchain: ${{ matrix.toolchain }} + override: true + + - name: Run cargo test in release mode + uses: actions-rs/cargo@v1 + with: + command: test + args: --all --release diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..457bfbf --- /dev/null +++ b/.gitignore @@ -0,0 +1,12 @@ +# Generated by Cargo +# will have compiled files and executables +/target/ + +# Remove Cargo.lock from gitignore if creating an executable, leave it for libraries +# More information here https://doc.rust-lang.org/cargo/guide/cargo-toml-vs-cargo-lock.html +Cargo.lock + +# These are backup files generated by rustfmt +**/*.rs.bk + +/build \ No newline at end of file diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..3edfbcc --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,21 @@ +[package] +authors = ["Remi Lafage "] +name = "slsqp" +version = "0.0.1" +edition = "2021" +license-file = "LICENSE.md" +description = "SLSQP optimizer for Rust" +readme = "README.md" +repository = "https://github.com/relf/slsqp/" +keywords = ["optimizer", "optimization", "constrained", "derivative-free"] +categories = ["algorithms", "api-bindings", "science"] +documentation = "https://docs.rs/slsqp" + +[dependencies] +libc = "0.2" +argmin = "0.8.0" +serde = "1" +instant = "0.1" + +[dev-dependencies] +approx = "0.5" diff --git a/README.md b/README.md new file mode 100644 index 0000000..6a152c8 --- /dev/null +++ b/README.md @@ -0,0 +1,12 @@ +# slsqp + +[![tests](https://github.com/relf/slsqp/workflows/tests/badge.svg)](https://github.com/relf/slsqp/actions?query=workflow%3Atests) +[![crates.io](https://img.shields.io/crates/v/slsqp)](https://crates.io/crates/slsqp) +[![docs](https://docs.rs/slsqp/badge.svg)](https://docs.rs/slsqp) + +## slsqp 0.1.x + +Rust wrapper for SLSQP optimizer (SLSQP stands for ). +slsqp Rust code was generated from NLopt 2.7.1 slsqp C code thanks to c2rust trabspiler then +manually edited to make it work. + diff --git a/slsqp/CMakeLists.txt b/slsqp/CMakeLists.txt new file mode 100644 index 0000000..c0e0d28 --- /dev/null +++ b/slsqp/CMakeLists.txt @@ -0,0 +1,44 @@ +cmake_minimum_required (VERSION 3.0.0) +project (slsqp) + +include (CheckIncludeFiles) +include (CheckFunctionExists) +include (CheckTypeSize) +include (CheckCCompilerFlag) +include (CheckCXXSymbolExists) +include (CheckCSourceCompiles) +include (CheckCXXCompilerFlag) +include (CheckLibraryExists) + +check_include_file (getopt.h HAVE_GETOPT_H) +check_include_file (unistd.h HAVE_UNISTD_H) +check_include_file (stdint.h HAVE_STDINT_H) +check_include_file (time.h HAVE_TIME_H) +check_include_file (sys/time.h HAVE_SYS_TIME_H) +if (HAVE_TIME_H AND HAVE_SYS_TIME_H) + set (TIME_WITH_SYS_TIME TRUE) +endif () +check_function_exists (getpid HAVE_GETPID) +check_function_exists (syscall HAVE_GETTID_SYSCALL) +check_function_exists (isinf HAVE_ISINF) +check_function_exists (isnan HAVE_ISNAN) +check_function_exists (gettimeofday HAVE_GETTIMEOFDAY) +check_function_exists (qsort_r HAVE_QSORT_R) +check_function_exists (time HAVE_TIME) +check_function_exists (copysign HAVE_COPYSIGN) +check_function_exists (getopt HAVE_GETOPT) +check_type_size ("uint32_t" SIZEOF_UINT32_T) +set (HAVE_UINT32_T ${SIZEOF_UINT32_T}) +check_type_size ("unsigned int" SIZEOF_UNSIGNED_INT) +check_type_size ("unsigned long" SIZEOF_UNSIGNED_LONG) + +check_library_exists ("m" sqrt "" HAVE_LIBM) +if (HAVE_LIBM) + set (M_LIBRARY m) +endif() + +configure_file (${CMAKE_CURRENT_SOURCE_DIR}/nlopt_config.h.in ${CMAKE_CURRENT_SOURCE_DIR}/nlopt_config.h IMMEDIATE) + +add_library (slsqp slsqp.c) + +target_include_directories (slsqp PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) \ No newline at end of file diff --git a/slsqp/README.md b/slsqp/README.md new file mode 100644 index 0000000..ae8c615 --- /dev/null +++ b/slsqp/README.md @@ -0,0 +1,11 @@ +To generate compile commands json file required by c2rust in top directory + +> mkdir build +> cd build +> cmake -DCMAKE_EXPORT_COMPILE_COMMANDS=1 ../slsqp +> cp compile_commands.json ../slsqp/compile_commands.json + +Then rust generation from C + +> cd ../slsqp +> c2rust transpile compile_commands.json \ No newline at end of file diff --git a/slsqp/compile_commands.json b/slsqp/compile_commands.json new file mode 100644 index 0000000..cf52991 --- /dev/null +++ b/slsqp/compile_commands.json @@ -0,0 +1,7 @@ +[ +{ + "directory": "/home/rlafage/workspace/slsqp/build", + "command": "/usr/bin/cc -I/home/rlafage/workspace/slsqp/slsqp -o CMakeFiles/slsqp.dir/slsqp.c.o -c /home/rlafage/workspace/slsqp/slsqp/slsqp.c", + "file": "/home/rlafage/workspace/slsqp/slsqp/slsqp.c" +} +] \ No newline at end of file diff --git a/slsqp/nlopt-util.h b/slsqp/nlopt-util.h new file mode 100644 index 0000000..135eb9f --- /dev/null +++ b/slsqp/nlopt-util.h @@ -0,0 +1,140 @@ +/* Copyright (c) 2007-2014 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#ifndef NLOPT_UTIL_H +#define NLOPT_UTIL_H + +#include +#include +#include +#include "nlopt_config.h" + +#include "nlopt.h" + +/* workaround for Solaris + gcc 3.4.x bug (see configure.ac) */ +#if defined(__GNUC__) && defined(REPLACEMENT_HUGE_VAL) +# undef HUGE_VAL +# define HUGE_VAL REPLACEMENT_HUGE_VAL +#endif + +#ifndef HAVE_COPYSIGN + /* not quite right for y == -0, but good enough for us */ +# define copysign(x, y) ((y) < 0 ? -fabs(x) : fabs(x)) +#elif !defined(__cplusplus) && __STDC_VERSION__ < 199901 +extern double copysign(double x, double y); /* may not be declared in C89 */ +#endif + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + + int nlopt_isinf(double x); + int nlopt_isfinite(double x); + int nlopt_istiny(double x); + int nlopt_isnan(double x); + +/* re-entrant qsort, uses the BSD convention */ + extern void nlopt_qsort_r(void *base_, size_t nmemb, size_t size, void *thunk, int (*compar) (void *, const void *, const void *)); + +/* seconds timer */ + extern double nlopt_seconds(void); + extern unsigned long nlopt_time_seed(void); + +/* pseudorandom number generation by Mersenne twister algorithm */ + extern void nlopt_init_genrand(unsigned long s); + extern double nlopt_urand(double a, double b); + extern int nlopt_iurand(int n); + extern double nlopt_nrand(double mean, double stddev); + +/* Sobol' low-discrepancy-sequence generation */ + typedef struct nlopt_soboldata_s *nlopt_sobol; + extern nlopt_sobol nlopt_sobol_create(unsigned sdim); + extern void nlopt_sobol_destroy(nlopt_sobol s); + extern void nlopt_sobol_next01(nlopt_sobol s, double *x); + extern void nlopt_sobol_next(nlopt_sobol s, double *x, const double *lb, const double *ub); + extern void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x); + +/* stopping criteria */ + typedef struct { + unsigned n; + double minf_max; + double ftol_rel; + double ftol_abs; + double xtol_rel; + const double *xtol_abs; + const double *x_weights; + int *nevals_p, maxeval; + double maxtime, start; + int *force_stop; + char **stop_msg; /* pointer to msg string to update */ + } nlopt_stopping; + extern int nlopt_stop_f(const nlopt_stopping * stop, double f, double oldf); + extern int nlopt_stop_ftol(const nlopt_stopping * stop, double f, double oldf); + extern int nlopt_stop_x(const nlopt_stopping * stop, const double *x, const double *oldx); + extern int nlopt_stop_dx(const nlopt_stopping * stop, const double *x, const double *dx); + extern int nlopt_stop_xs(const nlopt_stopping * stop, const double *xs, const double *oldxs, const double *scale_min, const double *scale_max); + extern int nlopt_stop_evals(const nlopt_stopping * stop); + extern int nlopt_stop_time_(double start, double maxtime); + extern int nlopt_stop_time(const nlopt_stopping * stop); + extern int nlopt_stop_evalstime(const nlopt_stopping * stop); + extern int nlopt_stop_forced(const nlopt_stopping * stop); + +/* like vsprintf, but reallocs p to whatever size is needed */ + extern char *nlopt_vsprintf(char *p, const char *format, va_list ap); + extern void nlopt_stop_msg(const nlopt_stopping * s, const char *format, ...) +#ifdef __GNUC__ + __attribute__ ((format(printf, 2, 3))) +#endif + ; + +/* for local optimizations, temporarily setting eval/time limits */ + extern nlopt_result nlopt_optimize_limited(nlopt_opt opt, double *x, double *minf, int maxevals, double maxtime); + +/* data structure for nonlinear inequality or equality constraint + (f <= 0 or f = 0, respectively). tol (>= 0) is a tolerance + that is used for stopping criteria -- the point is considered + "feasible" for purposes of stopping if the constraint is violated + by at most tol. */ + typedef struct { + unsigned m; /* dimensional of constraint: mf maps R^n -> R^m */ + nlopt_func f; /* one-dimensional constraint, requires m == 1 */ + nlopt_mfunc mf; + nlopt_precond pre; /* preconditioner for f (NULL if none or if mf) */ + void *f_data; + double *tol; + } nlopt_constraint; + + extern unsigned nlopt_count_constraints(unsigned p, const nlopt_constraint * c); + extern unsigned nlopt_max_constraint_dim(unsigned p, const nlopt_constraint * c); + extern void nlopt_eval_constraint(double *result, double *grad, const nlopt_constraint * c, unsigned n, const double *x); + +/* rescale.c: */ + double *nlopt_compute_rescaling(unsigned n, const double *dx); + double *nlopt_new_rescaled(unsigned n, const double *s, const double *x); + void nlopt_rescale(unsigned n, const double *s, const double *x, double *xs); + void nlopt_unscale(unsigned n, const double *s, const double *x, double *xs); + void nlopt_reorder_bounds(unsigned n, double *lb, double *ub); + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ +#endif diff --git a/slsqp/nlopt.h b/slsqp/nlopt.h new file mode 100644 index 0000000..fbfeed3 --- /dev/null +++ b/slsqp/nlopt.h @@ -0,0 +1,350 @@ +/* Copyright (c) 2007-2014 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#ifndef NLOPT_H +#define NLOPT_H + +#include /* for ptrdiff_t and size_t */ + +/* Change 0 to 1 to use stdcall convention under Win32 */ +#if 0 && (defined(_WIN32) || defined(__WIN32__)) +# if defined(__GNUC__) +# define NLOPT_STDCALL __attribute__((stdcall)) +# elif defined(_MSC_VER) || defined(_ICC) || defined(_STDCALL_SUPPORTED) +# define NLOPT_STDCALL __stdcall +# else +# define NLOPT_STDCALL +# endif +#else +# define NLOPT_STDCALL +#endif + +/* for Windows compilers, you should add a line + #define NLOPT_DLL + when using NLopt from a DLL, in order to do the proper + Windows importing nonsense. */ +#if defined(NLOPT_DLL) && (defined(_WIN32) || defined(__WIN32__)) && !defined(__LCC__) +/* annoying Windows syntax for calling functions in a DLL */ +# if defined(NLOPT_DLL_EXPORT) +# define NLOPT_EXTERN(T) extern __declspec(dllexport) T NLOPT_STDCALL +# else +# define NLOPT_EXTERN(T) extern __declspec(dllimport) T NLOPT_STDCALL +# endif +#else +# define NLOPT_EXTERN(T) extern T NLOPT_STDCALL +#endif + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +typedef double (*nlopt_func) (unsigned n, const double *x, + double *gradient, /* NULL if not needed */ + void *func_data); + +typedef void (*nlopt_mfunc) (unsigned m, double *result, unsigned n, const double *x, + double *gradient, /* NULL if not needed */ + void *func_data); + +/* A preconditioner, which preconditions v at x to return vpre. + (The meaning of "preconditioning" is algorithm-dependent.) */ +typedef void (*nlopt_precond) (unsigned n, const double *x, const double *v, double *vpre, void *data); + +typedef enum { + /* Naming conventions: + + NLOPT_{G/L}{D/N}_* + = global/local derivative/no-derivative optimization, + respectively + + *_RAND algorithms involve some randomization. + + *_NOSCAL algorithms are *not* scaled to a unit hypercube + (i.e. they are sensitive to the units of x) + */ + + NLOPT_GN_DIRECT = 0, + NLOPT_GN_DIRECT_L, + NLOPT_GN_DIRECT_L_RAND, + NLOPT_GN_DIRECT_NOSCAL, + NLOPT_GN_DIRECT_L_NOSCAL, + NLOPT_GN_DIRECT_L_RAND_NOSCAL, + + NLOPT_GN_ORIG_DIRECT, + NLOPT_GN_ORIG_DIRECT_L, + + NLOPT_GD_STOGO, + NLOPT_GD_STOGO_RAND, + + NLOPT_LD_LBFGS_NOCEDAL, + + NLOPT_LD_LBFGS, + + NLOPT_LN_PRAXIS, + + NLOPT_LD_VAR1, + NLOPT_LD_VAR2, + + NLOPT_LD_TNEWTON, + NLOPT_LD_TNEWTON_RESTART, + NLOPT_LD_TNEWTON_PRECOND, + NLOPT_LD_TNEWTON_PRECOND_RESTART, + + NLOPT_GN_CRS2_LM, + + NLOPT_GN_MLSL, + NLOPT_GD_MLSL, + NLOPT_GN_MLSL_LDS, + NLOPT_GD_MLSL_LDS, + + NLOPT_LD_MMA, + + NLOPT_LN_COBYLA, + + NLOPT_LN_NEWUOA, + NLOPT_LN_NEWUOA_BOUND, + + NLOPT_LN_NELDERMEAD, + NLOPT_LN_SBPLX, + + NLOPT_LN_AUGLAG, + NLOPT_LD_AUGLAG, + NLOPT_LN_AUGLAG_EQ, + NLOPT_LD_AUGLAG_EQ, + + NLOPT_LN_BOBYQA, + + NLOPT_GN_ISRES, + + /* new variants that require local_optimizer to be set, + not with older constants for backwards compatibility */ + NLOPT_AUGLAG, + NLOPT_AUGLAG_EQ, + NLOPT_G_MLSL, + NLOPT_G_MLSL_LDS, + + NLOPT_LD_SLSQP, + + NLOPT_LD_CCSAQ, + + NLOPT_GN_ESCH, + + NLOPT_GN_AGS, + + NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ +} nlopt_algorithm; + +NLOPT_EXTERN(const char *) nlopt_algorithm_name(nlopt_algorithm a); + +/* nlopt_algorithm enum <-> string conversion */ +NLOPT_EXTERN(const char *) nlopt_algorithm_to_string(nlopt_algorithm algorithm); +NLOPT_EXTERN(nlopt_algorithm) nlopt_algorithm_from_string(const char *name); + +typedef enum { + NLOPT_FAILURE = -1, /* generic failure code */ + NLOPT_INVALID_ARGS = -2, + NLOPT_OUT_OF_MEMORY = -3, + NLOPT_ROUNDOFF_LIMITED = -4, + NLOPT_FORCED_STOP = -5, + NLOPT_NUM_FAILURES = -6, /* not a result, just the number of possible failures */ + NLOPT_SUCCESS = 1, /* generic success code */ + NLOPT_STOPVAL_REACHED = 2, + NLOPT_FTOL_REACHED = 3, + NLOPT_XTOL_REACHED = 4, + NLOPT_MAXEVAL_REACHED = 5, + NLOPT_MAXTIME_REACHED = 6, + NLOPT_NUM_RESULTS /* not a result, just the number of possible successes */ +} nlopt_result; + +/* nlopt_result enum <-> string conversion */ +NLOPT_EXTERN(const char *) nlopt_result_to_string(nlopt_result algorithm); +NLOPT_EXTERN(nlopt_result) nlopt_result_from_string(const char *name); + +#define NLOPT_MINF_MAX_REACHED NLOPT_STOPVAL_REACHED + +NLOPT_EXTERN(void) nlopt_srand(unsigned long seed); +NLOPT_EXTERN(void) nlopt_srand_time(void); + +NLOPT_EXTERN(void) nlopt_version(int *major, int *minor, int *bugfix); + +/*************************** OBJECT-ORIENTED API **************************/ +/* The style here is that we create an nlopt_opt "object" (an opaque pointer), + then set various optimization parameters, and then execute the + algorithm. In this way, we can add more and more optimization parameters + (including algorithm-specific ones) without breaking backwards + compatibility, having functions with zillions of parameters, or + relying non-reentrantly on global variables.*/ + +struct nlopt_opt_s; /* opaque structure, defined internally */ +typedef struct nlopt_opt_s *nlopt_opt; + +/* the only immutable parameters of an optimization are the algorithm and + the dimension n of the problem, since changing either of these could + have side-effects on lots of other parameters */ +NLOPT_EXTERN(nlopt_opt) nlopt_create(nlopt_algorithm algorithm, unsigned n); +NLOPT_EXTERN(void) nlopt_destroy(nlopt_opt opt); +NLOPT_EXTERN(nlopt_opt) nlopt_copy(const nlopt_opt opt); + +NLOPT_EXTERN(nlopt_result) nlopt_optimize(nlopt_opt opt, double *x, double *opt_f); + +NLOPT_EXTERN(nlopt_result) nlopt_set_min_objective(nlopt_opt opt, nlopt_func f, void *f_data); +NLOPT_EXTERN(nlopt_result) nlopt_set_max_objective(nlopt_opt opt, nlopt_func f, void *f_data); + +NLOPT_EXTERN(nlopt_result) nlopt_set_precond_min_objective(nlopt_opt opt, nlopt_func f, nlopt_precond pre, void *f_data); +NLOPT_EXTERN(nlopt_result) nlopt_set_precond_max_objective(nlopt_opt opt, nlopt_func f, nlopt_precond pre, void *f_data); + +NLOPT_EXTERN(nlopt_algorithm) nlopt_get_algorithm(const nlopt_opt opt); +NLOPT_EXTERN(unsigned) nlopt_get_dimension(const nlopt_opt opt); + +NLOPT_EXTERN(const char *) nlopt_get_errmsg(nlopt_opt opt); + +/* generic algorithm parameters: */ +NLOPT_EXTERN(nlopt_result) nlopt_set_param(nlopt_opt opt, const char *name, double val); +NLOPT_EXTERN(double) nlopt_get_param(const nlopt_opt opt, const char *name, double defaultval); +NLOPT_EXTERN(int) nlopt_has_param(const nlopt_opt opt, const char *name); +NLOPT_EXTERN(unsigned) nlopt_num_params(const nlopt_opt opt); +NLOPT_EXTERN(const char *) nlopt_nth_param(const nlopt_opt opt, unsigned n); + +/* constraints: */ + +NLOPT_EXTERN(nlopt_result) nlopt_set_lower_bounds(nlopt_opt opt, const double *lb); +NLOPT_EXTERN(nlopt_result) nlopt_set_lower_bounds1(nlopt_opt opt, double lb); +NLOPT_EXTERN(nlopt_result) nlopt_set_lower_bound(nlopt_opt opt, int i, double lb); +NLOPT_EXTERN(nlopt_result) nlopt_get_lower_bounds(const nlopt_opt opt, double *lb); +NLOPT_EXTERN(nlopt_result) nlopt_set_upper_bounds(nlopt_opt opt, const double *ub); +NLOPT_EXTERN(nlopt_result) nlopt_set_upper_bounds1(nlopt_opt opt, double ub); +NLOPT_EXTERN(nlopt_result) nlopt_set_upper_bound(nlopt_opt opt, int i, double ub); +NLOPT_EXTERN(nlopt_result) nlopt_get_upper_bounds(const nlopt_opt opt, double *ub); + +NLOPT_EXTERN(nlopt_result) nlopt_remove_inequality_constraints(nlopt_opt opt); +NLOPT_EXTERN(nlopt_result) nlopt_add_inequality_constraint(nlopt_opt opt, nlopt_func fc, void *fc_data, double tol); +NLOPT_EXTERN(nlopt_result) nlopt_add_precond_inequality_constraint(nlopt_opt opt, nlopt_func fc, nlopt_precond pre, void *fc_data, double tol); +NLOPT_EXTERN(nlopt_result) nlopt_add_inequality_mconstraint(nlopt_opt opt, unsigned m, nlopt_mfunc fc, void *fc_data, const double *tol); + +NLOPT_EXTERN(nlopt_result) nlopt_remove_equality_constraints(nlopt_opt opt); +NLOPT_EXTERN(nlopt_result) nlopt_add_equality_constraint(nlopt_opt opt, nlopt_func h, void *h_data, double tol); +NLOPT_EXTERN(nlopt_result) nlopt_add_precond_equality_constraint(nlopt_opt opt, nlopt_func h, nlopt_precond pre, void *h_data, double tol); +NLOPT_EXTERN(nlopt_result) nlopt_add_equality_mconstraint(nlopt_opt opt, unsigned m, nlopt_mfunc h, void *h_data, const double *tol); + +/* stopping criteria: */ + +NLOPT_EXTERN(nlopt_result) nlopt_set_stopval(nlopt_opt opt, double stopval); +NLOPT_EXTERN(double) nlopt_get_stopval(const nlopt_opt opt); + +NLOPT_EXTERN(nlopt_result) nlopt_set_ftol_rel(nlopt_opt opt, double tol); +NLOPT_EXTERN(double) nlopt_get_ftol_rel(const nlopt_opt opt); +NLOPT_EXTERN(nlopt_result) nlopt_set_ftol_abs(nlopt_opt opt, double tol); +NLOPT_EXTERN(double) nlopt_get_ftol_abs(const nlopt_opt opt); + +NLOPT_EXTERN(nlopt_result) nlopt_set_xtol_rel(nlopt_opt opt, double tol); +NLOPT_EXTERN(double) nlopt_get_xtol_rel(const nlopt_opt opt); +NLOPT_EXTERN(nlopt_result) nlopt_set_xtol_abs1(nlopt_opt opt, double tol); +NLOPT_EXTERN(nlopt_result) nlopt_set_xtol_abs(nlopt_opt opt, const double *tol); +NLOPT_EXTERN(nlopt_result) nlopt_get_xtol_abs(const nlopt_opt opt, double *tol); +NLOPT_EXTERN(nlopt_result) nlopt_set_x_weights1(nlopt_opt opt, double w); +NLOPT_EXTERN(nlopt_result) nlopt_set_x_weights(nlopt_opt opt, const double *w); +NLOPT_EXTERN(nlopt_result) nlopt_get_x_weights(const nlopt_opt opt, double *w); + +NLOPT_EXTERN(nlopt_result) nlopt_set_maxeval(nlopt_opt opt, int maxeval); +NLOPT_EXTERN(int) nlopt_get_maxeval(const nlopt_opt opt); + +NLOPT_EXTERN(int) nlopt_get_numevals(const nlopt_opt opt); + +NLOPT_EXTERN(nlopt_result) nlopt_set_maxtime(nlopt_opt opt, double maxtime); +NLOPT_EXTERN(double) nlopt_get_maxtime(const nlopt_opt opt); + +NLOPT_EXTERN(nlopt_result) nlopt_force_stop(nlopt_opt opt); +NLOPT_EXTERN(nlopt_result) nlopt_set_force_stop(nlopt_opt opt, int val); +NLOPT_EXTERN(int) nlopt_get_force_stop(const nlopt_opt opt); + +/* more algorithm-specific parameters */ + +NLOPT_EXTERN(nlopt_result) nlopt_set_local_optimizer(nlopt_opt opt, const nlopt_opt local_opt); + +NLOPT_EXTERN(nlopt_result) nlopt_set_population(nlopt_opt opt, unsigned pop); +NLOPT_EXTERN(unsigned) nlopt_get_population(const nlopt_opt opt); + +NLOPT_EXTERN(nlopt_result) nlopt_set_vector_storage(nlopt_opt opt, unsigned dim); +NLOPT_EXTERN(unsigned) nlopt_get_vector_storage(const nlopt_opt opt); + +NLOPT_EXTERN(nlopt_result) nlopt_set_default_initial_step(nlopt_opt opt, const double *x); +NLOPT_EXTERN(nlopt_result) nlopt_set_initial_step(nlopt_opt opt, const double *dx); +NLOPT_EXTERN(nlopt_result) nlopt_set_initial_step1(nlopt_opt opt, double dx); +NLOPT_EXTERN(nlopt_result) nlopt_get_initial_step(const nlopt_opt opt, const double *x, double *dx); + +/* the following are functions mainly designed to be used internally + by the Fortran and SWIG wrappers, allow us to tel nlopt_destroy and + nlopt_copy to do something to the f_data pointers (e.g. free or + duplicate them, respectively) */ +typedef void *(*nlopt_munge) (void *p); +NLOPT_EXTERN(void) nlopt_set_munge(nlopt_opt opt, nlopt_munge munge_on_destroy, nlopt_munge munge_on_copy); +typedef void *(*nlopt_munge2) (void *p, void *data); +NLOPT_EXTERN(void) nlopt_munge_data(nlopt_opt opt, nlopt_munge2 munge, void *data); + +/*************************** DEPRECATED API **************************/ +/* The new "object-oriented" API is preferred, since it allows us to + gracefully add new features and algorithm-specific options in a + re-entrant way, and we can automatically assume reasonable defaults + for unspecified parameters. */ + +/* Where possible (e.g. for gcc >= 3.1), enable a compiler warning + for code that uses a deprecated function */ +#if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__==3 && __GNUC_MINOR__ > 0)) +# define NLOPT_DEPRECATED __attribute__((deprecated)) +#else +# define NLOPT_DEPRECATED +#endif + +typedef double (*nlopt_func_old) (int n, const double *x, double *gradient, /* NULL if not needed */ + void *func_data); + +NLOPT_EXTERN(nlopt_result) nlopt_minimize(nlopt_algorithm algorithm, int n, nlopt_func_old f, void *f_data, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, /* out: minimum */ + double minf_max, double ftol_rel, double ftol_abs, double xtol_rel, const double *xtol_abs, int maxeval, double maxtime) NLOPT_DEPRECATED; + +NLOPT_EXTERN(nlopt_result) nlopt_minimize_constrained(nlopt_algorithm algorithm, int n, nlopt_func_old f, void *f_data, int m, nlopt_func_old fc, void *fc_data, ptrdiff_t fc_datum_size, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, /* out: minimum */ + double minf_max, double ftol_rel, double ftol_abs, double xtol_rel, const double *xtol_abs, int maxeval, double maxtime) NLOPT_DEPRECATED; + +NLOPT_EXTERN(nlopt_result) nlopt_minimize_econstrained(nlopt_algorithm algorithm, int n, nlopt_func_old f, void *f_data, int m, nlopt_func_old fc, void *fc_data, ptrdiff_t fc_datum_size, int p, nlopt_func_old h, void *h_data, ptrdiff_t h_datum_size, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, /* out: minimum */ + double minf_max, double ftol_rel, double ftol_abs, + double xtol_rel, const double *xtol_abs, double htol_rel, double htol_abs, int maxeval, double maxtime) NLOPT_DEPRECATED; + +NLOPT_EXTERN(void) nlopt_get_local_search_algorithm(nlopt_algorithm * deriv, nlopt_algorithm * nonderiv, int *maxeval) NLOPT_DEPRECATED; +NLOPT_EXTERN(void) nlopt_set_local_search_algorithm(nlopt_algorithm deriv, nlopt_algorithm nonderiv, int maxeval) NLOPT_DEPRECATED; + +NLOPT_EXTERN(int) nlopt_get_stochastic_population(void) NLOPT_DEPRECATED; +NLOPT_EXTERN(void) nlopt_set_stochastic_population(int pop) NLOPT_DEPRECATED; + +/*********************************************************************/ + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ +#endif diff --git a/slsqp/nlopt_config.h b/slsqp/nlopt_config.h new file mode 100644 index 0000000..40dd1f2 --- /dev/null +++ b/slsqp/nlopt_config.h @@ -0,0 +1,139 @@ +/*============================================================================== +# NLOPT CMake configuration file +# +# NLopt is a free/open-source library for nonlinear optimization, providing +# a common interface for a number of different free optimization routines +# available online as well as original implementations of various other +# algorithms +# WEBSITE: http://ab-initio.mit.edu/wiki/index.php/NLopt +# AUTHOR: Steven G. Johnson +# +# This config.cmake.h.in file was created to compile NLOPT with the CMAKE utility. +# Benoit Scherrer, 2010 CRL, Harvard Medical School +# Copyright (c) 2008-2009 Children's Hospital Boston +# +# Minor changes to the source was applied to make possible the compilation with +# Cmake under Linux/Win32 +#============================================================================*/ + +/* Bugfix version number. */ +#define BUGFIX_VERSION + +/* Define to enable extra debugging code. */ +#undef DEBUG + +/* Define to 1 if you have the `BSDgettimeofday' function. */ +#undef HAVE_BSDGETTIMEOFDAY + +/* Define if the copysign function/macro is available. */ +#define HAVE_COPYSIGN + +/* Define if the fpclassify() function/macro is available. */ +/* #undef HAVE_FPCLASSIFY */ + +/* Define to 1 if you have the header file. */ +#define HAVE_GETOPT_H + +/* Define to 1 if you have the getopt function in your standard library. */ +#define HAVE_GETOPT + +/* Define to 1 if you have the `getpid' function. */ +#define HAVE_GETPID + +/* Define if syscall(SYS_gettid) available. */ +#undef HAVE_GETTID_SYSCALL + +/* Define to 1 if you have the `gettimeofday' function. */ +#define HAVE_GETTIMEOFDAY + +/* Define if the isinf() function/macro is available. */ +#define HAVE_ISINF + +/* Define if the isnan() function/macro is available. */ +#define HAVE_ISNAN + +/* Define to 1 if you have the `m' library (-lm). */ +#undef HAVE_LIBM + +/* Define to 1 if you have the `qsort_r' function. */ +#define HAVE_QSORT_R + +/* Define to 1 if you have the header file. */ +#define HAVE_STDINT_H + +/* Define to 1 if you have the header file. */ +#define HAVE_SYS_TIME_H + +/* Define to 1 if you have the `time' function. */ +#define HAVE_TIME + +/* Define to 1 if the system has the type `uint32_t'. */ +#define HAVE_UINT32_T + +/* Define to 1 if you have the header file. */ +#define HAVE_UNISTD_H + +/* Define to the sub-directory in which libtool stores uninstalled libraries. + */ +#undef LT_OBJDIR + +/* Major version number. */ +#define MAJOR_VERSION + +/* Minor version number. */ +#define MINOR_VERSION + +/* Name of package */ +#undef PACKAGE + +/* Define to the address where bug reports for this package should be sent. */ +#undef PACKAGE_BUGREPORT + +/* Define to the full name of this package. */ +#undef PACKAGE_NAME + +/* Define to the full name and version of this package. */ +#undef PACKAGE_STRING + +/* Define to the one symbol short name of this package. */ +#undef PACKAGE_TARNAME + +/* Define to the home page for this package. */ +#undef PACKAGE_URL + +/* Define to the version of this package. */ +#undef PACKAGE_VERSION + +/* replacement for broken HUGE_VAL macro, if needed */ +#undef REPLACEMENT_HUGE_VAL + +/* The size of `unsigned int', as computed by sizeof. */ +#define SIZEOF_UNSIGNED_INT 4 + +/* The size of `unsigned long', as computed by sizeof. */ +#define SIZEOF_UNSIGNED_LONG 8 + +/* Define to 1 if you have the ANSI C header files. */ +#undef STDC_HEADERS + +/* Define to C thread-local keyword, or to nothing if this is not supported in + your compiler. */ +#define THREADLOCAL + +/* Define to 1 if you can safely include both and . */ +#define TIME_WITH_SYS_TIME + +/* Version number of package */ +#undef VERSION + +/* Define if compiled including C++-based routines */ +/* #undef NLOPT_CXX */ + +/* Define to empty if `const' does not conform to ANSI C. */ +#undef const + +/* Define to `__inline__' or `__inline' if that's what the C compiler + calls it, or to nothing if 'inline' is not supported under any name. */ +#ifndef __cplusplus +#undef inline +#endif diff --git a/slsqp/nlopt_config.h.in b/slsqp/nlopt_config.h.in new file mode 100644 index 0000000..25daf4e --- /dev/null +++ b/slsqp/nlopt_config.h.in @@ -0,0 +1,139 @@ +/*============================================================================== +# NLOPT CMake configuration file +# +# NLopt is a free/open-source library for nonlinear optimization, providing +# a common interface for a number of different free optimization routines +# available online as well as original implementations of various other +# algorithms +# WEBSITE: http://ab-initio.mit.edu/wiki/index.php/NLopt +# AUTHOR: Steven G. Johnson +# +# This config.cmake.h.in file was created to compile NLOPT with the CMAKE utility. +# Benoit Scherrer, 2010 CRL, Harvard Medical School +# Copyright (c) 2008-2009 Children's Hospital Boston +# +# Minor changes to the source was applied to make possible the compilation with +# Cmake under Linux/Win32 +#============================================================================*/ + +/* Bugfix version number. */ +#define BUGFIX_VERSION @NLOPT_BUGFIX_VERSION@ + +/* Define to enable extra debugging code. */ +#undef DEBUG + +/* Define to 1 if you have the `BSDgettimeofday' function. */ +#undef HAVE_BSDGETTIMEOFDAY + +/* Define if the copysign function/macro is available. */ +#cmakedefine HAVE_COPYSIGN + +/* Define if the fpclassify() function/macro is available. */ +#cmakedefine HAVE_FPCLASSIFY + +/* Define to 1 if you have the header file. */ +#cmakedefine HAVE_GETOPT_H + +/* Define to 1 if you have the getopt function in your standard library. */ +#cmakedefine HAVE_GETOPT + +/* Define to 1 if you have the `getpid' function. */ +#cmakedefine HAVE_GETPID + +/* Define if syscall(SYS_gettid) available. */ +#undef HAVE_GETTID_SYSCALL + +/* Define to 1 if you have the `gettimeofday' function. */ +#cmakedefine HAVE_GETTIMEOFDAY + +/* Define if the isinf() function/macro is available. */ +#cmakedefine HAVE_ISINF + +/* Define if the isnan() function/macro is available. */ +#cmakedefine HAVE_ISNAN + +/* Define to 1 if you have the `m' library (-lm). */ +#undef HAVE_LIBM + +/* Define to 1 if you have the `qsort_r' function. */ +#cmakedefine HAVE_QSORT_R + +/* Define to 1 if you have the header file. */ +#cmakedefine HAVE_STDINT_H + +/* Define to 1 if you have the header file. */ +#cmakedefine HAVE_SYS_TIME_H + +/* Define to 1 if you have the `time' function. */ +#cmakedefine HAVE_TIME + +/* Define to 1 if the system has the type `uint32_t'. */ +#cmakedefine HAVE_UINT32_T + +/* Define to 1 if you have the header file. */ +#cmakedefine HAVE_UNISTD_H + +/* Define to the sub-directory in which libtool stores uninstalled libraries. + */ +#undef LT_OBJDIR + +/* Major version number. */ +#define MAJOR_VERSION @NLOPT_MAJOR_VERSION@ + +/* Minor version number. */ +#define MINOR_VERSION @NLOPT_MINOR_VERSION@ + +/* Name of package */ +#undef PACKAGE + +/* Define to the address where bug reports for this package should be sent. */ +#undef PACKAGE_BUGREPORT + +/* Define to the full name of this package. */ +#undef PACKAGE_NAME + +/* Define to the full name and version of this package. */ +#undef PACKAGE_STRING + +/* Define to the one symbol short name of this package. */ +#undef PACKAGE_TARNAME + +/* Define to the home page for this package. */ +#undef PACKAGE_URL + +/* Define to the version of this package. */ +#undef PACKAGE_VERSION + +/* replacement for broken HUGE_VAL macro, if needed */ +#undef REPLACEMENT_HUGE_VAL + +/* The size of `unsigned int', as computed by sizeof. */ +#define SIZEOF_UNSIGNED_INT @SIZEOF_UNSIGNED_INT@ + +/* The size of `unsigned long', as computed by sizeof. */ +#define SIZEOF_UNSIGNED_LONG @SIZEOF_UNSIGNED_LONG@ + +/* Define to 1 if you have the ANSI C header files. */ +#undef STDC_HEADERS + +/* Define to C thread-local keyword, or to nothing if this is not supported in + your compiler. */ +#define THREADLOCAL @THREADLOCAL@ + +/* Define to 1 if you can safely include both and . */ +#cmakedefine TIME_WITH_SYS_TIME + +/* Version number of package */ +#undef VERSION + +/* Define if compiled including C++-based routines */ +#cmakedefine NLOPT_CXX + +/* Define to empty if `const' does not conform to ANSI C. */ +#undef const + +/* Define to `__inline__' or `__inline' if that's what the C compiler + calls it, or to nothing if 'inline' is not supported under any name. */ +#ifndef __cplusplus +#undef inline +#endif diff --git a/slsqp/rescale.c b/slsqp/rescale.c new file mode 100644 index 0000000..d2f82b9 --- /dev/null +++ b/slsqp/rescale.c @@ -0,0 +1,96 @@ +/* Copyright (c) 2007-2014 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include "nlopt-util.h" + +/* Return a new array of length n (> 0) that gives a rescaling factor + for each dimension, or NULL if out of memory, with dx being the + array of nonzero initial steps in each dimension. */ +double *nlopt_compute_rescaling(unsigned n, const double *dx) +{ + double *s = (double *) malloc(sizeof(double) * n); + unsigned i; + + if (!s) + return NULL; + for (i = 0; i < n; ++i) + s[i] = 1.0; /* default: no rescaling */ + if (n == 1) + return s; + + for (i = 1; i < n && dx[i] == dx[i - 1]; ++i); + if (i < n) { /* unequal initial steps, rescale to make equal to dx[0] */ + for (i = 1; i < n; ++i) + s[i] = dx[i] / dx[0]; + } + return s; +} + +void nlopt_rescale(unsigned n, const double *s, const double *x, double *xs) +{ + unsigned i; + if (!s) { + for (i = 0; i < n; ++i) + xs[i] = x[i]; + } else { + for (i = 0; i < n; ++i) + xs[i] = x[i] / s[i]; + } +} + +void nlopt_unscale(unsigned n, const double *s, const double *x, double *xs) +{ + unsigned i; + if (!s) { + for (i = 0; i < n; ++i) + xs[i] = x[i]; + } else { + for (i = 0; i < n; ++i) + xs[i] = x[i] * s[i]; + } +} + +/* return a new array of length n equal to the original array + x divided by the scale factors s, or NULL on a memory error */ +double *nlopt_new_rescaled(unsigned n, const double *s, const double *x) +{ + double *xs = (double *) malloc(sizeof(double) * n); + if (!xs) + return NULL; + nlopt_rescale(n, s, x, xs); + return xs; +} + +/* since rescaling can flip the signs of the x components and the bounds, + we may have to re-order the bounds in order to ensure that they + remain in the correct order */ +void nlopt_reorder_bounds(unsigned n, double *lb, double *ub) +{ + unsigned i; + for (i = 0; i < n; ++i) + if (lb[i] > ub[i]) { + double t = lb[i]; + lb[i] = ub[i]; + ub[i] = t; + } +} diff --git a/slsqp/slsqp.c b/slsqp/slsqp.c new file mode 100644 index 0000000..ee8927a --- /dev/null +++ b/slsqp/slsqp.c @@ -0,0 +1,2638 @@ +/* SLSQP: Sequentional Least Squares Programming (aka sequential quadratic programming SQP) + method for nonlinearly constrained nonlinear optimization, by Dieter Kraft (1991). + Fortran released under a free (BSD) license by ACM to the SciPy project and used there. + C translation via f2c + hand-cleanup and incorporation into NLopt by S. G. Johnson (2009). */ + +/* Table of constant values */ + +#include +#include +#include + +#include "slsqp.h" + +/* ALGORITHM 733, COLLECTED ALGORITHMS FROM ACM. */ +/* TRANSACTIONS ON MATHEMATICAL SOFTWARE, */ +/* VOL. 20, NO. 3, SEPTEMBER, 1994, PP. 262-281. */ +/* http://doi.acm.org/10.1145/192115.192124 */ + + +/* http://permalink.gmane.org/gmane.comp.python.scientific.devel/6725 */ +/* ------ */ +/* From: Deborah Cotton */ +/* Date: Fri, 14 Sep 2007 12:35:55 -0500 */ +/* Subject: RE: Algorithm License requested */ +/* To: Alan Isaac */ + +/* Prof. Issac, */ + +/* In that case, then because the author consents to [the ACM] releasing */ +/* the code currently archived at http://www.netlib.org/toms/733 under the */ +/* BSD license, the ACM hereby releases this code under the BSD license. */ + +/* Regards, */ + +/* Deborah Cotton, Copyright & Permissions */ +/* ACM Publications */ +/* 2 Penn Plaza, Suite 701** */ +/* New York, NY 10121-0701 */ +/* permissions@acm.org */ +/* 212.869.7440 ext. 652 */ +/* Fax. 212.869.0481 */ +/* ------ */ + +/********************************* BLAS1 routines *************************/ + +/* COPIES A VECTOR, X, TO A VECTOR, Y, with the given increments */ +static void dcopy___(int *n_, const double *dx, int incx, + double *dy, int incy) +{ + int i, n = *n_; + + if (n <= 0) return; + if (incx == 1 && incy == 1) + memcpy(dy, dx, sizeof(double) * ((unsigned) n)); + else if (incx == 0 && incy == 1) { + double x = dx[0]; + for (i = 0; i < n; ++i) dy[i] = x; + } + else { + for (i = 0; i < n; ++i) dy[i*incy] = dx[i*incx]; + } +} /* dcopy___ */ + +/* CONSTANT TIMES A VECTOR PLUS A VECTOR. */ +static void daxpy_sl__(int *n_, const double *da_, const double *dx, + int incx, double *dy, int incy) +{ + int n = *n_, i; + double da = *da_; + + if (n <= 0 || da == 0) return; + for (i = 0; i < n; ++i) dy[i*incy] += da * dx[i*incx]; +} + +/* dot product dx dot dy. */ +static double ddot_sl__(int *n_, double *dx, int incx, double *dy, int incy) +{ + int n = *n_, i; + double sum = 0; + if (n <= 0) return 0; + for (i = 0; i < n; ++i) sum += dx[i*incx] * dy[i*incy]; + return (double) sum; +} + +/* compute the L2 norm of array DX of length N, stride INCX */ +static double dnrm2___(int *n_, double *dx, int incx) +{ + int i, n = *n_; + double xmax = 0, scale; + double sum = 0; + for (i = 0; i < n; ++i) { + double xabs = fabs(dx[incx*i]); + if (xmax < xabs) xmax = xabs; + } + if (xmax == 0) return 0; + scale = 1.0 / xmax; + for (i = 0; i < n; ++i) { + double xs = scale * dx[incx*i]; + sum += xs * xs; + } + return xmax * sqrt((double) sum); +} + +/* apply Givens rotation */ +static void dsrot_(int n, double *dx, int incx, + double *dy, int incy, double *c__, double *s_) +{ + int i; + double c = *c__, s = *s_; + + for (i = 0; i < n; ++i) { + double x = dx[incx*i], y = dy[incy*i]; + dx[incx*i] = c * x + s * y; + dy[incy*i] = c * y - s * x; + } +} + +/* construct Givens rotation */ +static void dsrotg_(double *da, double *db, double *c, double *s) +{ + double absa, absb, roe, scale; + + absa = fabs(*da); absb = fabs(*db); + if (absa > absb) { + roe = *da; + scale = absa; + } + else { + roe = *db; + scale = absb; + } + + if (scale != 0) { + double r, iscale = 1 / scale; + double tmpa = (*da) * iscale, tmpb = (*db) * iscale; + r = (roe < 0 ? -scale : scale) * sqrt((tmpa * tmpa) + (tmpb * tmpb)); + *c = *da / r; *s = *db / r; + *da = r; + if (*c != 0 && fabs(*c) <= *s) *db = 1 / *c; + else *db = *s; + } + else { + *c = 1; + *s = *da = *db = 0; + } +} + +/* scales vector X(n) by constant da */ +static void dscal_sl__(int *n_, const double *da, double *dx, int incx) +{ + int i, n = *n_; + double alpha = *da; + for (i = 0; i < n; ++i) dx[i*incx] *= alpha; +} + +/**************************************************************************/ + +static const int c__0 = 0; +static const int c__1 = 1; +static const int c__2 = 2; + +#define MIN2(a,b) ((a) <= (b) ? (a) : (b)) +#define MAX2(a,b) ((a) >= (b) ? (a) : (b)) + +static void h12_(const int *mode, int *lpivot, int *l1, + int *m, double *u, const int *iue, double *up, + double *c__, const int *ice, const int *icv, const int *ncv) +{ + /* Initialized data */ + + const double one = 1.; + + /* System generated locals */ + int u_dim1, u_offset, i__1, i__2; + double d__1; + + /* Local variables */ + double b; + int i__, j, i2, i3, i4; + double cl, sm; + int incr; + double clinv; + +/* C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 */ +/* TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 */ +/* CONSTRUCTION AND/OR APPLICATION OF A SINGLE */ +/* HOUSEHOLDER TRANSFORMATION Q = I + U*(U**T)/B */ +/* MODE = 1 OR 2 TO SELECT ALGORITHM H1 OR H2 . */ +/* LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */ +/* L1,M IF L1 <= M THE TRANSFORMATION WILL BE CONSTRUCTED TO */ +/* ZERO ELEMENTS INDEXED FROM L1 THROUGH M. */ +/* IF L1 > M THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */ +/* U(),IUE,UP */ +/* ON ENTRY TO H1 U() STORES THE PIVOT VECTOR. */ +/* IUE IS THE STORAGE INCREMENT BETWEEN ELEMENTS. */ +/* ON EXIT FROM H1 U() AND UP STORE QUANTITIES DEFINING */ +/* THE VECTOR U OF THE HOUSEHOLDER TRANSFORMATION. */ +/* ON ENTRY TO H2 U() AND UP */ +/* SHOULD STORE QUANTITIES PREVIOUSLY COMPUTED BY H1. */ +/* THESE WILL NOT BE MODIFIED BY H2. */ +/* C() ON ENTRY TO H1 OR H2 C() STORES A MATRIX WHICH WILL BE */ +/* REGARDED AS A SET OF VECTORS TO WHICH THE HOUSEHOLDER */ +/* TRANSFORMATION IS TO BE APPLIED. */ +/* ON EXIT C() STORES THE SET OF TRANSFORMED VECTORS. */ +/* ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */ +/* ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). */ +/* NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. */ +/* IF NCV <= 0 NO OPERATIONS WILL BE DONE ON C(). */ + /* Parameter adjustments */ + u_dim1 = *iue; + u_offset = 1 + u_dim1; + u -= u_offset; + --c__; + + /* Function Body */ + if (0 >= *lpivot || *lpivot >= *l1 || *l1 > *m) { + goto L80; + } + cl = (d__1 = u[*lpivot * u_dim1 + 1], fabs(d__1)); + if (*mode == 2) { + goto L30; + } +/* ****** CONSTRUCT THE TRANSFORMATION ****** */ + i__1 = *m; + for (j = *l1; j <= i__1; ++j) { + sm = (d__1 = u[j * u_dim1 + 1], fabs(d__1)); +/* L10: */ + cl = MAX2(sm,cl); + } + if (cl <= 0.0) { + goto L80; + } + clinv = one / cl; +/* Computing 2nd power */ + d__1 = u[*lpivot * u_dim1 + 1] * clinv; + sm = d__1 * d__1; + i__1 = *m; + for (j = *l1; j <= i__1; ++j) { +/* L20: */ +/* Computing 2nd power */ + d__1 = u[j * u_dim1 + 1] * clinv; + sm += d__1 * d__1; + } + cl *= sqrt(sm); + if (u[*lpivot * u_dim1 + 1] > 0.0) { + cl = -cl; + } + *up = u[*lpivot * u_dim1 + 1] - cl; + u[*lpivot * u_dim1 + 1] = cl; + goto L40; +/* ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C ****** */ +L30: + if (cl <= 0.0) { + goto L80; + } +L40: + if (*ncv <= 0) { + goto L80; + } + b = *up * u[*lpivot * u_dim1 + 1]; + if (b >= 0.0) { + goto L80; + } + b = one / b; + i2 = 1 - *icv + *ice * (*lpivot - 1); + incr = *ice * (*l1 - *lpivot); + i__1 = *ncv; + for (j = 1; j <= i__1; ++j) { + i2 += *icv; + i3 = i2 + incr; + i4 = i3; + sm = c__[i2] * *up; + i__2 = *m; + for (i__ = *l1; i__ <= i__2; ++i__) { + sm += c__[i3] * u[i__ * u_dim1 + 1]; +/* L50: */ + i3 += *ice; + } + if (sm == 0.0) { + goto L70; + } + sm *= b; + c__[i2] += sm * *up; + i__2 = *m; + for (i__ = *l1; i__ <= i__2; ++i__) { + c__[i4] += sm * u[i__ * u_dim1 + 1]; +/* L60: */ + i4 += *ice; + } +L70: + ; + } +L80: + return; +} /* h12_ */ + +static void nnls_(double *a, int *mda, int *m, int * + n, double *b, double *x, double *rnorm, double *w, + double *z__, int *indx, int *mode) +{ + /* Initialized data */ + + const double one = 1.; + const double factor = .01; + + /* System generated locals */ + int a_dim1, a_offset, i__1, i__2; + double d__1; + + /* Local variables */ + double c__; + int i__, j, k, l; + double s, t; + int ii, jj, ip, iz, jz; + double up; + int iz1, iz2, npp1, iter; + double wmax, alpha, asave; + int itmax, izmax = 0, nsetp; + double unorm; + +/* C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY: */ +/* 'SOLVING LEAST SQUARES PROBLEMS'. PRENTICE-HALL.1974 */ +/* ********** NONNEGATIVE LEAST SQUARES ********** */ +/* GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN */ +/* N-VECTOR, X, WHICH SOLVES THE LEAST SQUARES PROBLEM */ +/* A*X = B SUBJECT TO X >= 0 */ +/* A(),MDA,M,N */ +/* MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE ARRAY,A(). */ +/* ON ENTRY A() CONTAINS THE M BY N MATRIX,A. */ +/* ON EXIT A() CONTAINS THE PRODUCT Q*A, */ +/* WHERE Q IS AN M BY M ORTHOGONAL MATRIX GENERATED */ +/* IMPLICITLY BY THIS SUBROUTINE. */ +/* EITHER M>=N OR M iz2 || nsetp >= *m) { + goto L280; + } + i__1 = iz2; + for (iz = iz1; iz <= i__1; ++iz) { + j = indx[iz]; +/* L120: */ + i__2 = *m - nsetp; + w[j] = ddot_sl__(&i__2, &a[npp1 + j * a_dim1], 1, &b[npp1], 1) + ; + } +/* STEP THREE (TEST DUAL VARIABLES) */ +L130: + wmax = 0.0; + i__2 = iz2; + for (iz = iz1; iz <= i__2; ++iz) { + j = indx[iz]; + if (w[j] <= wmax) { + goto L140; + } + wmax = w[j]; + izmax = iz; +L140: + ; + } +/* .....EXIT LOOP A */ + if (wmax <= 0.0) { + goto L280; + } + iz = izmax; + j = indx[iz]; +/* STEP FOUR (TEST INDX J FOR LINEAR DEPENDENCY) */ + asave = a[npp1 + j * a_dim1]; + i__2 = npp1 + 1; + h12_(&c__1, &npp1, &i__2, m, &a[j * a_dim1 + 1], &c__1, &up, &z__[1], & + c__1, &c__1, &c__0); + unorm = dnrm2___(&nsetp, &a[j * a_dim1 + 1], 1); + t = factor * (d__1 = a[npp1 + j * a_dim1], fabs(d__1)); + d__1 = unorm + t; + if (d__1 - unorm <= 0.0) { + goto L150; + } + dcopy___(m, &b[1], 1, &z__[1], 1); + i__2 = npp1 + 1; + h12_(&c__2, &npp1, &i__2, m, &a[j * a_dim1 + 1], &c__1, &up, &z__[1], & + c__1, &c__1, &c__1); + if (z__[npp1] / a[npp1 + j * a_dim1] > 0.0) { + goto L160; + } +L150: + a[npp1 + j * a_dim1] = asave; + w[j] = 0.0; + goto L130; +/* STEP FIVE (ADD COLUMN) */ +L160: + dcopy___(m, &z__[1], 1, &b[1], 1); + indx[iz] = indx[iz1]; + indx[iz1] = j; + ++iz1; + nsetp = npp1; + ++npp1; + i__2 = iz2; + for (jz = iz1; jz <= i__2; ++jz) { + jj = indx[jz]; +/* L170: */ + h12_(&c__2, &nsetp, &npp1, m, &a[j * a_dim1 + 1], &c__1, &up, &a[jj * + a_dim1 + 1], &c__1, mda, &c__1); + } + k = MIN2(npp1,*mda); + w[j] = 0.0; + i__2 = *m - nsetp; + dcopy___(&i__2, &w[j], 0, &a[k + j * a_dim1], 1); +/* STEP SIX (SOLVE LEAST SQUARES SUB-PROBLEM) */ +/* .....ENTRY LOOP B */ +L180: + for (ip = nsetp; ip >= 1; --ip) { + if (ip == nsetp) { + goto L190; + } + d__1 = -z__[ip + 1]; + daxpy_sl__(&ip, &d__1, &a[jj * a_dim1 + 1], 1, &z__[1], 1); +L190: + jj = indx[ip]; +/* L200: */ + z__[ip] /= a[ip + jj * a_dim1]; + } + ++iter; + if (iter <= itmax) { + goto L220; + } +L210: + *mode = 3; + goto L280; +/* STEP SEVEN TO TEN (STEP LENGTH ALGORITHM) */ +L220: + alpha = one; + jj = 0; + i__2 = nsetp; + for (ip = 1; ip <= i__2; ++ip) { + if (z__[ip] > 0.0) { + goto L230; + } + l = indx[ip]; + t = -x[l] / (z__[ip] - x[l]); + if (alpha < t) { + goto L230; + } + alpha = t; + jj = ip; +L230: + ; + } + i__2 = nsetp; + for (ip = 1; ip <= i__2; ++ip) { + l = indx[ip]; +/* L240: */ + x[l] = (one - alpha) * x[l] + alpha * z__[ip]; + } +/* .....EXIT LOOP B */ + if (jj == 0) { + goto L110; + } +/* STEP ELEVEN (DELETE COLUMN) */ + i__ = indx[jj]; +L250: + x[i__] = 0.0; + ++jj; + i__2 = nsetp; + for (j = jj; j <= i__2; ++j) { + ii = indx[j]; + indx[j - 1] = ii; + dsrotg_(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], &c__, &s); + t = a[j - 1 + ii * a_dim1]; + dsrot_(*n, &a[j - 1 + a_dim1], *mda, &a[j + a_dim1], *mda, &c__, &s); + a[j - 1 + ii * a_dim1] = t; + a[j + ii * a_dim1] = 0.0; +/* L260: */ + dsrot_(1, &b[j - 1], 1, &b[j], 1, &c__, &s); + } + npp1 = nsetp; + --nsetp; + --iz1; + indx[iz1] = i__; + if (nsetp <= 0) { + goto L210; + } + i__2 = nsetp; + for (jj = 1; jj <= i__2; ++jj) { + i__ = indx[jj]; + if (x[i__] <= 0.0) { + goto L250; + } +/* L270: */ + } + dcopy___(m, &b[1], 1, &z__[1], 1); + goto L180; +/* STEP TWELVE (SOLUTION) */ +L280: + k = MIN2(npp1,*m); + i__2 = *m - nsetp; + *rnorm = dnrm2___(&i__2, &b[k], 1); + if (npp1 > *m) { + w[1] = 0.0; + dcopy___(n, &w[1], 0, &w[1], 1); + } +/* END OF SUBROUTINE NNLS */ +L290: + return; +} /* nnls_ */ + +static void ldp_(double *g, int *mg, int *m, int *n, + double *h__, double *x, double *xnorm, double *w, + int *indx, int *mode) +{ + /* Initialized data */ + + const double one = 1.; + + /* System generated locals */ + int g_dim1, g_offset, i__1, i__2; + double d__1; + + /* Local variables */ + int i__, j, n1, if__, iw, iy, iz; + double fac; + double rnorm; + int iwdual; + +/* T */ +/* MINIMIZE 1/2 X X SUBJECT TO G * X >= H. */ +/* C.L. LAWSON, R.J. HANSON: 'SOLVING LEAST SQUARES PROBLEMS' */ +/* PRENTICE HALL, ENGLEWOOD CLIFFS, NEW JERSEY, 1974. */ +/* PARAMETER DESCRIPTION: */ +/* G(),MG,M,N ON ENTRY G() STORES THE M BY N MATRIX OF */ +/* LINEAR INEQUALITY CONSTRAINTS. G() HAS FIRST */ +/* DIMENSIONING PARAMETER MG */ +/* H() ON ENTRY H() STORES THE M VECTOR H REPRESENTING */ +/* THE RIGHT SIDE OF THE INEQUALITY SYSTEM */ +/* REMARK: G(),H() WILL NOT BE CHANGED DURING CALCULATIONS BY LDP */ +/* X() ON ENTRY X() NEED NOT BE INITIALIZED. */ +/* ON EXIT X() STORES THE SOLUTION VECTOR X IF MODE=1. */ +/* XNORM ON EXIT XNORM STORES THE EUCLIDIAN NORM OF THE */ +/* SOLUTION VECTOR IF COMPUTATION IS SUCCESSFUL */ +/* W() W IS A ONE DIMENSIONAL WORKING SPACE, THE LENGTH */ +/* OF WHICH SHOULD BE AT LEAST (M+2)*(N+1) + 2*M */ +/* ON EXIT W() STORES THE LAGRANGE MULTIPLIERS */ +/* ASSOCIATED WITH THE CONSTRAINTS */ +/* AT THE SOLUTION OF PROBLEM LDP */ +/* INDX() INDX() IS A ONE DIMENSIONAL INT WORKING SPACE */ +/* OF LENGTH AT LEAST M */ +/* MODE MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */ +/* MEANINGS: */ +/* MODE=1: SUCCESSFUL COMPUTATION */ +/* 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N.LE.0) */ +/* 3: ITERATION COUNT EXCEEDED BY NNLS */ +/* 4: INEQUALITY CONSTRAINTS INCOMPATIBLE */ + /* Parameter adjustments */ + --indx; + --h__; + --x; + g_dim1 = *mg; + g_offset = 1 + g_dim1; + g -= g_offset; + --w; + + /* Function Body */ + *mode = 2; + if (*n <= 0) { + goto L50; + } +/* STATE DUAL PROBLEM */ + *mode = 1; + x[1] = 0.0; + dcopy___(n, &x[1], 0, &x[1], 1); + *xnorm = 0.0; + if (*m == 0) { + goto L50; + } + iw = 0; + i__1 = *m; + for (j = 1; j <= i__1; ++j) { + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { + ++iw; +/* L10: */ + w[iw] = g[j + i__ * g_dim1]; + } + ++iw; +/* L20: */ + w[iw] = h__[j]; + } + if__ = iw + 1; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + ++iw; +/* L30: */ + w[iw] = 0.0; + } + w[iw + 1] = one; + n1 = *n + 1; + iz = iw + 2; + iy = iz + n1; + iwdual = iy + *m; +/* SOLVE DUAL PROBLEM */ + nnls_(&w[1], &n1, &n1, m, &w[if__], &w[iy], &rnorm, &w[iwdual], &w[iz], & + indx[1], mode); + if (*mode != 1) { + goto L50; + } + *mode = 4; + if (rnorm <= 0.0) { + goto L50; + } +/* COMPUTE SOLUTION OF PRIMAL PROBLEM */ + fac = one - ddot_sl__(m, &h__[1], 1, &w[iy], 1); + d__1 = one + fac; + if (d__1 - one <= 0.0) { + goto L50; + } + *mode = 1; + fac = one / fac; + i__1 = *n; + for (j = 1; j <= i__1; ++j) { +/* L40: */ + x[j] = fac * ddot_sl__(m, &g[j * g_dim1 + 1], 1, &w[iy], 1); + } + *xnorm = dnrm2___(n, &x[1], 1); +/* COMPUTE LAGRANGE MULTIPLIERS FOR PRIMAL PROBLEM */ + w[1] = 0.0; + dcopy___(m, &w[1], 0, &w[1], 1); + daxpy_sl__(m, &fac, &w[iy], 1, &w[1], 1); +/* END OF SUBROUTINE LDP */ +L50: + return; +} /* ldp_ */ + +static void lsi_(double *e, double *f, double *g, + double *h__, int *le, int *me, int *lg, int *mg, + int *n, double *x, double *xnorm, double *w, int * + jw, int *mode) +{ + /* Initialized data */ + + const double epmach = 2.22e-16; + const double one = 1.; + + /* System generated locals */ + int e_dim1, e_offset, g_dim1, g_offset, i__1, i__2, i__3; + double d__1; + + /* Local variables */ + int i__, j; + double t; + +/* FOR MODE=1, THE SUBROUTINE RETURNS THE SOLUTION X OF */ +/* INEQUALITY CONSTRAINED LINEAR LEAST SQUARES PROBLEM: */ +/* MIN ||E*X-F|| */ +/* X */ +/* S.T. G*X >= H */ +/* THE ALGORITHM IS BASED ON QR DECOMPOSITION AS DESCRIBED IN */ +/* CHAPTER 23.5 OF LAWSON & HANSON: SOLVING LEAST SQUARES PROBLEMS */ +/* THE FOLLOWING DIMENSIONS OF THE ARRAYS DEFINING THE PROBLEM */ +/* ARE NECESSARY */ +/* DIM(E) : FORMAL (LE,N), ACTUAL (ME,N) */ +/* DIM(F) : FORMAL (LE ), ACTUAL (ME ) */ +/* DIM(G) : FORMAL (LG,N), ACTUAL (MG,N) */ +/* DIM(H) : FORMAL (LG ), ACTUAL (MG ) */ +/* DIM(X) : N */ +/* DIM(W) : (N+1)*(MG+2) + 2*MG */ +/* DIM(JW): LG */ +/* ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS E, F, G, AND H. */ +/* ON RETURN, ALL ARRAYS WILL BE CHANGED BY THE SUBROUTINE. */ +/* X STORES THE SOLUTION VECTOR */ +/* XNORM STORES THE RESIDUUM OF THE SOLUTION IN EUCLIDIAN NORM */ +/* W STORES THE VECTOR OF LAGRANGE MULTIPLIERS IN ITS FIRST */ +/* MG ELEMENTS */ +/* MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS: */ +/* MODE=1: SUCCESSFUL COMPUTATION */ +/* 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1) */ +/* 3: ITERATION COUNT EXCEEDED BY NNLS */ +/* 4: INEQUALITY CONSTRAINTS INCOMPATIBLE */ +/* 5: MATRIX E IS NOT OF FULL RANK */ +/* 03.01.1980, DIETER KRAFT: CODED */ +/* 20.03.1987, DIETER KRAFT: REVISED TO FORTRAN 77 */ + /* Parameter adjustments */ + --f; + --jw; + --h__; + --x; + g_dim1 = *lg; + g_offset = 1 + g_dim1; + g -= g_offset; + e_dim1 = *le; + e_offset = 1 + e_dim1; + e -= e_offset; + --w; + + /* Function Body */ +/* QR-FACTORS OF E AND APPLICATION TO F */ + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { +/* Computing MIN */ + i__2 = i__ + 1; + j = MIN2(i__2,*n); + i__2 = i__ + 1; + i__3 = *n - i__; + h12_(&c__1, &i__, &i__2, me, &e[i__ * e_dim1 + 1], &c__1, &t, &e[j * + e_dim1 + 1], &c__1, le, &i__3); +/* L10: */ + i__2 = i__ + 1; + h12_(&c__2, &i__, &i__2, me, &e[i__ * e_dim1 + 1], &c__1, &t, &f[1], & + c__1, &c__1, &c__1); + } +/* TRANSFORM G AND H TO GET LEAST DISTANCE PROBLEM */ + *mode = 5; + i__2 = *mg; + for (i__ = 1; i__ <= i__2; ++i__) { + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + if ((d__1 = e[j + j * e_dim1], fabs(d__1)) < epmach) { + goto L50; + } +/* L20: */ + i__3 = j - 1; + g[i__ + j * g_dim1] = (g[i__ + j * g_dim1] - ddot_sl__(&i__3, &g[ + i__ + g_dim1], *lg, &e[j * e_dim1 + 1], 1)) / e[j + j * + e_dim1]; + } +/* L30: */ + h__[i__] -= ddot_sl__(n, &g[i__ + g_dim1], *lg, &f[1], 1); + } +/* SOLVE LEAST DISTANCE PROBLEM */ + ldp_(&g[g_offset], lg, mg, n, &h__[1], &x[1], xnorm, &w[1], &jw[1], mode); + if (*mode != 1) { + goto L50; + } +/* SOLUTION OF ORIGINAL PROBLEM */ + daxpy_sl__(n, &one, &f[1], 1, &x[1], 1); + for (i__ = *n; i__ >= 1; --i__) { +/* Computing MIN */ + i__2 = i__ + 1; + j = MIN2(i__2,*n); +/* L40: */ + i__2 = *n - i__; + x[i__] = (x[i__] - ddot_sl__(&i__2, &e[i__ + j * e_dim1], *le, &x[j], 1)) + / e[i__ + i__ * e_dim1]; + } +/* Computing MIN */ + i__2 = *n + 1; + j = MIN2(i__2,*me); + i__2 = *me - *n; + t = dnrm2___(&i__2, &f[j], 1); + *xnorm = sqrt(*xnorm * *xnorm + t * t); +/* END OF SUBROUTINE LSI */ +L50: + return; +} /* lsi_ */ + +static void hfti_(double *a, int *mda, int *m, int * + n, double *b, int *mdb, const int *nb, double *tau, int + *krank, double *rnorm, double *h__, double *g, int * + ip) +{ + /* Initialized data */ + + const double factor = .001; + + /* System generated locals */ + int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3; + double d__1; + + /* Local variables */ + int i__, j, k, l; + int jb, kp1; + double tmp, hmax; + int lmax, ldiag; + +/* RANK-DEFICIENT LEAST SQUARES ALGORITHM AS DESCRIBED IN: */ +/* C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 */ +/* TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 */ +/* A(*,*),MDA,M,N THE ARRAY A INITIALLY CONTAINS THE M x N MATRIX A */ +/* OF THE LEAST SQUARES PROBLEM AX = B. */ +/* THE FIRST DIMENSIONING PARAMETER MDA MUST SATISFY */ +/* MDA >= M. EITHER M >= N OR M < N IS PERMITTED. */ +/* THERE IS NO RESTRICTION ON THE RANK OF A. */ +/* THE MATRIX A WILL BE MODIFIED BY THE SUBROUTINE. */ +/* B(*,*),MDB,NB IF NB = 0 THE SUBROUTINE WILL MAKE NO REFERENCE */ +/* TO THE ARRAY B. IF NB > 0 THE ARRAY B() MUST */ +/* INITIALLY CONTAIN THE M x NB MATRIX B OF THE */ +/* THE LEAST SQUARES PROBLEM AX = B AND ON RETURN */ +/* THE ARRAY B() WILL CONTAIN THE N x NB SOLUTION X. */ +/* IF NB>1 THE ARRAY B() MUST BE DOUBLE SUBSCRIPTED */ +/* WITH FIRST DIMENSIONING PARAMETER MDB>=MAX(M,N), */ +/* IF NB=1 THE ARRAY B() MAY BE EITHER SINGLE OR */ +/* DOUBLE SUBSCRIPTED. */ +/* TAU ABSOLUTE TOLERANCE PARAMETER FOR PSEUDORANK */ +/* DETERMINATION, PROVIDED BY THE USER. */ +/* KRANK PSEUDORANK OF A, SET BY THE SUBROUTINE. */ +/* RNORM ON EXIT, RNORM(J) WILL CONTAIN THE EUCLIDIAN */ +/* NORM OF THE RESIDUAL VECTOR FOR THE PROBLEM */ +/* DEFINED BY THE J-TH COLUMN VECTOR OF THE ARRAY B. */ +/* H(), G() ARRAYS OF WORKING SPACE OF LENGTH >= N. */ +/* IP() INT ARRAY OF WORKING SPACE OF LENGTH >= N */ +/* RECORDING PERMUTATION INDICES OF COLUMN VECTORS */ + /* Parameter adjustments */ + --ip; + --g; + --h__; + a_dim1 = *mda; + a_offset = 1 + a_dim1; + a -= a_offset; + --rnorm; + b_dim1 = *mdb; + b_offset = 1 + b_dim1; + b -= b_offset; + + /* Function Body */ + k = 0; + ldiag = MIN2(*m,*n); + if (ldiag <= 0) { + goto L270; + } +/* COMPUTE LMAX */ + i__1 = ldiag; + for (j = 1; j <= i__1; ++j) { + if (j == 1) { + goto L20; + } + lmax = j; + i__2 = *n; + for (l = j; l <= i__2; ++l) { +/* Computing 2nd power */ + d__1 = a[j - 1 + l * a_dim1]; + h__[l] -= d__1 * d__1; +/* L10: */ + if (h__[l] > h__[lmax]) { + lmax = l; + } + } + d__1 = hmax + factor * h__[lmax]; + if (d__1 - hmax > 0.0) { + goto L50; + } +L20: + lmax = j; + i__2 = *n; + for (l = j; l <= i__2; ++l) { + h__[l] = 0.0; + i__3 = *m; + for (i__ = j; i__ <= i__3; ++i__) { +/* L30: */ +/* Computing 2nd power */ + d__1 = a[i__ + l * a_dim1]; + h__[l] += d__1 * d__1; + } +/* L40: */ + if (h__[l] > h__[lmax]) { + lmax = l; + } + } + hmax = h__[lmax]; +/* COLUMN INTERCHANGES IF NEEDED */ +L50: + ip[j] = lmax; + if (ip[j] == j) { + goto L70; + } + i__2 = *m; + for (i__ = 1; i__ <= i__2; ++i__) { + tmp = a[i__ + j * a_dim1]; + a[i__ + j * a_dim1] = a[i__ + lmax * a_dim1]; +/* L60: */ + a[i__ + lmax * a_dim1] = tmp; + } + h__[lmax] = h__[j]; +/* J-TH TRANSFORMATION AND APPLICATION TO A AND B */ +L70: +/* Computing MIN */ + i__2 = j + 1; + i__ = MIN2(i__2,*n); + i__2 = j + 1; + i__3 = *n - j; + h12_(&c__1, &j, &i__2, m, &a[j * a_dim1 + 1], &c__1, &h__[j], &a[i__ * + a_dim1 + 1], &c__1, mda, &i__3); +/* L80: */ + i__2 = j + 1; + h12_(&c__2, &j, &i__2, m, &a[j * a_dim1 + 1], &c__1, &h__[j], &b[ + b_offset], &c__1, mdb, nb); + } +/* DETERMINE PSEUDORANK */ + i__2 = ldiag; + for (j = 1; j <= i__2; ++j) { +/* L90: */ + if ((d__1 = a[j + j * a_dim1], fabs(d__1)) <= *tau) { + goto L100; + } + } + k = ldiag; + goto L110; +L100: + k = j - 1; +L110: + kp1 = k + 1; +/* NORM OF RESIDUALS */ + i__2 = *nb; + for (jb = 1; jb <= i__2; ++jb) { +/* L130: */ + i__1 = *m - k; + rnorm[jb] = dnrm2___(&i__1, &b[kp1 + jb * b_dim1], 1); + } + if (k > 0) { + goto L160; + } + i__1 = *nb; + for (jb = 1; jb <= i__1; ++jb) { + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L150: */ + b[i__ + jb * b_dim1] = 0.0; + } + } + goto L270; +L160: + if (k == *n) { + goto L180; + } +/* HOUSEHOLDER DECOMPOSITION OF FIRST K ROWS */ + for (i__ = k; i__ >= 1; --i__) { +/* L170: */ + i__2 = i__ - 1; + h12_(&c__1, &i__, &kp1, n, &a[i__ + a_dim1], mda, &g[i__], &a[ + a_offset], mda, &c__1, &i__2); + } +L180: + i__2 = *nb; + for (jb = 1; jb <= i__2; ++jb) { +/* SOLVE K*K TRIANGULAR SYSTEM */ + for (i__ = k; i__ >= 1; --i__) { +/* Computing MIN */ + i__1 = i__ + 1; + j = MIN2(i__1,*n); +/* L210: */ + i__1 = k - i__; + b[i__ + jb * b_dim1] = (b[i__ + jb * b_dim1] - ddot_sl__(&i__1, & + a[i__ + j * a_dim1], *mda, &b[j + jb * b_dim1], 1)) / + a[i__ + i__ * a_dim1]; + } +/* COMPLETE SOLUTION VECTOR */ + if (k == *n) { + goto L240; + } + i__1 = *n; + for (j = kp1; j <= i__1; ++j) { +/* L220: */ + b[j + jb * b_dim1] = 0.0; + } + i__1 = k; + for (i__ = 1; i__ <= i__1; ++i__) { +/* L230: */ + h12_(&c__2, &i__, &kp1, n, &a[i__ + a_dim1], mda, &g[i__], &b[jb * + b_dim1 + 1], &c__1, mdb, &c__1); + } +/* REORDER SOLUTION ACCORDING TO PREVIOUS COLUMN INTERCHANGES */ +L240: + for (j = ldiag; j >= 1; --j) { + if (ip[j] == j) { + goto L250; + } + l = ip[j]; + tmp = b[l + jb * b_dim1]; + b[l + jb * b_dim1] = b[j + jb * b_dim1]; + b[j + jb * b_dim1] = tmp; +L250: + ; + } + } +L270: + *krank = k; +} /* hfti_ */ + +static void lsei_(double *c__, double *d__, double *e, + double *f, double *g, double *h__, int *lc, int * + mc, int *le, int *me, int *lg, int *mg, int *n, + double *x, double *xnrm, double *w, int *jw, int * + mode) +{ + /* Initialized data */ + + const double epmach = 2.22e-16; + + /* System generated locals */ + int c_dim1, c_offset, e_dim1, e_offset, g_dim1, g_offset, i__1, i__2, + i__3; + double d__1; + + /* Local variables */ + int i__, j, k, l; + double t; + int ie, if__, ig, iw, mc1; + int krank; + +/* FOR MODE=1, THE SUBROUTINE RETURNS THE SOLUTION X OF */ +/* EQUALITY & INEQUALITY CONSTRAINED LEAST SQUARES PROBLEM LSEI : */ +/* MIN ||E*X - F|| */ +/* X */ +/* S.T. C*X = D, */ +/* G*X >= H. */ +/* USING QR DECOMPOSITION & ORTHOGONAL BASIS OF NULLSPACE OF C */ +/* CHAPTER 23.6 OF LAWSON & HANSON: SOLVING LEAST SQUARES PROBLEMS. */ +/* THE FOLLOWING DIMENSIONS OF THE ARRAYS DEFINING THE PROBLEM */ +/* ARE NECESSARY */ +/* DIM(E) : FORMAL (LE,N), ACTUAL (ME,N) */ +/* DIM(F) : FORMAL (LE ), ACTUAL (ME ) */ +/* DIM(C) : FORMAL (LC,N), ACTUAL (MC,N) */ +/* DIM(D) : FORMAL (LC ), ACTUAL (MC ) */ +/* DIM(G) : FORMAL (LG,N), ACTUAL (MG,N) */ +/* DIM(H) : FORMAL (LG ), ACTUAL (MG ) */ +/* DIM(X) : FORMAL (N ), ACTUAL (N ) */ +/* DIM(W) : 2*MC+ME+(ME+MG)*(N-MC) for LSEI */ +/* +(N-MC+1)*(MG+2)+2*MG for LSI */ +/* DIM(JW): MAX(MG,L) */ +/* ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS C, D, E, F, G, AND H. */ +/* ON RETURN, ALL ARRAYS WILL BE CHANGED BY THE SUBROUTINE. */ +/* X STORES THE SOLUTION VECTOR */ +/* XNORM STORES THE RESIDUUM OF THE SOLUTION IN EUCLIDIAN NORM */ +/* W STORES THE VECTOR OF LAGRANGE MULTIPLIERS IN ITS FIRST */ +/* MC+MG ELEMENTS */ +/* MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS: */ +/* MODE=1: SUCCESSFUL COMPUTATION */ +/* 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1) */ +/* 3: ITERATION COUNT EXCEEDED BY NNLS */ +/* 4: INEQUALITY CONSTRAINTS INCOMPATIBLE */ +/* 5: MATRIX E IS NOT OF FULL RANK */ +/* 6: MATRIX C IS NOT OF FULL RANK */ +/* 7: RANK DEFECT IN HFTI */ +/* 18.5.1981, DIETER KRAFT, DFVLR OBERPFAFFENHOFEN */ +/* 20.3.1987, DIETER KRAFT, DFVLR OBERPFAFFENHOFEN */ + /* Parameter adjustments */ + --d__; + --f; + --h__; + --x; + g_dim1 = *lg; + g_offset = 1 + g_dim1; + g -= g_offset; + e_dim1 = *le; + e_offset = 1 + e_dim1; + e -= e_offset; + c_dim1 = *lc; + c_offset = 1 + c_dim1; + c__ -= c_offset; + --w; + --jw; + + /* Function Body */ + *mode = 2; + if (*mc > *n) { + goto L75; + } + l = *n - *mc; + mc1 = *mc + 1; + iw = (l + 1) * (*mg + 2) + (*mg << 1) + *mc; + ie = iw + *mc + 1; + if__ = ie + *me * l; + ig = if__ + *me; +/* TRIANGULARIZE C AND APPLY FACTORS TO E AND G */ + i__1 = *mc; + for (i__ = 1; i__ <= i__1; ++i__) { +/* Computing MIN */ + i__2 = i__ + 1; + j = MIN2(i__2,*lc); + i__2 = i__ + 1; + i__3 = *mc - i__; + h12_(&c__1, &i__, &i__2, n, &c__[i__ + c_dim1], lc, &w[iw + i__], & + c__[j + c_dim1], lc, &c__1, &i__3); + i__2 = i__ + 1; + h12_(&c__2, &i__, &i__2, n, &c__[i__ + c_dim1], lc, &w[iw + i__], &e[ + e_offset], le, &c__1, me); +/* L10: */ + i__2 = i__ + 1; + h12_(&c__2, &i__, &i__2, n, &c__[i__ + c_dim1], lc, &w[iw + i__], &g[ + g_offset], lg, &c__1, mg); + } +/* SOLVE C*X=D AND MODIFY F */ + *mode = 6; + i__2 = *mc; + for (i__ = 1; i__ <= i__2; ++i__) { + if ((d__1 = c__[i__ + i__ * c_dim1], fabs(d__1)) < epmach) { + goto L75; + } + i__1 = i__ - 1; + x[i__] = (d__[i__] - ddot_sl__(&i__1, &c__[i__ + c_dim1], *lc, &x[1], 1)) + / c__[i__ + i__ * c_dim1]; +/* L15: */ + } + *mode = 1; + w[mc1] = 0.0; + i__2 = *mg; /* BUGFIX for *mc == *n: changed from *mg - *mc, SGJ 2010 */ + dcopy___(&i__2, &w[mc1], 0, &w[mc1], 1); + if (*mc == *n) { + goto L50; + } + i__2 = *me; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L20: */ + w[if__ - 1 + i__] = f[i__] - ddot_sl__(mc, &e[i__ + e_dim1], *le, &x[1], 1); + } +/* STORE TRANSFORMED E & G */ + i__2 = *me; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L25: */ + dcopy___(&l, &e[i__ + mc1 * e_dim1], *le, &w[ie - 1 + i__], *me); + } + i__2 = *mg; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L30: */ + dcopy___(&l, &g[i__ + mc1 * g_dim1], *lg, &w[ig - 1 + i__], *mg); + } + if (*mg > 0) { + goto L40; + } +/* SOLVE LS WITHOUT INEQUALITY CONSTRAINTS */ + *mode = 7; + k = MAX2(*le,*n); + t = sqrt(epmach); + hfti_(&w[ie], me, me, &l, &w[if__], &k, &c__1, &t, &krank, xnrm, &w[1], & + w[l + 1], &jw[1]); + dcopy___(&l, &w[if__], 1, &x[mc1], 1); + if (krank != l) { + goto L75; + } + *mode = 1; + goto L50; +/* MODIFY H AND SOLVE INEQUALITY CONSTRAINED LS PROBLEM */ +L40: + i__2 = *mg; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L45: */ + h__[i__] -= ddot_sl__(mc, &g[i__ + g_dim1], *lg, &x[1], 1); + } + lsi_(&w[ie], &w[if__], &w[ig], &h__[1], me, me, mg, mg, &l, &x[mc1], xnrm, + &w[mc1], &jw[1], mode); + if (*mc == 0) { + goto L75; + } + t = dnrm2___(mc, &x[1], 1); + *xnrm = sqrt(*xnrm * *xnrm + t * t); + if (*mode != 1) { + goto L75; + } +/* SOLUTION OF ORIGINAL PROBLEM AND LAGRANGE MULTIPLIERS */ +L50: + i__2 = *me; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L55: */ + f[i__] = ddot_sl__(n, &e[i__ + e_dim1], *le, &x[1], 1) - f[i__]; + } + i__2 = *mc; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L60: */ + d__[i__] = ddot_sl__(me, &e[i__ * e_dim1 + 1], 1, &f[1], 1) - + ddot_sl__(mg, &g[i__ * g_dim1 + 1], 1, &w[mc1], 1); + } + for (i__ = *mc; i__ >= 1; --i__) { +/* L65: */ + i__2 = i__ + 1; + h12_(&c__2, &i__, &i__2, n, &c__[i__ + c_dim1], lc, &w[iw + i__], &x[ + 1], &c__1, &c__1, &c__1); + } + for (i__ = *mc; i__ >= 1; --i__) { +/* Computing MIN */ + i__2 = i__ + 1; + j = MIN2(i__2,*lc); + i__2 = *mc - i__; + w[i__] = (d__[i__] - ddot_sl__(&i__2, &c__[j + i__ * c_dim1], 1, & + w[j], 1)) / c__[i__ + i__ * c_dim1]; +/* L70: */ + } +/* END OF SUBROUTINE LSEI */ +L75: + return; +} /* lsei_ */ + +static void lsq_(int *m, int *meq, int *n, int *nl, + int *la, double *l, double *g, double *a, double * + b, const double *xl, const double *xu, double *x, double *y, + double *w, int *jw, int *mode) +{ + /* Initialized data */ + + const double one = 1.; + + /* System generated locals */ + int a_dim1, a_offset, i__1, i__2; + double d__1; + + /* Local variables */ + int i__, i1, i2, i3, i4, m1, n1, n2, n3, ic, id, ie, if__, ig, ih, il, + im, ip, iu, iw; + double diag; + int mineq; + double xnorm; + +/* MINIMIZE with respect to X */ +/* ||E*X - F|| */ +/* 1/2 T */ +/* WITH UPPER TRIANGULAR MATRIX E = +D *L , */ +/* -1/2 -1 */ +/* AND VECTOR F = -D *L *G, */ +/* WHERE THE UNIT LOWER TRIDIANGULAR MATRIX L IS STORED COLUMNWISE */ +/* DENSE IN THE N*(N+1)/2 ARRAY L WITH VECTOR D STORED IN ITS */ +/* 'DIAGONAL' THUS SUBSTITUTING THE ONE-ELEMENTS OF L */ +/* SUBJECT TO */ +/* A(J)*X - B(J) = 0 , J=1,...,MEQ, */ +/* A(J)*X - B(J) >=0, J=MEQ+1,...,M, */ +/* XL(I) <= X(I) <= XU(I), I=1,...,N, */ +/* ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS L, G, A, B, XL, XU. */ +/* WITH DIMENSIONS: L(N*(N+1)/2), G(N), A(LA,N), B(M), XL(N), XU(N) */ +/* THE WORKING ARRAY W MUST HAVE AT LEAST THE FOLLOWING DIMENSION: */ +/* DIM(W) = (3*N+M)*(N+1) for LSQ */ +/* +(N-MEQ+1)*(MINEQ+2) + 2*MINEQ for LSI */ +/* +(N+MINEQ)*(N-MEQ) + 2*MEQ + N for LSEI */ +/* with MINEQ = M - MEQ + 2*N */ +/* ON RETURN, NO ARRAY WILL BE CHANGED BY THE SUBROUTINE. */ +/* X STORES THE N-DIMENSIONAL SOLUTION VECTOR */ +/* Y STORES THE VECTOR OF LAGRANGE MULTIPLIERS OF DIMENSION */ +/* M+N+N (CONSTRAINTS+LOWER+UPPER BOUNDS) */ +/* MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS: */ +/* MODE=1: SUCCESSFUL COMPUTATION */ +/* 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1) */ +/* 3: ITERATION COUNT EXCEEDED BY NNLS */ +/* 4: INEQUALITY CONSTRAINTS INCOMPATIBLE */ +/* 5: MATRIX E IS NOT OF FULL RANK */ +/* 6: MATRIX C IS NOT OF FULL RANK */ +/* 7: RANK DEFECT IN HFTI */ +/* coded Dieter Kraft, april 1987 */ +/* revised march 1989 */ + /* Parameter adjustments */ + --y; + --x; + --xu; + --xl; + --g; + --l; + --b; + a_dim1 = *la; + a_offset = 1 + a_dim1; + a -= a_offset; + --w; + --jw; + + /* Function Body */ + n1 = *n + 1; + mineq = *m - *meq; + m1 = mineq + *n + *n; +/* determine whether to solve problem */ +/* with inconsistent linerarization (n2=1) */ +/* or not (n2=0) */ + n2 = n1 * *n / 2 + 1; + if (n2 == *nl) { + n2 = 0; + } else { + n2 = 1; + } + n3 = *n - n2; +/* RECOVER MATRIX E AND VECTOR F FROM L AND G */ + i2 = 1; + i3 = 1; + i4 = 1; + ie = 1; + if__ = *n * *n + 1; + i__1 = n3; + for (i__ = 1; i__ <= i__1; ++i__) { + i1 = n1 - i__; + diag = sqrt(l[i2]); + w[i3] = 0.0; + dcopy___(&i1, &w[i3], 0, &w[i3], 1); + i__2 = i1 - n2; + dcopy___(&i__2, &l[i2], 1, &w[i3], *n); + i__2 = i1 - n2; + dscal_sl__(&i__2, &diag, &w[i3], *n); + w[i3] = diag; + i__2 = i__ - 1; + w[if__ - 1 + i__] = (g[i__] - ddot_sl__(&i__2, &w[i4], 1, &w[if__] + , 1)) / diag; + i2 = i2 + i1 - n2; + i3 += n1; + i4 += *n; +/* L10: */ + } + if (n2 == 1) { + w[i3] = l[*nl]; + w[i4] = 0.0; + dcopy___(&n3, &w[i4], 0, &w[i4], 1); + w[if__ - 1 + *n] = 0.0; + } + d__1 = -one; + dscal_sl__(n, &d__1, &w[if__], 1); + ic = if__ + *n; + id = ic + *meq * *n; + if (*meq > 0) { +/* RECOVER MATRIX C FROM UPPER PART OF A */ + i__1 = *meq; + for (i__ = 1; i__ <= i__1; ++i__) { + dcopy___(n, &a[i__ + a_dim1], *la, &w[ic - 1 + i__], *meq); +/* L20: */ + } +/* RECOVER VECTOR D FROM UPPER PART OF B */ + dcopy___(meq, &b[1], 1, &w[id], 1); + d__1 = -one; + dscal_sl__(meq, &d__1, &w[id], 1); + } + ig = id + *meq; + if (mineq > 0) { +/* RECOVER MATRIX G FROM LOWER PART OF A */ + i__1 = mineq; + for (i__ = 1; i__ <= i__1; ++i__) { + dcopy___(n, &a[*meq + i__ + a_dim1], *la, &w[ig - 1 + i__], m1); +/* L30: */ + } + } +/* AUGMENT MATRIX G BY +I AND -I */ + ip = ig + mineq; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + w[ip - 1 + i__] = 0.0; + dcopy___(n, &w[ip - 1 + i__], 0, &w[ip - 1 + i__], m1); +/* L40: */ + } + i__1 = m1 + 1; + /* SGJ, 2010: skip constraints for infinite bounds */ + for (i__ = 1; i__ <= *n; ++i__) + if (!nlopt_isinf(xl[i__])) w[(ip - i__1) + i__ * i__1] = +1.0; + /* Old code: w[ip] = one; dcopy___(n, &w[ip], 0, &w[ip], i__1); */ + im = ip + *n; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + w[im - 1 + i__] = 0.0; + dcopy___(n, &w[im - 1 + i__], 0, &w[im - 1 + i__], m1); +/* L50: */ + } + i__1 = m1 + 1; + /* SGJ, 2010: skip constraints for infinite bounds */ + for (i__ = 1; i__ <= *n; ++i__) + if (!nlopt_isinf(xu[i__])) w[(im - i__1) + i__ * i__1] = -1.0; + /* Old code: w[im] = -one; dcopy___(n, &w[im], 0, &w[im], i__1); */ + ih = ig + m1 * *n; + if (mineq > 0) { +/* RECOVER H FROM LOWER PART OF B */ + dcopy___(&mineq, &b[*meq + 1], 1, &w[ih], 1); + d__1 = -one; + dscal_sl__(&mineq, &d__1, &w[ih], 1); + } +/* AUGMENT VECTOR H BY XL AND XU */ + il = ih + mineq; + iu = il + *n; + /* SGJ, 2010: skip constraints for infinite bounds */ + for (i__ = 1; i__ <= *n; ++i__) { + w[(il-1) + i__] = nlopt_isinf(xl[i__]) ? 0 : xl[i__]; + w[(iu-1) + i__] = nlopt_isinf(xu[i__]) ? 0 : -xu[i__]; + } + /* Old code: dcopy___(n, &xl[1], 1, &w[il], 1); + dcopy___(n, &xu[1], 1, &w[iu], 1); + d__1 = -one; dscal_sl__(n, &d__1, &w[iu], 1); */ + iw = iu + *n; + i__1 = MAX2(1,*meq); + lsei_(&w[ic], &w[id], &w[ie], &w[if__], &w[ig], &w[ih], &i__1, meq, n, n, + &m1, &m1, n, &x[1], &xnorm, &w[iw], &jw[1], mode); + if (*mode == 1) { +/* restore Lagrange multipliers */ + dcopy___(m, &w[iw], 1, &y[1], 1); + dcopy___(&n3, &w[iw + *m], 1, &y[*m + 1], 1); + dcopy___(&n3, &w[iw + *m + *n], 1, &y[*m + n3 + 1], 1); + + /* SGJ, 2010: make sure bound constraints are satisfied, since + roundoff error sometimes causes slight violations and + NLopt guarantees that bounds are strictly obeyed */ + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + if (x[i__] < xl[i__]) x[i__] = xl[i__]; + else if (x[i__] > xu[i__]) x[i__] = xu[i__]; + } + } +/* END OF SUBROUTINE LSQ */ +} /* lsq_ */ + +static void ldl_(int *n, double *a, double *z__, + double *sigma, double *w) +{ + /* Initialized data */ + + const double one = 1.; + const double four = 4.; + const double epmach = 2.22e-16; + + /* System generated locals */ + int i__1, i__2; + + /* Local variables */ + int i__, j; + double t, u, v; + int ij; + double tp, beta, gamma_, alpha, delta; + +/* LDL LDL' - RANK-ONE - UPDATE */ +/* PURPOSE: */ +/* UPDATES THE LDL' FACTORS OF MATRIX A BY RANK-ONE MATRIX */ +/* SIGMA*Z*Z' */ +/* INPUT ARGUMENTS: (* MEANS PARAMETERS ARE CHANGED DURING EXECUTION) */ +/* N : ORDER OF THE COEFFICIENT MATRIX A */ +/* * A : POSITIVE DEFINITE MATRIX OF DIMENSION N; */ +/* ONLY THE LOWER TRIANGLE IS USED AND IS STORED COLUMN BY */ +/* COLUMN AS ONE DIMENSIONAL ARRAY OF DIMENSION N*(N+1)/2. */ +/* * Z : VECTOR OF DIMENSION N OF UPDATING ELEMENTS */ +/* SIGMA : SCALAR FACTOR BY WHICH THE MODIFYING DYADE Z*Z' IS */ +/* MULTIPLIED */ +/* OUTPUT ARGUMENTS: */ +/* A : UPDATED LDL' FACTORS */ +/* WORKING ARRAY: */ +/* W : VECTOR OP DIMENSION N (USED ONLY IF SIGMA .LT. ZERO) */ +/* METHOD: */ +/* THAT OF FLETCHER AND POWELL AS DESCRIBED IN : */ +/* FLETCHER,R.,(1974) ON THE MODIFICATION OF LDL' FACTORIZATION. */ +/* POWELL,M.J.D. MATH.COMPUTATION 28, 1067-1078. */ +/* IMPLEMENTED BY: */ +/* KRAFT,D., DFVLR - INSTITUT FUER DYNAMIK DER FLUGSYSTEME */ +/* D-8031 OBERPFAFFENHOFEN */ +/* STATUS: 15. JANUARY 1980 */ +/* SUBROUTINES REQUIRED: NONE */ + /* Parameter adjustments */ + --w; + --z__; + --a; + + /* Function Body */ + if (*sigma == 0.0) { + goto L280; + } + ij = 1; + t = one / *sigma; + if (*sigma > 0.0) { + goto L220; + } +/* PREPARE NEGATIVE UPDATE */ + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { +/* L150: */ + w[i__] = z__[i__]; + } + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + v = w[i__]; + t += v * v / a[ij]; + i__2 = *n; + for (j = i__ + 1; j <= i__2; ++j) { + ++ij; +/* L160: */ + w[j] -= v * a[ij]; + } +/* L170: */ + ++ij; + } + if (t >= 0.0) { + t = epmach / *sigma; + } + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + j = *n + 1 - i__; + ij -= i__; + u = w[j]; + w[j] = t; +/* L210: */ + t -= u * u / a[ij]; + } +L220: +/* HERE UPDATING BEGINS */ + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + v = z__[i__]; + delta = v / a[ij]; + if (*sigma < 0.0) { + tp = w[i__]; + } + else /* if (*sigma > 0.0), since *sigma != 0 from above */ { + tp = t + delta * v; + } + alpha = tp / t; + a[ij] = alpha * a[ij]; + if (i__ == *n) { + goto L280; + } + beta = delta / tp; + if (alpha > four) { + goto L240; + } + i__2 = *n; + for (j = i__ + 1; j <= i__2; ++j) { + ++ij; + z__[j] -= v * a[ij]; +/* L230: */ + a[ij] += beta * z__[j]; + } + goto L260; +L240: + gamma_ = t / tp; + i__2 = *n; + for (j = i__ + 1; j <= i__2; ++j) { + ++ij; + u = a[ij]; + a[ij] = gamma_ * u + beta * z__[j]; +/* L250: */ + z__[j] -= v * u; + } +L260: + ++ij; +/* L270: */ + t = tp; + } +L280: + return; +/* END OF LDL */ +} /* ldl_ */ + +#if 0 +/* we don't want to use this linmin function, for two reasons: + 1) it was apparently written assuming an old version of Fortran where all variables + are saved by default, hence it was missing a "save" statement ... I would + need to go through and figure out which variables need to be declared "static" + (or, better yet, save them like I did in slsqp to make it re-entrant) + 2) it doesn't exploit gradient information, which is stupid since we have this info + 3) in the context of NLopt, it makes much more sence to use the local_opt algorithm + to do the line minimization recursively by whatever algorithm the user wants + (defaulting to a gradient-based method like LBFGS) */ +static double d_sign(double a, double s) { return s < 0 ? -a : a; } +static double linmin_(int *mode, const double *ax, const double *bx, double * + f, double *tol) +{ + /* Initialized data */ + + const double c__ = .381966011; + const double eps = 1.5e-8; + + /* System generated locals */ + double ret_val, d__1; + + /* Local variables */ + double a, b, d__, e, m, p, q, r__, u, v, w, x, fu, fv, fw, fx, tol1, + tol2; + +/* LINMIN LINESEARCH WITHOUT DERIVATIVES */ +/* PURPOSE: */ +/* TO FIND THE ARGUMENT LINMIN WHERE THE FUNCTION F TAKES ITS MINIMUM */ +/* ON THE INTERVAL AX, BX. */ +/* COMBINATION OF GOLDEN SECTION AND SUCCESSIVE QUADRATIC INTERPOLATION. */ +/* INPUT ARGUMENTS: (* MEANS PARAMETERS ARE CHANGED DURING EXECUTION) */ +/* *MODE SEE OUTPUT ARGUMENTS */ +/* AX LEFT ENDPOINT OF INITIAL INTERVAL */ +/* BX RIGHT ENDPOINT OF INITIAL INTERVAL */ +/* F FUNCTION VALUE AT LINMIN WHICH IS TO BE BROUGHT IN BY */ +/* REVERSE COMMUNICATION CONTROLLED BY MODE */ +/* TOL DESIRED LENGTH OF INTERVAL OF UNCERTAINTY OF FINAL RESULT */ +/* OUTPUT ARGUMENTS: */ +/* LINMIN ABSCISSA APPROXIMATING THE POINT WHERE F ATTAINS A MINIMUM */ +/* MODE CONTROLS REVERSE COMMUNICATION */ +/* MUST BE SET TO 0 INITIALLY, RETURNS WITH INTERMEDIATE */ +/* VALUES 1 AND 2 WHICH MUST NOT BE CHANGED BY THE USER, */ +/* ENDS WITH CONVERGENCE WITH VALUE 3. */ +/* WORKING ARRAY: */ +/* NONE */ +/* METHOD: */ +/* THIS FUNCTION SUBPROGRAM IS A SLIGHTLY MODIFIED VERSION OF THE */ +/* ALGOL 60 PROCEDURE LOCALMIN GIVEN IN */ +/* R.P. BRENT: ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES, */ +/* PRENTICE-HALL (1973). */ +/* IMPLEMENTED BY: */ +/* KRAFT, D., DFVLR - INSTITUT FUER DYNAMIK DER FLUGSYSTEME */ +/* D-8031 OBERPFAFFENHOFEN */ +/* STATUS: 31. AUGUST 1984 */ +/* SUBROUTINES REQUIRED: NONE */ +/* EPS = SQUARE - ROOT OF MACHINE PRECISION */ +/* C = GOLDEN SECTION RATIO = (3-SQRT(5))/2 */ + switch (*mode) { + case 1: goto L10; + case 2: goto L55; + } +/* INITIALIZATION */ + a = *ax; + b = *bx; + e = 0.0; + v = a + c__ * (b - a); + w = v; + x = w; + ret_val = x; + *mode = 1; + goto L100; +/* MAIN LOOP STARTS HERE */ +L10: + fx = *f; + fv = fx; + fw = fv; +L20: + m = (a + b) * .5; + tol1 = eps * fabs(x) + *tol; + tol2 = tol1 + tol1; +/* TEST CONVERGENCE */ + if ((d__1 = x - m, fabs(d__1)) <= tol2 - (b - a) * .5) { + goto L90; + } + r__ = 0.0; + q = r__; + p = q; + if (fabs(e) <= tol1) { + goto L30; + } +/* FIT PARABOLA */ + r__ = (x - w) * (fx - fv); + q = (x - v) * (fx - fw); + p = (x - v) * q - (x - w) * r__; + q -= r__; + q += q; + if (q > 0.0) { + p = -p; + } + if (q < 0.0) { + q = -q; + } + r__ = e; + e = d__; +/* IS PARABOLA ACCEPTABLE */ +L30: + if (fabs(p) >= (d__1 = q * r__, fabs(d__1)) * .5 || p <= q * (a - x) || p >= + q * (b - x)) { + goto L40; + } +/* PARABOLIC INTERPOLATION STEP */ + d__ = p / q; +/* F MUST NOT BE EVALUATED TOO CLOSE TO A OR B */ + if (u - a < tol2) { + d__1 = m - x; + d__ = d_sign(tol1, d__1); + } + if (b - u < tol2) { + d__1 = m - x; + d__ = d_sign(tol1, d__1); + } + goto L50; +/* GOLDEN SECTION STEP */ +L40: + if (x >= m) { + e = a - x; + } + if (x < m) { + e = b - x; + } + d__ = c__ * e; +/* F MUST NOT BE EVALUATED TOO CLOSE TO X */ +L50: + if (fabs(d__) < tol1) { + d__ = d_sign(tol1, d__); + } + u = x + d__; + ret_val = u; + *mode = 2; + goto L100; +L55: + fu = *f; +/* UPDATE A, B, V, W, AND X */ + if (fu > fx) { + goto L60; + } + if (u >= x) { + a = x; + } + if (u < x) { + b = x; + } + v = w; + fv = fw; + w = x; + fw = fx; + x = u; + fx = fu; + goto L85; +L60: + if (u < x) { + a = u; + } + if (u >= x) { + b = u; + } + if (fu <= fw || w == x) { + goto L70; + } + if (fu <= fv || v == x || v == w) { + goto L80; + } + goto L85; +L70: + v = w; + fv = fw; + w = u; + fw = fu; + goto L85; +L80: + v = u; + fv = fu; +L85: + goto L20; +/* END OF MAIN LOOP */ +L90: + ret_val = x; + *mode = 3; +L100: + return ret_val; +/* END OF LINMIN */ +} /* linmin_ */ +#endif + +// typedef struct { +// double t, f0, h1, h2, h3, h4; +// int n1, n2, n3; +// double t0, gs; +// double tol; +// int line; +// double alpha; +// int iexact; +// int incons, ireset, itermx; +// double *x0; +// } slsqpb_state; + +#define SS(var) state->var = var +#define SAVE_STATE \ + SS(t); SS(f0); SS(h1); SS(h2); SS(h3); SS(h4); \ + SS(n1); SS(n2); SS(n3); \ + SS(t0); SS(gs); \ + SS(tol); \ + SS(line); \ + SS(alpha); \ + SS(iexact); \ + SS(incons); SS(ireset); SS(itermx) + +#define RS(var) var = state->var +#define RESTORE_STATE \ + RS(t); RS(f0); RS(h1); RS(h2); RS(h3); RS(h4); \ + RS(n1); RS(n2); RS(n3); \ + RS(t0); RS(gs); \ + RS(tol); \ + RS(line); \ + RS(alpha); \ + RS(iexact); \ + RS(incons); RS(ireset); RS(itermx) + +static void slsqpb_(int *m, int *meq, int *la, int * + n, double *x, const double *xl, const double *xu, double *f, + double *c__, double *g, double *a, double *acc, + int *iter, int *mode, double *r__, double *l, + double *x0, double *mu, double *s, double *u, + double *v, double *w, int *iw, + slsqpb_state *state) +{ + /* Initialized data */ + + const double one = 1.; + const double alfmin = .1; + const double hun = 100.; + const double ten = 10.; + const double two = 2.; + + /* System generated locals */ + int a_dim1, a_offset, i__1, i__2; + double d__1, d__2; + + /* Local variables */ + int i__, j, k; + + /* saved state from one call to the next; + SGJ 2010: save/restore via state parameter, to make re-entrant. */ + double t, f0, h1, h2, h3, h4; + int n1, n2, n3; + double t0, gs; + double tol; + int line; + double alpha; + int iexact; + int incons, ireset, itermx; + RESTORE_STATE; + +/* NONLINEAR PROGRAMMING BY SOLVING SEQUENTIALLY QUADRATIC PROGRAMS */ +/* - L1 - LINE SEARCH, POSITIVE DEFINITE BFGS UPDATE - */ +/* BODY SUBROUTINE FOR SLSQP */ +/* dim(W) = N1*(N1+1) + MEQ*(N1+1) + MINEQ*(N1+1) for LSQ */ +/* +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ */ +/* +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 for LSEI */ +/* with MINEQ = M - MEQ + 2*N1 & N1 = N+1 */ + /* Parameter adjustments */ + --mu; + --c__; + --v; + --u; + --s; + --x0; + --l; + --r__; + a_dim1 = *la; + a_offset = 1 + a_dim1; + a -= a_offset; + --g; + --xu; + --xl; + --x; + --w; + --iw; + + /* Function Body */ + if (*mode == -1) { + goto L260; + } else if (*mode == 0) { + goto L100; + } else { + goto L220; + } +L100: + itermx = *iter; + if (*acc >= 0.0) { + iexact = 0; + } else { + iexact = 1; + } + *acc = fabs(*acc); + tol = ten * *acc; + *iter = 0; + ireset = 0; + n1 = *n + 1; + n2 = n1 * *n / 2; + n3 = n2 + 1; + s[1] = 0.0; + mu[1] = 0.0; + dcopy___(n, &s[1], 0, &s[1], 1); + dcopy___(m, &mu[1], 0, &mu[1], 1); +/* RESET BFGS MATRIX */ +L110: + ++ireset; + if (ireset > 5) { + goto L255; + } + l[1] = 0.0; + dcopy___(&n2, &l[1], 0, &l[1], 1); + j = 1; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + l[j] = one; + j = j + n1 - i__; +/* L120: */ + } +/* MAIN ITERATION : SEARCH DIRECTION, STEPLENGTH, LDL'-UPDATE */ +L130: + ++(*iter); + *mode = 9; + if (*iter > itermx && itermx > 0) { /* SGJ 2010: ignore if itermx <= 0 */ + goto L330; + } +/* SEARCH DIRECTION AS SOLUTION OF QP - SUBPROBLEM */ + dcopy___(n, &xl[1], 1, &u[1], 1); + dcopy___(n, &xu[1], 1, &v[1], 1); + d__1 = -one; + daxpy_sl__(n, &d__1, &x[1], 1, &u[1], 1); + d__1 = -one; + daxpy_sl__(n, &d__1, &x[1], 1, &v[1], 1); + h4 = one; + lsq_(m, meq, n, &n3, la, &l[1], &g[1], &a[a_offset], &c__[1], &u[1], &v[1] + , &s[1], &r__[1], &w[1], &iw[1], mode); + +/* AUGMENTED PROBLEM FOR INCONSISTENT LINEARIZATION */ + if (*mode == 6) { + if (*n == *meq) { + *mode = 4; + } + } + if (*mode == 4) { + i__1 = *m; + for (j = 1; j <= i__1; ++j) { + if (j <= *meq) { + a[j + n1 * a_dim1] = -c__[j]; + } else { +/* Computing MAX */ + d__1 = -c__[j]; + a[j + n1 * a_dim1] = MAX2(d__1,0.0); + } +/* L140: */ + } + s[1] = 0.0; + dcopy___(n, &s[1], 0, &s[1], 1); + h3 = 0.0; + g[n1] = 0.0; + l[n3] = hun; + s[n1] = one; + u[n1] = 0.0; + v[n1] = one; + incons = 0; +L150: + lsq_(m, meq, &n1, &n3, la, &l[1], &g[1], &a[a_offset], &c__[1], &u[1], + &v[1], &s[1], &r__[1], &w[1], &iw[1], mode); + h4 = one - s[n1]; + if (*mode == 4) { + l[n3] = ten * l[n3]; + ++incons; + if (incons > 5) { + goto L330; + } + goto L150; + } else if (*mode != 1) { + goto L330; + } + } else if (*mode != 1) { + goto L330; + } +/* UPDATE MULTIPLIERS FOR L1-TEST */ + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + v[i__] = g[i__] - ddot_sl__(m, &a[i__ * a_dim1 + 1], 1, &r__[1], 1); +/* L160: */ + } + f0 = *f; + dcopy___(n, &x[1], 1, &x0[1], 1); + gs = ddot_sl__(n, &g[1], 1, &s[1], 1); + h1 = fabs(gs); + h2 = 0.0; + i__1 = *m; + for (j = 1; j <= i__1; ++j) { + if (j <= *meq) { + h3 = c__[j]; + } else { + h3 = 0.0; + } +/* Computing MAX */ + d__1 = -c__[j]; + h2 += MAX2(d__1,h3); + h3 = (d__1 = r__[j], fabs(d__1)); +/* Computing MAX */ + d__1 = h3, d__2 = (mu[j] + h3) / two; + mu[j] = MAX2(d__1,d__2); + h1 += h3 * (d__1 = c__[j], fabs(d__1)); +/* L170: */ + } +/* CHECK CONVERGENCE */ + *mode = 0; + if (h1 < *acc && h2 < *acc) { + goto L330; + } + h1 = 0.0; + i__1 = *m; + for (j = 1; j <= i__1; ++j) { + if (j <= *meq) { + h3 = c__[j]; + } else { + h3 = 0.0; + } +/* Computing MAX */ + d__1 = -c__[j]; + h1 += mu[j] * MAX2(d__1,h3); +/* L180: */ + } + t0 = *f + h1; + h3 = gs - h1 * h4; + *mode = 8; + if (h3 >= 0.0) { + goto L110; + } +/* LINE SEARCH WITH AN L1-TESTFUNCTION */ + line = 0; + alpha = one; + if (iexact == 1) { + goto L210; + } +/* INEXACT LINESEARCH */ +L190: + ++line; + h3 = alpha * h3; + dscal_sl__(n, &alpha, &s[1], 1); + dcopy___(n, &x0[1], 1, &x[1], 1); + daxpy_sl__(n, &one, &s[1], 1, &x[1], 1); + + /* SGJ 2010: ensure roundoff doesn't push us past bound constraints */ + i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { + if (x[i__] < xl[i__]) x[i__] = xl[i__]; + else if (x[i__] > xu[i__]) x[i__] = xu[i__]; + } + + /* SGJ 2010: optimizing for the common case where the inexact line + search succeeds in one step, use special mode = -2 here to + eliminate a a subsequent unnecessary mode = -1 call, at the + expense of extra gradient evaluations when more than one inexact + line-search step is required */ + *mode = line == 1 ? -2 : 1; + goto L330; +L200: + if (nlopt_isfinite(h1)) { + if (h1 <= h3 / ten || line > 10) { + goto L240; + } + /* Computing MAX */ + d__1 = h3 / (two * (h3 - h1)); + alpha = MAX2(d__1,alfmin); + } else { + alpha = MAX2(alpha*.5,alfmin); + } + goto L190; +/* EXACT LINESEARCH */ +L210: +#if 0 + /* SGJ: see comments by linmin_ above: if we want to do an exact linesearch + (which usually we probably don't), we should call NLopt recursively */ + if (line != 3) { + alpha = linmin_(&line, &alfmin, &one, &t, &tol); + dcopy___(n, &x0[1], 1, &x[1], 1); + daxpy_sl__(n, &alpha, &s[1], 1, &x[1], 1); + *mode = 1; + goto L330; + } +#else + *mode = 9 /* will yield nlopt_failure */; return; +#endif + dscal_sl__(n, &alpha, &s[1], 1); + goto L240; +/* CALL FUNCTIONS AT CURRENT X */ +L220: + t = *f; + i__1 = *m; + for (j = 1; j <= i__1; ++j) { + if (j <= *meq) { + h1 = c__[j]; + } else { + h1 = 0.0; + } +/* Computing MAX */ + d__1 = -c__[j]; + t += mu[j] * MAX2(d__1,h1); +/* L230: */ + } + h1 = t - t0; + switch (iexact + 1) { + case 1: goto L200; + case 2: goto L210; + } +/* CHECK CONVERGENCE */ +L240: + h3 = 0.0; + i__1 = *m; + for (j = 1; j <= i__1; ++j) { + if (j <= *meq) { + h1 = c__[j]; + } else { + h1 = 0.0; + } +/* Computing MAX */ + d__1 = -c__[j]; + h3 += MAX2(d__1,h1); +/* L250: */ + } + if (((d__1 = *f - f0, fabs(d__1)) < *acc || dnrm2___(n, &s[1], 1) < * + acc) && h3 < *acc) { + *mode = 0; + } else { + *mode = -1; + } + goto L330; +/* CHECK relaxed CONVERGENCE in case of positive directional derivative */ +L255: + if (((d__1 = *f - f0, fabs(d__1)) < tol || dnrm2___(n, &s[1], 1) < tol) + && h3 < tol) { + *mode = 0; + } else { + *mode = 8; + } + goto L330; +/* CALL JACOBIAN AT CURRENT X */ +/* UPDATE CHOLESKY-FACTORS OF HESSIAN MATRIX BY MODIFIED BFGS FORMULA */ +L260: + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + u[i__] = g[i__] - ddot_sl__(m, &a[i__ * a_dim1 + 1], 1, &r__[1], 1) - v[i__]; +/* L270: */ + } +/* L'*S */ + k = 0; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + h1 = 0.0; + ++k; + i__2 = *n; + for (j = i__ + 1; j <= i__2; ++j) { + ++k; + h1 += l[k] * s[j]; +/* L280: */ + } + v[i__] = s[i__] + h1; +/* L290: */ + } +/* D*L'*S */ + k = 1; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + v[i__] = l[k] * v[i__]; + k = k + n1 - i__; +/* L300: */ + } +/* L*D*L'*S */ + for (i__ = *n; i__ >= 1; --i__) { + h1 = 0.0; + k = i__; + i__1 = i__ - 1; + for (j = 1; j <= i__1; ++j) { + h1 += l[k] * v[j]; + k = k + *n - j; +/* L310: */ + } + v[i__] += h1; +/* L320: */ + } + h1 = ddot_sl__(n, &s[1], 1, &u[1], 1); + h2 = ddot_sl__(n, &s[1], 1, &v[1], 1); + h3 = h2 * .2; + if (h1 < h3) { + h4 = (h2 - h3) / (h2 - h1); + h1 = h3; + dscal_sl__(n, &h4, &u[1], 1); + d__1 = one - h4; + daxpy_sl__(n, &d__1, &v[1], 1, &u[1], 1); + } + d__1 = one / h1; + ldl_(n, &l[1], &u[1], &d__1, &v[1]); + d__1 = -one / h2; + ldl_(n, &l[1], &v[1], &d__1, &u[1]); +/* END OF MAIN ITERATION */ + goto L130; +/* END OF SLSQPB */ +L330: + SAVE_STATE; +} /* slsqpb_ */ + +/* *********************************************************************** */ +/* optimizer * */ +/* *********************************************************************** */ +void slsqp(int *m, int *meq, int *la, int *n, + double *x, const double *xl, const double *xu, double *f, + double *c__, double *g, double *a, double *acc, + int *iter, int *mode, double *w, int *l_w__, int * + jw, int *l_jw__, slsqpb_state *state) +{ + /* System generated locals */ + int a_dim1, a_offset, i__1, i__2; + + /* Local variables */ + int n1, il, im, ir, is, iu, iv, iw, ix, mineq; + +/* SLSQP S EQUENTIAL L EAST SQ UARES P ROGRAMMING */ +/* TO SOLVE GENERAL NONLINEAR OPTIMIZATION PROBLEMS */ +/* *********************************************************************** */ +/* * * */ +/* * * */ +/* * A NONLINEAR PROGRAMMING METHOD WITH * */ +/* * QUADRATIC PROGRAMMING SUBPROBLEMS * */ +/* * * */ +/* * * */ +/* * THIS SUBROUTINE SOLVES THE GENERAL NONLINEAR PROGRAMMING PROBLEM * */ +/* * * */ +/* * MINIMIZE F(X) * */ +/* * * */ +/* * SUBJECT TO C (X) .EQ. 0 , J = 1,...,MEQ * */ +/* * J * */ +/* * * */ +/* * C (X) .GE. 0 , J = MEQ+1,...,M * */ +/* * J * */ +/* * * */ +/* * XL .LE. X .LE. XU , I = 1,...,N. * */ +/* * I I I * */ +/* * * */ +/* * THE ALGORITHM IMPLEMENTS THE METHOD OF HAN AND POWELL * */ +/* * WITH BFGS-UPDATE OF THE B-MATRIX AND L1-TEST FUNCTION * */ +/* * WITHIN THE STEPLENGTH ALGORITHM. * */ +/* * * */ +/* * PARAMETER DESCRIPTION: * */ +/* * ( * MEANS THIS PARAMETER WILL BE CHANGED DURING CALCULATION ) * */ +/* * * */ +/* * M IS THE TOTAL NUMBER OF CONSTRAINTS, M .GE. 0 * */ +/* * MEQ IS THE NUMBER OF EQUALITY CONSTRAINTS, MEQ .GE. 0 * */ +/* * LA SEE A, LA .GE. MAX(M,1) * */ +/* * N IS THE NUMBER OF VARIABLES, N .GE. 1 * */ +/* * * X() X() STORES THE CURRENT ITERATE OF THE N VECTOR X * */ +/* * ON ENTRY X() MUST BE INITIALIZED. ON EXIT X() * */ +/* * STORES THE SOLUTION VECTOR X IF MODE = 0. * */ +/* * XL() XL() STORES AN N VECTOR OF LOWER BOUNDS XL TO X. * */ +/* * XU() XU() STORES AN N VECTOR OF UPPER BOUNDS XU TO X. * */ +/* * F IS THE VALUE OF THE OBJECTIVE FUNCTION. * */ +/* * C() C() STORES THE M VECTOR C OF CONSTRAINTS, * */ +/* * EQUALITY CONSTRAINTS (IF ANY) FIRST. * */ +/* * DIMENSION OF C MUST BE GREATER OR EQUAL LA, * */ +/* * which must be GREATER OR EQUAL MAX(1,M). * */ +/* * G() G() STORES THE N VECTOR G OF PARTIALS OF THE * */ +/* * OBJECTIVE FUNCTION; DIMENSION OF G MUST BE * */ +/* * GREATER OR EQUAL N+1. * */ +/* * A(),LA,M,N THE LA BY N + 1 ARRAY A() STORES * */ +/* * THE M BY N MATRIX A OF CONSTRAINT NORMALS. * */ +/* * A() HAS FIRST DIMENSIONING PARAMETER LA, * */ +/* * WHICH MUST BE GREATER OR EQUAL MAX(1,M). * */ +/* * F,C,G,A MUST ALL BE SET BY THE USER BEFORE EACH CALL. * */ +/* * * ACC ABS(ACC) CONTROLS THE FINAL ACCURACY. * */ +/* * IF ACC .LT. ZERO AN EXACT LINESEARCH IS PERFORMED,* */ +/* * OTHERWISE AN ARMIJO-TYPE LINESEARCH IS USED. * */ +/* * * ITER PRESCRIBES THE MAXIMUM NUMBER OF ITERATIONS. * */ +/* * ON EXIT ITER INDICATES THE NUMBER OF ITERATIONS. * */ +/* * * MODE MODE CONTROLS CALCULATION: * */ +/* * REVERSE COMMUNICATION IS USED IN THE SENSE THAT * */ +/* * THE PROGRAM IS INITIALIZED BY MODE = 0; THEN IT IS* */ +/* * TO BE CALLED REPEATEDLY BY THE USER UNTIL A RETURN* */ +/* * WITH MODE .NE. IABS(1) TAKES PLACE. * */ +/* * IF MODE = -1 GRADIENTS HAVE TO BE CALCULATED, * */ +/* * WHILE WITH MODE = 1 FUNCTIONS HAVE TO BE CALCULATED */ +/* * MODE MUST NOT BE CHANGED BETWEEN SUBSEQUENT CALLS * */ +/* * OF SQP. * */ +/* * EVALUATION MODES: * */ +/* * MODE = -2,-1: GRADIENT EVALUATION, (G&A) * */ +/* * 0: ON ENTRY: INITIALIZATION, (F,G,C&A) * */ +/* * ON EXIT : REQUIRED ACCURACY FOR SOLUTION OBTAINED * */ +/* * 1: FUNCTION EVALUATION, (F&C) * */ +/* * * */ +/* * FAILURE MODES: * */ +/* * 2: NUMBER OF EQUALITY CONSTRAINTS LARGER THAN N * */ +/* * 3: MORE THAN 3*N ITERATIONS IN LSQ SUBPROBLEM * */ +/* * 4: INEQUALITY CONSTRAINTS INCOMPATIBLE * */ +/* * 5: SINGULAR MATRIX E IN LSQ SUBPROBLEM * */ +/* * 6: SINGULAR MATRIX C IN LSQ SUBPROBLEM * */ +/* * 7: RANK-DEFICIENT EQUALITY CONSTRAINT SUBPROBLEM HFTI* */ +/* * 8: POSITIVE DIRECTIONAL DERIVATIVE FOR LINESEARCH * */ +/* * 9: MORE THAN ITER ITERATIONS IN SQP * */ +/* * >=10: WORKING SPACE W OR JW TOO SMALL, * */ +/* * W SHOULD BE ENLARGED TO L_W=MODE/1000 * */ +/* * JW SHOULD BE ENLARGED TO L_JW=MODE-1000*L_W * */ +/* * * W(), L_W W() IS A ONE DIMENSIONAL WORKING SPACE, * */ +/* * THE LENGTH L_W OF WHICH SHOULD BE AT LEAST * */ +/* * (3*N1+M)*(N1+1) for LSQ * */ +/* * +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ for LSI * */ +/* * +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 for LSEI * */ +/* * + N1*N/2 + 2*M + 3*N + 3*N1 + 1 for SLSQPB * */ +/* * with MINEQ = M - MEQ + 2*N1 & N1 = N+1 * */ +/* * NOTICE: FOR PROPER DIMENSIONING OF W IT IS RECOMMENDED TO * */ +/* * COPY THE FOLLOWING STATEMENTS INTO THE HEAD OF * */ +/* * THE CALLING PROGRAM (AND REMOVE THE COMMENT C) * */ +/* ####################################################################### */ +/* INT LEN_W, LEN_JW, M, N, N1, MEQ, MINEQ */ +/* PARAMETER (M=... , MEQ=... , N=... ) */ +/* PARAMETER (N1= N+1, MINEQ= M-MEQ+N1+N1) */ +/* PARAMETER (LEN_W= */ +/* $ (3*N1+M)*(N1+1) */ +/* $ +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ */ +/* $ +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 */ +/* $ +(N+1)*N/2 + 2*M + 3*N + 3*N1 + 1, */ +/* $ LEN_JW=MINEQ) */ +/* DOUBLE PRECISION W(LEN_W) */ +/* INT JW(LEN_JW) */ +/* ####################################################################### */ +/* * THE FIRST M+N+N*N1/2 ELEMENTS OF W MUST NOT BE * */ +/* * CHANGED BETWEEN SUBSEQUENT CALLS OF SLSQP. * */ +/* * ON RETURN W(1) ... W(M) CONTAIN THE MULTIPLIERS * */ +/* * ASSOCIATED WITH THE GENERAL CONSTRAINTS, WHILE * */ +/* * W(M+1) ... W(M+N(N+1)/2) STORE THE CHOLESKY FACTOR* */ +/* * L*D*L(T) OF THE APPROXIMATE HESSIAN OF THE * */ +/* * LAGRANGIAN COLUMNWISE DENSE AS LOWER TRIANGULAR * */ +/* * UNIT MATRIX L WITH D IN ITS 'DIAGONAL' and * */ +/* * W(M+N(N+1)/2+N+2 ... W(M+N(N+1)/2+N+2+M+2N) * */ +/* * CONTAIN THE MULTIPLIERS ASSOCIATED WITH ALL * */ +/* * ALL CONSTRAINTS OF THE QUADRATIC PROGRAM FINDING * */ +/* * THE SEARCH DIRECTION TO THE SOLUTION X* * */ +/* * * JW(), L_JW JW() IS A ONE DIMENSIONAL INT WORKING SPACE * */ +/* * THE LENGTH L_JW OF WHICH SHOULD BE AT LEAST * */ +/* * MINEQ * */ +/* * with MINEQ = M - MEQ + 2*N1 & N1 = N+1 * */ +/* * * */ +/* * THE USER HAS TO PROVIDE THE FOLLOWING SUBROUTINES: * */ +/* * LDL(N,A,Z,SIG,W) : UPDATE OF THE LDL'-FACTORIZATION. * */ +/* * LINMIN(A,B,F,TOL) : LINESEARCH ALGORITHM IF EXACT = 1 * */ +/* * LSQ(M,MEQ,LA,N,NC,C,D,A,B,XL,XU,X,LAMBDA,W,....) : * */ +/* * * */ +/* * SOLUTION OF THE QUADRATIC PROGRAM * */ +/* * QPSOL IS RECOMMENDED: * */ +/* * PE GILL, W MURRAY, MA SAUNDERS, MH WRIGHT: * */ +/* * USER'S GUIDE FOR SOL/QPSOL: * */ +/* * A FORTRAN PACKAGE FOR QUADRATIC PROGRAMMING, * */ +/* * TECHNICAL REPORT SOL 83-7, JULY 1983 * */ +/* * DEPARTMENT OF OPERATIONS RESEARCH, STANFORD UNIVERSITY * */ +/* * STANFORD, CA 94305 * */ +/* * QPSOL IS THE MOST ROBUST AND EFFICIENT QP-SOLVER * */ +/* * AS IT ALLOWS WARM STARTS WITH PROPER WORKING SETS * */ +/* * * */ +/* * IF IT IS NOT AVAILABLE USE LSEI, A CONSTRAINT LINEAR LEAST * */ +/* * SQUARES SOLVER IMPLEMENTED USING THE SOFTWARE HFTI, LDP, NNLS * */ +/* * FROM C.L. LAWSON, R.J.HANSON: SOLVING LEAST SQUARES PROBLEMS, * */ +/* * PRENTICE HALL, ENGLEWOOD CLIFFS, 1974. * */ +/* * LSEI COMES WITH THIS PACKAGE, together with all necessary SR's. * */ +/* * * */ +/* * TOGETHER WITH A COUPLE OF SUBROUTINES FROM BLAS LEVEL 1 * */ +/* * * */ +/* * SQP IS HEAD SUBROUTINE FOR BODY SUBROUTINE SQPBDY * */ +/* * IN WHICH THE ALGORITHM HAS BEEN IMPLEMENTED. * */ +/* * * */ +/* * IMPLEMENTED BY: DIETER KRAFT, DFVLR OBERPFAFFENHOFEN * */ +/* * as described in Dieter Kraft: A Software Package for * */ +/* * Sequential Quadratic Programming * */ +/* * DFVLR-FB 88-28, 1988 * */ +/* * which should be referenced if the user publishes results of SLSQP * */ +/* * * */ +/* * DATE: APRIL - OCTOBER, 1981. * */ +/* * STATUS: DECEMBER, 31-ST, 1984. * */ +/* * STATUS: MARCH , 21-ST, 1987, REVISED TO FORTRAN 77 * */ +/* * STATUS: MARCH , 20-th, 1989, REVISED TO MS-FORTRAN * */ +/* * STATUS: APRIL , 14-th, 1989, HESSE in-line coded * */ +/* * STATUS: FEBRUARY, 28-th, 1991, FORTRAN/2 Version 1.04 * */ +/* * accepts Statement Functions * */ +/* * STATUS: MARCH , 1-st, 1991, tested with SALFORD * */ +/* * FTN77/386 COMPILER VERS 2.40* */ +/* * in protected mode * */ +/* * * */ +/* *********************************************************************** */ +/* * * */ +/* * Copyright 1991: Dieter Kraft, FHM * */ +/* * * */ +/* *********************************************************************** */ +/* dim(W) = N1*(N1+1) + MEQ*(N1+1) + MINEQ*(N1+1) for LSQ */ +/* +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ for LSI */ +/* +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 for LSEI */ +/* + N1*N/2 + 2*M + 3*N +3*N1 + 1 for SLSQPB */ +/* with MINEQ = M - MEQ + 2*N1 & N1 = N+1 */ +/* CHECK LENGTH OF WORKING ARRAYS */ + /* Parameter adjustments */ + --c__; + a_dim1 = *la; + a_offset = 1 + a_dim1; + a -= a_offset; + --g; + --xu; + --xl; + --x; + --w; + --jw; + + /* Function Body */ + n1 = *n + 1; + mineq = *m - *meq + n1 + n1; + il = (n1 * 3 + *m) * (n1 + 1) + (n1 - *meq + 1) * (mineq + 2) + (mineq << + 1) + (n1 + mineq) * (n1 - *meq) + (*meq << 1) + n1 * *n / 2 + (*m + << 1) + *n * 3 + (n1 << 2) + 1; +/* Computing MAX */ + i__1 = mineq, i__2 = n1 - *meq; + im = MAX2(i__1,i__2); + if (*l_w__ < il || *l_jw__ < im) { + *mode = MAX2(10,il) * 1000; + *mode += MAX2(10,im); + return; + } +/* PREPARE DATA FOR CALLING SQPBDY - INITIAL ADDRESSES IN W */ + im = 1; + il = im + MAX2(1,*m); + il = im + *la; + ix = il + n1 * *n / 2 + 1; + ir = ix + *n; + is = ir + *n + *n + MAX2(1,*m); + is = ir + *n + *n + *la; + iu = is + n1; + iv = iu + n1; + iw = iv + n1; + slsqpb_(m, meq, la, n, &x[1], &xl[1], &xu[1], f, &c__[1], &g[1], &a[ + a_offset], acc, iter, mode, &w[ir], &w[il], &w[ix], &w[im], &w[is] + , &w[iu], &w[iv], &w[iw], &jw[1], state); + state->x0 = &w[ix]; + return; +} /* slsqp_ */ + +static void length_work(int *LEN_W, int *LEN_JW, int M, int MEQ, int N) +{ + int N1 = N+1, MINEQ = M-MEQ+N1+N1; + *LEN_W = (3*N1+M)*(N1+1) + +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ + +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 + +(N+1)*N/2 + 2*M + 3*N + 3*N1 + 1; + *LEN_JW = MINEQ; +} + +// nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data, +// unsigned m, nlopt_constraint *fc, +// unsigned p, nlopt_constraint *h, +// const double *lb, const double *ub, +// double *x, double *minf, +// nlopt_stopping *stop) +// { +// slsqpb_state state = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,NULL}; +// unsigned mtot = nlopt_count_constraints(m, fc); +// unsigned ptot = nlopt_count_constraints(p, h); +// double *work, *cgrad, *c, *grad, *w, +// fcur, *xcur, fprev, *xprev, *cgradtmp; +// int mpi = (int) (mtot + ptot), pi = (int) ptot, ni = (int) n, mpi1 = mpi > 0 ? mpi : 1; +// int len_w, len_jw, *jw; +// int mode = 0, prev_mode = 0; +// double acc = 0; /* we do our own convergence tests below */ +// int iter = 0; /* tell sqsqp to ignore this check, since we check evaluation counts ourselves */ +// unsigned i, ii; +// nlopt_result ret = NLOPT_SUCCESS; +// int feasible, feasible_cur; +// double infeasibility = HUGE_VAL, infeasibility_cur = HUGE_VAL; +// unsigned max_cdim; +// int want_grad = 1; + +// max_cdim = MAX2(nlopt_max_constraint_dim(m, fc), +// nlopt_max_constraint_dim(p, h)); +// length_work(&len_w, &len_jw, mpi, pi, ni); + +// #define U(n) ((unsigned) (n)) +// work = (double *) malloc(sizeof(double) * (U(mpi1) * (n + 1) +// + U(mpi) +// + n+1 + n + n + max_cdim*n +// + U(len_w)) +// + sizeof(int) * U(len_jw)); +// if (!work) return NLOPT_OUT_OF_MEMORY; +// cgrad = work; +// c = cgrad + U(mpi1) * (n + 1); +// grad = c + mpi; +// xcur = grad + n+1; +// xprev = xcur + n; +// cgradtmp = xprev + n; +// w = cgradtmp + max_cdim*n; +// jw = (int *) (w + len_w); + +// memcpy(xcur, x, sizeof(double) * n); +// memcpy(xprev, x, sizeof(double) * n); +// fprev = fcur = *minf = HUGE_VAL; +// feasible = feasible_cur = 0; + +// goto eval_f_and_grad; /* eval before calling slsqp the first time */ + +// do { +// slsqp(&mpi, &pi, &mpi1, &ni, +// xcur, lb, ub, &fcur, +// c, grad, cgrad, +// &acc, &iter, &mode, +// w, &len_w, jw, &len_jw, +// &state); + +// switch (mode) { +// case -1: /* objective & gradient evaluation */ +// if (prev_mode == -2 && !want_grad) break; /* just evaluated this point */ +// /* fall through */ +// case -2: +// eval_f_and_grad: +// want_grad = 1; +// /* fall through */ +// case 1:{ /* don't need grad unless we don't have it yet */ +// double *newgrad = 0; +// double *newcgrad = 0; +// if (want_grad) { +// newgrad = grad; +// newcgrad = cgradtmp; +// } +// feasible_cur = 1; infeasibility_cur = 0; +// fcur = f(n, xcur, newgrad, f_data); +// ++ *(stop->nevals_p); +// if (nlopt_stop_forced(stop)) { +// fcur = HUGE_VAL; ret = NLOPT_FORCED_STOP; goto done; } +// if (nlopt_isfinite(fcur)) { +// want_grad = 0; +// ii = 0; +// for (i = 0; i < p; ++i) { +// unsigned j, k; +// nlopt_eval_constraint(c+ii, newcgrad, h+i, n, xcur); +// if (nlopt_stop_forced(stop)) { +// ret = NLOPT_FORCED_STOP; goto done; } +// for (k = 0; k < h[i].m; ++k, ++ii) { +// infeasibility_cur = +// MAX2(infeasibility_cur, fabs(c[ii])); +// feasible_cur = +// feasible_cur && fabs(c[ii]) <= h[i].tol[k]; +// if (newcgrad) { +// for (j = 0; j < n; ++ j) +// cgrad[j*U(mpi1) + ii] = cgradtmp[k*n + j]; +// } +// } +// } +// for (i = 0; i < m; ++i) { +// unsigned j, k; +// nlopt_eval_constraint(c+ii, newcgrad, fc+i, n, xcur); +// if (nlopt_stop_forced(stop)) { +// ret = NLOPT_FORCED_STOP; goto done; } +// for (k = 0; k < fc[i].m; ++k, ++ii) { +// infeasibility_cur = +// MAX2(infeasibility_cur, c[ii]); +// feasible_cur = +// feasible_cur && c[ii] <= fc[i].tol[k]; +// if (newcgrad) { +// for (j = 0; j < n; ++ j) +// cgrad[j*U(mpi1) + ii] = -cgradtmp[k*n + j]; +// } +// c[ii] = -c[ii]; /* slsqp sign convention */ +// } +// } +// } +// break;} +// case 0: /* required accuracy for solution obtained */ +// goto done; +// case 8: /* positive directional derivative for linesearch */ +// /* relaxed convergence check for a feasible_cur point, +// as in the SLSQP code (except xtol as well as ftol) */ +// ret = NLOPT_ROUNDOFF_LIMITED; /* usually why deriv>0 */ +// if (feasible_cur) { +// double save_ftol_rel = stop->ftol_rel; +// double save_xtol_rel = stop->xtol_rel; +// double save_ftol_abs = stop->ftol_abs; +// stop->ftol_rel *= 10; +// stop->ftol_abs *= 10; +// stop->xtol_rel *= 10; +// if (nlopt_stop_ftol(stop, fcur, state.f0)) +// ret = NLOPT_FTOL_REACHED; +// else if (nlopt_stop_x(stop, xcur, state.x0)) +// ret = NLOPT_XTOL_REACHED; +// stop->ftol_rel = save_ftol_rel; +// stop->ftol_abs = save_ftol_abs; +// stop->xtol_rel = save_xtol_rel; +// } +// goto done; +// case 5: /* singular matrix E in LSQ subproblem */ +// case 6: /* singular matrix C in LSQ subproblem */ +// case 7: /* rank-deficient equality constraint subproblem HFTI */ +// ret = NLOPT_ROUNDOFF_LIMITED; +// goto done; +// case 4: /* inequality constraints incompatible */ +// case 3: /* more than 3*n iterations in LSQ subproblem */ +// case 9: /* more than iter iterations in SQP */ +// nlopt_stop_msg(stop, "bug: more than iter SQP iterations"); +// ret = NLOPT_FAILURE; +// goto done; +// case 2: /* number of equality constraints > n */ +// default: /* >= 10: working space w or jw too small */ +// nlopt_stop_msg(stop, "bug: workspace is too small"); +// ret = NLOPT_INVALID_ARGS; +// goto done; +// } +// prev_mode = mode; + +// /* update best point so far */ +// if (nlopt_isfinite(fcur) && ((fcur < *minf && (feasible_cur || !feasible)) +// || (!feasible && infeasibility_cur < infeasibility))) { +// *minf = fcur; +// feasible = feasible_cur; +// infeasibility = infeasibility_cur; +// memcpy(x, xcur, sizeof(double)*n); +// } + +// /* note: mode == -1 corresponds to the completion of a line search, +// and is the only time we should check convergence (as in original slsqp code) */ +// if (mode == -1) { +// if (!nlopt_isinf(fprev) && feasible) { +// if (nlopt_stop_ftol(stop, fcur, fprev)) +// ret = NLOPT_FTOL_REACHED; +// else if (nlopt_stop_x(stop, xcur, xprev)) +// ret = NLOPT_XTOL_REACHED; +// } +// fprev = fcur; +// memcpy(xprev, xcur, sizeof(double)*n); +// } + +// /* do some additional termination tests */ +// if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED; +// else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED; +// else if (feasible && *minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED; +// } while (ret == NLOPT_SUCCESS); + +// done: +// if (nlopt_isinf(*minf)) { /* didn't find any feasible points, just return last point evaluated */ +// if (nlopt_isinf(fcur)) { /* invalid cur. point, use previous pt. */ +// *minf = fprev; +// memcpy(x, xprev, sizeof(double)*n); +// } +// else { +// *minf = fcur; +// memcpy(x, xcur, sizeof(double)*n); +// } +// } + +// free(work); +// return ret; +// } diff --git a/slsqp/slsqp.h b/slsqp/slsqp.h new file mode 100644 index 0000000..c85b853 --- /dev/null +++ b/slsqp/slsqp.h @@ -0,0 +1,62 @@ +#ifndef SLSQP_H +#define SLSQP_H + +#include +#include +// #include "nlopt.h" +// #include "nlopt-util.h" + +// #ifdef __cplusplus +// extern "C" +// { +// #endif /* __cplusplus */ + +// nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data, +// unsigned m, nlopt_constraint *fc, +// unsigned p, nlopt_constraint *h, +// const double *lb, const double *ub, +// double *x, double *minf, +// nlopt_stopping *stop); +// #ifdef __cplusplus +// } /* extern "C" */ +// #endif /* __cplusplus */ + + +int nlopt_isnan(double x) +{ + return isnan(x); +} + +int nlopt_isinf(double x) +{ + return (fabs(x) >= HUGE_VAL * 0.99) + || isinf(x) + ; +} + +int nlopt_isfinite(double x) +{ + return (fabs(x) <= DBL_MAX) + || isfinite(x) + ; +} + +typedef struct { + double t, f0, h1, h2, h3, h4; + int n1, n2, n3; + double t0, gs; + double tol; + int line; + double alpha; + int iexact; + int incons, ireset, itermx; + double *x0; +} slsqpb_state; + +void slsqp(int *m, int *meq, int *la, int *n, + double *x, const double *xl, const double *xu, double *f, + double *c__, double *g, double *a, double *acc, + int *iter, int *mode, double *w, int *l_w__, int * + jw, int *l_jw__, slsqpb_state *state); + +#endif diff --git a/slsqp/slsqp.rs b/slsqp/slsqp.rs new file mode 100644 index 0000000..a58fb11 --- /dev/null +++ b/slsqp/slsqp.rs @@ -0,0 +1,2868 @@ +#![allow(dead_code, mutable_transmutes, non_camel_case_types, non_snake_case, non_upper_case_globals, unused_assignments, unused_mut)] +#![register_tool(c2rust)] +#![feature(register_tool)] +extern "C" { + fn sqrt(_: libc::c_double) -> libc::c_double; + fn fabs(_: libc::c_double) -> libc::c_double; + fn memcpy( + _: *mut libc::c_void, + _: *const libc::c_void, + _: libc::c_ulong, + ) -> *mut libc::c_void; +} +#[derive(Copy, Clone)] +#[repr(C)] +pub struct slsqpb_state { + pub t: libc::c_double, + pub f0: libc::c_double, + pub h1: libc::c_double, + pub h2: libc::c_double, + pub h3: libc::c_double, + pub h4: libc::c_double, + pub n1: libc::c_int, + pub n2: libc::c_int, + pub n3: libc::c_int, + pub t0: libc::c_double, + pub gs: libc::c_double, + pub tol: libc::c_double, + pub line: libc::c_int, + pub alpha: libc::c_double, + pub iexact: libc::c_int, + pub incons: libc::c_int, + pub ireset: libc::c_int, + pub itermx: libc::c_int, + pub x0: *mut libc::c_double, +} +#[no_mangle] +pub unsafe extern "C" fn nlopt_isnan(mut x: libc::c_double) -> libc::c_int { + return x.is_nan() as i32; +} +#[no_mangle] +pub unsafe extern "C" fn nlopt_isinf(mut x: libc::c_double) -> libc::c_int { + return (fabs(x) >= ::std::f64::INFINITY * 0.99f64 + || if x.is_infinite() { if x.is_sign_positive() { 1 } else { -1 } } else { 0 } + != 0) as libc::c_int; +} +#[no_mangle] +pub unsafe extern "C" fn nlopt_isfinite(mut x: libc::c_double) -> libc::c_int { + return (fabs(x) <= 1.7976931348623157e+308f64 || x.is_finite() as i32 != 0) + as libc::c_int; +} +unsafe extern "C" fn dcopy___( + mut n_: *mut libc::c_int, + mut dx: *const libc::c_double, + mut incx: libc::c_int, + mut dy: *mut libc::c_double, + mut incy: libc::c_int, +) { + let mut i: libc::c_int = 0; + let mut n: libc::c_int = *n_; + if n <= 0 as libc::c_int { + return; + } + if incx == 1 as libc::c_int && incy == 1 as libc::c_int { + memcpy( + dy as *mut libc::c_void, + dx as *const libc::c_void, + (::std::mem::size_of::() as libc::c_ulong) + .wrapping_mul(n as libc::c_uint as libc::c_ulong), + ); + } else if incx == 0 as libc::c_int && incy == 1 as libc::c_int { + let mut x: libc::c_double = *dx.offset(0 as libc::c_int as isize); + i = 0 as libc::c_int; + while i < n { + *dy.offset(i as isize) = x; + i += 1; + } + } else { + i = 0 as libc::c_int; + while i < n { + *dy.offset((i * incy) as isize) = *dx.offset((i * incx) as isize); + i += 1; + } + }; +} +unsafe extern "C" fn daxpy_sl__( + mut n_: *mut libc::c_int, + mut da_: *const libc::c_double, + mut dx: *const libc::c_double, + mut incx: libc::c_int, + mut dy: *mut libc::c_double, + mut incy: libc::c_int, +) { + let mut n: libc::c_int = *n_; + let mut i: libc::c_int = 0; + let mut da: libc::c_double = *da_; + if n <= 0 as libc::c_int || da == 0 as libc::c_int as libc::c_double { + return; + } + i = 0 as libc::c_int; + while i < n { + *dy.offset((i * incy) as isize) += da * *dx.offset((i * incx) as isize); + i += 1; + } +} +unsafe extern "C" fn ddot_sl__( + mut n_: *mut libc::c_int, + mut dx: *mut libc::c_double, + mut incx: libc::c_int, + mut dy: *mut libc::c_double, + mut incy: libc::c_int, +) -> libc::c_double { + let mut n: libc::c_int = *n_; + let mut i: libc::c_int = 0; + let mut sum: libc::c_double = 0 as libc::c_int as libc::c_double; + if n <= 0 as libc::c_int { + return 0 as libc::c_int as libc::c_double; + } + i = 0 as libc::c_int; + while i < n { + sum += *dx.offset((i * incx) as isize) * *dy.offset((i * incy) as isize); + i += 1; + } + return sum; +} +unsafe extern "C" fn dnrm2___( + mut n_: *mut libc::c_int, + mut dx: *mut libc::c_double, + mut incx: libc::c_int, +) -> libc::c_double { + let mut i: libc::c_int = 0; + let mut n: libc::c_int = *n_; + let mut xmax: libc::c_double = 0 as libc::c_int as libc::c_double; + let mut scale: libc::c_double = 0.; + let mut sum: libc::c_double = 0 as libc::c_int as libc::c_double; + i = 0 as libc::c_int; + while i < n { + let mut xabs: libc::c_double = fabs(*dx.offset((incx * i) as isize)); + if xmax < xabs { + xmax = xabs; + } + i += 1; + } + if xmax == 0 as libc::c_int as libc::c_double { + return 0 as libc::c_int as libc::c_double; + } + scale = 1.0f64 / xmax; + i = 0 as libc::c_int; + while i < n { + let mut xs: libc::c_double = scale * *dx.offset((incx * i) as isize); + sum += xs * xs; + i += 1; + } + return xmax * sqrt(sum); +} +unsafe extern "C" fn dsrot_( + mut n: libc::c_int, + mut dx: *mut libc::c_double, + mut incx: libc::c_int, + mut dy: *mut libc::c_double, + mut incy: libc::c_int, + mut c__: *mut libc::c_double, + mut s_: *mut libc::c_double, +) { + let mut i: libc::c_int = 0; + let mut c: libc::c_double = *c__; + let mut s: libc::c_double = *s_; + i = 0 as libc::c_int; + while i < n { + let mut x: libc::c_double = *dx.offset((incx * i) as isize); + let mut y: libc::c_double = *dy.offset((incy * i) as isize); + *dx.offset((incx * i) as isize) = c * x + s * y; + *dy.offset((incy * i) as isize) = c * y - s * x; + i += 1; + } +} +unsafe extern "C" fn dsrotg_( + mut da: *mut libc::c_double, + mut db: *mut libc::c_double, + mut c: *mut libc::c_double, + mut s: *mut libc::c_double, +) { + let mut absa: libc::c_double = 0.; + let mut absb: libc::c_double = 0.; + let mut roe: libc::c_double = 0.; + let mut scale: libc::c_double = 0.; + absa = fabs(*da); + absb = fabs(*db); + if absa > absb { + roe = *da; + scale = absa; + } else { + roe = *db; + scale = absb; + } + if scale != 0 as libc::c_int as libc::c_double { + let mut r: libc::c_double = 0.; + let mut iscale: libc::c_double = 1 as libc::c_int as libc::c_double / scale; + let mut tmpa: libc::c_double = *da * iscale; + let mut tmpb: libc::c_double = *db * iscale; + r = (if roe < 0 as libc::c_int as libc::c_double { -scale } else { scale }) + * sqrt(tmpa * tmpa + tmpb * tmpb); + *c = *da / r; + *s = *db / r; + *da = r; + if *c != 0 as libc::c_int as libc::c_double && fabs(*c) <= *s { + *db = 1 as libc::c_int as libc::c_double / *c; + } else { + *db = *s; + } + } else { + *c = 1 as libc::c_int as libc::c_double; + *db = 0 as libc::c_int as libc::c_double; + *da = *db; + *s = *da; + }; +} +unsafe extern "C" fn dscal_sl__( + mut n_: *mut libc::c_int, + mut da: *const libc::c_double, + mut dx: *mut libc::c_double, + mut incx: libc::c_int, +) { + let mut i: libc::c_int = 0; + let mut n: libc::c_int = *n_; + let mut alpha: libc::c_double = *da; + i = 0 as libc::c_int; + while i < n { + *dx.offset((i * incx) as isize) *= alpha; + i += 1; + } +} +static mut c__0: libc::c_int = 0 as libc::c_int; +static mut c__1: libc::c_int = 1 as libc::c_int; +static mut c__2: libc::c_int = 2 as libc::c_int; +unsafe extern "C" fn h12_( + mut mode: *const libc::c_int, + mut lpivot: *mut libc::c_int, + mut l1: *mut libc::c_int, + mut m: *mut libc::c_int, + mut u: *mut libc::c_double, + mut iue: *const libc::c_int, + mut up: *mut libc::c_double, + mut c__: *mut libc::c_double, + mut ice: *const libc::c_int, + mut icv: *const libc::c_int, + mut ncv: *const libc::c_int, +) { + let mut current_block: u64; + let one: libc::c_double = 1.0f64; + let mut u_dim1: libc::c_int = 0; + let mut u_offset: libc::c_int = 0; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut d__1: libc::c_double = 0.; + let mut b: libc::c_double = 0.; + let mut i__: libc::c_int = 0; + let mut j: libc::c_int = 0; + let mut i2: libc::c_int = 0; + let mut i3: libc::c_int = 0; + let mut i4: libc::c_int = 0; + let mut cl: libc::c_double = 0.; + let mut sm: libc::c_double = 0.; + let mut incr: libc::c_int = 0; + let mut clinv: libc::c_double = 0.; + u_dim1 = *iue; + u_offset = 1 as libc::c_int + u_dim1; + u = u.offset(-(u_offset as isize)); + c__ = c__.offset(-1); + if !(0 as libc::c_int >= *lpivot || *lpivot >= *l1 || *l1 > *m) { + d__1 = *u.offset((*lpivot * u_dim1 + 1 as libc::c_int) as isize); + cl = fabs(d__1); + if *mode == 2 as libc::c_int { + if cl <= 0.0f64 { + current_block = 5221028069996397600; + } else { + current_block = 15706997569754958587; + } + } else { + i__1 = *m; + j = *l1; + while j <= i__1 { + d__1 = *u.offset((j * u_dim1 + 1 as libc::c_int) as isize); + sm = fabs(d__1); + cl = if sm >= cl { sm } else { cl }; + j += 1; + } + if cl <= 0.0f64 { + current_block = 5221028069996397600; + } else { + clinv = one / cl; + d__1 = *u.offset((*lpivot * u_dim1 + 1 as libc::c_int) as isize) * clinv; + sm = d__1 * d__1; + i__1 = *m; + j = *l1; + while j <= i__1 { + d__1 = *u.offset((j * u_dim1 + 1 as libc::c_int) as isize) * clinv; + sm += d__1 * d__1; + j += 1; + } + cl *= sqrt(sm); + if *u.offset((*lpivot * u_dim1 + 1 as libc::c_int) as isize) > 0.0f64 { + cl = -cl; + } + *up = *u.offset((*lpivot * u_dim1 + 1 as libc::c_int) as isize) - cl; + *u.offset((*lpivot * u_dim1 + 1 as libc::c_int) as isize) = cl; + current_block = 15706997569754958587; + } + } + match current_block { + 5221028069996397600 => {} + _ => { + if !(*ncv <= 0 as libc::c_int) { + b = *up * *u.offset((*lpivot * u_dim1 + 1 as libc::c_int) as isize); + if !(b >= 0.0f64) { + b = one / b; + i2 = 1 as libc::c_int - *icv + + *ice * (*lpivot - 1 as libc::c_int); + incr = *ice * (*l1 - *lpivot); + i__1 = *ncv; + j = 1 as libc::c_int; + while j <= i__1 { + i2 += *icv; + i3 = i2 + incr; + i4 = i3; + sm = *c__.offset(i2 as isize) * *up; + i__2 = *m; + i__ = *l1; + while i__ <= i__2 { + sm + += *c__.offset(i3 as isize) + * *u.offset((i__ * u_dim1 + 1 as libc::c_int) as isize); + i3 += *ice; + i__ += 1; + } + if !(sm == 0.0f64) { + sm *= b; + *c__.offset(i2 as isize) += sm * *up; + i__2 = *m; + i__ = *l1; + while i__ <= i__2 { + *c__.offset(i4 as isize) + += sm + * *u.offset((i__ * u_dim1 + 1 as libc::c_int) as isize); + i4 += *ice; + i__ += 1; + } + } + j += 1; + } + } + } + } + } + } +} +unsafe extern "C" fn nnls_( + mut a: *mut libc::c_double, + mut mda: *mut libc::c_int, + mut m: *mut libc::c_int, + mut n: *mut libc::c_int, + mut b: *mut libc::c_double, + mut x: *mut libc::c_double, + mut rnorm: *mut libc::c_double, + mut w: *mut libc::c_double, + mut z__: *mut libc::c_double, + mut indx: *mut libc::c_int, + mut mode: *mut libc::c_int, +) { + let mut current_block: u64; + let one: libc::c_double = 1.0f64; + let factor: libc::c_double = 0.01f64; + let mut a_dim1: libc::c_int = 0; + let mut a_offset: libc::c_int = 0; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut d__1: libc::c_double = 0.; + let mut c__: libc::c_double = 0.; + let mut i__: libc::c_int = 0; + let mut j: libc::c_int = 0; + let mut k: libc::c_int = 0; + let mut l: libc::c_int = 0; + let mut s: libc::c_double = 0.; + let mut t: libc::c_double = 0.; + let mut ii: libc::c_int = 0; + let mut jj: libc::c_int = 0; + let mut ip: libc::c_int = 0; + let mut iz: libc::c_int = 0; + let mut jz: libc::c_int = 0; + let mut up: libc::c_double = 0.; + let mut iz1: libc::c_int = 0; + let mut iz2: libc::c_int = 0; + let mut npp1: libc::c_int = 0; + let mut iter: libc::c_int = 0; + let mut wmax: libc::c_double = 0.; + let mut alpha: libc::c_double = 0.; + let mut asave: libc::c_double = 0.; + let mut itmax: libc::c_int = 0; + let mut izmax: libc::c_int = 0 as libc::c_int; + let mut nsetp: libc::c_int = 0; + let mut unorm: libc::c_double = 0.; + z__ = z__.offset(-1); + b = b.offset(-1); + indx = indx.offset(-1); + w = w.offset(-1); + x = x.offset(-1); + a_dim1 = *mda; + a_offset = 1 as libc::c_int + a_dim1; + a = a.offset(-(a_offset as isize)); + *mode = 2 as libc::c_int; + if !(*m <= 0 as libc::c_int || *n <= 0 as libc::c_int) { + *mode = 1 as libc::c_int; + iter = 0 as libc::c_int; + itmax = *n * 3 as libc::c_int; + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + *indx.offset(i__ as isize) = i__; + i__ += 1; + } + iz1 = 1 as libc::c_int; + iz2 = *n; + nsetp = 0 as libc::c_int; + npp1 = 1 as libc::c_int; + *x.offset(1 as libc::c_int as isize) = 0.0f64; + dcopy___( + n, + &mut *x.offset(1 as libc::c_int as isize), + 0 as libc::c_int, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + 'c_7799: loop { + if iz1 > iz2 || nsetp >= *m { + current_block = 12767425942910239388; + break; + } + i__1 = iz2; + iz = iz1; + while iz <= i__1 { + j = *indx.offset(iz as isize); + i__2 = *m - nsetp; + *w + .offset( + j as isize, + ) = ddot_sl__( + &mut i__2, + &mut *a.offset((npp1 + j * a_dim1) as isize), + 1 as libc::c_int, + &mut *b.offset(npp1 as isize), + 1 as libc::c_int, + ); + iz += 1; + } + loop { + wmax = 0.0f64; + i__2 = iz2; + iz = iz1; + while iz <= i__2 { + j = *indx.offset(iz as isize); + if !(*w.offset(j as isize) <= wmax) { + wmax = *w.offset(j as isize); + izmax = iz; + } + iz += 1; + } + if wmax <= 0.0f64 { + current_block = 12767425942910239388; + break 'c_7799; + } + iz = izmax; + j = *indx.offset(iz as isize); + asave = *a.offset((npp1 + j * a_dim1) as isize); + i__2 = npp1 + 1 as libc::c_int; + h12_( + &c__1, + &mut npp1, + &mut i__2, + m, + &mut *a.offset((j * a_dim1 + 1 as libc::c_int) as isize), + &c__1, + &mut up, + &mut *z__.offset(1 as libc::c_int as isize), + &c__1, + &c__1, + &c__0, + ); + unorm = dnrm2___( + &mut nsetp, + &mut *a.offset((j * a_dim1 + 1 as libc::c_int) as isize), + 1 as libc::c_int, + ); + d__1 = *a.offset((npp1 + j * a_dim1) as isize); + t = factor * fabs(d__1); + d__1 = unorm + t; + if !(d__1 - unorm <= 0.0f64) { + dcopy___( + m, + &mut *b.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *z__.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + i__2 = npp1 + 1 as libc::c_int; + h12_( + &c__2, + &mut npp1, + &mut i__2, + m, + &mut *a.offset((j * a_dim1 + 1 as libc::c_int) as isize), + &c__1, + &mut up, + &mut *z__.offset(1 as libc::c_int as isize), + &c__1, + &c__1, + &c__1, + ); + if *z__.offset(npp1 as isize) + / *a.offset((npp1 + j * a_dim1) as isize) > 0.0f64 + { + break; + } + } + *a.offset((npp1 + j * a_dim1) as isize) = asave; + *w.offset(j as isize) = 0.0f64; + } + dcopy___( + m, + &mut *z__.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *b.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + *indx.offset(iz as isize) = *indx.offset(iz1 as isize); + *indx.offset(iz1 as isize) = j; + iz1 += 1; + nsetp = npp1; + npp1 += 1; + i__2 = iz2; + jz = iz1; + while jz <= i__2 { + jj = *indx.offset(jz as isize); + h12_( + &c__2, + &mut nsetp, + &mut npp1, + m, + &mut *a.offset((j * a_dim1 + 1 as libc::c_int) as isize), + &c__1, + &mut up, + &mut *a.offset((jj * a_dim1 + 1 as libc::c_int) as isize), + &c__1, + mda, + &c__1, + ); + jz += 1; + } + k = if npp1 <= *mda { npp1 } else { *mda }; + *w.offset(j as isize) = 0.0f64; + i__2 = *m - nsetp; + dcopy___( + &mut i__2, + &mut *w.offset(j as isize), + 0 as libc::c_int, + &mut *a.offset((k + j * a_dim1) as isize), + 1 as libc::c_int, + ); + loop { + ip = nsetp; + while ip >= 1 as libc::c_int { + if !(ip == nsetp) { + d__1 = -*z__.offset((ip + 1 as libc::c_int) as isize); + daxpy_sl__( + &mut ip, + &mut d__1, + &mut *a.offset((jj * a_dim1 + 1 as libc::c_int) as isize), + 1 as libc::c_int, + &mut *z__.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + } + jj = *indx.offset(ip as isize); + *z__.offset(ip as isize) /= *a.offset((ip + jj * a_dim1) as isize); + ip -= 1; + } + iter += 1; + if !(iter <= itmax) { + current_block = 14151449179633537743; + break 'c_7799; + } + alpha = one; + jj = 0 as libc::c_int; + i__2 = nsetp; + ip = 1 as libc::c_int; + while ip <= i__2 { + if !(*z__.offset(ip as isize) > 0.0f64) { + l = *indx.offset(ip as isize); + t = -*x.offset(l as isize) + / (*z__.offset(ip as isize) - *x.offset(l as isize)); + if !(alpha < t) { + alpha = t; + jj = ip; + } + } + ip += 1; + } + i__2 = nsetp; + ip = 1 as libc::c_int; + while ip <= i__2 { + l = *indx.offset(ip as isize); + *x + .offset( + l as isize, + ) = (one - alpha) * *x.offset(l as isize) + + alpha * *z__.offset(ip as isize); + ip += 1; + } + if jj == 0 as libc::c_int { + break; + } + i__ = *indx.offset(jj as isize); + 'c_7847: loop { + *x.offset(i__ as isize) = 0.0f64; + jj += 1; + i__2 = nsetp; + j = jj; + while j <= i__2 { + ii = *indx.offset(j as isize); + *indx.offset((j - 1 as libc::c_int) as isize) = ii; + dsrotg_( + &mut *a + .offset((j - 1 as libc::c_int + ii * a_dim1) as isize), + &mut *a.offset((j + ii * a_dim1) as isize), + &mut c__, + &mut s, + ); + t = *a.offset((j - 1 as libc::c_int + ii * a_dim1) as isize); + dsrot_( + *n, + &mut *a.offset((j - 1 as libc::c_int + a_dim1) as isize), + *mda, + &mut *a.offset((j + a_dim1) as isize), + *mda, + &mut c__, + &mut s, + ); + *a.offset((j - 1 as libc::c_int + ii * a_dim1) as isize) = t; + *a.offset((j + ii * a_dim1) as isize) = 0.0f64; + dsrot_( + 1 as libc::c_int, + &mut *b.offset((j - 1 as libc::c_int) as isize), + 1 as libc::c_int, + &mut *b.offset(j as isize), + 1 as libc::c_int, + &mut c__, + &mut s, + ); + j += 1; + } + npp1 = nsetp; + nsetp -= 1; + iz1 -= 1; + *indx.offset(iz1 as isize) = i__; + if nsetp <= 0 as libc::c_int { + current_block = 14151449179633537743; + break 'c_7799; + } + i__2 = nsetp; + jj = 1 as libc::c_int; + loop { + if !(jj <= i__2) { + break 'c_7847; + } + i__ = *indx.offset(jj as isize); + if *x.offset(i__ as isize) <= 0.0f64 { + break; + } + jj += 1; + } + } + dcopy___( + m, + &mut *b.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *z__.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + } + } + match current_block { + 14151449179633537743 => { + *mode = 3 as libc::c_int; + } + _ => {} + } + k = if npp1 <= *m { npp1 } else { *m }; + i__2 = *m - nsetp; + *rnorm = dnrm2___(&mut i__2, &mut *b.offset(k as isize), 1 as libc::c_int); + if npp1 > *m { + *w.offset(1 as libc::c_int as isize) = 0.0f64; + dcopy___( + n, + &mut *w.offset(1 as libc::c_int as isize), + 0 as libc::c_int, + &mut *w.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + } + } +} +unsafe extern "C" fn ldp_( + mut g: *mut libc::c_double, + mut mg: *mut libc::c_int, + mut m: *mut libc::c_int, + mut n: *mut libc::c_int, + mut h__: *mut libc::c_double, + mut x: *mut libc::c_double, + mut xnorm: *mut libc::c_double, + mut w: *mut libc::c_double, + mut indx: *mut libc::c_int, + mut mode: *mut libc::c_int, +) { + let one: libc::c_double = 1.0f64; + let mut g_dim1: libc::c_int = 0; + let mut g_offset: libc::c_int = 0; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut d__1: libc::c_double = 0.; + let mut i__: libc::c_int = 0; + let mut j: libc::c_int = 0; + let mut n1: libc::c_int = 0; + let mut if__: libc::c_int = 0; + let mut iw: libc::c_int = 0; + let mut iy: libc::c_int = 0; + let mut iz: libc::c_int = 0; + let mut fac: libc::c_double = 0.; + let mut rnorm: libc::c_double = 0.; + let mut iwdual: libc::c_int = 0; + indx = indx.offset(-1); + h__ = h__.offset(-1); + x = x.offset(-1); + g_dim1 = *mg; + g_offset = 1 as libc::c_int + g_dim1; + g = g.offset(-(g_offset as isize)); + w = w.offset(-1); + *mode = 2 as libc::c_int; + if !(*n <= 0 as libc::c_int) { + *mode = 1 as libc::c_int; + *x.offset(1 as libc::c_int as isize) = 0.0f64; + dcopy___( + n, + &mut *x.offset(1 as libc::c_int as isize), + 0 as libc::c_int, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + *xnorm = 0.0f64; + if !(*m == 0 as libc::c_int) { + iw = 0 as libc::c_int; + i__1 = *m; + j = 1 as libc::c_int; + while j <= i__1 { + i__2 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__2 { + iw += 1; + *w.offset(iw as isize) = *g.offset((j + i__ * g_dim1) as isize); + i__ += 1; + } + iw += 1; + *w.offset(iw as isize) = *h__.offset(j as isize); + j += 1; + } + if__ = iw + 1 as libc::c_int; + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + iw += 1; + *w.offset(iw as isize) = 0.0f64; + i__ += 1; + } + *w.offset((iw + 1 as libc::c_int) as isize) = one; + n1 = *n + 1 as libc::c_int; + iz = iw + 2 as libc::c_int; + iy = iz + n1; + iwdual = iy + *m; + nnls_( + &mut *w.offset(1 as libc::c_int as isize), + &mut n1, + &mut n1, + m, + &mut *w.offset(if__ as isize), + &mut *w.offset(iy as isize), + &mut rnorm, + &mut *w.offset(iwdual as isize), + &mut *w.offset(iz as isize), + &mut *indx.offset(1 as libc::c_int as isize), + mode, + ); + if !(*mode != 1 as libc::c_int) { + *mode = 4 as libc::c_int; + if !(rnorm <= 0.0f64) { + fac = one + - ddot_sl__( + m, + &mut *h__.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *w.offset(iy as isize), + 1 as libc::c_int, + ); + d__1 = one + fac; + if !(d__1 - one <= 0.0f64) { + *mode = 1 as libc::c_int; + fac = one / fac; + i__1 = *n; + j = 1 as libc::c_int; + while j <= i__1 { + *x + .offset( + j as isize, + ) = fac + * ddot_sl__( + m, + &mut *g.offset((j * g_dim1 + 1 as libc::c_int) as isize), + 1 as libc::c_int, + &mut *w.offset(iy as isize), + 1 as libc::c_int, + ); + j += 1; + } + *xnorm = dnrm2___( + n, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + *w.offset(1 as libc::c_int as isize) = 0.0f64; + dcopy___( + m, + &mut *w.offset(1 as libc::c_int as isize), + 0 as libc::c_int, + &mut *w.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + daxpy_sl__( + m, + &mut fac, + &mut *w.offset(iy as isize), + 1 as libc::c_int, + &mut *w.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + } + } + } + } + } +} +unsafe extern "C" fn lsi_( + mut e: *mut libc::c_double, + mut f: *mut libc::c_double, + mut g: *mut libc::c_double, + mut h__: *mut libc::c_double, + mut le: *mut libc::c_int, + mut me: *mut libc::c_int, + mut lg: *mut libc::c_int, + mut mg: *mut libc::c_int, + mut n: *mut libc::c_int, + mut x: *mut libc::c_double, + mut xnorm: *mut libc::c_double, + mut w: *mut libc::c_double, + mut jw: *mut libc::c_int, + mut mode: *mut libc::c_int, +) { + let mut current_block: u64; + let epmach: libc::c_double = 2.22e-16f64; + let one: libc::c_double = 1.0f64; + let mut e_dim1: libc::c_int = 0; + let mut e_offset: libc::c_int = 0; + let mut g_dim1: libc::c_int = 0; + let mut g_offset: libc::c_int = 0; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut i__3: libc::c_int = 0; + let mut d__1: libc::c_double = 0.; + let mut i__: libc::c_int = 0; + let mut j: libc::c_int = 0; + let mut t: libc::c_double = 0.; + f = f.offset(-1); + jw = jw.offset(-1); + h__ = h__.offset(-1); + x = x.offset(-1); + g_dim1 = *lg; + g_offset = 1 as libc::c_int + g_dim1; + g = g.offset(-(g_offset as isize)); + e_dim1 = *le; + e_offset = 1 as libc::c_int + e_dim1; + e = e.offset(-(e_offset as isize)); + w = w.offset(-1); + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + i__2 = i__ + 1 as libc::c_int; + j = if i__2 <= *n { i__2 } else { *n }; + i__2 = i__ + 1 as libc::c_int; + i__3 = *n - i__; + h12_( + &c__1, + &mut i__, + &mut i__2, + me, + &mut *e.offset((i__ * e_dim1 + 1 as libc::c_int) as isize), + &c__1, + &mut t, + &mut *e.offset((j * e_dim1 + 1 as libc::c_int) as isize), + &c__1, + le, + &mut i__3, + ); + i__2 = i__ + 1 as libc::c_int; + h12_( + &c__2, + &mut i__, + &mut i__2, + me, + &mut *e.offset((i__ * e_dim1 + 1 as libc::c_int) as isize), + &c__1, + &mut t, + &mut *f.offset(1 as libc::c_int as isize), + &c__1, + &c__1, + &c__1, + ); + i__ += 1; + } + *mode = 5 as libc::c_int; + i__2 = *mg; + i__ = 1 as libc::c_int; + 's_121: loop { + if !(i__ <= i__2) { + current_block = 14434620278749266018; + break; + } + i__1 = *n; + j = 1 as libc::c_int; + while j <= i__1 { + d__1 = *e.offset((j + j * e_dim1) as isize); + if fabs(d__1) < epmach { + current_block = 16236312898130099128; + break 's_121; + } + i__3 = j - 1 as libc::c_int; + *g + .offset( + (i__ + j * g_dim1) as isize, + ) = (*g.offset((i__ + j * g_dim1) as isize) + - ddot_sl__( + &mut i__3, + &mut *g.offset((i__ + g_dim1) as isize), + *lg, + &mut *e.offset((j * e_dim1 + 1 as libc::c_int) as isize), + 1 as libc::c_int, + )) / *e.offset((j + j * e_dim1) as isize); + j += 1; + } + *h__.offset(i__ as isize) + -= ddot_sl__( + n, + &mut *g.offset((i__ + g_dim1) as isize), + *lg, + &mut *f.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + i__ += 1; + } + match current_block { + 14434620278749266018 => { + ldp_( + &mut *g.offset(g_offset as isize), + lg, + mg, + n, + &mut *h__.offset(1 as libc::c_int as isize), + &mut *x.offset(1 as libc::c_int as isize), + xnorm, + &mut *w.offset(1 as libc::c_int as isize), + &mut *jw.offset(1 as libc::c_int as isize), + mode, + ); + if !(*mode != 1 as libc::c_int) { + daxpy_sl__( + n, + &one, + &mut *f.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + i__ = *n; + while i__ >= 1 as libc::c_int { + i__2 = i__ + 1 as libc::c_int; + j = if i__2 <= *n { i__2 } else { *n }; + i__2 = *n - i__; + *x + .offset( + i__ as isize, + ) = (*x.offset(i__ as isize) + - ddot_sl__( + &mut i__2, + &mut *e.offset((i__ + j * e_dim1) as isize), + *le, + &mut *x.offset(j as isize), + 1 as libc::c_int, + )) / *e.offset((i__ + i__ * e_dim1) as isize); + i__ -= 1; + } + i__2 = *n + 1 as libc::c_int; + j = if i__2 <= *me { i__2 } else { *me }; + i__2 = *me - *n; + t = dnrm2___(&mut i__2, &mut *f.offset(j as isize), 1 as libc::c_int); + *xnorm = sqrt(*xnorm * *xnorm + t * t); + } + } + _ => {} + }; +} +unsafe extern "C" fn hfti_( + mut a: *mut libc::c_double, + mut mda: *mut libc::c_int, + mut m: *mut libc::c_int, + mut n: *mut libc::c_int, + mut b: *mut libc::c_double, + mut mdb: *mut libc::c_int, + mut nb: *const libc::c_int, + mut tau: *mut libc::c_double, + mut krank: *mut libc::c_int, + mut rnorm: *mut libc::c_double, + mut h__: *mut libc::c_double, + mut g: *mut libc::c_double, + mut ip: *mut libc::c_int, +) { + let mut current_block: u64; + let factor: libc::c_double = 0.001f64; + let mut a_dim1: libc::c_int = 0; + let mut a_offset: libc::c_int = 0; + let mut b_dim1: libc::c_int = 0; + let mut b_offset: libc::c_int = 0; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut i__3: libc::c_int = 0; + let mut d__1: libc::c_double = 0.; + let mut i__: libc::c_int = 0; + let mut j: libc::c_int = 0; + let mut k: libc::c_int = 0; + let mut l: libc::c_int = 0; + let mut jb: libc::c_int = 0; + let mut kp1: libc::c_int = 0; + let mut tmp: libc::c_double = 0.; + let mut hmax: libc::c_double = 0.; + let mut lmax: libc::c_int = 0; + let mut ldiag: libc::c_int = 0; + ip = ip.offset(-1); + g = g.offset(-1); + h__ = h__.offset(-1); + a_dim1 = *mda; + a_offset = 1 as libc::c_int + a_dim1; + a = a.offset(-(a_offset as isize)); + rnorm = rnorm.offset(-1); + b_dim1 = *mdb; + b_offset = 1 as libc::c_int + b_dim1; + b = b.offset(-(b_offset as isize)); + k = 0 as libc::c_int; + ldiag = if *m <= *n { *m } else { *n }; + if !(ldiag <= 0 as libc::c_int) { + i__1 = ldiag; + j = 1 as libc::c_int; + while j <= i__1 { + let mut current_block_56: u64; + if j == 1 as libc::c_int { + current_block_56 = 14274665232496457393; + } else { + lmax = j; + i__2 = *n; + l = j; + while l <= i__2 { + d__1 = *a.offset((j - 1 as libc::c_int + l * a_dim1) as isize); + *h__.offset(l as isize) -= d__1 * d__1; + if *h__.offset(l as isize) > *h__.offset(lmax as isize) { + lmax = l; + } + l += 1; + } + d__1 = hmax + factor * *h__.offset(lmax as isize); + if d__1 - hmax > 0.0f64 { + current_block_56 = 8942514791948061562; + } else { + current_block_56 = 14274665232496457393; + } + } + match current_block_56 { + 14274665232496457393 => { + lmax = j; + i__2 = *n; + l = j; + while l <= i__2 { + *h__.offset(l as isize) = 0.0f64; + i__3 = *m; + i__ = j; + while i__ <= i__3 { + d__1 = *a.offset((i__ + l * a_dim1) as isize); + *h__.offset(l as isize) += d__1 * d__1; + i__ += 1; + } + if *h__.offset(l as isize) > *h__.offset(lmax as isize) { + lmax = l; + } + l += 1; + } + hmax = *h__.offset(lmax as isize); + } + _ => {} + } + *ip.offset(j as isize) = lmax; + if !(*ip.offset(j as isize) == j) { + i__2 = *m; + i__ = 1 as libc::c_int; + while i__ <= i__2 { + tmp = *a.offset((i__ + j * a_dim1) as isize); + *a + .offset( + (i__ + j * a_dim1) as isize, + ) = *a.offset((i__ + lmax * a_dim1) as isize); + *a.offset((i__ + lmax * a_dim1) as isize) = tmp; + i__ += 1; + } + *h__.offset(lmax as isize) = *h__.offset(j as isize); + } + i__2 = j + 1 as libc::c_int; + i__ = if i__2 <= *n { i__2 } else { *n }; + i__2 = j + 1 as libc::c_int; + i__3 = *n - j; + h12_( + &c__1, + &mut j, + &mut i__2, + m, + &mut *a.offset((j * a_dim1 + 1 as libc::c_int) as isize), + &c__1, + &mut *h__.offset(j as isize), + &mut *a.offset((i__ * a_dim1 + 1 as libc::c_int) as isize), + &c__1, + mda, + &mut i__3, + ); + i__2 = j + 1 as libc::c_int; + h12_( + &c__2, + &mut j, + &mut i__2, + m, + &mut *a.offset((j * a_dim1 + 1 as libc::c_int) as isize), + &c__1, + &mut *h__.offset(j as isize), + &mut *b.offset(b_offset as isize), + &c__1, + mdb, + nb, + ); + j += 1; + } + i__2 = ldiag; + j = 1 as libc::c_int; + loop { + if !(j <= i__2) { + current_block = 18038362259723567392; + break; + } + d__1 = *a.offset((j + j * a_dim1) as isize); + if fabs(d__1) <= *tau { + current_block = 8916208649952454714; + break; + } + j += 1; + } + match current_block { + 18038362259723567392 => { + k = ldiag; + } + _ => { + k = j - 1 as libc::c_int; + } + } + kp1 = k + 1 as libc::c_int; + i__2 = *nb; + jb = 1 as libc::c_int; + while jb <= i__2 { + i__1 = *m - k; + *rnorm + .offset( + jb as isize, + ) = dnrm2___( + &mut i__1, + &mut *b.offset((kp1 + jb * b_dim1) as isize), + 1 as libc::c_int, + ); + jb += 1; + } + if k > 0 as libc::c_int { + if !(k == *n) { + i__ = k; + while i__ >= 1 as libc::c_int { + i__2 = i__ - 1 as libc::c_int; + h12_( + &c__1, + &mut i__, + &mut kp1, + n, + &mut *a.offset((i__ + a_dim1) as isize), + mda, + &mut *g.offset(i__ as isize), + &mut *a.offset(a_offset as isize), + mda, + &c__1, + &mut i__2, + ); + i__ -= 1; + } + } + i__2 = *nb; + jb = 1 as libc::c_int; + while jb <= i__2 { + i__ = k; + while i__ >= 1 as libc::c_int { + i__1 = i__ + 1 as libc::c_int; + j = if i__1 <= *n { i__1 } else { *n }; + i__1 = k - i__; + *b + .offset( + (i__ + jb * b_dim1) as isize, + ) = (*b.offset((i__ + jb * b_dim1) as isize) + - ddot_sl__( + &mut i__1, + &mut *a.offset((i__ + j * a_dim1) as isize), + *mda, + &mut *b.offset((j + jb * b_dim1) as isize), + 1 as libc::c_int, + )) / *a.offset((i__ + i__ * a_dim1) as isize); + i__ -= 1; + } + if !(k == *n) { + i__1 = *n; + j = kp1; + while j <= i__1 { + *b.offset((j + jb * b_dim1) as isize) = 0.0f64; + j += 1; + } + i__1 = k; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + h12_( + &c__2, + &mut i__, + &mut kp1, + n, + &mut *a.offset((i__ + a_dim1) as isize), + mda, + &mut *g.offset(i__ as isize), + &mut *b.offset((jb * b_dim1 + 1 as libc::c_int) as isize), + &c__1, + mdb, + &c__1, + ); + i__ += 1; + } + } + j = ldiag; + while j >= 1 as libc::c_int { + if !(*ip.offset(j as isize) == j) { + l = *ip.offset(j as isize); + tmp = *b.offset((l + jb * b_dim1) as isize); + *b + .offset( + (l + jb * b_dim1) as isize, + ) = *b.offset((j + jb * b_dim1) as isize); + *b.offset((j + jb * b_dim1) as isize) = tmp; + } + j -= 1; + } + jb += 1; + } + } else { + i__1 = *nb; + jb = 1 as libc::c_int; + while jb <= i__1 { + i__2 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__2 { + *b.offset((i__ + jb * b_dim1) as isize) = 0.0f64; + i__ += 1; + } + jb += 1; + } + } + } + *krank = k; +} +unsafe extern "C" fn lsei_( + mut c__: *mut libc::c_double, + mut d__: *mut libc::c_double, + mut e: *mut libc::c_double, + mut f: *mut libc::c_double, + mut g: *mut libc::c_double, + mut h__: *mut libc::c_double, + mut lc: *mut libc::c_int, + mut mc: *mut libc::c_int, + mut le: *mut libc::c_int, + mut me: *mut libc::c_int, + mut lg: *mut libc::c_int, + mut mg: *mut libc::c_int, + mut n: *mut libc::c_int, + mut x: *mut libc::c_double, + mut xnrm: *mut libc::c_double, + mut w: *mut libc::c_double, + mut jw: *mut libc::c_int, + mut mode: *mut libc::c_int, +) { + let mut current_block: u64; + let epmach: libc::c_double = 2.22e-16f64; + let mut c_dim1: libc::c_int = 0; + let mut c_offset: libc::c_int = 0; + let mut e_dim1: libc::c_int = 0; + let mut e_offset: libc::c_int = 0; + let mut g_dim1: libc::c_int = 0; + let mut g_offset: libc::c_int = 0; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut i__3: libc::c_int = 0; + let mut d__1: libc::c_double = 0.; + let mut i__: libc::c_int = 0; + let mut j: libc::c_int = 0; + let mut k: libc::c_int = 0; + let mut l: libc::c_int = 0; + let mut t: libc::c_double = 0.; + let mut ie: libc::c_int = 0; + let mut if__: libc::c_int = 0; + let mut ig: libc::c_int = 0; + let mut iw: libc::c_int = 0; + let mut mc1: libc::c_int = 0; + let mut krank: libc::c_int = 0; + d__ = d__.offset(-1); + f = f.offset(-1); + h__ = h__.offset(-1); + x = x.offset(-1); + g_dim1 = *lg; + g_offset = 1 as libc::c_int + g_dim1; + g = g.offset(-(g_offset as isize)); + e_dim1 = *le; + e_offset = 1 as libc::c_int + e_dim1; + e = e.offset(-(e_offset as isize)); + c_dim1 = *lc; + c_offset = 1 as libc::c_int + c_dim1; + c__ = c__.offset(-(c_offset as isize)); + w = w.offset(-1); + jw = jw.offset(-1); + *mode = 2 as libc::c_int; + if !(*mc > *n) { + l = *n - *mc; + mc1 = *mc + 1 as libc::c_int; + iw = (l + 1 as libc::c_int) * (*mg + 2 as libc::c_int) + + (*mg << 1 as libc::c_int) + *mc; + ie = iw + *mc + 1 as libc::c_int; + if__ = ie + *me * l; + ig = if__ + *me; + i__1 = *mc; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + i__2 = i__ + 1 as libc::c_int; + j = if i__2 <= *lc { i__2 } else { *lc }; + i__2 = i__ + 1 as libc::c_int; + i__3 = *mc - i__; + h12_( + &c__1, + &mut i__, + &mut i__2, + n, + &mut *c__.offset((i__ + c_dim1) as isize), + lc, + &mut *w.offset((iw + i__) as isize), + &mut *c__.offset((j + c_dim1) as isize), + lc, + &c__1, + &mut i__3, + ); + i__2 = i__ + 1 as libc::c_int; + h12_( + &c__2, + &mut i__, + &mut i__2, + n, + &mut *c__.offset((i__ + c_dim1) as isize), + lc, + &mut *w.offset((iw + i__) as isize), + &mut *e.offset(e_offset as isize), + le, + &c__1, + me, + ); + i__2 = i__ + 1 as libc::c_int; + h12_( + &c__2, + &mut i__, + &mut i__2, + n, + &mut *c__.offset((i__ + c_dim1) as isize), + lc, + &mut *w.offset((iw + i__) as isize), + &mut *g.offset(g_offset as isize), + lg, + &c__1, + mg, + ); + i__ += 1; + } + *mode = 6 as libc::c_int; + i__2 = *mc; + i__ = 1 as libc::c_int; + loop { + if !(i__ <= i__2) { + current_block = 3222590281903869779; + break; + } + d__1 = *c__.offset((i__ + i__ * c_dim1) as isize); + if fabs(d__1) < epmach { + current_block = 7757720978182437014; + break; + } + i__1 = i__ - 1 as libc::c_int; + *x + .offset( + i__ as isize, + ) = (*d__.offset(i__ as isize) + - ddot_sl__( + &mut i__1, + &mut *c__.offset((i__ + c_dim1) as isize), + *lc, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + )) / *c__.offset((i__ + i__ * c_dim1) as isize); + i__ += 1; + } + match current_block { + 7757720978182437014 => {} + _ => { + *mode = 1 as libc::c_int; + *w.offset(mc1 as isize) = 0.0f64; + i__2 = *mg; + dcopy___( + &mut i__2, + &mut *w.offset(mc1 as isize), + 0 as libc::c_int, + &mut *w.offset(mc1 as isize), + 1 as libc::c_int, + ); + if *mc == *n { + current_block = 13396851843167794670; + } else { + i__2 = *me; + i__ = 1 as libc::c_int; + while i__ <= i__2 { + *w + .offset( + (if__ - 1 as libc::c_int + i__) as isize, + ) = *f.offset(i__ as isize) + - ddot_sl__( + mc, + &mut *e.offset((i__ + e_dim1) as isize), + *le, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + i__ += 1; + } + i__2 = *me; + i__ = 1 as libc::c_int; + while i__ <= i__2 { + dcopy___( + &mut l, + &mut *e.offset((i__ + mc1 * e_dim1) as isize), + *le, + &mut *w.offset((ie - 1 as libc::c_int + i__) as isize), + *me, + ); + i__ += 1; + } + i__2 = *mg; + i__ = 1 as libc::c_int; + while i__ <= i__2 { + dcopy___( + &mut l, + &mut *g.offset((i__ + mc1 * g_dim1) as isize), + *lg, + &mut *w.offset((ig - 1 as libc::c_int + i__) as isize), + *mg, + ); + i__ += 1; + } + if *mg > 0 as libc::c_int { + i__2 = *mg; + i__ = 1 as libc::c_int; + while i__ <= i__2 { + *h__.offset(i__ as isize) + -= ddot_sl__( + mc, + &mut *g.offset((i__ + g_dim1) as isize), + *lg, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + i__ += 1; + } + lsi_( + &mut *w.offset(ie as isize), + &mut *w.offset(if__ as isize), + &mut *w.offset(ig as isize), + &mut *h__.offset(1 as libc::c_int as isize), + me, + me, + mg, + mg, + &mut l, + &mut *x.offset(mc1 as isize), + xnrm, + &mut *w.offset(mc1 as isize), + &mut *jw.offset(1 as libc::c_int as isize), + mode, + ); + if *mc == 0 as libc::c_int { + current_block = 7757720978182437014; + } else { + t = dnrm2___( + mc, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + *xnrm = sqrt(*xnrm * *xnrm + t * t); + if *mode != 1 as libc::c_int { + current_block = 7757720978182437014; + } else { + current_block = 13396851843167794670; + } + } + } else { + *mode = 7 as libc::c_int; + k = if *le >= *n { *le } else { *n }; + t = sqrt(epmach); + hfti_( + &mut *w.offset(ie as isize), + me, + me, + &mut l, + &mut *w.offset(if__ as isize), + &mut k, + &c__1, + &mut t, + &mut krank, + xnrm, + &mut *w.offset(1 as libc::c_int as isize), + &mut *w.offset((l + 1 as libc::c_int) as isize), + &mut *jw.offset(1 as libc::c_int as isize), + ); + dcopy___( + &mut l, + &mut *w.offset(if__ as isize), + 1 as libc::c_int, + &mut *x.offset(mc1 as isize), + 1 as libc::c_int, + ); + if krank != l { + current_block = 7757720978182437014; + } else { + *mode = 1 as libc::c_int; + current_block = 13396851843167794670; + } + } + } + match current_block { + 7757720978182437014 => {} + _ => { + i__2 = *me; + i__ = 1 as libc::c_int; + while i__ <= i__2 { + *f + .offset( + i__ as isize, + ) = ddot_sl__( + n, + &mut *e.offset((i__ + e_dim1) as isize), + *le, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ) - *f.offset(i__ as isize); + i__ += 1; + } + i__2 = *mc; + i__ = 1 as libc::c_int; + while i__ <= i__2 { + *d__ + .offset( + i__ as isize, + ) = ddot_sl__( + me, + &mut *e.offset((i__ * e_dim1 + 1 as libc::c_int) as isize), + 1 as libc::c_int, + &mut *f.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ) + - ddot_sl__( + mg, + &mut *g.offset((i__ * g_dim1 + 1 as libc::c_int) as isize), + 1 as libc::c_int, + &mut *w.offset(mc1 as isize), + 1 as libc::c_int, + ); + i__ += 1; + } + i__ = *mc; + while i__ >= 1 as libc::c_int { + i__2 = i__ + 1 as libc::c_int; + h12_( + &c__2, + &mut i__, + &mut i__2, + n, + &mut *c__.offset((i__ + c_dim1) as isize), + lc, + &mut *w.offset((iw + i__) as isize), + &mut *x.offset(1 as libc::c_int as isize), + &c__1, + &c__1, + &c__1, + ); + i__ -= 1; + } + i__ = *mc; + while i__ >= 1 as libc::c_int { + i__2 = i__ + 1 as libc::c_int; + j = if i__2 <= *lc { i__2 } else { *lc }; + i__2 = *mc - i__; + *w + .offset( + i__ as isize, + ) = (*d__.offset(i__ as isize) + - ddot_sl__( + &mut i__2, + &mut *c__.offset((j + i__ * c_dim1) as isize), + 1 as libc::c_int, + &mut *w.offset(j as isize), + 1 as libc::c_int, + )) / *c__.offset((i__ + i__ * c_dim1) as isize); + i__ -= 1; + } + } + } + } + } + } +} +unsafe extern "C" fn lsq_( + mut m: *mut libc::c_int, + mut meq: *mut libc::c_int, + mut n: *mut libc::c_int, + mut nl: *mut libc::c_int, + mut la: *mut libc::c_int, + mut l: *mut libc::c_double, + mut g: *mut libc::c_double, + mut a: *mut libc::c_double, + mut b: *mut libc::c_double, + mut xl: *const libc::c_double, + mut xu: *const libc::c_double, + mut x: *mut libc::c_double, + mut y: *mut libc::c_double, + mut w: *mut libc::c_double, + mut jw: *mut libc::c_int, + mut mode: *mut libc::c_int, +) { + let one: libc::c_double = 1.0f64; + let mut a_dim1: libc::c_int = 0; + let mut a_offset: libc::c_int = 0; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut d__1: libc::c_double = 0.; + let mut i__: libc::c_int = 0; + let mut i1: libc::c_int = 0; + let mut i2: libc::c_int = 0; + let mut i3: libc::c_int = 0; + let mut i4: libc::c_int = 0; + let mut m1: libc::c_int = 0; + let mut n1: libc::c_int = 0; + let mut n2: libc::c_int = 0; + let mut n3: libc::c_int = 0; + let mut ic: libc::c_int = 0; + let mut id: libc::c_int = 0; + let mut ie: libc::c_int = 0; + let mut if__: libc::c_int = 0; + let mut ig: libc::c_int = 0; + let mut ih: libc::c_int = 0; + let mut il: libc::c_int = 0; + let mut im: libc::c_int = 0; + let mut ip: libc::c_int = 0; + let mut iu: libc::c_int = 0; + let mut iw: libc::c_int = 0; + let mut diag: libc::c_double = 0.; + let mut mineq: libc::c_int = 0; + let mut xnorm: libc::c_double = 0.; + y = y.offset(-1); + x = x.offset(-1); + xu = xu.offset(-1); + xl = xl.offset(-1); + g = g.offset(-1); + l = l.offset(-1); + b = b.offset(-1); + a_dim1 = *la; + a_offset = 1 as libc::c_int + a_dim1; + a = a.offset(-(a_offset as isize)); + w = w.offset(-1); + jw = jw.offset(-1); + n1 = *n + 1 as libc::c_int; + mineq = *m - *meq; + m1 = mineq + *n + *n; + n2 = n1 * *n / 2 as libc::c_int + 1 as libc::c_int; + if n2 == *nl { + n2 = 0 as libc::c_int; + } else { + n2 = 1 as libc::c_int; + } + n3 = *n - n2; + i2 = 1 as libc::c_int; + i3 = 1 as libc::c_int; + i4 = 1 as libc::c_int; + ie = 1 as libc::c_int; + if__ = *n * *n + 1 as libc::c_int; + i__1 = n3; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + i1 = n1 - i__; + diag = sqrt(*l.offset(i2 as isize)); + *w.offset(i3 as isize) = 0.0f64; + dcopy___( + &mut i1, + &mut *w.offset(i3 as isize), + 0 as libc::c_int, + &mut *w.offset(i3 as isize), + 1 as libc::c_int, + ); + i__2 = i1 - n2; + dcopy___( + &mut i__2, + &mut *l.offset(i2 as isize), + 1 as libc::c_int, + &mut *w.offset(i3 as isize), + *n, + ); + i__2 = i1 - n2; + dscal_sl__(&mut i__2, &mut diag, &mut *w.offset(i3 as isize), *n); + *w.offset(i3 as isize) = diag; + i__2 = i__ - 1 as libc::c_int; + *w + .offset( + (if__ - 1 as libc::c_int + i__) as isize, + ) = (*g.offset(i__ as isize) + - ddot_sl__( + &mut i__2, + &mut *w.offset(i4 as isize), + 1 as libc::c_int, + &mut *w.offset(if__ as isize), + 1 as libc::c_int, + )) / diag; + i2 = i2 + i1 - n2; + i3 += n1; + i4 += *n; + i__ += 1; + } + if n2 == 1 as libc::c_int { + *w.offset(i3 as isize) = *l.offset(*nl as isize); + *w.offset(i4 as isize) = 0.0f64; + dcopy___( + &mut n3, + &mut *w.offset(i4 as isize), + 0 as libc::c_int, + &mut *w.offset(i4 as isize), + 1 as libc::c_int, + ); + *w.offset((if__ - 1 as libc::c_int + *n) as isize) = 0.0f64; + } + d__1 = -one; + dscal_sl__(n, &mut d__1, &mut *w.offset(if__ as isize), 1 as libc::c_int); + ic = if__ + *n; + id = ic + *meq * *n; + if *meq > 0 as libc::c_int { + i__1 = *meq; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + dcopy___( + n, + &mut *a.offset((i__ + a_dim1) as isize), + *la, + &mut *w.offset((ic - 1 as libc::c_int + i__) as isize), + *meq, + ); + i__ += 1; + } + dcopy___( + meq, + &mut *b.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *w.offset(id as isize), + 1 as libc::c_int, + ); + d__1 = -one; + dscal_sl__(meq, &mut d__1, &mut *w.offset(id as isize), 1 as libc::c_int); + } + ig = id + *meq; + if mineq > 0 as libc::c_int { + i__1 = mineq; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + dcopy___( + n, + &mut *a.offset((*meq + i__ + a_dim1) as isize), + *la, + &mut *w.offset((ig - 1 as libc::c_int + i__) as isize), + m1, + ); + i__ += 1; + } + } + ip = ig + mineq; + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + *w.offset((ip - 1 as libc::c_int + i__) as isize) = 0.0f64; + dcopy___( + n, + &mut *w.offset((ip - 1 as libc::c_int + i__) as isize), + 0 as libc::c_int, + &mut *w.offset((ip - 1 as libc::c_int + i__) as isize), + m1, + ); + i__ += 1; + } + i__1 = m1 + 1 as libc::c_int; + i__ = 1 as libc::c_int; + while i__ <= *n { + if nlopt_isinf(*xl.offset(i__ as isize)) == 0 { + *w.offset((ip - i__1 + i__ * i__1) as isize) = 1.0f64; + } + i__ += 1; + } + im = ip + *n; + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + *w.offset((im - 1 as libc::c_int + i__) as isize) = 0.0f64; + dcopy___( + n, + &mut *w.offset((im - 1 as libc::c_int + i__) as isize), + 0 as libc::c_int, + &mut *w.offset((im - 1 as libc::c_int + i__) as isize), + m1, + ); + i__ += 1; + } + i__1 = m1 + 1 as libc::c_int; + i__ = 1 as libc::c_int; + while i__ <= *n { + if nlopt_isinf(*xu.offset(i__ as isize)) == 0 { + *w.offset((im - i__1 + i__ * i__1) as isize) = -1.0f64; + } + i__ += 1; + } + ih = ig + m1 * *n; + if mineq > 0 as libc::c_int { + dcopy___( + &mut mineq, + &mut *b.offset((*meq + 1 as libc::c_int) as isize), + 1 as libc::c_int, + &mut *w.offset(ih as isize), + 1 as libc::c_int, + ); + d__1 = -one; + dscal_sl__(&mut mineq, &mut d__1, &mut *w.offset(ih as isize), 1 as libc::c_int); + } + il = ih + mineq; + iu = il + *n; + i__ = 1 as libc::c_int; + while i__ <= *n { + *w + .offset( + (il - 1 as libc::c_int + i__) as isize, + ) = if nlopt_isinf(*xl.offset(i__ as isize)) != 0 { + 0 as libc::c_int as libc::c_double + } else { + *xl.offset(i__ as isize) + }; + *w + .offset( + (iu - 1 as libc::c_int + i__) as isize, + ) = if nlopt_isinf(*xu.offset(i__ as isize)) != 0 { + 0 as libc::c_int as libc::c_double + } else { + -*xu.offset(i__ as isize) + }; + i__ += 1; + } + iw = iu + *n; + i__1 = if 1 as libc::c_int >= *meq { 1 as libc::c_int } else { *meq }; + lsei_( + &mut *w.offset(ic as isize), + &mut *w.offset(id as isize), + &mut *w.offset(ie as isize), + &mut *w.offset(if__ as isize), + &mut *w.offset(ig as isize), + &mut *w.offset(ih as isize), + &mut i__1, + meq, + n, + n, + &mut m1, + &mut m1, + n, + &mut *x.offset(1 as libc::c_int as isize), + &mut xnorm, + &mut *w.offset(iw as isize), + &mut *jw.offset(1 as libc::c_int as isize), + mode, + ); + if *mode == 1 as libc::c_int { + dcopy___( + m, + &mut *w.offset(iw as isize), + 1 as libc::c_int, + &mut *y.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + dcopy___( + &mut n3, + &mut *w.offset((iw + *m) as isize), + 1 as libc::c_int, + &mut *y.offset((*m + 1 as libc::c_int) as isize), + 1 as libc::c_int, + ); + dcopy___( + &mut n3, + &mut *w.offset((iw + *m + *n) as isize), + 1 as libc::c_int, + &mut *y.offset((*m + n3 + 1 as libc::c_int) as isize), + 1 as libc::c_int, + ); + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + if *x.offset(i__ as isize) < *xl.offset(i__ as isize) { + *x.offset(i__ as isize) = *xl.offset(i__ as isize); + } else if *x.offset(i__ as isize) > *xu.offset(i__ as isize) { + *x.offset(i__ as isize) = *xu.offset(i__ as isize); + } + i__ += 1; + } + } +} +unsafe extern "C" fn ldl_( + mut n: *mut libc::c_int, + mut a: *mut libc::c_double, + mut z__: *mut libc::c_double, + mut sigma: *mut libc::c_double, + mut w: *mut libc::c_double, +) { + let one: libc::c_double = 1.0f64; + let four: libc::c_double = 4.0f64; + let epmach: libc::c_double = 2.22e-16f64; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut i__: libc::c_int = 0; + let mut j: libc::c_int = 0; + let mut t: libc::c_double = 0.; + let mut u: libc::c_double = 0.; + let mut v: libc::c_double = 0.; + let mut ij: libc::c_int = 0; + let mut tp: libc::c_double = 0.; + let mut beta: libc::c_double = 0.; + let mut gamma_: libc::c_double = 0.; + let mut alpha: libc::c_double = 0.; + let mut delta: libc::c_double = 0.; + w = w.offset(-1); + z__ = z__.offset(-1); + a = a.offset(-1); + if !(*sigma == 0.0f64) { + ij = 1 as libc::c_int; + t = one / *sigma; + if !(*sigma > 0.0f64) { + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + *w.offset(i__ as isize) = *z__.offset(i__ as isize); + i__ += 1; + } + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + v = *w.offset(i__ as isize); + t += v * v / *a.offset(ij as isize); + i__2 = *n; + j = i__ + 1 as libc::c_int; + while j <= i__2 { + ij += 1; + *w.offset(j as isize) -= v * *a.offset(ij as isize); + j += 1; + } + ij += 1; + i__ += 1; + } + if t >= 0.0f64 { + t = epmach / *sigma; + } + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + j = *n + 1 as libc::c_int - i__; + ij -= i__; + u = *w.offset(j as isize); + *w.offset(j as isize) = t; + t -= u * u / *a.offset(ij as isize); + i__ += 1; + } + } + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + v = *z__.offset(i__ as isize); + delta = v / *a.offset(ij as isize); + if *sigma < 0.0f64 { + tp = *w.offset(i__ as isize); + } else { + tp = t + delta * v; + } + alpha = tp / t; + *a.offset(ij as isize) = alpha * *a.offset(ij as isize); + if i__ == *n { + break; + } + beta = delta / tp; + if alpha > four { + gamma_ = t / tp; + i__2 = *n; + j = i__ + 1 as libc::c_int; + while j <= i__2 { + ij += 1; + u = *a.offset(ij as isize); + *a.offset(ij as isize) = gamma_ * u + beta * *z__.offset(j as isize); + *z__.offset(j as isize) -= v * u; + j += 1; + } + } else { + i__2 = *n; + j = i__ + 1 as libc::c_int; + while j <= i__2 { + ij += 1; + *z__.offset(j as isize) -= v * *a.offset(ij as isize); + *a.offset(ij as isize) += beta * *z__.offset(j as isize); + j += 1; + } + } + ij += 1; + t = tp; + i__ += 1; + } + } +} +unsafe extern "C" fn slsqpb_( + mut m: *mut libc::c_int, + mut meq: *mut libc::c_int, + mut la: *mut libc::c_int, + mut n: *mut libc::c_int, + mut x: *mut libc::c_double, + mut xl: *const libc::c_double, + mut xu: *const libc::c_double, + mut f: *mut libc::c_double, + mut c__: *mut libc::c_double, + mut g: *mut libc::c_double, + mut a: *mut libc::c_double, + mut acc: *mut libc::c_double, + mut iter: *mut libc::c_int, + mut mode: *mut libc::c_int, + mut r__: *mut libc::c_double, + mut l: *mut libc::c_double, + mut x0: *mut libc::c_double, + mut mu: *mut libc::c_double, + mut s: *mut libc::c_double, + mut u: *mut libc::c_double, + mut v: *mut libc::c_double, + mut w: *mut libc::c_double, + mut iw: *mut libc::c_int, + mut state: *mut slsqpb_state, +) { + let mut current_block: u64; + let one: libc::c_double = 1.0f64; + let alfmin: libc::c_double = 0.1f64; + let hun: libc::c_double = 100.0f64; + let ten: libc::c_double = 10.0f64; + let two: libc::c_double = 2.0f64; + let mut a_dim1: libc::c_int = 0; + let mut a_offset: libc::c_int = 0; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut d__1: libc::c_double = 0.; + let mut d__2: libc::c_double = 0.; + let mut i__: libc::c_int = 0; + let mut j: libc::c_int = 0; + let mut k: libc::c_int = 0; + let mut t: libc::c_double = 0.; + let mut f0: libc::c_double = 0.; + let mut h1: libc::c_double = 0.; + let mut h2: libc::c_double = 0.; + let mut h3: libc::c_double = 0.; + let mut h4: libc::c_double = 0.; + let mut n1: libc::c_int = 0; + let mut n2: libc::c_int = 0; + let mut n3: libc::c_int = 0; + let mut t0: libc::c_double = 0.; + let mut gs: libc::c_double = 0.; + let mut tol: libc::c_double = 0.; + let mut line: libc::c_int = 0; + let mut alpha: libc::c_double = 0.; + let mut iexact: libc::c_int = 0; + let mut incons: libc::c_int = 0; + let mut ireset: libc::c_int = 0; + let mut itermx: libc::c_int = 0; + t = (*state).t; + f0 = (*state).f0; + h1 = (*state).h1; + h2 = (*state).h2; + h3 = (*state).h3; + h4 = (*state).h4; + n1 = (*state).n1; + n2 = (*state).n2; + n3 = (*state).n3; + t0 = (*state).t0; + gs = (*state).gs; + tol = (*state).tol; + line = (*state).line; + alpha = (*state).alpha; + iexact = (*state).iexact; + incons = (*state).incons; + ireset = (*state).ireset; + itermx = (*state).itermx; + mu = mu.offset(-1); + c__ = c__.offset(-1); + v = v.offset(-1); + u = u.offset(-1); + s = s.offset(-1); + x0 = x0.offset(-1); + l = l.offset(-1); + r__ = r__.offset(-1); + a_dim1 = *la; + a_offset = 1 as libc::c_int + a_dim1; + a = a.offset(-(a_offset as isize)); + g = g.offset(-1); + xu = xu.offset(-1); + xl = xl.offset(-1); + x = x.offset(-1); + w = w.offset(-1); + iw = iw.offset(-1); + if *mode == -(1 as libc::c_int) { + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + *u + .offset( + i__ as isize, + ) = *g.offset(i__ as isize) + - ddot_sl__( + m, + &mut *a.offset((i__ * a_dim1 + 1 as libc::c_int) as isize), + 1 as libc::c_int, + &mut *r__.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ) - *v.offset(i__ as isize); + i__ += 1; + } + k = 0 as libc::c_int; + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + h1 = 0.0f64; + k += 1; + i__2 = *n; + j = i__ + 1 as libc::c_int; + while j <= i__2 { + k += 1; + h1 += *l.offset(k as isize) * *s.offset(j as isize); + j += 1; + } + *v.offset(i__ as isize) = *s.offset(i__ as isize) + h1; + i__ += 1; + } + k = 1 as libc::c_int; + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + *v.offset(i__ as isize) = *l.offset(k as isize) * *v.offset(i__ as isize); + k = k + n1 - i__; + i__ += 1; + } + i__ = *n; + while i__ >= 1 as libc::c_int { + h1 = 0.0f64; + k = i__; + i__1 = i__ - 1 as libc::c_int; + j = 1 as libc::c_int; + while j <= i__1 { + h1 += *l.offset(k as isize) * *v.offset(j as isize); + k = k + *n - j; + j += 1; + } + *v.offset(i__ as isize) += h1; + i__ -= 1; + } + h1 = ddot_sl__( + n, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *u.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + h2 = ddot_sl__( + n, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *v.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + h3 = h2 * 0.2f64; + if h1 < h3 { + h4 = (h2 - h3) / (h2 - h1); + h1 = h3; + dscal_sl__( + n, + &mut h4, + &mut *u.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + d__1 = one - h4; + daxpy_sl__( + n, + &mut d__1, + &mut *v.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *u.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + } + d__1 = one / h1; + ldl_( + n, + &mut *l.offset(1 as libc::c_int as isize), + &mut *u.offset(1 as libc::c_int as isize), + &mut d__1, + &mut *v.offset(1 as libc::c_int as isize), + ); + d__1 = -one / h2; + ldl_( + n, + &mut *l.offset(1 as libc::c_int as isize), + &mut *v.offset(1 as libc::c_int as isize), + &mut d__1, + &mut *u.offset(1 as libc::c_int as isize), + ); + current_block = 17386100730026577537; + } else if *mode == 0 as libc::c_int { + itermx = *iter; + if *acc >= 0.0f64 { + iexact = 0 as libc::c_int; + } else { + iexact = 1 as libc::c_int; + } + *acc = fabs(*acc); + tol = ten * *acc; + *iter = 0 as libc::c_int; + ireset = 0 as libc::c_int; + n1 = *n + 1 as libc::c_int; + n2 = n1 * *n / 2 as libc::c_int; + n3 = n2 + 1 as libc::c_int; + *s.offset(1 as libc::c_int as isize) = 0.0f64; + *mu.offset(1 as libc::c_int as isize) = 0.0f64; + dcopy___( + n, + &mut *s.offset(1 as libc::c_int as isize), + 0 as libc::c_int, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + dcopy___( + m, + &mut *mu.offset(1 as libc::c_int as isize), + 0 as libc::c_int, + &mut *mu.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + current_block = 3938304504335184834; + } else { + t = *f; + i__1 = *m; + j = 1 as libc::c_int; + while j <= i__1 { + if j <= *meq { + h1 = *c__.offset(j as isize); + } else { + h1 = 0.0f64; + } + d__1 = -*c__.offset(j as isize); + t += *mu.offset(j as isize) * (if d__1 >= h1 { d__1 } else { h1 }); + j += 1; + } + h1 = t - t0; + match iexact + 1 as libc::c_int { + 1 => { + current_block = 13298064082901038618; + match current_block { + 13298064082901038618 => { + if nlopt_isfinite(h1) != 0 { + if h1 <= h3 / ten || line > 10 as libc::c_int { + current_block = 778573858929849943; + } else { + d__1 = h3 / (two * (h3 - h1)); + alpha = if d__1 >= alfmin { d__1 } else { alfmin }; + current_block = 2990001543308207002; + } + } else { + alpha = if alpha * 0.5f64 >= alfmin { + alpha * 0.5f64 + } else { + alfmin + }; + current_block = 2990001543308207002; + } + } + _ => {} + } + match current_block { + 2990001543308207002 => {} + _ => { + h3 = 0.0f64; + i__1 = *m; + j = 1 as libc::c_int; + while j <= i__1 { + if j <= *meq { + h1 = *c__.offset(j as isize); + } else { + h1 = 0.0f64; + } + d__1 = -*c__.offset(j as isize); + h3 += if d__1 >= h1 { d__1 } else { h1 }; + j += 1; + } + d__1 = *f - f0; + if (fabs(d__1) < *acc + || dnrm2___( + n, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ) < *acc) && h3 < *acc + { + *mode = 0 as libc::c_int; + } else { + *mode = -(1 as libc::c_int); + } + current_block = 5266944821211112966; + } + } + } + 2 => { + current_block = 17680132173918469312; + } + _ => { + current_block = 778573858929849943; + match current_block { + 13298064082901038618 => { + if nlopt_isfinite(h1) != 0 { + if h1 <= h3 / ten || line > 10 as libc::c_int { + current_block = 778573858929849943; + } else { + d__1 = h3 / (two * (h3 - h1)); + alpha = if d__1 >= alfmin { d__1 } else { alfmin }; + current_block = 2990001543308207002; + } + } else { + alpha = if alpha * 0.5f64 >= alfmin { + alpha * 0.5f64 + } else { + alfmin + }; + current_block = 2990001543308207002; + } + } + _ => {} + } + match current_block { + 2990001543308207002 => {} + _ => { + h3 = 0.0f64; + i__1 = *m; + j = 1 as libc::c_int; + while j <= i__1 { + if j <= *meq { + h1 = *c__.offset(j as isize); + } else { + h1 = 0.0f64; + } + d__1 = -*c__.offset(j as isize); + h3 += if d__1 >= h1 { d__1 } else { h1 }; + j += 1; + } + d__1 = *f - f0; + if (fabs(d__1) < *acc + || dnrm2___( + n, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ) < *acc) && h3 < *acc + { + *mode = 0 as libc::c_int; + } else { + *mode = -(1 as libc::c_int); + } + current_block = 5266944821211112966; + } + } + } + } + } + 'c_2994: loop { + match current_block { + 17680132173918469312 => { + *mode = 9 as libc::c_int; + return; + } + 2990001543308207002 => { + line += 1; + h3 = alpha * h3; + dscal_sl__( + n, + &mut alpha, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + dcopy___( + n, + &mut *x0.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + daxpy_sl__( + n, + &one, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + if *x.offset(i__ as isize) < *xl.offset(i__ as isize) { + *x.offset(i__ as isize) = *xl.offset(i__ as isize); + } else if *x.offset(i__ as isize) > *xu.offset(i__ as isize) { + *x.offset(i__ as isize) = *xu.offset(i__ as isize); + } + i__ += 1; + } + *mode = if line == 1 as libc::c_int { + -(2 as libc::c_int) + } else { + 1 as libc::c_int + }; + current_block = 5266944821211112966; + } + 17386100730026577537 => { + *iter += 1; + *mode = 9 as libc::c_int; + if *iter > itermx && itermx > 0 as libc::c_int { + current_block = 5266944821211112966; + } else { + dcopy___( + n, + &*xl.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *u.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + dcopy___( + n, + &*xu.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *v.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + d__1 = -one; + daxpy_sl__( + n, + &mut d__1, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *u.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + d__1 = -one; + daxpy_sl__( + n, + &mut d__1, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *v.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + h4 = one; + lsq_( + m, + meq, + n, + &mut n3, + la, + &mut *l.offset(1 as libc::c_int as isize), + &mut *g.offset(1 as libc::c_int as isize), + &mut *a.offset(a_offset as isize), + &mut *c__.offset(1 as libc::c_int as isize), + &mut *u.offset(1 as libc::c_int as isize), + &mut *v.offset(1 as libc::c_int as isize), + &mut *s.offset(1 as libc::c_int as isize), + &mut *r__.offset(1 as libc::c_int as isize), + &mut *w.offset(1 as libc::c_int as isize), + &mut *iw.offset(1 as libc::c_int as isize), + mode, + ); + if *mode == 6 as libc::c_int { + if *n == *meq { + *mode = 4 as libc::c_int; + } + } + if *mode == 4 as libc::c_int { + i__1 = *m; + j = 1 as libc::c_int; + while j <= i__1 { + if j <= *meq { + *a + .offset( + (j + n1 * a_dim1) as isize, + ) = -*c__.offset(j as isize); + } else { + d__1 = -*c__.offset(j as isize); + *a + .offset( + (j + n1 * a_dim1) as isize, + ) = if d__1 >= 0.0f64 { d__1 } else { 0.0f64 }; + } + j += 1; + } + *s.offset(1 as libc::c_int as isize) = 0.0f64; + dcopy___( + n, + &mut *s.offset(1 as libc::c_int as isize), + 0 as libc::c_int, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + h3 = 0.0f64; + *g.offset(n1 as isize) = 0.0f64; + *l.offset(n3 as isize) = hun; + *s.offset(n1 as isize) = one; + *u.offset(n1 as isize) = 0.0f64; + *v.offset(n1 as isize) = one; + incons = 0 as libc::c_int; + loop { + lsq_( + m, + meq, + &mut n1, + &mut n3, + la, + &mut *l.offset(1 as libc::c_int as isize), + &mut *g.offset(1 as libc::c_int as isize), + &mut *a.offset(a_offset as isize), + &mut *c__.offset(1 as libc::c_int as isize), + &mut *u.offset(1 as libc::c_int as isize), + &mut *v.offset(1 as libc::c_int as isize), + &mut *s.offset(1 as libc::c_int as isize), + &mut *r__.offset(1 as libc::c_int as isize), + &mut *w.offset(1 as libc::c_int as isize), + &mut *iw.offset(1 as libc::c_int as isize), + mode, + ); + h4 = one - *s.offset(n1 as isize); + if !(*mode == 4 as libc::c_int) { + break; + } + *l.offset(n3 as isize) = ten * *l.offset(n3 as isize); + incons += 1; + if incons > 5 as libc::c_int { + current_block = 5266944821211112966; + continue 'c_2994; + } + } + if *mode != 1 as libc::c_int { + current_block = 5266944821211112966; + continue; + } + } else if *mode != 1 as libc::c_int { + current_block = 5266944821211112966; + continue; + } + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + *v + .offset( + i__ as isize, + ) = *g.offset(i__ as isize) + - ddot_sl__( + m, + &mut *a.offset((i__ * a_dim1 + 1 as libc::c_int) as isize), + 1 as libc::c_int, + &mut *r__.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + i__ += 1; + } + f0 = *f; + dcopy___( + n, + &mut *x.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *x0.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + gs = ddot_sl__( + n, + &mut *g.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + h1 = fabs(gs); + h2 = 0.0f64; + i__1 = *m; + j = 1 as libc::c_int; + while j <= i__1 { + if j <= *meq { + h3 = *c__.offset(j as isize); + } else { + h3 = 0.0f64; + } + d__1 = -*c__.offset(j as isize); + h2 += if d__1 >= h3 { d__1 } else { h3 }; + d__1 = *r__.offset(j as isize); + h3 = fabs(d__1); + d__1 = h3; + d__2 = (*mu.offset(j as isize) + h3) / two; + *mu.offset(j as isize) = if d__1 >= d__2 { d__1 } else { d__2 }; + d__1 = *c__.offset(j as isize); + h1 += h3 * fabs(d__1); + j += 1; + } + *mode = 0 as libc::c_int; + if h1 < *acc && h2 < *acc { + current_block = 5266944821211112966; + continue; + } + h1 = 0.0f64; + i__1 = *m; + j = 1 as libc::c_int; + while j <= i__1 { + if j <= *meq { + h3 = *c__.offset(j as isize); + } else { + h3 = 0.0f64; + } + d__1 = -*c__.offset(j as isize); + h1 + += *mu.offset(j as isize) + * (if d__1 >= h3 { d__1 } else { h3 }); + j += 1; + } + t0 = *f + h1; + h3 = gs - h1 * h4; + *mode = 8 as libc::c_int; + if h3 >= 0.0f64 { + current_block = 3938304504335184834; + continue; + } + line = 0 as libc::c_int; + alpha = one; + if iexact == 1 as libc::c_int { + current_block = 17680132173918469312; + } else { + current_block = 2990001543308207002; + } + } + } + 3938304504335184834 => { + ireset += 1; + if ireset > 5 as libc::c_int { + d__1 = *f - f0; + if (fabs(d__1) < tol + || dnrm2___( + n, + &mut *s.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ) < tol) && h3 < tol + { + *mode = 0 as libc::c_int; + } else { + *mode = 8 as libc::c_int; + } + current_block = 5266944821211112966; + } else { + *l.offset(1 as libc::c_int as isize) = 0.0f64; + dcopy___( + &mut n2, + &mut *l.offset(1 as libc::c_int as isize), + 0 as libc::c_int, + &mut *l.offset(1 as libc::c_int as isize), + 1 as libc::c_int, + ); + j = 1 as libc::c_int; + i__1 = *n; + i__ = 1 as libc::c_int; + while i__ <= i__1 { + *l.offset(j as isize) = one; + j = j + n1 - i__; + i__ += 1; + } + current_block = 17386100730026577537; + } + } + _ => { + (*state).t = t; + (*state).f0 = f0; + (*state).h1 = h1; + (*state).h2 = h2; + (*state).h3 = h3; + (*state).h4 = h4; + (*state).n1 = n1; + (*state).n2 = n2; + (*state).n3 = n3; + (*state).t0 = t0; + (*state).gs = gs; + (*state).tol = tol; + (*state).line = line; + (*state).alpha = alpha; + (*state).iexact = iexact; + (*state).incons = incons; + (*state).ireset = ireset; + (*state).itermx = itermx; + return; + } + } + }; +} +#[no_mangle] +pub unsafe extern "C" fn slsqp( + mut m: *mut libc::c_int, + mut meq: *mut libc::c_int, + mut la: *mut libc::c_int, + mut n: *mut libc::c_int, + mut x: *mut libc::c_double, + mut xl: *const libc::c_double, + mut xu: *const libc::c_double, + mut f: *mut libc::c_double, + mut c__: *mut libc::c_double, + mut g: *mut libc::c_double, + mut a: *mut libc::c_double, + mut acc: *mut libc::c_double, + mut iter: *mut libc::c_int, + mut mode: *mut libc::c_int, + mut w: *mut libc::c_double, + mut l_w__: *mut libc::c_int, + mut jw: *mut libc::c_int, + mut l_jw__: *mut libc::c_int, + mut state: *mut slsqpb_state, +) { + let mut a_dim1: libc::c_int = 0; + let mut a_offset: libc::c_int = 0; + let mut i__1: libc::c_int = 0; + let mut i__2: libc::c_int = 0; + let mut n1: libc::c_int = 0; + let mut il: libc::c_int = 0; + let mut im: libc::c_int = 0; + let mut ir: libc::c_int = 0; + let mut is: libc::c_int = 0; + let mut iu: libc::c_int = 0; + let mut iv: libc::c_int = 0; + let mut iw: libc::c_int = 0; + let mut ix: libc::c_int = 0; + let mut mineq: libc::c_int = 0; + c__ = c__.offset(-1); + a_dim1 = *la; + a_offset = 1 as libc::c_int + a_dim1; + a = a.offset(-(a_offset as isize)); + g = g.offset(-1); + xu = xu.offset(-1); + xl = xl.offset(-1); + x = x.offset(-1); + w = w.offset(-1); + jw = jw.offset(-1); + n1 = *n + 1 as libc::c_int; + mineq = *m - *meq + n1 + n1; + il = (n1 * 3 as libc::c_int + *m) * (n1 + 1 as libc::c_int) + + (n1 - *meq + 1 as libc::c_int) * (mineq + 2 as libc::c_int) + + (mineq << 1 as libc::c_int) + (n1 + mineq) * (n1 - *meq) + + (*meq << 1 as libc::c_int) + n1 * *n / 2 as libc::c_int + + (*m << 1 as libc::c_int) + *n * 3 as libc::c_int + (n1 << 2 as libc::c_int) + + 1 as libc::c_int; + i__1 = mineq; + i__2 = n1 - *meq; + im = if i__1 >= i__2 { i__1 } else { i__2 }; + if *l_w__ < il || *l_jw__ < im { + *mode = (if 10 as libc::c_int >= il { 10 as libc::c_int } else { il }) + * 1000 as libc::c_int; + *mode += if 10 as libc::c_int >= im { 10 as libc::c_int } else { im }; + return; + } + im = 1 as libc::c_int; + il = im + (if 1 as libc::c_int >= *m { 1 as libc::c_int } else { *m }); + il = im + *la; + ix = il + n1 * *n / 2 as libc::c_int + 1 as libc::c_int; + ir = ix + *n; + is = ir + *n + *n + (if 1 as libc::c_int >= *m { 1 as libc::c_int } else { *m }); + is = ir + *n + *n + *la; + iu = is + n1; + iv = iu + n1; + iw = iv + n1; + slsqpb_( + m, + meq, + la, + n, + &mut *x.offset(1 as libc::c_int as isize), + &*xl.offset(1 as libc::c_int as isize), + &*xu.offset(1 as libc::c_int as isize), + f, + &mut *c__.offset(1 as libc::c_int as isize), + &mut *g.offset(1 as libc::c_int as isize), + &mut *a.offset(a_offset as isize), + acc, + iter, + mode, + &mut *w.offset(ir as isize), + &mut *w.offset(il as isize), + &mut *w.offset(ix as isize), + &mut *w.offset(im as isize), + &mut *w.offset(is as isize), + &mut *w.offset(iu as isize), + &mut *w.offset(iv as isize), + &mut *w.offset(iw as isize), + &mut *jw.offset(1 as libc::c_int as isize), + state, + ); + let ref mut fresh0 = (*state).x0; + *fresh0 = &mut *w.offset(ix as isize) as *mut libc::c_double; +} diff --git a/slsqp/stop.c b/slsqp/stop.c new file mode 100644 index 0000000..94b31bc --- /dev/null +++ b/slsqp/stop.c @@ -0,0 +1,266 @@ +/* Copyright (c) 2007-2014 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include +#include +#include "nlopt-util.h" + +/* utility routines to implement the various stopping criteria */ + +static double sc(double x, double smin, double smax) +{ + return smin + x * (smax - smin); +} + +static double vector_norm(unsigned n, const double *vec, const double *w, const double *scale_min, const double *scale_max) +{ + unsigned i; + double ret = 0; + if (scale_min && scale_max) { + if (w) + for (i = 0; i < n; i++) + ret += w[i] * fabs(sc(vec[i], scale_min[i], scale_max[i])); + else + for (i = 0; i < n; i++) + ret += fabs(sc(vec[i], scale_min[i], scale_max[i])); + } else { + if (w) + for (i = 0; i < n; i++) + ret += w[i] * fabs(vec[i]); + else + for (i = 0; i < n; i++) + ret += fabs(vec[i]); + } + return ret; +} + +static double diff_norm(unsigned n, const double *x, const double *oldx, const double *w, const double *scale_min, const double *scale_max) +{ + unsigned i; + double ret = 0; + if (scale_min && scale_max) { + if (w) + for (i = 0; i < n; i++) + ret += w[i] * fabs(sc(x[i], scale_min[i], scale_max[i]) - sc(oldx[i], scale_min[i], scale_max[i])); + else + for (i = 0; i < n; i++) + ret += fabs(sc(x[i], scale_min[i], scale_max[i]) - sc(oldx[i], scale_min[i], scale_max[i])); + } else { + if (w) + for (i = 0; i < n; i++) + ret += w[i] * fabs(x[i] - oldx[i]); + else + for (i = 0; i < n; i++) + ret += fabs(x[i] - oldx[i]); + } + return ret; +} + +static int relstop(double vold, double vnew, double reltol, double abstol) +{ + if (nlopt_isinf(vold)) + return 0; + return (fabs(vnew - vold) < abstol || fabs(vnew - vold) < reltol * (fabs(vnew) + fabs(vold)) * 0.5 || (reltol > 0 && vnew == vold)); /* catch vnew == vold == 0 */ +} + +int nlopt_stop_ftol(const nlopt_stopping * s, double f, double oldf) +{ + return (relstop(oldf, f, s->ftol_rel, s->ftol_abs)); +} + +int nlopt_stop_f(const nlopt_stopping * s, double f, double oldf) +{ + return (f <= s->minf_max || nlopt_stop_ftol(s, f, oldf)); +} + +int nlopt_stop_x(const nlopt_stopping * s, const double *x, const double *oldx) +{ + unsigned i; + if (diff_norm(s->n, x, oldx, s->x_weights, NULL, NULL) < s->xtol_rel * vector_norm(s->n, x, s->x_weights, NULL, NULL)) + return 1; + if (!s->xtol_abs) return 0; + for (i = 0; i < s->n; ++i) + if (fabs(x[i] - oldx[i]) >= s->xtol_abs[i]) + return 0; + return 1; +} + +int nlopt_stop_dx(const nlopt_stopping * s, const double *x, const double *dx) +{ + unsigned i; + if (vector_norm(s->n, dx, s->x_weights, NULL, NULL) < s->xtol_rel * vector_norm(s->n, x, s->x_weights, NULL, NULL)) + return 1; + if (!s->xtol_abs) return 0; + for (i = 0; i < s->n; ++i) + if (fabs(dx[i]) >= s->xtol_abs[i]) + return 0; + return 1; +} + +/* some of the algorithms rescale x to a unit hypercube, so we need to + scale back before we can compare to the tolerances */ +int nlopt_stop_xs(const nlopt_stopping * s, const double *xs, const double *oldxs, const double *scale_min, const double *scale_max) +{ + unsigned i; + if (diff_norm(s->n, xs, oldxs, s->x_weights, scale_min, scale_max) < s->xtol_rel * vector_norm(s->n, xs, s->x_weights, scale_min, scale_max)) + return 1; + if (!s->xtol_abs) return 0; + for (i = 0; i < s->n; ++i) + if (fabs(sc(xs[i], scale_min[i], scale_max[i]) - sc(oldxs[i], scale_min[i], scale_max[i])) >= s->xtol_abs[i]) + return 0; + return 1; +} + +int nlopt_stop_evals(const nlopt_stopping * s) +{ + return (s->maxeval > 0 && *(s->nevals_p) >= s->maxeval); +} + +int nlopt_stop_time_(double start, double maxtime) +{ + return (maxtime > 0 && nlopt_seconds() - start >= maxtime); +} + +int nlopt_stop_time(const nlopt_stopping * s) +{ + return nlopt_stop_time_(s->start, s->maxtime); +} + +int nlopt_stop_evalstime(const nlopt_stopping * stop) +{ + return nlopt_stop_evals(stop) || nlopt_stop_time(stop); +} + +int nlopt_stop_forced(const nlopt_stopping * stop) +{ + return stop->force_stop && *(stop->force_stop); +} + +unsigned nlopt_count_constraints(unsigned p, const nlopt_constraint * c) +{ + unsigned i, count = 0; + for (i = 0; i < p; ++i) + count += c[i].m; + return count; +} + +unsigned nlopt_max_constraint_dim(unsigned p, const nlopt_constraint * c) +{ + unsigned i, max_dim = 0; + for (i = 0; i < p; ++i) + if (c[i].m > max_dim) + max_dim = c[i].m; + return max_dim; +} + +void nlopt_eval_constraint(double *result, double *grad, const nlopt_constraint * c, unsigned n, const double *x) +{ + if (c->f) + result[0] = c->f(n, x, grad, c->f_data); + else + c->mf(c->m, result, n, x, grad, c->f_data); +} + +char *nlopt_vsprintf(char *p, const char *format, va_list ap) +{ + size_t len = strlen(format) + 128; + int ret; + + p = (char *) realloc(p, len); + if (!p) + abort(); + + /* TODO: check HAVE_VSNPRINTF, and fallback to vsprintf otherwise */ + while ((ret = vsnprintf(p, len, format, ap)) < 0 || (size_t) ret >= len) { + /* C99 vsnprintf returns the required number of bytes (excluding \0) + if the buffer is too small; older versions (e.g. MS) return -1 */ + len = ret >= 0 ? (size_t) (ret + 1) : (len * 3) >> 1; + p = (char *) realloc(p, len); + if (!p) + abort(); + } + return p; +} + +void nlopt_stop_msg(const nlopt_stopping * s, const char *format, ...) +{ + va_list ap; + if (s->stop_msg) { + va_start(ap, format); + *(s->stop_msg) = nlopt_vsprintf(*(s->stop_msg), format, ap); + va_end(ap); + } +} + +/*************************************************************************/ + +int nlopt_isinf(double x) +{ + return (fabs(x) >= HUGE_VAL * 0.99) +#if defined(HAVE_ISINF) + || isinf(x) +#else + || (!nlopt_isnan(x) && nlopt_isnan(x - x)) +#endif + ; +} + +int nlopt_isfinite(double x) +{ + return (fabs(x) <= DBL_MAX) +#if defined(HAVE_ISFINITE) + || isfinite(x) +#elif defined(_WIN32) + || _finite(x) +#endif + ; +} + +int nlopt_istiny(double x) +{ + if (x == 0.0) + return 1; + else { +#if defined(HAVE_FPCLASSIFY) + return fpclassify(x) == FP_SUBNORMAL; +#elif defined(_WIN32) + int c = _fpclass(x); + return c == _FPCLASS_ND || c == _FPCLASS_PD; +#else + return fabs(x) < 2.2250738585072014e-308; /* assume IEEE 754 double */ +#endif + } +} + +int nlopt_isnan(double x) +{ +#if defined(HAVE_ISNAN) + return isnan(x); +#elif defined(_WIN32) + return _isnan(x); +#else + return (x != x); /* might fail with aggressive optimization */ +#endif +} diff --git a/slsqp/timer.c b/slsqp/timer.c new file mode 100644 index 0000000..727d98e --- /dev/null +++ b/slsqp/timer.c @@ -0,0 +1,93 @@ +/* Copyright (c) 2007-2014 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include "nlopt-util.h" + +#ifdef TIME_WITH_SYS_TIME +# include +# include +#else +# ifdef HAVE_SYS_TIME_H +# include +# else +# include +# endif +#endif + +#if defined(_WIN32) || defined(__WIN32__) +# include +#endif + +/* return time in seconds since some arbitrary point in the past */ +double nlopt_seconds(void) +{ + static THREADLOCAL int start_inited = 0; /* whether start time has been initialized */ +#if defined(HAVE_GETTIMEOFDAY) + static THREADLOCAL struct timeval start; + struct timeval tv; + if (!start_inited) { + start_inited = 1; + gettimeofday(&start, NULL); + } + gettimeofday(&tv, NULL); + return (tv.tv_sec - start.tv_sec) + 1.e-6 * (tv.tv_usec - start.tv_usec); +#elif defined(HAVE_TIME) + return time(NULL); +#elif defined(_WIN32) || defined(__WIN32__) + static THREADLOCAL ULONGLONG start; + FILETIME ft; + if (!start_inited) { + start_inited = 1; + GetSystemTimeAsFileTime(&ft); + start = (((ULONGLONG) ft.dwHighDateTime) << 32) + ft.dwLowDateTime; + } + GetSystemTimeAsFileTime(&ft); + return 100e-9 * (((((ULONGLONG) ft.dwHighDateTime) << 32) + ft.dwLowDateTime) - start); +#else + /* use clock() as a fallback... this is somewhat annoying + because clock() may wrap around with a fairly short period */ + static THREADLOCAL clock_t start; + if (!start_inited) { + start_inited = 1; + start = clock(); + } + return (clock() - start) * 1.0 / CLOCKS_PER_SEC; +#endif +} + +/* number based on time for use as random seed */ +unsigned long nlopt_time_seed(void) +{ +#if defined(HAVE_GETTIMEOFDAY) + struct timeval tv; + gettimeofday(&tv, NULL); + return (tv.tv_sec ^ tv.tv_usec); +#elif defined(HAVE_TIME) + return time(NULL); +#elif defined(_WIN32) || defined(__WIN32__) + FILETIME ft; + GetSystemTimeAsFileTime(&ft); + return ft.dwHighDateTime ^ ft.dwLowDateTime; +#else + return clock(); +#endif +} diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..e7a11a9 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,3 @@ +fn main() { + println!("Hello, world!"); +}