upgrade ast 8.7.1

11 files changed, 2504 insertions, 0 deletions

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 <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
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 <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
new file mode 100644
index 0000000..81c6fcd
--- /dev/null
+++ b/ast/cminpack/dpmpar.c
@@ -0,0 +1,201 @@
+#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
new file mode 100644
index 0000000..ad10824
--- /dev/null
+++ b/ast/cminpack/enorm.c
@@ -0,0 +1,157 @@
+#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
new file mode 100644
index 0000000..7f57428
--- /dev/null
+++ b/ast/cminpack/lmder.c
@@ -0,0 +1,526 @@
+#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
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 <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
new file mode 100644
index 0000000..1573772
--- /dev/null
+++ b/ast/cminpack/qrfac.c
@@ -0,0 +1,285 @@
+#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
new file mode 100644
index 0000000..6ab9e98
--- /dev/null
+++ b/ast/cminpack/qrsolv.c
@@ -0,0 +1,218 @@
+#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_ */
+