From facb9273c571b6bef430e71b31ff229abd6d1a46 Mon Sep 17 00:00:00 2001

From: Mark Abraham <mark.j.abraham@gmail.com>

Date: Mon, 11 Jun 2018 08:30:19 +0200

Subject: [PATCH] Bump lmfit support to version 7

The breaking API change means that distributions will not be able to

build GROMACS reliably with support for lmfit of arbitrary versions.

Nor does lmfit provide any support for querying the version. Tools

that call such fitting will now issue a fatal error.

Updated the FindLmfit and other cmake code to make better use of

modern cmake idioms for managing build targets. lmfit support is now

handled by a tri-state GMX_USE_LMFIT cmake variable, which defaults to

using the bundled version.

Reorganized aspects of gmx_lmcurve to better encapsulate the

dependency, now that we have to support the absence of lmfit.

Updated install guide and COPYING file.

Refs #2533

Change-Id: I6b14e08be216f803326b0ad9215b8d89bb0962c1

---

diff --git a/COPYING b/COPYING

index 4cc52a2..ba2710e 100644

--- a/COPYING

+++ b/COPYING

@@ -1249,7 +1249,7 @@

 --

   Copyright (c) 1980-1999 University of Chicago,

                           as operator of Argonne National Laboratory

-  Copyright (c) 2004-2013 Joachim Wuttke, Forschungszentrum Juelich GmbH

+  Copyright (c) 2004-2015 Joachim Wuttke, Forschungszentrum Juelich GmbH

   All rights reserved.

diff --git a/cmake/FindLmfit.cmake b/cmake/FindLmfit.cmake

index 5bf489a..9dd23f1 100644

--- a/cmake/FindLmfit.cmake

+++ b/cmake/FindLmfit.cmake

@@ -1,7 +1,7 @@

 #

 # This file is part of the GROMACS molecular simulation package.

 #

-# Copyright (c) 2016, by the GROMACS development team, led by

+# Copyright (c) 2016,2018, by the GROMACS development team, led by

 # Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,

 # and including many others, as listed in the AUTHORS file in the

 # top-level source directory and at http://www.gromacs.org.

@@ -32,19 +32,20 @@

 # To help us fund GROMACS development, we humbly ask that you cite

 # the research papers on the package. Check out http://www.gromacs.org.

-# This package tries to find an external lmfit library. It is intended

-# to work with pkg-config, because that is the mechanism supported in

-# lmfit. Upon exit, the following variables may be set:

+# This package tries to find an external lmfit library, version

+# 7. Note that pkg-config support was removed before version 7, so is

+# no longer supported here. Upon exit, the following variables may be

+# set:

 #

 # LMFIT_FOUND       - lmfit was found

 # LMFIT_INCLUDE_DIR - lmfit include directory

-# LMFIT_LIBRARIES   - lmfit libraries

+# LMFIT_LIBRARY     - lmfit library

 # LMFIT_LINKS_OK    - lmfit libraries link correctly

 # LMFIT_VERSION     - lmfit version string as "major.minor"

 #

-# If you cannot use pkg-config for some reason, then setting

-# LMFIT_INCLUDE_DIRS and LMFIT_LIBRARY_DIRS on the cmake command line

-# to suitable values will work.

+# CMake will search the CMAKE_PREFIX_PATH in the usual way, but if you

+# need more control then setting LMFIT_INCLUDE_DIR and LMFIT_LIBRARY

+# on the cmake command line to suitable values will work.

 include(CMakePushCheckState)

 cmake_push_check_state()

@@ -55,12 +56,12 @@

         # lmfit doesn't support CMake-based find_package version

         # checking in 6.1, so this code does nothing.

         if(LMFIT_FIND_VERSION_EXACT)

-            pkg_check_modules(PC_LMFIT lmfit=${LMFIT_FIND_VERSION})

+            pkg_check_modules(PC_LMFIT QUIET lmfit=${LMFIT_FIND_VERSION})

         else()

-            pkg_check_modules(PC_LMFIT lmfit>=${LMFIT_FIND_VERSION})

+            pkg_check_modules(PC_LMFIT QUIET lmfit>=${LMFIT_FIND_VERSION})

         endif()

     else()

-        pkg_check_modules(PC_LMFIT lmfit)

+        pkg_check_modules(PC_LMFIT QUIET lmfit)

         if (PC_LMFIT_VERSION)

             string(REGEX REPLACE "^([0-9]+):([0-9]+)" "\\1.\\2" LMFIT_VERSION "${PC_LMFIT_VERSION}")

         endif()

@@ -68,20 +69,15 @@

 endif()

 # Try to find lmfit, perhaps with help from pkg-config

-find_path(LMFIT_INCLUDE_DIRS lmcurve.h HINTS "${PC_LMFIT_INCLUDE_DIRS}" PATH_SUFFIXES include)

-find_library(LMFIT_LIBRARY_DIRS NAMES lmfit HINTS "${PC_LMFIT_LIBRARY_DIRS}" PATH_SUFFIXES lib64 lib)

+find_path(LMFIT_INCLUDE_DIR lmcurve.h HINTS "${PC_LMFIT_INCLUDE_DIRS}" PATH_SUFFIXES include)

+find_library(LMFIT_LIBRARY NAMES lmfit HINTS "${PC_LMFIT_LIBRARY_DIRS}" PATH_SUFFIXES lib64 lib)

 # Make sure we can also link, so that cross-compilation is properly supported

-if (LMFIT_INCLUDE_DIRS AND LMFIT_LIBRARY_DIRS)

+if (LMFIT_INCLUDE_DIR AND LMFIT_LIBRARY)

     include(CheckCXXSourceCompiles)

-    set(CMAKE_REQUIRED_INCLUDES ${LMFIT_INCLUDE_DIRS})

-    set(CMAKE_REQUIRED_LIBRARIES ${LMFIT_LIBRARY_DIRS})

+    set(CMAKE_REQUIRED_INCLUDES ${LMFIT_INCLUDE_DIR})

+    set(CMAKE_REQUIRED_LIBRARIES ${LMFIT_LIBRARY})

     check_cxx_source_compiles("#include <lmcurve.h>\nint main(){lmcurve(0,0,0,0,0,0,0,0);}" LMFIT_LINKS_OK)

-endif()

-

-if (LMFIT_LINKS_OK)

-    set(LMFIT_INCLUDE_DIR ${LMFIT_INCLUDE_DIRS})

-    set(LMFIT_LIBRARIES ${LMFIT_LIBRARY_DIRS})

 endif()

 include(FindPackageHandleStandardArgs)

@@ -90,11 +86,19 @@

     LMFIT_FOUND

     REQUIRED_VARS

     LMFIT_INCLUDE_DIR

-    LMFIT_LIBRARIES

+    LMFIT_LIBRARY

     LMFIT_LINKS_OK

     VERSION_VAR

     LMFIT_VERSION)

-mark_as_advanced(LMFIT_INCLUDE_DIRS LMFIT_LIBRARY_DIRS)

+mark_as_advanced(LMFIT_INCLUDE_DIR LMFIT_LIBRARY)

+

+if (LMFIT_FOUND)

+    add_library(lmfit INTERFACE IMPORTED)

+    set_target_properties(lmfit PROPERTIES

+        INTERFACE_INCLUDE_DIRECTORIES "${LMFIT_INCLUDE_DIR}"

+        INTERFACE_LINK_LIBRARIES "${LMFIT_LIBRARY}"

+        )

+endif()

 cmake_pop_check_state()

diff --git a/cmake/gmxManageLmfit.cmake b/cmake/gmxManageLmfit.cmake

index 182928f..5953a19 100644

--- a/cmake/gmxManageLmfit.cmake

+++ b/cmake/gmxManageLmfit.cmake

@@ -1,7 +1,7 @@

 #

 # This file is part of the GROMACS molecular simulation package.

 #

-# Copyright (c) 2016, by the GROMACS development team, led by

+# Copyright (c) 2016,2018, by the GROMACS development team, led by

 # Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,

 # and including many others, as listed in the AUTHORS file in the

 # top-level source directory and at http://www.gromacs.org.

@@ -32,21 +32,16 @@

 # To help us fund GROMACS development, we humbly ask that you cite

 # the research papers on the package. Check out http://www.gromacs.org.

-set(GMX_LMFIT_MINIMUM_REQUIRED_VERSION "6.1")

-set(GMX_BUNDLED_LMFIT_DIR "${CMAKE_SOURCE_DIR}/src/external/lmfit")

+# Note that lmfit does not have a stable API, so GROMACS only supports

+# the same version that it bundles.

+set(GMX_LMFIT_REQUIRED_VERSION "7.0")

-option(GMX_EXTERNAL_LMFIT "Use external lmfit instead of compiling the version bundled with GROMACS." OFF)

-mark_as_advanced(GMX_EXTERNAL_LMFIT)

+include(gmxOptionUtilities)

-macro(manage_lmfit)

-    if(GMX_EXTERNAL_LMFIT)

-        # Find an external lmfit library.

-        find_package(Lmfit ${GMX_LMFIT_MINIMUM_REQUIRED_VERSION})

-        if(NOT LMFIT_FOUND)

-            message(FATAL_ERROR "External lmfit could not be found, please adjust your pkg-config path to include the lmfit.pc file")

-        endif()

-    endif()

-endmacro()

+gmx_option_multichoice(GMX_USE_LMFIT

+    "Use external lmfit instead of compiling the version bundled with GROMACS."

+    INTERNAL

+    INTERNAL EXTERNAL NONE)

 macro(get_lmfit_properties LMFIT_SOURCES_VAR LMFIT_LIBRARIES_VAR LMFIT_INCLUDE_DIR_VAR LMFIT_INCLUDE_DIR_ORDER_VAR)

     if (GMX_EXTERNAL_LMFIT)

@@ -62,5 +57,29 @@

     endif()

 endmacro()

-manage_lmfit()

+function(manage_lmfit)

+    if(GMX_USE_LMFIT STREQUAL "INTERNAL")

+        set(BUNDLED_LMFIT_DIR "${CMAKE_SOURCE_DIR}/src/external/lmfit")

+        file(GLOB LMFIT_SOURCES ${BUNDLED_LMFIT_DIR}/*.cpp)

+        # Create a fake library for lmfit for libgromacs to depend on

+        add_library(lmfit INTERFACE)

+        target_sources(lmfit INTERFACE ${LMFIT_SOURCES})

+        target_include_directories(lmfit INTERFACE

+            $<BUILD_INTERFACE:${BUNDLED_LMFIT_DIR}>)

+        target_link_libraries(libgromacs PRIVATE lmfit)

+

+        set(HAVE_LMFIT_VALUE TRUE)

+    elseif(GMX_USE_LMFIT STREQUAL "EXTERNAL")

+        # Find an external lmfit library.

+        find_package(Lmfit ${GMX_LMFIT_MINIMUM_REQUIRED_VERSION})

+        if(NOT LMFIT_FOUND)

+            message(FATAL_ERROR "External lmfit could not be found, please adjust your pkg-config path to include the lmfit.pc file")

+        endif()

+        target_link_libraries(libgromacs PRIVATE lmfit)

+        set(HAVE_LMFIT_VALUE TRUE)

+    else()

+        set(HAVE_LMFIT_VALUE FALSE)

+    endif()

+    set(HAVE_LMFIT ${HAVE_LMFIT_VALUE} CACHE BOOL "Whether lmfit library support is enabled")

+endfunction()

diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt

index 9b38d75..3f050b5 100644

--- a/docs/CMakeLists.txt

+++ b/docs/CMakeLists.txt

@@ -366,7 +366,7 @@

             REQUIRED_CUDA_COMPUTE_CAPABILITY REGRESSIONTEST_VERSION

             SOURCE_MD5SUM REGRESSIONTEST_MD5SUM_STRING

             GMX_TNG_MINIMUM_REQUIRED_VERSION

-            GMX_LMFIT_MINIMUM_REQUIRED_VERSION

+            GMX_LMFIT_REQUIRED_VERSION

         COMMENT "Configuring Sphinx configuration file")

     gmx_add_sphinx_input_file(${SPHINX_CONFIG_VARS_FILE})

     gmx_add_sphinx_source_files(FILES ${SPHINX_SOURCE_FILES})

diff --git a/docs/conf-vars.py.cmakein b/docs/conf-vars.py.cmakein

index f4c3e02..038014a 100644

--- a/docs/conf-vars.py.cmakein

+++ b/docs/conf-vars.py.cmakein

@@ -1,7 +1,7 @@

 #

 # This file is part of the GROMACS molecular simulation package.

 #

-# Copyright (c) 2015,2016,2017, by the GROMACS development team, led by

+# Copyright (c) 2015,2016,2017,2018, by the GROMACS development team, led by

 # Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,

 # and including many others, as listed in the AUTHORS file in the

 # top-level source directory and at http://www.gromacs.org.

@@ -48,5 +48,5 @@

         ('SOURCE_MD5SUM', '@SOURCE_MD5SUM@'),

         ('REGRESSIONTEST_MD5SUM', '@REGRESSIONTEST_MD5SUM_STRING@'),

         ('GMX_TNG_MINIMUM_REQUIRED_VERSION', '@GMX_TNG_MINIMUM_REQUIRED_VERSION@'),

-        ('GMX_LMFIT_MINIMUM_REQUIRED_VERSION', '@GMX_LMFIT_MINIMUM_REQUIRED_VERSION@')

+        ('GMX_LMFIT_REQUIRED_VERSION', '@GMX_LMFIT_REQUIRED_VERSION@')

     ]

diff --git a/docs/install-guide/index.rst b/docs/install-guide/index.rst

index 48fe2c3..5324629 100644

--- a/docs/install-guide/index.rst

+++ b/docs/install-guide/index.rst

@@ -340,10 +340,14 @@

   by setting ``-DGMX_EXTERNAL_TNG=yes``, but TNG

   |GMX_TNG_MINIMUM_REQUIRED_VERSION| is bundled in the |Gromacs|

   source already.

-* An external lmfit library for Levenberg-Marquardt curve fitting

-  can be used by setting ``-DGMX_EXTERNAL_LMFIT=yes``, but lmfit

-  |GMX_LMFIT_MINIMUM_REQUIRED_VERSION| is bundled in the |Gromacs|

-  source already.

+* The lmfit library for Levenberg-Marquardt curve fitting is used in

+  |Gromacs|. Only lmfit |GMX_LMFIT_REQUIRED_VERSION| is supported.  A

+  reduced version of that library is bundled in the |Gromacs|

+  distribution, and the default build uses it. That default may be

+  explicitly enabled with ``-DGMX_USE_LMFIT=internal``. To use an

+  external lmfit library, set ``-DGMX_USE_LMFIT=external``, and adjust

+  ``CMAKE_PREFIX_PATH`` as needed.  lmfit support can be disabled with

+  ``-DGMX_USE_LMFIT=none``.

 * zlib is used by TNG for compressing some kinds of trajectory data

 * Building the |Gromacs| documentation is optional, and requires

   ImageMagick, pdflatex, bibtex, doxygen, python 2.7, sphinx 

diff --git a/src/config.h.cmakein b/src/config.h.cmakein

index b7e5964..d0c3c35 100644

--- a/src/config.h.cmakein

+++ b/src/config.h.cmakein

@@ -377,6 +377,9 @@

 /* Define if we have feenableexcept */

 #cmakedefine01 HAVE_FEENABLEEXCEPT

+/* Define if we have lmfit support */

+#cmakedefine01 HAVE_LMFIT

+

 /*! \endcond */

 #endif

diff --git a/src/external/lmfit/README b/src/external/lmfit/README

index 3c0557b..4b8d7c6 100644

--- a/src/external/lmfit/README

+++ b/src/external/lmfit/README

@@ -1,5 +1,5 @@

 This directory contains only the barebones version of the

-lmfit package, version 6.1. Full version is available from

+lmfit package, version 7.0. Full version is available from

 http://apps.jcns.fz-juelich.de/doku/sc/lmfit

 The license under which lmfit is bundled may be found in the GROMACS

diff --git a/src/external/lmfit/lmmin.cpp b/src/external/lmfit/lmmin.cpp

index eb2d298..813fe5d 100644

--- a/src/external/lmfit/lmmin.cpp

+++ b/src/external/lmfit/lmmin.cpp

@@ -16,40 +16,41 @@

 #include <assert.h>

 #include <stdlib.h>

 #include <stdio.h>

-#include <cmath>

+#include <math.h>

 #include <float.h>

 #include "lmmin.h"

-using namespace std;

-

-static double lm_enorm(int n, const double* x);

-

-#define MIN(a, b) (((a) <= (b)) ? (a) : (b))

-#define MAX(a, b) (((a) >= (b)) ? (a) : (b))

-#define SQR(x) (x) * (x)

+#define MIN(a,b) (((a)<=(b)) ? (a) : (b))

+#define MAX(a,b) (((a)>=(b)) ? (a) : (b))

+#define SQR(x)   (x)*(x)

 /* Declare functions that do the heavy numerics.

    Implementions are in this source file, below lmmin.

-   Dependences: lmmin calls lmpar, which calls qrfac and qrsolv. */

-static void lm_lmpar(const int n, double* r, const int ldr, const int* Pivot,

-                     const double* diag, const double* qtb, const double delta,

-                     double* par, double* x, double* Sdiag, double* aux, double* xdi);

-static void lm_qrfac(const int m, const int n, double* A, int* Pivot, double* Rdiag,

-                     double* Acnorm, double* W);

-static void lm_qrsolv(const int n, double* r, const int ldr, const int* Pivot,

-                      const double* diag, const double* qtb, double* x,

-                      double* Sdiag, double* W);

+   Dependences: lmmin calls qrfac and lmpar; lmpar calls qrsolv. */

+void lm_lmpar(

+    const int n, double *const r, const int ldr, int *const ipvt,

+    double *const diag, double *const qtb, double delta, double *const par,

+    double *const x,

+    double *const sdiag, double *const aux, double *const xdi);

+void lm_qrfac(

+    const int m, const int n, double *const a, int *const ipvt,

+    double *const rdiag, double *const acnorm, double *const wa );

+void lm_qrsolv(

+    const int n, double *const r, const int ldr, int *const ipvt,

+    double *const diag, double *const qtb, double *const x,

+    double *const sdiag, double *const wa);

-/******************************************************************************/

-/*  Numeric constants                                                         */

-/******************************************************************************/

-/* Set machine-dependent constants to values from float.h. */

-#define LM_MACHEP DBL_EPSILON       /* resolution of arithmetic */

-#define LM_DWARF DBL_MIN            /* smallest nonzero number */

+/*****************************************************************************/

+/*  Numeric constants                                                        */

+/*****************************************************************************/

+

+/* machine-dependent constants from float.h */

+#define LM_MACHEP     DBL_EPSILON   /* resolution of arithmetic */

+#define LM_DWARF      DBL_MIN       /* smallest nonzero number */

 #define LM_SQRT_DWARF sqrt(DBL_MIN) /* square should not underflow */

 #define LM_SQRT_GIANT sqrt(DBL_MAX) /* square should not overflow */

-#define LM_USERTOL 30 * LM_MACHEP   /* users are recommended to require this */

+#define LM_USERTOL    30*LM_MACHEP  /* users are recommended to require this */

 /* If the above values do not work, the following seem good for an x86:

  LM_MACHEP     .555e-16

@@ -65,19 +66,19 @@

  LM_USER_TOL   1.e-14

 */

-/* Predefined control parameter sets (msgfile=NULL means stdout). */

 const lm_control_struct lm_control_double = {

-    LM_USERTOL, LM_USERTOL, LM_USERTOL, LM_USERTOL,

-    100., 100, 1, NULL, 0, -1, -1};

+    LM_USERTOL, LM_USERTOL, LM_USERTOL, LM_USERTOL, 100., 100, 1,

+    NULL, 0, -1, -1 };

 const lm_control_struct lm_control_float = {

-    1.e-7, 1.e-7, 1.e-7, 1.e-7,

-    100., 100, 1, NULL, 0, -1, -1};

+    1.e-7,      1.e-7,      1.e-7,      1.e-7,      100., 100, 1,

+    NULL, 0, -1, -1 };

-/******************************************************************************/

-/*  Message texts (indexed by status.info)                                    */

-/******************************************************************************/

-const char* lm_infmsg[] = {

+/*****************************************************************************/

+/*  Message texts (indexed by status.info)                                   */

+/*****************************************************************************/

+

+const char *lm_infmsg[] = {

     "found zero (sum of squares below underflow limit)",

     "converged  (the relative error in the sum of squares is at most tol)",

     "converged  (the relative error of the parameter vector is at most tol)",

@@ -90,66 +91,86 @@

     "crashed    (not enough memory)",

     "exploded   (fatal coding error: improper input parameters)",

     "stopped    (break requested within function evaluation)",

-    "found nan  (function value is not-a-number or infinite)"};

+    "found nan  (function value is not-a-number or infinite)",

+    "won't fit  (no free parameter)"

+};

-/******************************************************************************/

-/*  Monitoring auxiliaries.                                                   */

-/******************************************************************************/

+const char *lm_shortmsg[] = {

+    "found zero",      //  0

+    "converged (f)",   //  1

+    "converged (p)",   //  2

+    "converged (2)",   //  3

+    "degenerate",      //  4

+    "call limit",      //  5

+    "failed (f)",      //  6

+    "failed (p)",      //  7

+    "failed (o)",      //  8

+    "no memory",       //  9

+    "invalid input",   // 10

+    "user break",      // 11

+    "found nan",       // 12

+    "no free par"      // 13

+};

-static void lm_print_pars(int nout, const double* par, FILE* fout)

+

+/*****************************************************************************/

+/*  Monitoring auxiliaries.                                                  */

+/*****************************************************************************/

+

+static void lm_print_pars(const int nout, const double *par, FILE* fout)

 {

-    int i;

-    for (i = 0; i < nout; ++i)

-        fprintf(fout, " %16.9g", par[i]);

+    fprintf(fout, "  pars:");

+    for (int i = 0; i < nout; ++i)

+        fprintf(fout, " %23.16g", par[i]);

     fprintf(fout, "\n");

 }

-/******************************************************************************/

-/*  lmmin (main minimization routine)                                         */

-/******************************************************************************/

-void lmmin(const int n, double* x, const int m, const void* data,

-           void (*evaluate)(const double* par, const int m_dat,

-                            const void* data, double* fvec, int* userbreak),

-           const lm_control_struct* C, lm_status_struct* S)

+/*****************************************************************************/

+/*  lmmin (main minimization routine)                                        */

+/*****************************************************************************/

+

+void lmmin(

+    const int n, double *const x, const int m, const double* y,

+    const void *const data,

+    void (*const evaluate)(

+        const double *const par, const int m_dat, const void *const data,

+        double *const fvec, int *const userbreak),

+    const lm_control_struct *const C, lm_status_struct *const S)

 {

     int j, i;

-    double actred, dirder, fnorm, fnorm1, gnorm, pnorm, prered, ratio, step,

-        sum, temp, temp1, temp2, temp3;

-

-    /***  Initialize internal variables.  ***/

+    double actred, dirder, fnorm, fnorm1, gnorm, pnorm,

+        prered, ratio, step, sum, temp, temp1, temp2, temp3;

+    static double p1 = 0.1, p0001 = 1.0e-4;

     int maxfev = C->patience * (n+1);

-    int inner_success; /* flag for loop control */

-    double lmpar = 0;  /* Levenberg-Marquardt parameter */

+    int    inner_success; /* flag for loop control */

+    double lmpar = 0;     /* Levenberg-Marquardt parameter */

     double delta = 0;

     double xnorm = 0;

     double eps = sqrt(MAX(C->epsilon, LM_MACHEP)); /* for forward differences */

-    int nout = C->n_maxpri == -1 ? n : MIN(C->n_maxpri, n);

+    int nout = C->n_maxpri==-1 ? n : MIN(C->n_maxpri, n);

-    /* Reinterpret C->msgfile=NULL as stdout (which is unavailable for

-       compile-time initialization of lm_control_double and similar). */

+    /* The workaround msgfile=NULL is needed for default initialization */

     FILE* msgfile = C->msgfile ? C->msgfile : stdout;

-    /***  Default status info; must be set before first return statement.  ***/

-

-    S->outcome = 0; /* status code */

+    /* Default status info; must be set ahead of first return statements */

+    S->outcome = 0;      /* status code */

     S->userbreak = 0;

-    S->nfev = 0; /* function evaluation counter */

+    S->nfev = 0;      /* function evaluation counter */

-    /***  Check input parameters for errors.  ***/

+/***  Check input parameters for errors.  ***/

-    if (n <= 0) {

+    if ( n < 0 ) {

         fprintf(stderr, "lmmin: invalid number of parameters %i\n", n);

-        S->outcome = 10;

+        S->outcome = 10; /* invalid parameter */

         return;

     }

     if (m < n) {

         fprintf(stderr, "lmmin: number of data points (%i) "

-                        "smaller than number of parameters (%i)\n",

-                m, n);

+                "smaller than number of parameters (%i)\n", m, n);

         S->outcome = 10;

         return;

     }

@@ -173,93 +194,99 @@

     }

     if (C->scale_diag != 0 && C->scale_diag != 1) {

         fprintf(stderr, "lmmin: logical variable scale_diag=%i, "

-                        "should be 0 or 1\n",

-                C->scale_diag);

+                "should be 0 or 1\n", C->scale_diag);

         S->outcome = 10;

         return;

     }

-    /***  Allocate work space.  ***/

+/***  Allocate work space.  ***/

     /* Allocate total workspace with just one system call */

-    char* ws;

-    if ((ws = (char *)malloc((2*m + 5*n + m*n) * sizeof(double) +

-                             n * sizeof(int))) == NULL) {

+    char *ws;

+    if ( ( ws = static_cast<char *>(malloc(

+               (2*m+5*n+m*n)*sizeof(double) + n*sizeof(int)) ) ) == NULL ) {

         S->outcome = 9;

         return;

     }

     /* Assign workspace segments. */

-    char* pws = ws;

-    double* fvec = (double*)pws;

-    pws += m * sizeof(double) / sizeof(char);

-    double* diag = (double*)pws;

-    pws += n * sizeof(double) / sizeof(char);

-    double* qtf = (double*)pws;

-    pws += n * sizeof(double) / sizeof(char);

-    double* fjac = (double*)pws;

-    pws += n * m * sizeof(double) / sizeof(char);

-    double* wa1 = (double*)pws;

-    pws += n * sizeof(double) / sizeof(char);

-    double* wa2 = (double*)pws;

-    pws += n * sizeof(double) / sizeof(char);

-    double* wa3 = (double*)pws;

-    pws += n * sizeof(double) / sizeof(char);

-    double* wf = (double*)pws;

-    pws += m * sizeof(double) / sizeof(char);

-    int* Pivot = (int*)pws;

-    //pws += n * sizeof(int) / sizeof(char);

+    char *pws = ws;

+    double *fvec = (double*) pws; pws += m * sizeof(double)/sizeof(char);

+    double *diag = (double*) pws; pws += n * sizeof(double)/sizeof(char);

+    double *qtf  = (double*) pws; pws += n * sizeof(double)/sizeof(char);

+    double *fjac = (double*) pws; pws += n*m*sizeof(double)/sizeof(char);

+    double *wa1  = (double*) pws; pws += n * sizeof(double)/sizeof(char);

+    double *wa2  = (double*) pws; pws += n * sizeof(double)/sizeof(char);

+    double *wa3  = (double*) pws; pws += n * sizeof(double)/sizeof(char);

+    double *wf   = (double*) pws; pws += m * sizeof(double)/sizeof(char);

+    int    *ipvt = (int*)    pws; pws += n * sizeof(int)   /sizeof(char);

-    /* Initialize diag. */

-    if (!C->scale_diag)

+    /* Initialize diag */ // TODO: check whether this is still needed

+    if (!C->scale_diag) {

         for (j = 0; j < n; j++)

-            diag[j] = 1;

-

-    /***  Evaluate function at starting point and calculate norm.  ***/

-

-    if (C->verbosity) {

-        fprintf(msgfile, "lmmin start ");

-        lm_print_pars(nout, x, msgfile);

+            diag[j] = 1.;

     }

-    (*evaluate)(x, m, data, fvec, &(S->userbreak));

-    if (C->verbosity > 4)

-        for (i = 0; i < m; ++i)

-            fprintf(msgfile, "    fvec[%4i] = %18.8g\n", i, fvec[i]);

-    S->nfev = 1;

-    if (S->userbreak)

-        goto terminate;

-    fnorm = lm_enorm(m, fvec);

-    if (C->verbosity)

-        fprintf(msgfile, "  fnorm = %18.8g\n", fnorm);

-    if (!std::isfinite(fnorm)) {

+/***  Evaluate function at starting point and calculate norm.  ***/

+

+    if( C->verbosity&1 )

+        fprintf(msgfile, "lmmin start (ftol=%g gtol=%g xtol=%g)\n",

+                C->ftol, C->gtol, C->xtol);

+    if( C->verbosity&2 )

+        lm_print_pars(nout, x, msgfile);

+    (*evaluate)(x, m, data, fvec, &(S->userbreak));

+    if( C->verbosity&8 )

+    {

+        if (y) {

+            for( i=0; i

+                fprintf(msgfile, "    i, f, y-f: %4i %18.8g %18.8g\n",

+                        i, fvec[i], y[i]-fvec[i]);

+        } else {

+            for( i=0; i

+                fprintf(msgfile, "    i, f: %4i %18.8g\n", i, fvec[i]);

+        }

+    }

+    S->nfev = 1;

+    if ( S->userbreak )

+        goto terminate;

+    if ( n == 0 ) {

+        S->outcome = 13; /* won't fit */

+        goto terminate;

+    }

+    fnorm = lm_fnorm(m, fvec, y);

+    if( C->verbosity&2 )

+        fprintf(msgfile, "  fnorm = %24.16g\n", fnorm);

+    if( !isfinite(fnorm) ){

+        if( C->verbosity )

+            fprintf(msgfile, "nan case 1\n");

         S->outcome = 12; /* nan */

         goto terminate;

-    } else if (fnorm <= LM_DWARF) {

+    } else if( fnorm <= LM_DWARF ){

         S->outcome = 0; /* sum of squares almost zero, nothing to do */

         goto terminate;

     }

-    /***  The outer loop: compute gradient, then descend.  ***/

+/***  The outer loop: compute gradient, then descend.  ***/

-    for (int outer = 0;; ++outer) {

+    for( int outer=0; ; ++outer ) {

-        /** Calculate the Jacobian. **/

+/***  [outer]  Calculate the Jacobian.  ***/

+

         for (j = 0; j < n; j++) {

             temp = x[j];

-            step = MAX(eps * eps, eps * fabs(temp));

+            step = MAX(eps*eps, eps * fabs(temp));

             x[j] += step; /* replace temporarily */

             (*evaluate)(x, m, data, wf, &(S->userbreak));

             ++(S->nfev);

-            if (S->userbreak)

+            if ( S->userbreak )

                 goto terminate;

             for (i = 0; i < m; i++)

                 fjac[j*m+i] = (wf[i] - fvec[i]) / step;

             x[j] = temp; /* restore */

         }

-        if (C->verbosity >= 10) {

+        if ( C->verbosity&16 ) {

             /* print the entire matrix */

-            printf("\nlmmin Jacobian\n");

+            printf("Jacobian\n");

             for (i = 0; i < m; i++) {

                 printf("  ");

                 for (j = 0; j < n; j++)

@@ -268,34 +295,39 @@

             }

         }

-        /** Compute the QR factorization of the Jacobian. **/

+/***  [outer]  Compute the QR factorization of the Jacobian.  ***/

-        /* fjac is an m by n array. The upper n by n submatrix of fjac is made

-         *   to contain an upper triangular matrix R with diagonal elements of

-         *   nonincreasing magnitude such that

-         *

-         *         P^T*(J^T*J)*P = R^T*R

-         *

-         *         (NOTE: ^T stands for matrix transposition),

-         *

-         *   where P is a permutation matrix and J is the final calculated

-         *   Jacobian. Column j of P is column Pivot(j) of the identity matrix.

-         *   The lower trapezoidal part of fjac contains information generated

-         *   during the computation of R.

-         *

-         * Pivot is an integer array of length n. It 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 Pivot(j) of the identity matrix.

-         */

+/*      fjac is an m by n array. The upper n by n submatrix of fjac

+ *        is made to contain an upper triangular matrix R with diagonal

+ *        elements of nonincreasing magnitude such that

+ *

+ *              P^T*(J^T*J)*P = R^T*R

+ *

+ *              (NOTE: ^T stands for matrix transposition),

+ *

+ *        where P is a permutation matrix and J is the final calculated

+ *        Jacobian. Column j of P is column ipvt(j) of the identity matrix.

+ *        The lower trapezoidal part of fjac contains information generated

+ *        during the computation of R.

+ *

+ *      ipvt is an integer array of length n. It 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.

+ */