From 691a7ae55fc82b1e253e323d6741c7daa1e91ebf Mon Sep 17 00:00:00 2001 From: William Joye Date: Fri, 8 Dec 2017 14:39:24 -0500 Subject: upgrade AST --- ast/cminpack/CopyrightMINPACK.txt | 52 ++++ ast/cminpack/README.md | 128 ++++++++++ ast/cminpack/cminpack.h | 370 +++++++++++++++++++++++++++ ast/cminpack/cminpackP.h | 62 +++++ ast/cminpack/dpmpar.c | 201 +++++++++++++++ ast/cminpack/enorm.c | 157 ++++++++++++ ast/cminpack/lmder.c | 526 ++++++++++++++++++++++++++++++++++++++ ast/cminpack/lmder1.c | 167 ++++++++++++ ast/cminpack/lmpar.c | 338 ++++++++++++++++++++++++ ast/cminpack/qrfac.c | 285 +++++++++++++++++++++ ast/cminpack/qrsolv.c | 218 ++++++++++++++++ 11 files changed, 2504 insertions(+) create mode 100644 ast/cminpack/CopyrightMINPACK.txt create mode 100644 ast/cminpack/README.md create mode 100644 ast/cminpack/cminpack.h create mode 100644 ast/cminpack/cminpackP.h create mode 100644 ast/cminpack/dpmpar.c create mode 100644 ast/cminpack/enorm.c create mode 100644 ast/cminpack/lmder.c create mode 100644 ast/cminpack/lmder1.c create mode 100644 ast/cminpack/lmpar.c create mode 100644 ast/cminpack/qrfac.c create mode 100644 ast/cminpack/qrsolv.c diff --git a/ast/cminpack/CopyrightMINPACK.txt b/ast/cminpack/CopyrightMINPACK.txt new file mode 100644 index 0000000..ae7984d --- /dev/null +++ b/ast/cminpack/CopyrightMINPACK.txt @@ -0,0 +1,52 @@ +Minpack Copyright Notice (1999) University of Chicago. All rights reserved + +Redistribution and use in source and binary forms, with or +without modification, are permitted provided that the +following conditions are met: + +1. Redistributions of source code must retain the above +copyright notice, this list of conditions and the following +disclaimer. + +2. Redistributions in binary form must reproduce the above +copyright notice, this list of conditions and the following +disclaimer in the documentation and/or other materials +provided with the distribution. + +3. The end-user documentation included with the +redistribution, if any, must include the following +acknowledgment: + + "This product includes software developed by the + University of Chicago, as Operator of Argonne National + Laboratory. + +Alternately, this acknowledgment may appear in the software +itself, if and wherever such third-party acknowledgments +normally appear. + +4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" +WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE +UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND +THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES +OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE +OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY +OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR +USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF +THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) +DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION +UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL +BE CORRECTED. + +5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT +HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF +ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, +INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF +ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF +PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER +SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT +(INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, +EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE +POSSIBILITY OF SUCH LOSS OR DAMAGES. + diff --git a/ast/cminpack/README.md b/ast/cminpack/README.md new file mode 100644 index 0000000..66bbc6f --- /dev/null +++ b/ast/cminpack/README.md @@ -0,0 +1,128 @@ +C/C++ Minpack [![Build Status](https://api.travis-ci.org/devernay/cminpack.png?branch=master)](https://travis-ci.org/devernay/cminpack) +========== + +This is a C version of the minpack minimization package. +It has been derived from the fortran code using f2c and +some limited manual editing. Note that you need to link +against libf2c to use this version of minpack. Extern "C" +linkage permits the package routines to be called from C++. +Check ftp://netlib.bell-labs.com/netlib/f2c for the latest +f2c version. For general minpack info and test programs, see +the accompanying readme.txt and http://www.netlib.org/minpack/. + +Type `make` to compile and `make install` to install in /usr/local +or modify the makefile to suit your needs. + +This software has been tested on a RedHat 7.3 Linux machine - +usual 'use at your own risk' warnings apply. + +Manolis Lourakis -- lourakis at ics forth gr, July 2002 + Institute of Computer Science, + Foundation for Research and Technology - Hellas + Heraklion, Crete, Greece + +Repackaging by Frederic Devernay -- frederic dot devernay at m4x dot org + +The project home page is at http://devernay.free.fr/hacks/cminpack/ + +History +------ + +* version 1.3.3 (04/02/2014):: + - Add documentation and examples abouts how to add box constraints to the variables. + - continuous integration https://travis-ci.org/devernay/cminpack + +* version 1.3.2 (27/10/2013): + - Minor change in the CMake build: also set SOVERSION. + +* version 1.3.1 (02/10/2013): + - Fix CUDA examples compilation, and remove non-free files. + +* version 1.3.0 (09/06/2012): + - Optionally use LAPACK and CBLAS in lmpar, qrfac, and qrsolv. Added + "make lapack" to build the LAPACK-based cminpack and "make + checklapack" to test it (results of the test may depend on the + underlying LAPACK and BLAS implementations). + On 64-bits architectures, the preprocessor symbol __LP64__ must be + defined (see cminpackP.h) if the LAPACK library uses the LP64 + interface (i.e. 32-bits integer, vhereas the ILP interface uses 64 + bits integers). + +* version 1.2.2 (16/05/2012): + - Update Makefiles and documentation (see "Using CMinpack" above) for + easier building and testing. + +* version 1.2.1 (15/05/2012): + - The library can now be built as double, float or half + versions. Standard tests in the "examples" directory can now be + lauched using "make check" (to run common tests, including against + the float version), "make checkhalf" (to test the half version) and + "make checkfail" (to run all the tests, even those that fail). + +* version 1.2.0 (14/05/2012): +- Added original FORTRAN sources for better testing (type "make" in + directory fortran, then "make" in examples and follow the + instructions). Added driver tests lmsdrv, chkdrv, hyjdrv, + hybdrv. Typing "make alltest" in the examples directory will run all + possible test combinations (make sure you have gfortran installed). + +* version 1.1.5 (04/05/2012): + - cminpack now works in CUDA, thanks to Jordi Bataller Mascarell, type + "make" in the "cuda" subdir (be careful, though: this is a + straightforward port from C, and each problem is solved using a + single thread). cminpack can now also be compiled with + single-precision floating point computation (define + __cminpack_real__ to float when compiling and using the + library). Fix cmake support for CMINPACK_LIB_INSTALL_DIR. Update the + reference files for tests. + +* version 1.1.4 (30/10/2011): + - Translated all the Levenberg-Marquardt code (lmder, lmdif, lmstr, + lmder1, lmdif1, lmstr1, lmpar, qrfac, qrsolv, fdjac2, chkder) to use + C-style indices. + +* version 1.1.3 (16/03/2011): + - Minor fix: Change non-standard strnstr() to strstr() in + genf77tests.c. + +* version 1.1.2 (07/01/2011): + - Fix Windows DLL building (David Graeff) and document covar in + cminpack.h. + +* version 1.1.1 (04/12/2010): + - Complete rewrite of the C functions (without trailing underscore in + the function name). Using the original FORTRAN code, the original + algorithms structure was recovered, and many goto's were converted + to if...then...else. The code should now be both more readable and + easier to optimize, both for humans and for compilers. Added lmddrv + and lmfdrv test drivers, which test a lot of difficult functions + (these functions are explained in Testing Unconstrained Optimization + Software by Moré et al.). Also added the pkg-config files to the + cmake build, as well as an "uninstall" target, contributed by + Geoffrey Biggs. + +* version 1.0.4 (18/10/2010): + - Support for shared library building using CMake, thanks to Goeffrey + Biggs and Radu Bogdan Rusu from Willow Garage. Shared libraries can be + enabled using cmake options, as in; + cmake -DUSE_FPIC=ON -DSHARED_LIBS=ON -DBUILD_EXAMPLES=OFF path_to_sources + +* version 1.0.3 (18/03/2010): + - Added CMake support. + - XCode build is now Universal. + - Added tfdjac2_ and tfdjac2c examples, which test the accuracy of a + finite-differences approximation of the Jacobian. + - Bug fix in tlmstr1 (signaled by Thomas Capricelli). + +* version 1.0.2 (27/02/2009): + - Added Xcode and Visual Studio project files + +* version 1.0.1 (17/12/2007): + - bug fix in covar() and covar_(), the computation of tolr caused a + segfault (signaled by Timo Hartmann). + +* version 1.0.0 (24/04/2007): + - Added fortran and C examples + - Added documentation from Debian man pages + - Wrote pure C version + - Added covar() and covar_(), and use it in tlmdef/tlmdif diff --git a/ast/cminpack/cminpack.h b/ast/cminpack/cminpack.h new file mode 100644 index 0000000..6d3f757 --- /dev/null +++ b/ast/cminpack/cminpack.h @@ -0,0 +1,370 @@ +/* Header file for cminpack, by Frederic Devernay. + The documentation for all functions can be found in the file + minpack-documentation.txt from the distribution, or in the source + code of each function. */ + +#ifndef __CMINPACK_H__ +#define __CMINPACK_H__ + +/* The default floating-point type is "double" for C/C++ and "float" for CUDA, + but you can change this by defining one of the following symbols when + compiling the library, and before including cminpack.h when using it: + __cminpack_double__ for double + __cminpack_float__ for float + __cminpack_half__ for half from the OpenEXR library (in this case, you must + compile cminpack with a C++ compiler) +*/ +#ifdef __cminpack_double__ +#define __cminpack_real__ double +#endif + +#ifdef __cminpack_float__ +#define __cminpack_real__ float +#endif + +#ifdef __cminpack_half__ +#include +#define __cminpack_real__ half +#endif + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +/* Cmake will define cminpack_EXPORTS on Windows when it +configures to build a shared library. If you are going to use +another build system on windows or create the visual studio +projects by hand you need to define cminpack_EXPORTS when +building a DLL on windows. +*/ +#if defined (__GNUC__) +#define CMINPACK_DECLSPEC_EXPORT __declspec(__dllexport__) +#define CMINPACK_DECLSPEC_IMPORT __declspec(__dllimport__) +#endif +#if defined (_MSC_VER) || defined (__BORLANDC__) +#define CMINPACK_DECLSPEC_EXPORT __declspec(dllexport) +#define CMINPACK_DECLSPEC_IMPORT __declspec(dllimport) +#endif +#ifdef __WATCOMC__ +#define CMINPACK_DECLSPEC_EXPORT __export +#define CMINPACK_DECLSPEC_IMPORT __import +#endif +#ifdef __IBMC__ +#define CMINPACK_DECLSPEC_EXPORT _Export +#define CMINPACK_DECLSPEC_IMPORT _Import +#endif + +#if !defined(CMINPACK_NO_DLL) && (defined(__WIN32__) || defined(WIN32) || defined (_WIN32)) +#if defined(cminpack_EXPORTS) || defined(CMINPACK_EXPORTS) || defined(CMINPACK_DLL_EXPORTS) + #define CMINPACK_EXPORT CMINPACK_DECLSPEC_EXPORT + #else + #define CMINPACK_EXPORT CMINPACK_DECLSPEC_IMPORT + #endif /* cminpack_EXPORTS */ +#else /* defined (_WIN32) */ + #define CMINPACK_EXPORT +#endif + +#if defined(__CUDA_ARCH__) || defined(__CUDACC__) +#define __cminpack_attr__ __device__ +#ifndef __cminpack_real__ +#define __cminpack_float__ +#define __cminpack_real__ float +#endif +#define __cminpack_type_fcn_nn__ __cminpack_attr__ int fcn_nn +#define __cminpack_type_fcnder_nn__ __cminpack_attr__ int fcnder_nn +#define __cminpack_type_fcn_mn__ __cminpack_attr__ int fcn_mn +#define __cminpack_type_fcnder_mn__ __cminpack_attr__ int fcnder_mn +#define __cminpack_type_fcnderstr_mn__ __cminpack_attr__ int fcnderstr_mn +#define __cminpack_decl_fcn_nn__ +#define __cminpack_decl_fcnder_nn__ +#define __cminpack_decl_fcn_mn__ +#define __cminpack_decl_fcnder_mn__ +#define __cminpack_decl_fcnderstr_mn__ +#define __cminpack_param_fcn_nn__ +#define __cminpack_param_fcnder_nn__ +#define __cminpack_param_fcn_mn__ +#define __cminpack_param_fcnder_mn__ +#define __cminpack_param_fcnderstr_mn__ +#else +#define __cminpack_attr__ +#ifndef __cminpack_real__ +#define __cminpack_double__ +#define __cminpack_real__ double +#endif +#define __cminpack_type_fcn_nn__ typedef int (*cminpack_func_nn) +#define __cminpack_type_fcnder_nn__ typedef int (*cminpack_funcder_nn) +#define __cminpack_type_fcn_mn__ typedef int (*cminpack_func_mn) +#define __cminpack_type_fcnder_mn__ typedef int (*cminpack_funcder_mn) +#define __cminpack_type_fcnderstr_mn__ typedef int (*cminpack_funcderstr_mn) +#define __cminpack_decl_fcn_nn__ cminpack_func_nn fcn_nn, +#define __cminpack_decl_fcnder_nn__ cminpack_funcder_nn fcnder_nn, +#define __cminpack_decl_fcn_mn__ cminpack_func_mn fcn_mn, +#define __cminpack_decl_fcnder_mn__ cminpack_funcder_mn fcnder_mn, +#define __cminpack_decl_fcnderstr_mn__ cminpack_funcderstr_mn fcnderstr_mn, +#define __cminpack_param_fcn_nn__ fcn_nn, +#define __cminpack_param_fcnder_nn__ fcnder_nn, +#define __cminpack_param_fcn_mn__ fcn_mn, +#define __cminpack_param_fcnder_mn__ fcnder_mn, +#define __cminpack_param_fcnderstr_mn__ fcnderstr_mn, +#endif + +#ifdef __cminpack_double__ +#define __cminpack_func__(func) func +#endif + +#ifdef __cminpack_float__ +#define __cminpack_func__(func) s ## func +#endif + +#ifdef __cminpack_half__ +#define __cminpack_func__(func) h ## func +#endif + +/* Declarations for minpack */ + +/* Function types: */ +/* The first argument can be used to store extra function parameters, thus */ +/* avoiding the use of global variables. */ +/* the iflag parameter is input-only (with respect to the FORTRAN */ +/* version), the output iflag value is the return value of the function. */ +/* If iflag=0, the function shoulkd just print the current values (see */ +/* the nprint parameters below). */ + +/* for hybrd1 and hybrd: */ +/* calculate the functions at x and */ +/* return this vector in fvec. */ +/* return a negative value to terminate hybrd1/hybrd */ +__cminpack_type_fcn_nn__(void *p, int n, const __cminpack_real__ *x, __cminpack_real__ *fvec, int iflag ); + +/* for hybrj1 and hybrj */ +/* if iflag = 1 calculate the functions at x and */ +/* return this vector in fvec. do not alter fjac. */ +/* if iflag = 2 calculate the jacobian at x and */ +/* return this matrix in fjac. do not alter fvec. */ +/* return a negative value to terminate hybrj1/hybrj */ +__cminpack_type_fcnder_nn__(void *p, int n, const __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ *fjac, + int ldfjac, int iflag ); + +/* for lmdif1 and lmdif */ +/* calculate the functions at x and */ +/* return this vector in fvec. */ +/* if iflag = 1 the result is used to compute the residuals. */ +/* if iflag = 2 the result is used to compute the Jacobian by finite differences. */ +/* Jacobian computation requires exactly n function calls with iflag = 2. */ +/* return a negative value to terminate lmdif1/lmdif */ +__cminpack_type_fcn_mn__(void *p, int m, int n, const __cminpack_real__ *x, __cminpack_real__ *fvec, + int iflag ); + +/* for lmder1 and lmder */ +/* if iflag = 1 calculate the functions at x and */ +/* return this vector in fvec. do not alter fjac. */ +/* if iflag = 2 calculate the jacobian at x and */ +/* return this matrix in fjac. do not alter fvec. */ +/* return a negative value to terminate lmder1/lmder */ +__cminpack_type_fcnder_mn__(void *p, int m, int n, const __cminpack_real__ *x, __cminpack_real__ *fvec, + __cminpack_real__ *fjac, int ldfjac, int iflag ); + +/* for lmstr1 and lmstr */ +/* if iflag = 1 calculate the functions at x and */ +/* return this vector in fvec. */ +/* if iflag = i calculate the (i-1)-st row of the */ +/* jacobian at x and return this vector in fjrow. */ +/* return a negative value to terminate lmstr1/lmstr */ +__cminpack_type_fcnderstr_mn__(void *p, int m, int n, const __cminpack_real__ *x, __cminpack_real__ *fvec, + __cminpack_real__ *fjrow, int iflag ); + + + + + + +/* MINPACK functions: */ +/* the info parameter was removed from most functions: the return */ +/* value of the function is used instead. */ +/* The argument 'p' can be used to store extra function parameters, thus */ +/* avoiding the use of global variables. You can also think of it as a */ +/* 'this' pointer a la C++. */ + +/* find a zero of a system of N nonlinear functions in N variables by + a modification of the Powell hybrid method (Jacobian calculated by + a forward-difference approximation) */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(hybrd1)( __cminpack_decl_fcn_nn__ + void *p, int n, __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ tol, + __cminpack_real__ *wa, int lwa ); + +/* find a zero of a system of N nonlinear functions in N variables by + a modification of the Powell hybrid method (Jacobian calculated by + a forward-difference approximation, more general). */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(hybrd)( __cminpack_decl_fcn_nn__ + void *p, int n, __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ xtol, int maxfev, + int ml, int mu, __cminpack_real__ epsfcn, __cminpack_real__ *diag, int mode, + __cminpack_real__ factor, int nprint, int *nfev, + __cminpack_real__ *fjac, int ldfjac, __cminpack_real__ *r, int lr, __cminpack_real__ *qtf, + __cminpack_real__ *wa1, __cminpack_real__ *wa2, __cminpack_real__ *wa3, __cminpack_real__ *wa4); + +/* find a zero of a system of N nonlinear functions in N variables by + a modification of the Powell hybrid method (user-supplied Jacobian) */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(hybrj1)( __cminpack_decl_fcnder_nn__ void *p, int n, __cminpack_real__ *x, + __cminpack_real__ *fvec, __cminpack_real__ *fjac, int ldfjac, __cminpack_real__ tol, + __cminpack_real__ *wa, int lwa ); + +/* find a zero of a system of N nonlinear functions in N variables by + a modification of the Powell hybrid method (user-supplied Jacobian, + more general) */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(hybrj)( __cminpack_decl_fcnder_nn__ void *p, int n, __cminpack_real__ *x, + __cminpack_real__ *fvec, __cminpack_real__ *fjac, int ldfjac, __cminpack_real__ xtol, + int maxfev, __cminpack_real__ *diag, int mode, __cminpack_real__ factor, + int nprint, int *nfev, int *njev, __cminpack_real__ *r, + int lr, __cminpack_real__ *qtf, __cminpack_real__ *wa1, __cminpack_real__ *wa2, + __cminpack_real__ *wa3, __cminpack_real__ *wa4 ); + +/* minimize the sum of the squares of nonlinear functions in N + variables by a modification of the Levenberg-Marquardt algorithm + (Jacobian calculated by a forward-difference approximation) */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(lmdif1)( __cminpack_decl_fcn_mn__ + void *p, int m, int n, __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ tol, + int *iwa, __cminpack_real__ *wa, int lwa ); + +/* minimize the sum of the squares of nonlinear functions in N + variables by a modification of the Levenberg-Marquardt algorithm + (Jacobian calculated by a forward-difference approximation, more + general) */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(lmdif)( __cminpack_decl_fcn_mn__ + void *p, int m, int n, __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ ftol, + __cminpack_real__ xtol, __cminpack_real__ gtol, int maxfev, __cminpack_real__ epsfcn, + __cminpack_real__ *diag, int mode, __cminpack_real__ factor, int nprint, + int *nfev, __cminpack_real__ *fjac, int ldfjac, int *ipvt, + __cminpack_real__ *qtf, __cminpack_real__ *wa1, __cminpack_real__ *wa2, __cminpack_real__ *wa3, + __cminpack_real__ *wa4 ); + +/* minimize the sum of the squares of nonlinear functions in N + variables by a modification of the Levenberg-Marquardt algorithm + (user-supplied Jacobian) */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(lmder1)( __cminpack_decl_fcnder_mn__ + void *p, int m, int n, __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ *fjac, + int ldfjac, __cminpack_real__ tol, int *ipvt, + __cminpack_real__ *wa, int lwa ); + +/* minimize the sum of the squares of nonlinear functions in N + variables by a modification of the Levenberg-Marquardt algorithm + (user-supplied Jacobian, more general) */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(lmder)( __cminpack_decl_fcnder_mn__ + void *p, int m, int n, __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ *fjac, + int ldfjac, __cminpack_real__ ftol, __cminpack_real__ xtol, __cminpack_real__ gtol, + int maxfev, __cminpack_real__ *diag, int mode, __cminpack_real__ factor, + int nprint, int *nfev, int *njev, int *ipvt, + __cminpack_real__ *qtf, __cminpack_real__ *wa1, __cminpack_real__ *wa2, __cminpack_real__ *wa3, + __cminpack_real__ *wa4 ); + +/* minimize the sum of the squares of nonlinear functions in N + variables by a modification of the Levenberg-Marquardt algorithm + (user-supplied Jacobian, minimal storage) */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(lmstr1)( __cminpack_decl_fcnderstr_mn__ void *p, int m, int n, + __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ *fjac, int ldfjac, + __cminpack_real__ tol, int *ipvt, __cminpack_real__ *wa, int lwa ); + +/* minimize the sum of the squares of nonlinear functions in N + variables by a modification of the Levenberg-Marquardt algorithm + (user-supplied Jacobian, minimal storage, more general) */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(lmstr)( __cminpack_decl_fcnderstr_mn__ void *p, int m, + int n, __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ *fjac, + int ldfjac, __cminpack_real__ ftol, __cminpack_real__ xtol, __cminpack_real__ gtol, + int maxfev, __cminpack_real__ *diag, int mode, __cminpack_real__ factor, + int nprint, int *nfev, int *njev, int *ipvt, + __cminpack_real__ *qtf, __cminpack_real__ *wa1, __cminpack_real__ *wa2, __cminpack_real__ *wa3, + __cminpack_real__ *wa4 ); + +__cminpack_attr__ +void CMINPACK_EXPORT __cminpack_func__(chkder)( int m, int n, const __cminpack_real__ *x, __cminpack_real__ *fvec, __cminpack_real__ *fjac, + int ldfjac, __cminpack_real__ *xp, __cminpack_real__ *fvecp, int mode, + __cminpack_real__ *err ); + +__cminpack_attr__ +__cminpack_real__ CMINPACK_EXPORT __cminpack_func__(dpmpar)( int i ); + +__cminpack_attr__ +__cminpack_real__ CMINPACK_EXPORT __cminpack_func__(enorm)( int n, const __cminpack_real__ *x ); + +/* compute a forward-difference approximation to the m by n jacobian + matrix associated with a specified problem of m functions in n + variables. */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(fdjac2)(__cminpack_decl_fcn_mn__ + void *p, int m, int n, __cminpack_real__ *x, const __cminpack_real__ *fvec, __cminpack_real__ *fjac, + int ldfjac, __cminpack_real__ epsfcn, __cminpack_real__ *wa); + +/* compute a forward-difference approximation to the n by n jacobian + matrix associated with a specified problem of n functions in n + variables. if the jacobian has a banded form, then function + evaluations are saved by only approximating the nonzero terms. */ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(fdjac1)(__cminpack_decl_fcn_nn__ + void *p, int n, __cminpack_real__ *x, const __cminpack_real__ *fvec, __cminpack_real__ *fjac, int ldfjac, + int ml, int mu, __cminpack_real__ epsfcn, __cminpack_real__ *wa1, + __cminpack_real__ *wa2); + +/* compute inverse(JtJ) after a run of lmdif or lmder. The covariance matrix is obtained + by scaling the result by enorm(y)**2/(m-n). If JtJ is singular and k = rank(J), the + pseudo-inverse is computed, and the result has to be scaled by enorm(y)**2/(m-k). */ +__cminpack_attr__ +void CMINPACK_EXPORT __cminpack_func__(covar)(int n, __cminpack_real__ *r, int ldr, + const int *ipvt, __cminpack_real__ tol, __cminpack_real__ *wa); + +/* covar1 estimates the variance-covariance matrix: + C = sigma**2 (JtJ)**+ + where (JtJ)**+ is the inverse of JtJ or the pseudo-inverse of JtJ (in case J does not have full rank), + and sigma**2 = fsumsq / (m - k) + where fsumsq is the residual sum of squares and k is the rank of J. + The function returns 0 if J has full rank, else the rank of J. +*/ +__cminpack_attr__ +int CMINPACK_EXPORT __cminpack_func__(covar1)(int m, int n, __cminpack_real__ fsumsq, __cminpack_real__ *r, int ldr, + const int *ipvt, __cminpack_real__ tol, __cminpack_real__ *wa); + +/* internal MINPACK subroutines */ +__cminpack_attr__ +void __cminpack_func__(dogleg)(int n, const __cminpack_real__ *r, int lr, + const __cminpack_real__ *diag, const __cminpack_real__ *qtb, __cminpack_real__ delta, __cminpack_real__ *x, + __cminpack_real__ *wa1, __cminpack_real__ *wa2); +__cminpack_attr__ +void __cminpack_func__(qrfac)(int m, int n, __cminpack_real__ *a, int + lda, int pivot, int *ipvt, int lipvt, __cminpack_real__ *rdiag, + __cminpack_real__ *acnorm, __cminpack_real__ *wa); +__cminpack_attr__ +void __cminpack_func__(qrsolv)(int n, __cminpack_real__ *r, int ldr, + const int *ipvt, const __cminpack_real__ *diag, const __cminpack_real__ *qtb, __cminpack_real__ *x, + __cminpack_real__ *sdiag, __cminpack_real__ *wa); +__cminpack_attr__ +void __cminpack_func__(qform)(int m, int n, __cminpack_real__ *q, int + ldq, __cminpack_real__ *wa); +__cminpack_attr__ +void __cminpack_func__(r1updt)(int m, int n, __cminpack_real__ *s, int + ls, const __cminpack_real__ *u, __cminpack_real__ *v, __cminpack_real__ *w, int *sing); +__cminpack_attr__ +void __cminpack_func__(r1mpyq)(int m, int n, __cminpack_real__ *a, int + lda, const __cminpack_real__ *v, const __cminpack_real__ *w); +__cminpack_attr__ +void __cminpack_func__(lmpar)(int n, __cminpack_real__ *r, int ldr, + const int *ipvt, const __cminpack_real__ *diag, const __cminpack_real__ *qtb, __cminpack_real__ delta, + __cminpack_real__ *par, __cminpack_real__ *x, __cminpack_real__ *sdiag, __cminpack_real__ *wa1, + __cminpack_real__ *wa2); +__cminpack_attr__ +void __cminpack_func__(rwupdt)(int n, __cminpack_real__ *r, int ldr, + const __cminpack_real__ *w, __cminpack_real__ *b, __cminpack_real__ *alpha, __cminpack_real__ *cos, + __cminpack_real__ *sin); +#ifdef __cplusplus +} +#endif /* __cplusplus */ + + +#endif /* __CMINPACK_H__ */ diff --git a/ast/cminpack/cminpackP.h b/ast/cminpack/cminpackP.h new file mode 100644 index 0000000..4e8ba7b --- /dev/null +++ b/ast/cminpack/cminpackP.h @@ -0,0 +1,62 @@ +/* Internal header file for cminpack, by Frederic Devernay. */ +#ifndef __CMINPACKP_H__ +#define __CMINPACKP_H__ + +#ifndef __CMINPACK_H__ +#error "cminpackP.h in an internal cminpack header, and must be included after all other headers (including cminpack.h)" +#endif + +#if (defined (USE_CBLAS) || defined (USE_LAPACK)) && !defined (__cminpack_double__) +#error "cminpack can use cblas and lapack only in double precision mode" +#endif + +#ifdef USE_CBLAS +#ifdef __APPLE__ +#include +#else +#include +#endif +#define __cminpack_enorm__(n,x) cblas_dnrm2(n,x,1) +#else +#define __cminpack_enorm__(n,x) __cminpack_func__(enorm)(n,x) +#endif + +#ifdef USE_LAPACK +#ifdef __APPLE__ +#include +#else +#if defined(__LP64__) /* In LP64 match sizes with the 32 bit ABI */ +typedef int __CLPK_integer; +typedef int __CLPK_logical; +typedef float __CLPK_real; +typedef double __CLPK_doublereal; +typedef __CLPK_logical (*__CLPK_L_fp)(); +typedef int __CLPK_ftnlen; +#else +typedef long int __CLPK_integer; +typedef long int __CLPK_logical; +typedef float __CLPK_real; +typedef double __CLPK_doublereal; +typedef __CLPK_logical (*__CLPK_L_fp)(); +typedef long int __CLPK_ftnlen; +#endif +//extern void dlartg_(double *f, double *g, double *cs, double *sn, double *r__); +int dlartg_(__CLPK_doublereal *f, __CLPK_doublereal *g, __CLPK_doublereal *cs, + __CLPK_doublereal *sn, __CLPK_doublereal *r__) +//extern void dgeqp3_(int *m, int *n, double *a, int *lda, int *jpvt, double *tau, double *work, int *lwork, int *info); +int dgeqp3_(__CLPK_integer *m, __CLPK_integer *n, __CLPK_doublereal *a, __CLPK_integer * + lda, __CLPK_integer *jpvt, __CLPK_doublereal *tau, __CLPK_doublereal *work, __CLPK_integer *lwork, + __CLPK_integer *info) +//extern void dgeqrf_(int *m, int *n, double *a, int *lda, double *tau, double *work, int *lwork, int *info); +int dgeqrf_(__CLPK_integer *m, __CLPK_integer *n, __CLPK_doublereal *a, __CLPK_integer * + lda, __CLPK_doublereal *tau, __CLPK_doublereal *work, __CLPK_integer *lwork, __CLPK_integer *info) +#endif +#endif + +#define real __cminpack_real__ +#define min(a,b) ((a) <= (b) ? (a) : (b)) +#define max(a,b) ((a) >= (b) ? (a) : (b)) +#define TRUE_ (1) +#define FALSE_ (0) + +#endif /* !__CMINPACKP_H__ */ diff --git a/ast/cminpack/dpmpar.c b/ast/cminpack/dpmpar.c new file mode 100644 index 0000000..81c6fcd --- /dev/null +++ b/ast/cminpack/dpmpar.c @@ -0,0 +1,201 @@ +#include "cminpack.h" +#include +#include "cminpackP.h" + +#define double_EPSILON DBL_EPSILON +#define double_MIN DBL_MIN +#define double_MAX DBL_MAX +#define float_EPSILON FLT_EPSILON +#define float_MIN FLT_MIN +#define float_MAX FLT_MAX +#define half_EPSILON HALF_EPSILON +#define half_MIN HALF_NRM_MIN +#define half_MAX HALF_MAX + +#define DPMPAR(type,X) _DPMPAR(type,X) +#define _DPMPAR(type,X) type ## _ ## X + +__cminpack_attr__ +real __cminpack_func__(dpmpar)(int i) +{ +/* ********** */ + +/* Function dpmpar */ + +/* This function provides double precision machine parameters */ +/* when the appropriate set of data statements is activated (by */ +/* removing the c from column 1) and all other data statements are */ +/* rendered inactive. Most of the parameter values were obtained */ +/* from the corresponding Bell Laboratories Port Library function. */ + +/* The function statement is */ + +/* double precision function dpmpar(i) */ + +/* where */ + +/* i is an integer input variable set to 1, 2, or 3 which */ +/* selects the desired machine parameter. If the machine has */ +/* t base b digits and its smallest and largest exponents are */ +/* emin and emax, respectively, then these parameters are */ + +/* dpmpar(1) = b**(1 - t), the machine precision, */ + +/* dpmpar(2) = b**(emin - 1), the smallest magnitude, */ + +/* dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude. */ + +/* Argonne National Laboratory. MINPACK Project. November 1996. */ +/* Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More' */ + +/* ********** */ + +/* Machine constants for the IBM 360/370 series, */ +/* the Amdahl 470/V6, the ICL 2900, the Itel AS/6, */ +/* the Xerox Sigma 5/7/9 and the Sel systems 85/86. */ + +/* data mcheps(1),mcheps(2) / z34100000, z00000000 / */ +/* data minmag(1),minmag(2) / z00100000, z00000000 / */ +/* data maxmag(1),maxmag(2) / z7fffffff, zffffffff / */ + +/* Machine constants for the Honeywell 600/6000 series. */ + +/* data mcheps(1),mcheps(2) / o606400000000, o000000000000 / */ +/* data minmag(1),minmag(2) / o402400000000, o000000000000 / */ +/* data maxmag(1),maxmag(2) / o376777777777, o777777777777 / */ + +/* Machine constants for the CDC 6000/7000 series. */ + +/* data mcheps(1) / 15614000000000000000b / */ +/* data mcheps(2) / 15010000000000000000b / */ + +/* data minmag(1) / 00604000000000000000b / */ +/* data minmag(2) / 00000000000000000000b / */ + +/* data maxmag(1) / 37767777777777777777b / */ +/* data maxmag(2) / 37167777777777777777b / */ + +/* Machine constants for the PDP-10 (KA processor). */ + +/* data mcheps(1),mcheps(2) / "114400000000, "000000000000 / */ +/* data minmag(1),minmag(2) / "033400000000, "000000000000 / */ +/* data maxmag(1),maxmag(2) / "377777777777, "344777777777 / */ + +/* Machine constants for the PDP-10 (KI processor). */ + +/* data mcheps(1),mcheps(2) / "104400000000, "000000000000 / */ +/* data minmag(1),minmag(2) / "000400000000, "000000000000 / */ +/* data maxmag(1),maxmag(2) / "377777777777, "377777777777 / */ + +/* Machine constants for the PDP-11. */ + +/* data mcheps(1),mcheps(2) / 9472, 0 / */ +/* data mcheps(3),mcheps(4) / 0, 0 / */ + +/* data minmag(1),minmag(2) / 128, 0 / */ +/* data minmag(3),minmag(4) / 0, 0 / */ + +/* data maxmag(1),maxmag(2) / 32767, -1 / */ +/* data maxmag(3),maxmag(4) / -1, -1 / */ + +/* Machine constants for the Burroughs 6700/7700 systems. */ + +/* data mcheps(1) / o1451000000000000 / */ +/* data mcheps(2) / o0000000000000000 / */ + +/* data minmag(1) / o1771000000000000 / */ +/* data minmag(2) / o7770000000000000 / */ + +/* data maxmag(1) / o0777777777777777 / */ +/* data maxmag(2) / o7777777777777777 / */ + +/* Machine constants for the Burroughs 5700 system. */ + +/* data mcheps(1) / o1451000000000000 / */ +/* data mcheps(2) / o0000000000000000 / */ + +/* data minmag(1) / o1771000000000000 / */ +/* data minmag(2) / o0000000000000000 / */ + +/* data maxmag(1) / o0777777777777777 / */ +/* data maxmag(2) / o0007777777777777 / */ + +/* Machine constants for the Burroughs 1700 system. */ + +/* data mcheps(1) / zcc6800000 / */ +/* data mcheps(2) / z000000000 / */ + +/* data minmag(1) / zc00800000 / */ +/* data minmag(2) / z000000000 / */ + +/* data maxmag(1) / zdffffffff / */ +/* data maxmag(2) / zfffffffff / */ + +/* Machine constants for the Univac 1100 series. */ + +/* data mcheps(1),mcheps(2) / o170640000000, o000000000000 / */ +/* data minmag(1),minmag(2) / o000040000000, o000000000000 / */ +/* data maxmag(1),maxmag(2) / o377777777777, o777777777777 / */ + +/* Machine constants for the Data General Eclipse S/200. */ + +/* Note - it may be appropriate to include the following card - */ +/* static dmach(3) */ + +/* data minmag/20k,3*0/,maxmag/77777k,3*177777k/ */ +/* data mcheps/32020k,3*0/ */ + +/* Machine constants for the Harris 220. */ + +/* data mcheps(1),mcheps(2) / '20000000, '00000334 / */ +/* data minmag(1),minmag(2) / '20000000, '00000201 / */ +/* data maxmag(1),maxmag(2) / '37777777, '37777577 / */ + +/* Machine constants for the Cray-1. */ + +/* data mcheps(1) / 0376424000000000000000b / */ +/* data mcheps(2) / 0000000000000000000000b / */ + +/* data minmag(1) / 0200034000000000000000b / */ +/* data minmag(2) / 0000000000000000000000b / */ + +/* data maxmag(1) / 0577777777777777777777b / */ +/* data maxmag(2) / 0000007777777777777776b / */ + +/* Machine constants for the Prime 400. */ + +/* data mcheps(1),mcheps(2) / :10000000000, :00000000123 / */ +/* data minmag(1),minmag(2) / :10000000000, :00000100000 / */ +/* data maxmag(1),maxmag(2) / :17777777777, :37777677776 / */ + +/* Machine constants for the VAX-11. */ + +/* data mcheps(1),mcheps(2) / 9472, 0 / */ +/* data minmag(1),minmag(2) / 128, 0 / */ +/* data maxmag(1),maxmag(2) / -32769, -1 / */ + +/* Machine constants for IEEE machines. */ + +/* data dmach(1) /2.22044604926d-16/ */ +/* data dmach(2) /2.22507385852d-308/ */ +/* data dmach(3) /1.79769313485d+308/ */ + + switch(i) { + case 1: + return DPMPAR(real,EPSILON); /* 2.2204460492503131e-16 | 1.19209290e-07F */ + case 2: + return DPMPAR(real,MIN); /* 2.2250738585072014e-308 | 1.17549435e-38F */ + default: + return DPMPAR(real,MAX); /* 1.7976931348623157e+308 | 3.40282347e+38F */ + } + +/* Last card of function dpmpar. */ + +} /* dpmpar_ */ + +#undef mcheps +#undef maxmag +#undef minmag +#undef dmach + + diff --git a/ast/cminpack/enorm.c b/ast/cminpack/enorm.c new file mode 100644 index 0000000..ad10824 --- /dev/null +++ b/ast/cminpack/enorm.c @@ -0,0 +1,157 @@ +#include "cminpack.h" +#include +#include "cminpackP.h" + +/* + About the values for rdwarf and rgiant. + + The original values, both in signe-precision FORTRAN source code and in double-precision code were: +#define rdwarf 3.834e-20 +#define rgiant 1.304e19 + See for example: + http://www.netlib.org/slatec/src/denorm.f + http://www.netlib.org/slatec/src/enorm.f + However, rdwarf is smaller than sqrt(FLT_MIN) = 1.0842021724855044e-19, so that rdwarf**2 will + underflow. This contradicts the constraints expressed in the comments below. + + We changed these constants to be sqrt(dpmpar(2))*0.9 and sqrt(dpmpar(3))*0.9, as proposed by the + implementation found in MPFIT http://cow.physics.wisc.edu/~craigm/idl/fitting.html +*/ + +#define double_dwarf (1.4916681462400413e-154*0.9) +#define double_giant (1.3407807929942596e+154*0.9) +#define float_dwarf (1.0842021724855044e-19f*0.9f) +#define float_giant (1.8446743523953730e+19f*0.9f) +#define half_dwarf (2.4414062505039999e-4f*0.9f) +#define half_giant (255.93749236874225497222f*0.9f) + +#define dwarf(type) _dwarf(type) +#define _dwarf(type) type ## _dwarf +#define giant(type) _giant(type) +#define _giant(type) type ## _giant + +#define rdwarf dwarf(real) +#define rgiant giant(real) + +__cminpack_attr__ +real __cminpack_func__(enorm)(int n, const real *x) +{ +#ifdef USE_CBLAS + return cblas_dnrm2(n, x, 1); +#else /* !USE_CBLAS */ + /* System generated locals */ + real ret_val, d1; + + /* Local variables */ + int i; + real s1, s2, s3, xabs, x1max, x3max, agiant, floatn; + +/* ********** */ + +/* function enorm */ + +/* given an n-vector x, this function calculates the */ +/* euclidean norm of x. */ + +/* the euclidean norm is computed by accumulating the sum of */ +/* squares in three different sums. the sums of squares for the */ +/* small and large components are scaled so that no overflows */ +/* occur. non-destructive underflows are permitted. underflows */ +/* and overflows do not occur in the computation of the unscaled */ +/* sum of squares for the intermediate components. */ +/* the definitions of small, intermediate and large components */ +/* depend on two constants, rdwarf and rgiant. the main */ +/* restrictions on these constants are that rdwarf**2 not */ +/* underflow and rgiant**2 not overflow. the constants */ +/* given here are suitable for every known computer. */ + +/* the function statement is */ + +/* double precision function enorm(n,x) */ + +/* where */ + +/* n is a positive integer input variable. */ + +/* x is an input array of length n. */ + +/* subprograms called */ + +/* fortran-supplied ... dabs,dsqrt */ + +/* argonne national laboratory. minpack project. march 1980. */ +/* burton s. garbow, kenneth e. hillstrom, jorge j. more */ + +/* ********** */ + + s1 = 0.; + s2 = 0.; + s3 = 0.; + x1max = 0.; + x3max = 0.; + floatn = (real) (n); + agiant = rgiant / floatn; + for (i = 0; i < n; ++i) { + xabs = fabs(x[i]); + if (xabs <= rdwarf || xabs >= agiant) { + if (xabs > rdwarf) { + +/* sum for large components. */ + + if (xabs > x1max) { + /* Computing 2nd power */ + d1 = x1max / xabs; + s1 = 1. + s1 * (d1 * d1); + x1max = xabs; + } else { + /* Computing 2nd power */ + d1 = xabs / x1max; + s1 += d1 * d1; + } + } else { + +/* sum for small components. */ + + if (xabs > x3max) { + /* Computing 2nd power */ + d1 = x3max / xabs; + s3 = 1. + s3 * (d1 * d1); + x3max = xabs; + } else { + if (xabs != 0.) { + /* Computing 2nd power */ + d1 = xabs / x3max; + s3 += d1 * d1; + } + } + } + } else { + +/* sum for intermediate components. */ + + /* Computing 2nd power */ + s2 += xabs * xabs; + } + } + +/* calculation of norm. */ + + if (s1 != 0.) { + ret_val = x1max * sqrt(s1 + (s2 / x1max) / x1max); + } else { + if (s2 != 0.) { + if (s2 >= x3max) { + ret_val = sqrt(s2 * (1. + (x3max / s2) * (x3max * s3))); + } else { + ret_val = sqrt(x3max * ((s2 / x3max) + (x3max * s3))); + } + } else { + ret_val = x3max * sqrt(s3); + } + } + return ret_val; + +/* last card of function enorm. */ +#endif /* !USE_CBLAS */ +} /* enorm_ */ + diff --git a/ast/cminpack/lmder.c b/ast/cminpack/lmder.c new file mode 100644 index 0000000..7f57428 --- /dev/null +++ b/ast/cminpack/lmder.c @@ -0,0 +1,526 @@ +#include "cminpack.h" +#include +#include "cminpackP.h" + +__cminpack_attr__ +int __cminpack_func__(lmder)(__cminpack_decl_fcnder_mn__ void *p, int m, int n, real *x, + real *fvec, real *fjac, int ldfjac, real ftol, + real xtol, real gtol, int maxfev, real * + diag, int mode, real factor, int nprint, + int *nfev, int *njev, int *ipvt, real *qtf, + real *wa1, real *wa2, real *wa3, real *wa4) +{ + /* Initialized data */ + +#define p1 .1 +#define p5 .5 +#define p25 .25 +#define p75 .75 +#define p0001 1e-4 + + /* System generated locals */ + real d1, d2; + + /* Local variables */ + int i, j, l; + real par, sum; + int iter; + real temp, temp1, temp2; + int iflag; + real delta = 0.; + real ratio; + real fnorm, gnorm, pnorm, xnorm = 0., fnorm1, actred, dirder, + epsmch, prered; + int info; + +/* ********** */ + +/* subroutine lmder */ + +/* the purpose of lmder is to minimize the sum of the squares of */ +/* m nonlinear functions in n variables by a modification of */ +/* the levenberg-marquardt algorithm. the user must provide a */ +/* subroutine which calculates the functions and the jacobian. */ + +/* the subroutine statement is */ + +/* subroutine lmder(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol, */ +/* maxfev,diag,mode,factor,nprint,info,nfev, */ +/* njev,ipvt,qtf,wa1,wa2,wa3,wa4) */ + +/* where */ + +/* fcn is the name of the user-supplied subroutine which */ +/* calculates the functions and the jacobian. fcn must */ +/* be declared in an external statement in the user */ +/* calling program, and should be written as follows. */ + +/* subroutine fcn(m,n,x,fvec,fjac,ldfjac,iflag) */ +/* integer m,n,ldfjac,iflag */ +/* double precision x(n),fvec(m),fjac(ldfjac,n) */ +/* ---------- */ +/* if iflag = 1 calculate the functions at x and */ +/* return this vector in fvec. do not alter fjac. */ +/* if iflag = 2 calculate the jacobian at x and */ +/* return this matrix in fjac. do not alter fvec. */ +/* ---------- */ +/* return */ +/* end */ + +/* the value of iflag should not be changed by fcn unless */ +/* the user wants to terminate execution of lmder. */ +/* in this case set iflag to a negative integer. */ + +/* m is a positive integer input variable set to the number */ +/* of functions. */ + +/* n is a positive integer input variable set to the number */ +/* of variables. n must not exceed m. */ + +/* x is an array of length n. on input x must contain */ +/* an initial estimate of the solution vector. on output x */ +/* contains the final estimate of the solution vector. */ + +/* fvec is an output array of length m which contains */ +/* the functions evaluated at the output x. */ + +/* fjac is an output m by n array. the upper n by n submatrix */ +/* of fjac contains an upper triangular matrix r with */ +/* diagonal elements of nonincreasing magnitude such that */ + +/* t t t */ +/* p *(jac *jac)*p = r *r, */ + +/* where p is a permutation matrix and jac is the final */ +/* calculated jacobian. column j of p is column ipvt(j) */ +/* (see below) of the identity matrix. the lower trapezoidal */ +/* part of fjac contains information generated during */ +/* the computation of r. */ + +/* ldfjac is a positive integer input variable not less than m */ +/* which specifies the leading dimension of the array fjac. */ + +/* ftol is a nonnegative input variable. termination */ +/* occurs when both the actual and predicted relative */ +/* reductions in the sum of squares are at most ftol. */ +/* therefore, ftol measures the relative error desired */ +/* in the sum of squares. */ + +/* xtol is a nonnegative input variable. termination */ +/* occurs when the relative error between two consecutive */ +/* iterates is at most xtol. therefore, xtol measures the */ +/* relative error desired in the approximate solution. */ + +/* gtol is a nonnegative input variable. termination */ +/* occurs when the cosine of the angle between fvec and */ +/* any column of the jacobian is at most gtol in absolute */ +/* value. therefore, gtol measures the orthogonality */ +/* desired between the function vector and the columns */ +/* of the jacobian. */ + +/* maxfev is a positive integer input variable. termination */ +/* occurs when the number of calls to fcn with iflag = 1 */ +/* has reached maxfev. */ + +/* diag is an array of length n. if mode = 1 (see */ +/* below), diag is internally set. if mode = 2, diag */ +/* must contain positive entries that serve as */ +/* multiplicative scale factors for the variables. */ + +/* mode is an integer input variable. if mode = 1, the */ +/* variables will be scaled internally. if mode = 2, */ +/* the scaling is specified by the input diag. other */ +/* values of mode are equivalent to mode = 1. */ + +/* factor is a positive input variable used in determining the */ +/* initial step bound. this bound is set to the product of */ +/* factor and the euclidean norm of diag*x if nonzero, or else */ +/* to factor itself. in most cases factor should lie in the */ +/* interval (.1,100.).100. is a generally recommended value. */ + +/* nprint is an integer input variable that enables controlled */ +/* printing of iterates if it is positive. in this case, */ +/* fcn is called with iflag = 0 at the beginning of the first */ +/* iteration and every nprint iterations thereafter and */ +/* immediately prior to return, with x, fvec, and fjac */ +/* available for printing. fvec and fjac should not be */ +/* altered. if nprint is not positive, no special calls */ +/* of fcn with iflag = 0 are made. */ + +/* info is an integer output variable. if the user has */ +/* terminated execution, info is set to the (negative) */ +/* value of iflag. see description of fcn. otherwise, */ +/* info is set as follows. */ + +/* info = 0 improper input parameters. */ + +/* info = 1 both actual and predicted relative reductions */ +/* in the sum of squares are at most ftol. */ + +/* info = 2 relative error between two consecutive iterates */ +/* is at most xtol. */ + +/* info = 3 conditions for info = 1 and info = 2 both hold. */ + +/* info = 4 the cosine of the angle between fvec and any */ +/* column of the jacobian is at most gtol in */ +/* absolute value. */ + +/* info = 5 number of calls to fcn with iflag = 1 has */ +/* reached maxfev. */ + +/* info = 6 ftol is too small. no further reduction in */ +/* the sum of squares is possible. */ + +/* info = 7 xtol is too small. no further improvement in */ +/* the approximate solution x is possible. */ + +/* info = 8 gtol is too small. fvec is orthogonal to the */ +/* columns of the jacobian to machine precision. */ + +/* nfev is an integer output variable set to the number of */ +/* calls to fcn with iflag = 1. */ + +/* njev is an integer output variable set to the number of */ +/* calls to fcn with iflag = 2. */ + +/* ipvt is an integer output array of length n. ipvt */ +/* defines a permutation matrix p such that jac*p = q*r, */ +/* where jac is the final calculated jacobian, q is */ +/* orthogonal (not stored), and r is upper triangular */ +/* with diagonal elements of nonincreasing magnitude. */ +/* column j of p is column ipvt(j) of the identity matrix. */ + +/* qtf is an output array of length n which contains */ +/* the first n elements of the vector (q transpose)*fvec. */ + +/* wa1, wa2, and wa3 are work arrays of length n. */ + +/* wa4 is a work array of length m. */ + +/* subprograms called */ + +/* user-supplied ...... fcn */ + +/* minpack-supplied ... dpmpar,enorm,lmpar,qrfac */ + +/* fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod */ + +/* argonne national laboratory. minpack project. march 1980. */ +/* burton s. garbow, kenneth e. hillstrom, jorge j. more */ + +/* ********** */ + +/* epsmch is the machine precision. */ + + epsmch = __cminpack_func__(dpmpar)(1); + + info = 0; + iflag = 0; + *nfev = 0; + *njev = 0; + +/* check the input parameters for errors. */ + + if (n <= 0 || m < n || ldfjac < m || ftol < 0. || xtol < 0. || + gtol < 0. || maxfev <= 0 || factor <= 0.) { + goto TERMINATE; + } + if (mode == 2) { + for (j = 0; j < n; ++j) { + if (diag[j] <= 0.) { + goto TERMINATE; + } + } + } + +/* evaluate the function at the starting point */ +/* and calculate its norm. */ + + iflag = fcnder_mn(p, m, n, x, fvec, fjac, ldfjac, 1); + *nfev = 1; + if (iflag < 0) { + goto TERMINATE; + } + fnorm = __cminpack_enorm__(m, fvec); + +/* initialize levenberg-marquardt parameter and iteration counter. */ + + par = 0.; + iter = 1; + +/* beginning of the outer loop. */ + + for (;;) { + +/* calculate the jacobian matrix. */ + + iflag = fcnder_mn(p, m, n, x, fvec, fjac, ldfjac, 2); + ++(*njev); + if (iflag < 0) { + goto TERMINATE; + } + +/* if requested, call fcn to enable printing of iterates. */ + + if (nprint > 0) { + iflag = 0; + if ((iter - 1) % nprint == 0) { + iflag = fcnder_mn(p, m, n, x, fvec, fjac, ldfjac, 0); + } + if (iflag < 0) { + goto TERMINATE; + } + } + +/* compute the qr factorization of the jacobian. */ + + __cminpack_func__(qrfac)(m, n, fjac, ldfjac, TRUE_, ipvt, n, + wa1, wa2, wa3); + +/* on the first iteration and if mode is 1, scale according */ +/* to the norms of the columns of the initial jacobian. */ + + if (iter == 1) { + if (mode != 2) { + for (j = 0; j < n; ++j) { + diag[j] = wa2[j]; + if (wa2[j] == 0.) { + diag[j] = 1.; + } + } + } + +/* on the first iteration, calculate the norm of the scaled x */ +/* and initialize the step bound delta. */ + + for (j = 0; j < n; ++j) { + wa3[j] = diag[j] * x[j]; + } + xnorm = __cminpack_enorm__(n, wa3); + delta = factor * xnorm; + if (delta == 0.) { + delta = factor; + } + } + +/* form (q transpose)*fvec and store the first n components in */ +/* qtf. */ + + for (i = 0; i < m; ++i) { + wa4[i] = fvec[i]; + } + for (j = 0; j < n; ++j) { + if (fjac[j + j * ldfjac] != 0.) { + sum = 0.; + for (i = j; i < m; ++i) { + sum += fjac[i + j * ldfjac] * wa4[i]; + } + temp = -sum / fjac[j + j * ldfjac]; + for (i = j; i < m; ++i) { + wa4[i] += fjac[i + j * ldfjac] * temp; + } + } + fjac[j + j * ldfjac] = wa1[j]; + qtf[j] = wa4[j]; + } + +/* compute the norm of the scaled gradient. */ + + gnorm = 0.; + if (fnorm != 0.) { + for (j = 0; j < n; ++j) { + l = ipvt[j]-1; + if (wa2[l] != 0.) { + sum = 0.; + for (i = 0; i <= j; ++i) { + sum += fjac[i + j * ldfjac] * (qtf[i] / fnorm); + } + /* Computing MAX */ + d1 = fabs(sum / wa2[l]); + gnorm = max(gnorm,d1); + } + } + } + +/* test for convergence of the gradient norm. */ + + if (gnorm <= gtol) { + info = 4; + } + if (info != 0) { + goto TERMINATE; + } + +/* rescale if necessary. */ + + if (mode != 2) { + for (j = 0; j < n; ++j) { + /* Computing MAX */ + d1 = diag[j], d2 = wa2[j]; + diag[j] = max(d1,d2); + } + } + +/* beginning of the inner loop. */ + + do { + +/* determine the levenberg-marquardt parameter. */ + + __cminpack_func__(lmpar)(n, fjac, ldfjac, ipvt, diag, qtf, delta, + &par, wa1, wa2, wa3, wa4); + +/* store the direction p and x + p. calculate the norm of p. */ + + for (j = 0; j < n; ++j) { + wa1[j] = -wa1[j]; + wa2[j] = x[j] + wa1[j]; + wa3[j] = diag[j] * wa1[j]; + } + pnorm = __cminpack_enorm__(n, wa3); + +/* on the first iteration, adjust the initial step bound. */ + + if (iter == 1) { + delta = min(delta,pnorm); + } + +/* evaluate the function at x + p and calculate its norm. */ + + iflag = fcnder_mn(p, m, n, wa2, wa4, fjac, ldfjac, 1); + ++(*nfev); + if (iflag < 0) { + goto TERMINATE; + } + fnorm1 = __cminpack_enorm__(m, wa4); + +/* compute the scaled actual reduction. */ + + actred = -1.; + if (p1 * fnorm1 < fnorm) { + /* Computing 2nd power */ + d1 = fnorm1 / fnorm; + actred = 1. - d1 * d1; + } + +/* compute the scaled predicted reduction and */ +/* the scaled directional derivative. */ + + for (j = 0; j < n; ++j) { + wa3[j] = 0.; + l = ipvt[j]-1; + temp = wa1[l]; + for (i = 0; i <= j; ++i) { + wa3[i] += fjac[i + j * ldfjac] * temp; + } + } + temp1 = __cminpack_enorm__(n, wa3) / fnorm; + temp2 = (sqrt(par) * pnorm) / fnorm; + prered = temp1 * temp1 + temp2 * temp2 / p5; + dirder = -(temp1 * temp1 + temp2 * temp2); + +/* compute the ratio of the actual to the predicted */ +/* reduction. */ + + ratio = 0.; + if (prered != 0.) { + ratio = actred / prered; + } + +/* update the step bound. */ + + if (ratio <= p25) { + if (actred >= 0.) { + temp = p5; + } else { + temp = p5 * dirder / (dirder + p5 * actred); + } + if (p1 * fnorm1 >= fnorm || temp < p1) { + temp = p1; + } + /* Computing MIN */ + d1 = pnorm / p1; + delta = temp * min(delta,d1); + par /= temp; + } else { + if (par == 0. || ratio >= p75) { + delta = pnorm / p5; + par = p5 * par; + } + } + +/* test for successful iteration. */ + + if (ratio >= p0001) { + +/* successful iteration. update x, fvec, and their norms. */ + + for (j = 0; j < n; ++j) { + x[j] = wa2[j]; + wa2[j] = diag[j] * x[j]; + } + for (i = 0; i < m; ++i) { + fvec[i] = wa4[i]; + } + xnorm = __cminpack_enorm__(n, wa2); + fnorm = fnorm1; + ++iter; + } + +/* tests for convergence. */ + + if (fabs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) { + info = 1; + } + if (delta <= xtol * xnorm) { + info = 2; + } + if (fabs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info == 2) { + info = 3; + } + if (info != 0) { + goto TERMINATE; + } + +/* tests for termination and stringent tolerances. */ + + if (*nfev >= maxfev) { + info = 5; + } + if (fabs(actred) <= epsmch && prered <= epsmch && p5 * ratio <= 1.) { + info = 6; + } + if (delta <= epsmch * xnorm) { + info = 7; + } + if (gnorm <= epsmch) { + info = 8; + } + if (info != 0) { + goto TERMINATE; + } + +/* end of the inner loop. repeat if iteration unsuccessful. */ + + } while (ratio < p0001); + +/* end of the outer loop. */ + + } +TERMINATE: + +/* termination, either normal or user imposed. */ + + if (iflag < 0) { + info = iflag; + } + if (nprint > 0) { + fcnder_mn(p, m, n, x, fvec, fjac, ldfjac, 0); + } + return info; + +/* last card of subroutine lmder. */ + +} /* lmder_ */ + diff --git a/ast/cminpack/lmder1.c b/ast/cminpack/lmder1.c new file mode 100644 index 0000000..581462e --- /dev/null +++ b/ast/cminpack/lmder1.c @@ -0,0 +1,167 @@ +#include "cminpack.h" +#include "cminpackP.h" + +__cminpack_attr__ +int __cminpack_func__(lmder1)(__cminpack_decl_fcnder_mn__ void *p, int m, int n, real *x, + real *fvec, real *fjac, int ldfjac, real tol, + int *ipvt, real *wa, int lwa) +{ + /* Initialized data */ + + const real factor = 100.; + + /* Local variables */ + int mode, nfev, njev; + real ftol, gtol, xtol; + int maxfev, nprint; + int info; + +/* ********** */ + +/* subroutine lmder1 */ + +/* the purpose of lmder1 is to minimize the sum of the squares of */ +/* m nonlinear functions in n variables by a modification of the */ +/* levenberg-marquardt algorithm. this is done by using the more */ +/* general least-squares solver lmder. the user must provide a */ +/* subroutine which calculates the functions and the jacobian. */ + +/* the subroutine statement is */ + +/* subroutine lmder1(fcn,m,n,x,fvec,fjac,ldfjac,tol,info, */ +/* ipvt,wa,lwa) */ + +/* where */ + +/* fcn is the name of the user-supplied subroutine which */ +/* calculates the functions and the jacobian. fcn must */ +/* be declared in an external statement in the user */ +/* calling program, and should be written as follows. */ + +/* subroutine fcn(m,n,x,fvec,fjac,ldfjac,iflag) */ +/* integer m,n,ldfjac,iflag */ +/* double precision x(n),fvec(m),fjac(ldfjac,n) */ +/* ---------- */ +/* if iflag = 1 calculate the functions at x and */ +/* return this vector in fvec. do not alter fjac. */ +/* if iflag = 2 calculate the jacobian at x and */ +/* return this matrix in fjac. do not alter fvec. */ +/* ---------- */ +/* return */ +/* end */ + +/* the value of iflag should not be changed by fcn unless */ +/* the user wants to terminate execution of lmder1. */ +/* in this case set iflag to a negative integer. */ + +/* m is a positive integer input variable set to the number */ +/* of functions. */ + +/* n is a positive integer input variable set to the number */ +/* of variables. n must not exceed m. */ + +/* x is an array of length n. on input x must contain */ +/* an initial estimate of the solution vector. on output x */ +/* contains the final estimate of the solution vector. */ + +/* fvec is an output array of length m which contains */ +/* the functions evaluated at the output x. */ + +/* fjac is an output m by n array. the upper n by n submatrix */ +/* of fjac contains an upper triangular matrix r with */ +/* diagonal elements of nonincreasing magnitude such that */ + +/* t t t */ +/* p *(jac *jac)*p = r *r, */ + +/* where p is a permutation matrix and jac is the final */ +/* calculated jacobian. column j of p is column ipvt(j) */ +/* (see below) of the identity matrix. the lower trapezoidal */ +/* part of fjac contains information generated during */ +/* the computation of r. */ + +/* ldfjac is a positive integer input variable not less than m */ +/* which specifies the leading dimension of the array fjac. */ + +/* tol is a nonnegative input variable. termination occurs */ +/* when the algorithm estimates either that the relative */ +/* error in the sum of squares is at most tol or that */ +/* the relative error between x and the solution is at */ +/* most tol. */ + +/* info is an integer output variable. if the user has */ +/* terminated execution, info is set to the (negative) */ +/* value of iflag. see description of fcn. otherwise, */ +/* info is set as follows. */ + +/* info = 0 improper input parameters. */ + +/* info = 1 algorithm estimates that the relative error */ +/* in the sum of squares is at most tol. */ + +/* info = 2 algorithm estimates that the relative error */ +/* between x and the solution is at most tol. */ + +/* info = 3 conditions for info = 1 and info = 2 both hold. */ + +/* info = 4 fvec is orthogonal to the columns of the */ +/* jacobian to machine precision. */ + +/* info = 5 number of calls to fcn with iflag = 1 has */ +/* reached 100*(n+1). */ + +/* info = 6 tol is too small. no further reduction in */ +/* the sum of squares is possible. */ + +/* info = 7 tol is too small. no further improvement in */ +/* the approximate solution x is possible. */ + +/* ipvt is an integer output array of length n. ipvt */ +/* defines a permutation matrix p such that jac*p = q*r, */ +/* where jac is the final calculated jacobian, q is */ +/* orthogonal (not stored), and r is upper triangular */ +/* with diagonal elements of nonincreasing magnitude. */ +/* column j of p is column ipvt(j) of the identity matrix. */ + +/* wa is a work array of length lwa. */ + +/* lwa is a positive integer input variable not less than 5*n+m. */ + +/* subprograms called */ + +/* user-supplied ...... fcn */ + +/* minpack-supplied ... lmder */ + +/* argonne national laboratory. minpack project. march 1980. */ +/* burton s. garbow, kenneth e. hillstrom, jorge j. more */ + +/* ********** */ + +/* check the input parameters for errors. */ + + if (n <= 0 || m < n || ldfjac < m || tol < 0. || lwa < n * 5 + m) { + return 0; + } + +/* call lmder. */ + + maxfev = (n + 1) * 100; + ftol = tol; + xtol = tol; + gtol = 0.; + mode = 1; + nprint = 0; + info = __cminpack_func__(lmder)(__cminpack_param_fcnder_mn__ p, m, n, x, fvec, fjac, ldfjac, + ftol, xtol, gtol, maxfev, wa, mode, factor, nprint, + &nfev, &njev, ipvt, &wa[n], &wa[(n << 1)], & + wa[n * 3], &wa[(n << 2)], &wa[n * 5]); + if (info == 8) { + info = 4; + } + return info; + +/* last card of subroutine lmder1. */ + +} /* lmder1_ */ + diff --git a/ast/cminpack/lmpar.c b/ast/cminpack/lmpar.c new file mode 100644 index 0000000..108e687 --- /dev/null +++ b/ast/cminpack/lmpar.c @@ -0,0 +1,338 @@ +/* lmpar.f -- translated by f2c (version 20020621). + You must link the resulting object file with the libraries: + -lf2c -lm (in that order) +*/ + +#include "cminpack.h" +#include +#include "cminpackP.h" + +__cminpack_attr__ +void __cminpack_func__(lmpar)(int n, real *r, int ldr, + const int *ipvt, const real *diag, const real *qtb, real delta, + real *par, real *x, real *sdiag, real *wa1, + real *wa2) +{ + /* Initialized data */ + +#define p1 .1 +#define p001 .001 + + /* System generated locals */ + real d1, d2; + + /* Local variables */ + int j, l; + real fp; + real parc, parl; + int iter; + real temp, paru, dwarf; + int nsing; + real gnorm; + real dxnorm; + +/* ********** */ + +/* subroutine lmpar */ + +/* given an m by n matrix a, an n by n nonsingular diagonal */ +/* matrix d, an m-vector b, and a positive number delta, */ +/* the problem is to determine a value for the parameter */ +/* par such that if x solves the system */ + +/* a*x = b , sqrt(par)*d*x = 0 , */ + +/* in the least squares sense, and dxnorm is the euclidean */ +/* norm of d*x, then either par is zero and */ + +/* (dxnorm-delta) .le. 0.1*delta , */ + +/* or par is positive and */ + +/* abs(dxnorm-delta) .le. 0.1*delta . */ + +/* this subroutine completes the solution of the problem */ +/* if it is provided with the necessary information from the */ +/* qr factorization, with column pivoting, of a. that is, if */ +/* a*p = q*r, where p is a permutation matrix, q has orthogonal */ +/* columns, and r is an upper triangular matrix with diagonal */ +/* elements of nonincreasing magnitude, then lmpar expects */ +/* the full upper triangle of r, the permutation matrix p, */ +/* and the first n components of (q transpose)*b. on output */ +/* lmpar also provides an upper triangular matrix s such that */ + +/* t t t */ +/* p *(a *a + par*d*d)*p = s *s . */ + +/* s is employed within lmpar and may be of separate interest. */ + +/* only a few iterations are generally needed for convergence */ +/* of the algorithm. if, however, the limit of 10 iterations */ +/* is reached, then the output par will contain the best */ +/* value obtained so far. */ + +/* the subroutine statement is */ + +/* subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, */ +/* wa1,wa2) */ + +/* where */ + +/* n is a positive integer input variable set to the order of r. */ + +/* r is an n by n array. on input the full upper triangle */ +/* must contain the full upper triangle of the matrix r. */ +/* on output the full upper triangle is unaltered, and the */ +/* strict lower triangle contains the strict upper triangle */ +/* (transposed) of the upper triangular matrix s. */ + +/* ldr is a positive integer input variable not less than n */ +/* which specifies the leading dimension of the array r. */ + +/* ipvt is an integer input array of length n which defines the */ +/* permutation matrix p such that a*p = q*r. column j of p */ +/* is column ipvt(j) of the identity matrix. */ + +/* diag is an input array of length n which must contain the */ +/* diagonal elements of the matrix d. */ + +/* qtb is an input array of length n which must contain the first */ +/* n elements of the vector (q transpose)*b. */ + +/* delta is a positive input variable which specifies an upper */ +/* bound on the euclidean norm of d*x. */ + +/* par is a nonnegative variable. on input par contains an */ +/* initial estimate of the levenberg-marquardt parameter. */ +/* on output par contains the final estimate. */ + +/* x is an output array of length n which contains the least */ +/* squares solution of the system a*x = b, sqrt(par)*d*x = 0, */ +/* for the output par. */ + +/* sdiag is an output array of length n which contains the */ +/* diagonal elements of the upper triangular matrix s. */ + +/* wa1 and wa2 are work arrays of length n. */ + +/* subprograms called */ + +/* minpack-supplied ... dpmpar,enorm,qrsolv */ + +/* fortran-supplied ... dabs,dmax1,dmin1,dsqrt */ + +/* argonne national laboratory. minpack project. march 1980. */ +/* burton s. garbow, kenneth e. hillstrom, jorge j. more */ + +/* ********** */ + +/* dwarf is the smallest positive magnitude. */ + + dwarf = __cminpack_func__(dpmpar)(2); + +/* compute and store in x the gauss-newton direction. if the */ +/* jacobian is rank-deficient, obtain a least squares solution. */ + + nsing = n; + for (j = 0; j < n; ++j) { + wa1[j] = qtb[j]; + if (r[j + j * ldr] == 0. && nsing == n) { + nsing = j; + } + if (nsing < n) { + wa1[j] = 0.; + } + } +# ifdef USE_CBLAS + cblas_dtrsv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, nsing, r, ldr, wa1, 1); +# else + if (nsing >= 1) { + int k; + for (k = 1; k <= nsing; ++k) { + j = nsing - k; + wa1[j] /= r[j + j * ldr]; + temp = wa1[j]; + if (j >= 1) { + int i; + for (i = 0; i < j; ++i) { + wa1[i] -= r[i + j * ldr] * temp; + } + } + } + } +# endif + for (j = 0; j < n; ++j) { + l = ipvt[j]-1; + x[l] = wa1[j]; + } + +/* initialize the iteration counter. */ +/* evaluate the function at the origin, and test */ +/* for acceptance of the gauss-newton direction. */ + + iter = 0; + for (j = 0; j < n; ++j) { + wa2[j] = diag[j] * x[j]; + } + dxnorm = __cminpack_enorm__(n, wa2); + fp = dxnorm - delta; + if (fp <= p1 * delta) { + goto TERMINATE; + } + +/* if the jacobian is not rank deficient, the newton */ +/* step provides a lower bound, parl, for the zero of */ +/* the function. otherwise set this bound to zero. */ + + parl = 0.; + if (nsing >= n) { + for (j = 0; j < n; ++j) { + l = ipvt[j]-1; + wa1[j] = diag[l] * (wa2[l] / dxnorm); + } +# ifdef USE_CBLAS + cblas_dtrsv(CblasColMajor, CblasUpper, CblasTrans, CblasNonUnit, n, r, ldr, wa1, 1); +# else + for (j = 0; j < n; ++j) { + real sum = 0.; + if (j >= 1) { + int i; + for (i = 0; i < j; ++i) { + sum += r[i + j * ldr] * wa1[i]; + } + } + wa1[j] = (wa1[j] - sum) / r[j + j * ldr]; + } +# endif + temp = __cminpack_enorm__(n, wa1); + parl = fp / delta / temp / temp; + } + +/* calculate an upper bound, paru, for the zero of the function. */ + + for (j = 0; j < n; ++j) { + real sum; +# ifdef USE_CBLAS + sum = cblas_ddot(j+1, &r[j*ldr], 1, qtb, 1); +# else + int i; + sum = 0.; + for (i = 0; i <= j; ++i) { + sum += r[i + j * ldr] * qtb[i]; + } +# endif + l = ipvt[j]-1; + wa1[j] = sum / diag[l]; + } + gnorm = __cminpack_enorm__(n, wa1); + paru = gnorm / delta; + if (paru == 0.) { + paru = dwarf / min(delta,(real)p1) /* / p001 ??? */; + } + +/* if the input par lies outside of the interval (parl,paru), */ +/* set par to the closer endpoint. */ + + *par = max(*par,parl); + *par = min(*par,paru); + if (*par == 0.) { + *par = gnorm / dxnorm; + } + +/* beginning of an iteration. */ + + for (;;) { + ++iter; + +/* evaluate the function at the current value of par. */ + + if (*par == 0.) { + /* Computing MAX */ + d1 = dwarf, d2 = p001 * paru; + *par = max(d1,d2); + } + temp = sqrt(*par); + for (j = 0; j < n; ++j) { + wa1[j] = temp * diag[j]; + } + __cminpack_func__(qrsolv)(n, r, ldr, ipvt, wa1, qtb, x, sdiag, wa2); + for (j = 0; j < n; ++j) { + wa2[j] = diag[j] * x[j]; + } + dxnorm = __cminpack_enorm__(n, wa2); + temp = fp; + fp = dxnorm - delta; + +/* if the function is small enough, accept the current value */ +/* of par. also test for the exceptional cases where parl */ +/* is zero or the number of iterations has reached 10. */ + + if (fabs(fp) <= p1 * delta || (parl == 0. && fp <= temp && temp < 0.) || iter == 10) { + goto TERMINATE; + } + +/* compute the newton correction. */ + +# ifdef USE_CBLAS + for (j = 0; j < nsing; ++j) { + l = ipvt[j]-1; + wa1[j] = diag[l] * (wa2[l] / dxnorm); + } + for (j = nsing; j < n; ++j) { + wa1[j] = 0.; + } + /* exchange the diagonal of r with sdiag */ + cblas_dswap(n, r, ldr+1, sdiag, 1); + /* solve lower(r).x = wa1, result id put in wa1 */ + cblas_dtrsv(CblasColMajor, CblasLower, CblasNoTrans, CblasNonUnit, nsing, r, ldr, wa1, 1); + /* exchange the diagonal of r with sdiag */ + cblas_dswap( n, r, ldr+1, sdiag, 1); +# else /* !USE_CBLAS */ + for (j = 0; j < n; ++j) { + l = ipvt[j]-1; + wa1[j] = diag[l] * (wa2[l] / dxnorm); + } + for (j = 0; j < n; ++j) { + wa1[j] /= sdiag[j]; + temp = wa1[j]; + if (n > j+1) { + int i; + for (i = j+1; i < n; ++i) { + wa1[i] -= r[i + j * ldr] * temp; + } + } + } +# endif /* !USE_CBLAS */ + temp = __cminpack_enorm__(n, wa1); + parc = fp / delta / temp / temp; + +/* depending on the sign of the function, update parl or paru. */ + + if (fp > 0.) { + parl = max(parl,*par); + } + if (fp < 0.) { + paru = min(paru,*par); + } + +/* compute an improved estimate for par. */ + + /* Computing MAX */ + d1 = parl, d2 = *par + parc; + *par = max(d1,d2); + +/* end of an iteration. */ + + } +TERMINATE: + +/* termination. */ + + if (iter == 0) { + *par = 0.; + } + +/* last card of subroutine lmpar. */ + +} /* lmpar_ */ + diff --git a/ast/cminpack/qrfac.c b/ast/cminpack/qrfac.c new file mode 100644 index 0000000..1573772 --- /dev/null +++ b/ast/cminpack/qrfac.c @@ -0,0 +1,285 @@ +#include "cminpack.h" +#include +#ifdef USE_LAPACK +#include +#include +#include +#endif +#include "cminpackP.h" + +__cminpack_attr__ +void __cminpack_func__(qrfac)(int m, int n, real *a, int + lda, int pivot, int *ipvt, int lipvt, real *rdiag, + real *acnorm, real *wa) +{ +#ifdef USE_LAPACK + __CLPK_integer m_ = m; + __CLPK_integer n_ = n; + __CLPK_integer lda_ = lda; + __CLPK_integer *jpvt; + + int i, j, k; + double t; + double* tau = wa; + const __CLPK_integer ltau = m > n ? n : m; + __CLPK_integer lwork = -1; + __CLPK_integer info = 0; + double* work; + + if (pivot) { + assert( lipvt >= n ); + if (sizeof(__CLPK_integer) != sizeof(ipvt[0])) { + jpvt = malloc(n*sizeof(__CLPK_integer)); + } else { + /* __CLPK_integer is actually an int, just do a cast */ + jpvt = (__CLPK_integer *)ipvt; + } + /* set all columns free */ + memset(jpvt, 0, sizeof(int)*n); + } + + /* query optimal size of work */ + lwork = -1; + if (pivot) { + dgeqp3_(&m_,&n_,a,&lda_,jpvt,tau,tau,&lwork,&info); + lwork = (int)tau[0]; + assert( lwork >= 3*n+1 ); + } else { + dgeqrf_(&m_,&n_,a,&lda_,tau,tau,&lwork,&info); + lwork = (int)tau[0]; + assert( lwork >= 1 && lwork >= n ); + } + + assert( info == 0 ); + + /* alloc work area */ + work = (double *)malloc(sizeof(double)*lwork); + assert(work != NULL); + + /* set acnorm first (from the doc of qrfac, acnorm may point to the same area as rdiag) */ + if (acnorm != rdiag) { + for (j = 0; j < n; ++j) { + acnorm[j] = __cminpack_enorm__(m, &a[j * lda]); + } + } + + /* QR decomposition */ + if (pivot) { + dgeqp3_(&m_,&n_,a,&lda_,jpvt,tau,work,&lwork,&info); + } else { + dgeqrf_(&m_,&n_,a,&lda_,tau,work,&lwork,&info); + } + assert(info == 0); + + /* set rdiag, before the diagonal is replaced */ + memset(rdiag, 0, sizeof(double)*n); + for(i=0 ; i rdiag[kmax]) { + kmax = k; + } + } + if (kmax != j) { + for (i = 0; i < m; ++i) { + temp = a[i + j * lda]; + a[i + j * lda] = a[i + kmax * lda]; + a[i + kmax * lda] = temp; + } + rdiag[kmax] = rdiag[j]; + wa[kmax] = wa[j]; + k = ipvt[j]; + ipvt[j] = ipvt[kmax]; + ipvt[kmax] = k; + } + } + +/* compute the householder transformation to reduce the */ +/* j-th column of a to a multiple of the j-th unit vector. */ + + ajnorm = __cminpack_enorm__(m - (j+1) + 1, &a[j + j * lda]); + if (ajnorm != 0.) { + if (a[j + j * lda] < 0.) { + ajnorm = -ajnorm; + } + for (i = j; i < m; ++i) { + a[i + j * lda] /= ajnorm; + } + a[j + j * lda] += 1.; + +/* apply the transformation to the remaining columns */ +/* and update the norms. */ + + jp1 = j + 1; + if (n > jp1) { + for (k = jp1; k < n; ++k) { + sum = 0.; + for (i = j; i < m; ++i) { + sum += a[i + j * lda] * a[i + k * lda]; + } + temp = sum / a[j + j * lda]; + for (i = j; i < m; ++i) { + a[i + k * lda] -= temp * a[i + j * lda]; + } + if (pivot && rdiag[k] != 0.) { + temp = a[j + k * lda] / rdiag[k]; + /* Computing MAX */ + d1 = 1. - temp * temp; + rdiag[k] *= sqrt((max((real)0.,d1))); + /* Computing 2nd power */ + d1 = rdiag[k] / wa[k]; + if (p05 * (d1 * d1) <= epsmch) { + rdiag[k] = __cminpack_enorm__(m - (j+1), &a[jp1 + k * lda]); + wa[k] = rdiag[k]; + } + } + } + } + } + rdiag[j] = -ajnorm; + } + +/* last card of subroutine qrfac. */ +#endif /* !USE_LAPACK */ +} /* qrfac_ */ + diff --git a/ast/cminpack/qrsolv.c b/ast/cminpack/qrsolv.c new file mode 100644 index 0000000..6ab9e98 --- /dev/null +++ b/ast/cminpack/qrsolv.c @@ -0,0 +1,218 @@ +#include "cminpack.h" +#include +#include "cminpackP.h" + +__cminpack_attr__ +void __cminpack_func__(qrsolv)(int n, real *r, int ldr, + const int *ipvt, const real *diag, const real *qtb, real *x, + real *sdiag, real *wa) +{ + /* Initialized data */ + +#define p5 .5 +#define p25 .25 + + /* Local variables */ + int i, j, k, l; + real cos, sin, sum, temp; + int nsing; + real qtbpj; + +/* ********** */ + +/* subroutine qrsolv */ + +/* given an m by n matrix a, an n by n diagonal matrix d, */ +/* and an m-vector b, the problem is to determine an x which */ +/* solves the system */ + +/* a*x = b , d*x = 0 , */ + +/* in the least squares sense. */ + +/* this subroutine completes the solution of the problem */ +/* if it is provided with the necessary information from the */ +/* qr factorization, with column pivoting, of a. that is, if */ +/* a*p = q*r, where p is a permutation matrix, q has orthogonal */ +/* columns, and r is an upper triangular matrix with diagonal */ +/* elements of nonincreasing magnitude, then qrsolv expects */ +/* the full upper triangle of r, the permutation matrix p, */ +/* and the first n components of (q transpose)*b. the system */ +/* a*x = b, d*x = 0, is then equivalent to */ + +/* t t */ +/* r*z = q *b , p *d*p*z = 0 , */ + +/* where x = p*z. if this system does not have full rank, */ +/* then a least squares solution is obtained. on output qrsolv */ +/* also provides an upper triangular matrix s such that */ + +/* t t t */ +/* p *(a *a + d*d)*p = s *s . */ + +/* s is computed within qrsolv and may be of separate interest. */ + +/* the subroutine statement is */ + +/* subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa) */ + +/* where */ + +/* n is a positive integer input variable set to the order of r. */ + +/* r is an n by n array. on input the full upper triangle */ +/* must contain the full upper triangle of the matrix r. */ +/* on output the full upper triangle is unaltered, and the */ +/* strict lower triangle contains the strict upper triangle */ +/* (transposed) of the upper triangular matrix s. */ + +/* ldr is a positive integer input variable not less than n */ +/* which specifies the leading dimension of the array r. */ + +/* ipvt is an integer input array of length n which defines the */ +/* permutation matrix p such that a*p = q*r. column j of p */ +/* is column ipvt(j) of the identity matrix. */ + +/* diag is an input array of length n which must contain the */ +/* diagonal elements of the matrix d. */ + +/* qtb is an input array of length n which must contain the first */ +/* n elements of the vector (q transpose)*b. */ + +/* x is an output array of length n which contains the least */ +/* squares solution of the system a*x = b, d*x = 0. */ + +/* sdiag is an output array of length n which contains the */ +/* diagonal elements of the upper triangular matrix s. */ + +/* wa is a work array of length n. */ + +/* subprograms called */ + +/* fortran-supplied ... dabs,dsqrt */ + +/* argonne national laboratory. minpack project. march 1980. */ +/* burton s. garbow, kenneth e. hillstrom, jorge j. more */ + +/* ********** */ + +/* copy r and (q transpose)*b to preserve input and initialize s. */ +/* in particular, save the diagonal elements of r in x. */ + + for (j = 0; j < n; ++j) { + for (i = j; i < n; ++i) { + r[i + j * ldr] = r[j + i * ldr]; + } + x[j] = r[j + j * ldr]; + wa[j] = qtb[j]; + } + +/* eliminate the diagonal matrix d using a givens rotation. */ + + for (j = 0; j < n; ++j) { + +/* prepare the row of d to be eliminated, locating the */ +/* diagonal element using p from the qr factorization. */ + + l = ipvt[j]-1; + if (diag[l] != 0.) { + for (k = j; k < n; ++k) { + sdiag[k] = 0.; + } + sdiag[j] = diag[l]; + +/* the transformations to eliminate the row of d */ +/* modify only a single element of (q transpose)*b */ +/* beyond the first n, which is initially zero. */ + + qtbpj = 0.; + for (k = j; k < n; ++k) { + +/* determine a givens rotation which eliminates the */ +/* appropriate element in the current row of d. */ + + if (sdiag[k] != 0.) { +# ifdef USE_LAPACK + dlartg_( &r[k + k * ldr], &sdiag[k], &cos, &sin, &temp ); +# else /* !USE_LAPACK */ + if (fabs(r[k + k * ldr]) < fabs(sdiag[k])) { + real cotan; + cotan = r[k + k * ldr] / sdiag[k]; + sin = p5 / sqrt(p25 + p25 * (cotan * cotan)); + cos = sin * cotan; + } else { + real tan; + tan = sdiag[k] / r[k + k * ldr]; + cos = p5 / sqrt(p25 + p25 * (tan * tan)); + sin = cos * tan; + } + +/* compute the modified diagonal element of r and */ +/* the modified element of ((q transpose)*b,0). */ + +# endif /* !USE_LAPACK */ + temp = cos * wa[k] + sin * qtbpj; + qtbpj = -sin * wa[k] + cos * qtbpj; + wa[k] = temp; + +/* accumulate the tranformation in the row of s. */ +# ifdef USE_CBLAS + cblas_drot( n-k, &r[k + k * ldr], 1, &sdiag[k], 1, cos, sin ); +# else /* !USE_CBLAS */ + r[k + k * ldr] = cos * r[k + k * ldr] + sin * sdiag[k]; + if (n > k+1) { + for (i = k+1; i < n; ++i) { + temp = cos * r[i + k * ldr] + sin * sdiag[i]; + sdiag[i] = -sin * r[i + k * ldr] + cos * sdiag[i]; + r[i + k * ldr] = temp; + } + } +# endif /* !USE_CBLAS */ + } + } + } + +/* store the diagonal element of s and restore */ +/* the corresponding diagonal element of r. */ + + sdiag[j] = r[j + j * ldr]; + r[j + j * ldr] = x[j]; + } + +/* solve the triangular system for z. if the system is */ +/* singular, then obtain a least squares solution. */ + + nsing = n; + for (j = 0; j < n; ++j) { + if (sdiag[j] == 0. && nsing == n) { + nsing = j; + } + if (nsing < n) { + wa[j] = 0.; + } + } + if (nsing >= 1) { + for (k = 1; k <= nsing; ++k) { + j = nsing - k; + sum = 0.; + if (nsing > j+1) { + for (i = j+1; i < nsing; ++i) { + sum += r[i + j * ldr] * wa[i]; + } + } + wa[j] = (wa[j] - sum) / sdiag[j]; + } + } + +/* permute the components of z back to components of x. */ + + for (j = 0; j < n; ++j) { + l = ipvt[j]-1; + x[l] = wa[j]; + } + return; + +/* last card of subroutine qrsolv. */ + +} /* qrsolv_ */ + -- cgit v0.12