diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2018-01-09 19:26:44 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2018-01-09 19:26:44 (GMT) |
commit | 1332d38f2805d986ea130e43218c0d2e870b4dc1 (patch) | |
tree | aa72853cb8d0d8fcd53a6f5eddf196a374226706 /ast/cminpack | |
parent | 5e545ec8058cc5238dc870468b34b5d4617f307f (diff) | |
download | blt-1332d38f2805d986ea130e43218c0d2e870b4dc1.zip blt-1332d38f2805d986ea130e43218c0d2e870b4dc1.tar.gz blt-1332d38f2805d986ea130e43218c0d2e870b4dc1.tar.bz2 |
update ast 8.6.2
Diffstat (limited to 'ast/cminpack')
-rw-r--r-- | ast/cminpack/CopyrightMINPACK.txt | 52 | ||||
-rw-r--r-- | ast/cminpack/README.md | 128 | ||||
-rw-r--r-- | ast/cminpack/cminpack.h | 370 | ||||
-rw-r--r-- | ast/cminpack/cminpackP.h | 62 | ||||
-rw-r--r-- | ast/cminpack/dpmpar.c | 201 | ||||
-rw-r--r-- | ast/cminpack/enorm.c | 157 | ||||
-rw-r--r-- | ast/cminpack/lmder.c | 526 | ||||
-rw-r--r-- | ast/cminpack/lmder1.c | 167 | ||||
-rw-r--r-- | ast/cminpack/lmpar.c | 338 | ||||
-rw-r--r-- | ast/cminpack/qrfac.c | 285 | ||||
-rw-r--r-- | ast/cminpack/qrsolv.c | 218 |
11 files changed, 0 insertions, 2504 deletions
diff --git a/ast/cminpack/CopyrightMINPACK.txt b/ast/cminpack/CopyrightMINPACK.txt deleted file mode 100644 index ae7984d..0000000 --- a/ast/cminpack/CopyrightMINPACK.txt +++ /dev/null @@ -1,52 +0,0 @@ -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 deleted file mode 100644 index 66bbc6f..0000000 --- a/ast/cminpack/README.md +++ /dev/null @@ -1,128 +0,0 @@ -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 deleted file mode 100644 index 6d3f757..0000000 --- a/ast/cminpack/cminpack.h +++ /dev/null @@ -1,370 +0,0 @@ -/* 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 <OpenEXR/half.h> -#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 deleted file mode 100644 index 4e8ba7b..0000000 --- a/ast/cminpack/cminpackP.h +++ /dev/null @@ -1,62 +0,0 @@ -/* 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 <Accelerate/Accelerate.h> -#else -#include <cblas.h> -#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 <Accelerate/Accelerate.h> -#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 deleted file mode 100644 index 81c6fcd..0000000 --- a/ast/cminpack/dpmpar.c +++ /dev/null @@ -1,201 +0,0 @@ -#include "cminpack.h" -#include <float.h> -#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 deleted file mode 100644 index ad10824..0000000 --- a/ast/cminpack/enorm.c +++ /dev/null @@ -1,157 +0,0 @@ -#include "cminpack.h" -#include <math.h> -#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 deleted file mode 100644 index 7f57428..0000000 --- a/ast/cminpack/lmder.c +++ /dev/null @@ -1,526 +0,0 @@ -#include "cminpack.h" -#include <math.h> -#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 deleted file mode 100644 index 581462e..0000000 --- a/ast/cminpack/lmder1.c +++ /dev/null @@ -1,167 +0,0 @@ -#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 deleted file mode 100644 index 108e687..0000000 --- a/ast/cminpack/lmpar.c +++ /dev/null @@ -1,338 +0,0 @@ -/* 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 <math.h> -#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 deleted file mode 100644 index 1573772..0000000 --- a/ast/cminpack/qrfac.c +++ /dev/null @@ -1,285 +0,0 @@ -#include "cminpack.h" -#include <math.h> -#ifdef USE_LAPACK -#include <stdlib.h> -#include <string.h> -#include <assert.h> -#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<n ; ++i) { - rdiag[i] = a[i*lda+i]; - } - - /* modify lower trinagular part to look like qrfac's output */ - for(i=0 ; i<ltau ; ++i) { - k = i*lda+i; - t = tau[i]; - a[k] = t; - for(j=i+1 ; j<m ; j++) { - k++; - a[k] *= t; - } - } - - free(work); - if (pivot) { - /* convert back jpvt to ipvt */ - if (sizeof(__CLPK_integer) != sizeof(ipvt[0])) { - for(i=0; i<n; ++i) { - ipvt[i] = jpvt[i]; - } - free(jpvt); - } - } -#else /* !USE_LAPACK */ - /* Initialized data */ - -#define p05 .05 - - /* System generated locals */ - real d1; - - /* Local variables */ - int i, j, k, jp1; - real sum; - real temp; - int minmn; - real epsmch; - real ajnorm; - -/* ********** */ - -/* subroutine qrfac */ - -/* this subroutine uses householder transformations with column */ -/* pivoting (optional) to compute a qr factorization of the */ -/* m by n matrix a. that is, qrfac determines an orthogonal */ -/* matrix q, a permutation matrix p, and an upper trapezoidal */ -/* matrix r with diagonal elements of nonincreasing magnitude, */ -/* such that a*p = q*r. the householder transformation for */ -/* column k, k = 1,2,...,min(m,n), is of the form */ - -/* t */ -/* i - (1/u(k))*u*u */ - -/* where u has zeros in the first k-1 positions. the form of */ -/* this transformation and the method of pivoting first */ -/* appeared in the corresponding linpack subroutine. */ - -/* the subroutine statement is */ - -/* subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa) */ - -/* where */ - -/* m is a positive integer input variable set to the number */ -/* of rows of a. */ - -/* n is a positive integer input variable set to the number */ -/* of columns of a. */ - -/* a is an m by n array. on input a contains the matrix for */ -/* which the qr factorization is to be computed. on output */ -/* the strict upper trapezoidal part of a contains the strict */ -/* upper trapezoidal part of r, and the lower trapezoidal */ -/* part of a contains a factored form of q (the non-trivial */ -/* elements of the u vectors described above). */ - -/* lda is a positive integer input variable not less than m */ -/* which specifies the leading dimension of the array a. */ - -/* pivot is a logical input variable. if pivot is set true, */ -/* then column pivoting is enforced. if pivot is set false, */ -/* then no column pivoting is done. */ - -/* ipvt is an integer output array of length lipvt. ipvt */ -/* defines the permutation matrix p such that a*p = q*r. */ -/* column j of p is column ipvt(j) of the identity matrix. */ -/* if pivot is false, ipvt is not referenced. */ - -/* lipvt is a positive integer input variable. if pivot is false, */ -/* then lipvt may be as small as 1. if pivot is true, then */ -/* lipvt must be at least n. */ - -/* rdiag is an output array of length n which contains the */ -/* diagonal elements of r. */ - -/* acnorm is an output array of length n which contains the */ -/* norms of the corresponding columns of the input matrix a. */ -/* if this information is not needed, then acnorm can coincide */ -/* with rdiag. */ - -/* wa is a work array of length n. if pivot is false, then wa */ -/* can coincide with rdiag. */ - -/* subprograms called */ - -/* minpack-supplied ... dpmpar,enorm */ - -/* fortran-supplied ... dmax1,dsqrt,min0 */ - -/* argonne national laboratory. minpack project. march 1980. */ -/* burton s. garbow, kenneth e. hillstrom, jorge j. more */ - -/* ********** */ - (void)lipvt; - -/* epsmch is the machine precision. */ - - epsmch = __cminpack_func__(dpmpar)(1); - -/* compute the initial column norms and initialize several arrays. */ - - for (j = 0; j < n; ++j) { - acnorm[j] = __cminpack_enorm__(m, &a[j * lda + 0]); - rdiag[j] = acnorm[j]; - wa[j] = rdiag[j]; - if (pivot) { - ipvt[j] = j+1; - } - } - -/* reduce a to r with householder transformations. */ - - minmn = min(m,n); - for (j = 0; j < minmn; ++j) { - if (pivot) { - -/* bring the column of largest norm into the pivot position. */ - - int kmax = j; - for (k = j; k < n; ++k) { - if (rdiag[k] > 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 deleted file mode 100644 index 6ab9e98..0000000 --- a/ast/cminpack/qrsolv.c +++ /dev/null @@ -1,218 +0,0 @@ -#include "cminpack.h" -#include <math.h> -#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_ */ - |