diff --git a/.gitignore b/.gitignore index 93f6241..16e1189 100644 --- a/.gitignore +++ b/.gitignore @@ -57,3 +57,5 @@ manual-4.0.pdf /gromacs-2018.1.tar.gz /manual-2018.1.pdf /regressiontests-2018.1.tar.gz +/gromacs-2018.2.tar.gz +/manual-2018.2.pdf diff --git a/facb927.diff b/facb927.diff new file mode 100644 index 0000000..0764a91 --- /dev/null +++ b/facb927.diff @@ -0,0 +1,2421 @@ +From facb9273c571b6bef430e71b31ff229abd6d1a46 Mon Sep 17 00:00:00 2001 +From: Mark Abraham +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 \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 ++ $) ++ 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 + #include + #include +-#include ++#include + #include + #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(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; infev = 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. ++ */ + +- lm_qrfac(m, n, fjac, Pivot, wa1, wa2, wa3); +- /* return values are Pivot, wa1=rdiag, wa2=acnorm */ ++ lm_qrfac(m, n, fjac, ipvt, wa1, wa2, wa3); ++ /* return values are ipvt, wa1=rdiag, wa2=acnorm */ + +- /** Form Q^T * fvec, and store first n components in qtf. **/ +- for (i = 0; i < m; i++) +- wf[i] = fvec[i]; ++/*** [outer] Form Q^T * fvec, and store first n components in qtf. ***/ ++ ++ if (y) ++ for (i = 0; i < m; i++) ++ wf[i] = fvec[i] - y[i]; ++ else ++ for (i = 0; i < m; i++) ++ wf[i] = fvec[i]; + + for (j = 0; j < n; j++) { + temp3 = fjac[j*m+j]; +@@ -311,15 +343,16 @@ + qtf[j] = wf[j]; + } + +- /** Compute norm of scaled gradient and detect degeneracy. **/ ++/*** [outer] Compute norm of scaled gradient and detect degeneracy. ***/ ++ + gnorm = 0; + for (j = 0; j < n; j++) { +- if (wa2[Pivot[j]] == 0) ++ if (wa2[ipvt[j]] == 0) + continue; + sum = 0; + for (i = 0; i <= j; i++) + sum += fjac[j*m+i] * qtf[i]; +- gnorm = MAX(gnorm, fabs(sum / wa2[Pivot[j]] / fnorm)); ++ gnorm = MAX(gnorm, fabs( sum / wa2[ipvt[j]] / fnorm )); + } + + if (gnorm <= C->gtol) { +@@ -327,8 +360,10 @@ + goto terminate; + } + +- /** Initialize or update diag and delta. **/ +- if (!outer) { /* first iteration only */ ++/*** [outer] Initialize / update diag and delta. ***/ ++ ++ if ( !outer ) { ++ /* first iteration only */ + if (C->scale_diag) { + /* diag := norms of the columns of the initial Jacobian */ + for (j = 0; j < n; j++) +@@ -337,129 +372,140 @@ + for (j = 0; j < n; j++) + wa3[j] = diag[j] * x[j]; + xnorm = lm_enorm(n, wa3); +- if (C->verbosity >= 2) { +- fprintf(msgfile, "lmmin diag "); +- lm_print_pars(nout, x, msgfile); // xnorm +- fprintf(msgfile, " xnorm = %18.8g\n", xnorm); +- } +- /* Only now print the header for the loop table. */ +- if (C->verbosity >= 3) { +- fprintf(msgfile, " o i lmpar prered" +- " ratio dirder delta" +- " pnorm fnorm"); +- for (i = 0; i < nout; ++i) +- fprintf(msgfile, " p%i", i); +- fprintf(msgfile, "\n"); +- } + } else { + xnorm = lm_enorm(n, x); + } +- if (!std::isfinite(xnorm)) { ++ if( !isfinite(xnorm) ){ ++ if( C->verbosity ) ++ fprintf(msgfile, "nan case 2\n"); + S->outcome = 12; /* nan */ + goto terminate; + } +- /* Initialize the step bound delta. */ +- if (xnorm) ++ /* initialize the step bound delta. */ ++ if ( xnorm ) + delta = C->stepbound * xnorm; + else + delta = C->stepbound; ++ /* only now print the header for the loop table */ ++ if( C->verbosity&2 ) { ++ fprintf(msgfile, " #o #i lmpar prered actred" ++ " ratio dirder delta" ++ " pnorm fnorm"); ++ for (i = 0; i < nout; ++i) ++ fprintf(msgfile, " p%i", i); ++ fprintf(msgfile, "\n"); ++ } + } else { + if (C->scale_diag) { + for (j = 0; j < n; j++) +- diag[j] = MAX(diag[j], wa2[j]); ++ diag[j] = MAX( diag[j], wa2[j] ); + } + } + +- /** The inner loop. **/ ++/*** The inner loop. ***/ + int inner = 0; + do { + +- /** Determine the Levenberg-Marquardt parameter. **/ +- lm_lmpar(n, fjac, m, Pivot, diag, qtf, delta, &lmpar, ++/*** [inner] Determine the Levenberg-Marquardt parameter. ***/ ++ ++ lm_lmpar(n, fjac, m, ipvt, diag, qtf, delta, &lmpar, + wa1, wa2, wf, wa3); + /* used return values are fjac (partly), lmpar, wa1=x, wa3=diag*x */ + +- /* Predict scaled reduction. */ ++ /* predict scaled reduction */ + pnorm = lm_enorm(n, wa3); +- if (!std::isfinite(pnorm)) { ++ if( !isfinite(pnorm) ){ ++ if( C->verbosity ) ++ fprintf(msgfile, "nan case 3\n"); + S->outcome = 12; /* nan */ + goto terminate; + } +- temp2 = lmpar * SQR(pnorm / fnorm); ++ temp2 = lmpar * SQR( pnorm / fnorm ); + for (j = 0; j < n; j++) { + wa3[j] = 0; + for (i = 0; i <= j; i++) +- wa3[i] -= fjac[j*m+i] * wa1[Pivot[j]]; ++ wa3[i] -= fjac[j*m+i] * wa1[ipvt[j]]; + } +- temp1 = SQR(lm_enorm(n, wa3) / fnorm); +- if (!std::isfinite(temp1)) { ++ temp1 = SQR( lm_enorm(n, wa3) / fnorm ); ++ if( !isfinite(temp1) ){ ++ if( C->verbosity ) ++ fprintf(msgfile, "nan case 4\n"); + S->outcome = 12; /* nan */ + goto terminate; + } +- prered = temp1 + 2*temp2; ++ prered = temp1 + 2 * temp2; + dirder = -temp1 + temp2; /* scaled directional derivative */ + +- /* At first call, adjust the initial step bound. */ +- if (!outer && pnorm < delta) ++ /* at first call, adjust the initial step bound. */ ++ if ( !outer && !inner && pnorm < delta ) + delta = pnorm; + +- /** Evaluate the function at x + p. **/ ++/*** [inner] Evaluate the function at x + p. ***/ ++ + for (j = 0; j < n; j++) + wa2[j] = x[j] - wa1[j]; +- (*evaluate)(wa2, m, data, wf, &(S->userbreak)); ++ ++ (*evaluate)( wa2, m, data, wf, &(S->userbreak) ); + ++(S->nfev); +- if (S->userbreak) ++ if ( S->userbreak ) + goto terminate; +- fnorm1 = lm_enorm(m, wf); +- if (!std::isfinite(fnorm1)) { +- S->outcome = 12; /* nan */ +- goto terminate; ++ fnorm1 = lm_fnorm(m, wf, y); ++ // exceptionally, for this norm we do not test for infinity ++ // because we can deal with it without terminating. ++ ++/*** [inner] Evaluate the scaled reduction. ***/ ++ ++ /* actual scaled reduction (supports even the case fnorm1=infty) */ ++ if (p1 * fnorm1 < fnorm) ++ actred = 1 - SQR(fnorm1 / fnorm); ++ else ++ actred = -1; ++ ++ /* ratio of actual to predicted reduction */ ++ ratio = prered ? actred/prered : 0; ++ ++ if( C->verbosity&32 ) { ++ if (y) { ++ for( i=0; iverbosity == 2) { +- fprintf(msgfile, "lmmin (%i:%i) ", outer, inner); +- lm_print_pars(nout, wa2, msgfile); // fnorm1, +- } else if (C->verbosity >= 3) { +- printf("%3i %2i %9.2g %9.2g %14.6g" ++ if( C->verbosity&2 ) { ++ printf("%3i %2i %9.2g %9.2g %9.2g %14.6g" + " %9.2g %10.3e %10.3e %21.15e", +- outer, inner, lmpar, prered, ratio, ++ outer, inner, lmpar, prered, actred, ratio, + dirder, delta, pnorm, fnorm1); + for (i = 0; i < nout; ++i) + fprintf(msgfile, " %16.9g", wa2[i]); + fprintf(msgfile, "\n"); + } + +- /* Update the step bound. */ +- if (ratio <= 0.25) { +- if (actred >= 0) +- temp = 0.5; +- else if (actred > -99) /* -99 = 1-1/0.1^2 */ +- temp = MAX(dirder / (2*dirder + actred), 0.1); +- else +- temp = 0.1; +- delta = temp * MIN(delta, pnorm / 0.1); +- lmpar /= temp; +- } else if (ratio >= 0.75) { +- delta = 2 * pnorm; +- lmpar *= 0.5; +- } else if (!lmpar) { +- delta = 2 * pnorm; +- } ++ /* update the step bound */ ++ if (ratio <= 0.25) { ++ if (actred >= 0) ++ temp = 0.5; ++ else ++ temp = 0.5 * dirder / (dirder + 0.5 * actred); ++ if (p1 * fnorm1 >= fnorm || temp < p1) ++ temp = p1; ++ delta = temp * MIN(delta, pnorm / p1); ++ lmpar /= temp; ++ } else if (lmpar == 0 || ratio >= 0.75) { ++ delta = 2 * pnorm; ++ lmpar *= 0.5; ++ } + +- /** On success, update solution, and test for convergence. **/ ++/*** [inner] On success, update solution, and test for convergence. ***/ + +- inner_success = ratio >= 1e-4; +- if (inner_success) { ++ inner_success = ratio >= p0001; ++ if ( inner_success ) { + +- /* Update x, fvec, and their norms. */ ++ /* update x, fvec, and their norms */ + if (C->scale_diag) { + for (j = 0; j < n; j++) { + x[j] = wa2[j]; +@@ -472,172 +518,193 @@ + for (i = 0; i < m; i++) + fvec[i] = wf[i]; + xnorm = lm_enorm(n, wa2); +- if (!std::isfinite(xnorm)) { ++ if( !isfinite(xnorm) ){ ++ if( C->verbosity ) ++ fprintf(msgfile, "nan case 6\n"); + S->outcome = 12; /* nan */ + goto terminate; + } + fnorm = fnorm1; + } + +- /* Convergence tests. */ ++ /* convergence tests */ + S->outcome = 0; +- if (fnorm <= LM_DWARF) +- goto terminate; /* success: sum of squares almost zero */ +- /* Test two criteria (both may be fulfilled). */ ++ if( fnorm<=LM_DWARF ) ++ goto terminate; /* success: sum of squares almost zero */ ++ /* test two criteria (both may be fulfilled) */ + if (fabs(actred) <= C->ftol && prered <= C->ftol && ratio <= 2) +- S->outcome = 1; /* success: x almost stable */ ++ S->outcome = 1; /* success: x almost stable */ + if (delta <= C->xtol * xnorm) + S->outcome += 2; /* success: sum of squares almost stable */ + if (S->outcome != 0) { + goto terminate; + } + +- /** Tests for termination and stringent tolerances. **/ +- if (S->nfev >= maxfev) { ++/*** [inner] Tests for termination and stringent tolerances. ***/ ++ ++ if ( S->nfev >= maxfev ){ + S->outcome = 5; + goto terminate; + } +- if (fabs(actred) <= LM_MACHEP && prered <= LM_MACHEP && +- ratio <= 2) { ++ if ( fabs(actred) <= LM_MACHEP && ++ prered <= LM_MACHEP && ratio <= 2 ){ + S->outcome = 6; + goto terminate; + } +- if (delta <= LM_MACHEP * xnorm) { ++ if ( delta <= LM_MACHEP*xnorm ){ + S->outcome = 7; + goto terminate; + } +- if (gnorm <= LM_MACHEP) { ++ if ( gnorm <= LM_MACHEP ){ + S->outcome = 8; + goto terminate; + } + +- /** End of the inner loop. Repeat if iteration unsuccessful. **/ +- ++inner; +- } while (!inner_success); ++/*** [inner] End of the loop. Repeat if iteration unsuccessful. ***/ + +- }; /*** End of the outer loop. ***/ ++ ++inner; ++ } while ( !inner_success ); ++ ++/*** [outer] End of the loop. ***/ ++ ++ }; + + terminate: +- S->fnorm = lm_enorm(m, fvec); +- if (C->verbosity >= 2) +- printf("lmmin outcome (%i) xnorm %g ftol %g xtol %g\n", S->outcome, +- xnorm, C->ftol, C->xtol); +- if (C->verbosity & 1) { +- fprintf(msgfile, "lmmin final "); +- lm_print_pars(nout, x, msgfile); // S->fnorm, +- fprintf(msgfile, " fnorm = %18.8g\n", S->fnorm); ++ S->fnorm = lm_fnorm(m, fvec, y); ++ if( C->verbosity&1 ) ++ fprintf(msgfile, "lmmin terminates with outcome %i\n", S->outcome); ++ if( C->verbosity&2 ) ++ lm_print_pars(nout, x, msgfile); ++ if( C->verbosity&8 ) { ++ if (y) { ++ for( i=0; iuserbreak) /* user-requested break */ ++ if( C->verbosity&2 ) ++ fprintf(msgfile, " fnorm=%24.16g xnorm=%24.16g\n", S->fnorm, xnorm); ++ if ( S->userbreak ) /* user-requested break */ + S->outcome = 11; + +- /*** Deallocate the workspace. ***/ ++/*** Deallocate the workspace. ***/ + free(ws); + + } /*** lmmin. ***/ + +-/******************************************************************************/ +-/* lm_lmpar (determine Levenberg-Marquardt parameter) */ +-/******************************************************************************/ + +-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) ++/*****************************************************************************/ ++/* lm_lmpar (determine Levenberg-Marquardt parameter) */ ++/*****************************************************************************/ ++ ++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) ++{ + /* 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 parameter value par such that if x solves the system + * + * A*x = b and sqrt(par)*D*x = 0 + * +- * in the least squares sense, and dxnorm is the Euclidean norm of D*x, +- * then either par=0 and (dxnorm-delta) < 0.1*delta, or par>0 and +- * abs(dxnorm-delta) < 0.1*delta. ++ * in the least squares sense, and dxnorm is the euclidean ++ * norm of D*x, then either par=0 and (dxnorm-delta) < 0.1*delta, ++ * or par>0 and abs(dxnorm-delta) < 0.1*delta. + * + * Using lm_qrsolv, 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^T*b. On output +- * lmpar also provides an upper triangular matrix S such that ++ * 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^T*b. On output lmpar also provides an ++ * upper triangular matrix S such that + * + * P^T*(A^T*A + par*D*D)*P = S^T*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. ++ * 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. + * + * Parameters: + * + * 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. ++ * 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. ++ * ldr is a positive integer INPUT variable not less than n ++ * which specifies the leading dimension of the array R. + * +- * Pivot 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 +- * Pivot(j) of the identity matrix. ++ * 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. ++ * 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^T*b. + * +- * delta is a positive INPUT variable which specifies an upper bound +- * on the Euclidean norm of D*x. ++ * 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. ++ * 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. ++ * 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 array of length n needed as workspace; on OUTPUT it +- * contains the diagonal elements of the upper triangular matrix S. ++ * sdiag is an array of length n needed as workspace; on OUTPUT ++ * it contains the diagonal elements of the upper triangular ++ * matrix S. + * + * aux is a multi-purpose work array of length n. + * + * xdi is a work array of length n. On OUTPUT: diag[j] * x[j]. + * + */ +-{ + int i, iter, j, nsing; + double dxnorm, fp, fp_old, gnorm, parc, parl, paru; + double sum, temp; + static double p1 = 0.1; + +- /*** Compute and store in x the Gauss-Newton direction. If the Jacobian +- is rank-deficient, obtain a least-squares solution. ***/ ++/*** lmpar: 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++) { + aux[j] = qtb[j]; +- if (r[j*ldr+j] == 0 && nsing == n) ++ if (r[j * ldr + j] == 0 && nsing == n) + nsing = j; + if (nsing < n) + aux[j] = 0; + } +- for (j = nsing-1; j >= 0; j--) { +- aux[j] = aux[j] / r[j+ldr*j]; ++ for (j = nsing - 1; j >= 0; j--) { ++ aux[j] = aux[j] / r[j + ldr * j]; + temp = aux[j]; + for (i = 0; i < j; i++) +- aux[i] -= r[j*ldr+i] * temp; ++ aux[i] -= r[j * ldr + i] * temp; + } + + for (j = 0; j < n; j++) +- x[Pivot[j]] = aux[j]; ++ x[ipvt[j]] = aux[j]; + +- /*** Initialize the iteration counter, evaluate the function at the origin, +- and test for acceptance of the Gauss-Newton direction. ***/ ++/*** lmpar: initialize the iteration counter, evaluate the function at the ++ origin, and test for acceptance of the gauss-newton direction. ***/ + + for (j = 0; j < n; j++) + xdi[j] = diag[j] * x[j]; +@@ -645,65 +712,67 @@ + fp = dxnorm - delta; + if (fp <= p1 * delta) { + #ifdef LMFIT_DEBUG_MESSAGES +- printf("debug lmpar nsing=%d, n=%d, terminate[fp<=p1*del]\n", nsing, n); ++ printf("debug lmpar nsing %d n %d, terminate (fp= n) { + for (j = 0; j < n; j++) +- aux[j] = diag[Pivot[j]] * xdi[Pivot[j]] / dxnorm; ++ aux[j] = diag[ipvt[j]] * xdi[ipvt[j]] / dxnorm; + + for (j = 0; j < n; j++) { + sum = 0; + for (i = 0; i < j; i++) +- sum += r[j*ldr+i] * aux[i]; +- aux[j] = (aux[j] - sum) / r[j+ldr*j]; ++ sum += r[j * ldr + i] * aux[i]; ++ aux[j] = (aux[j] - sum) / r[j + ldr * j]; + } + temp = lm_enorm(n, aux); + parl = fp / delta / temp / temp; + } + +- /*** Calculate an upper bound, paru, for the zero of the function. ***/ ++/*** lmpar: calculate an upper bound, paru, for the zero of the function. ***/ + + for (j = 0; j < n; j++) { + sum = 0; + for (i = 0; i <= j; i++) +- sum += r[j*ldr+i] * qtb[i]; +- aux[j] = sum / diag[Pivot[j]]; ++ sum += r[j * ldr + i] * qtb[i]; ++ aux[j] = sum / diag[ipvt[j]]; + } + gnorm = lm_enorm(n, aux); + paru = gnorm / delta; + if (paru == 0) + paru = LM_DWARF / MIN(delta, p1); + +- /*** If the input par lies outside of the interval (parl,paru), +- set par to the closer endpoint. ***/ ++/*** lmpar: 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; + +- /*** Iterate. ***/ ++/*** lmpar: iterate. ***/ + +- for (iter = 0;; iter++) { ++ for (iter=0; ; iter++) { + +- /** Evaluate the function at the current value of par. **/ ++ /** evaluate the function at the current value of par. **/ ++ + if (*par == 0) + *par = MAX(LM_DWARF, 0.001 * paru); + temp = sqrt(*par); + for (j = 0; j < n; j++) + aux[j] = temp * diag[j]; + +- lm_qrsolv(n, r, ldr, Pivot, aux, qtb, x, Sdiag, xdi); +- /* return values are r, x, Sdiag */ ++ lm_qrsolv(n, r, ldr, ipvt, aux, qtb, x, sdiag, xdi); ++ /* return values are r, x, sdiag */ + + for (j = 0; j < n; j++) + xdi[j] = diag[j] * x[j]; /* used as output */ +@@ -711,49 +780,58 @@ + fp_old = fp; + fp = dxnorm - delta; + +- /** If the function is small enough, accept the current value ++ /** 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 <= fp_old && fp_old < 0) || iter == 10) { ++ ++ if (fabs(fp) <= p1 * delta ++ || (parl == 0 && fp <= fp_old && fp_old < 0) ++ || iter == 10) { + #ifdef LMFIT_DEBUG_MESSAGES +- printf("debug lmpar nsing=%d, iter=%d, " +- "par=%.4e [%.4e %.4e], delta=%.4e, fp=%.4e\n", ++ printf("debug lmpar nsing %d iter %d " ++ "par %.4e [%.4e %.4e] delta %.4e fp %.4e\n", + nsing, iter, *par, parl, paru, delta, fp); + #endif + break; /* the only exit from the iteration. */ + } + +- /** Compute the Newton correction. **/ ++ /** compute the Newton correction. **/ ++ + for (j = 0; j < n; j++) +- aux[j] = diag[Pivot[j]] * xdi[Pivot[j]] / dxnorm; ++ aux[j] = diag[ipvt[j]] * xdi[ipvt[j]] / dxnorm; + + for (j = 0; j < n; j++) { +- aux[j] = aux[j] / Sdiag[j]; +- for (i = j+1; i < n; i++) +- aux[i] -= r[j*ldr+i] * aux[j]; ++ aux[j] = aux[j] / sdiag[j]; ++ for (i = j + 1; i < n; i++) ++ aux[i] -= r[j * ldr + i] * aux[j]; + } + temp = lm_enorm(n, aux); + parc = fp / delta / temp / temp; + +- /** Depending on the sign of the function, update parl or paru. **/ ++ /** depending on the sign of the function, update parl or paru. **/ ++ + if (fp > 0) + parl = MAX(parl, *par); +- else /* fp < 0 [the case fp==0 is precluded by the break condition] */ ++ else if (fp < 0) + paru = MIN(paru, *par); ++ /* the case fp==0 is precluded by the break condition */ + +- /** Compute an improved estimate for par. **/ ++ /** compute an improved estimate for par. **/ ++ + *par = MAX(parl, *par + parc); ++ + } + + } /*** lm_lmpar. ***/ + +-/******************************************************************************/ +-/* lm_qrfac (QR factorization, from lapack) */ +-/******************************************************************************/ ++/*****************************************************************************/ ++/* lm_qrfac (QR factorization, from lapack) */ ++/*****************************************************************************/ + +-static void lm_qrfac(const int m, const int n, double* A, int* Pivot, double* Rdiag, +- double* Acnorm, double* W) ++void lm_qrfac( ++ const int m, const int n, double *const A, int *const Pivot, ++ double *const Rdiag, double *const Acnorm, double *const W) ++{ + /* + * This subroutine uses Householder transformations with column pivoting + * to compute a QR factorization of the m by n matrix A. That is, qrfac +@@ -772,27 +850,27 @@ + * + * n is an INPUT parameter 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 vectors w described above). ++ * 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 vectors w described above). + * + * Pivot is an integer OUTPUT array of length n that describes the +- * permutation matrix P. Column j of P is column Pivot(j) of the +- * identity matrix. ++ * permutation matrix P: ++ * Column j of P is column ipvt(j) of the identity matrix. + * +- * Rdiag is an OUTPUT array of length n which contains the diagonal +- * elements of R. ++ * 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 share storage with Rdiag. ++ * 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 share storage with Rdiag. + * + * W is a work array of length n. + * + */ +-{ + int i, j, k, kmax; + double ajnorm, sum, temp; + +@@ -802,16 +880,19 @@ + + /** Compute initial column norms; + initialize Pivot with identity permutation. ***/ ++ + for (j = 0; j < n; j++) { + W[j] = Rdiag[j] = Acnorm[j] = lm_enorm(m, &A[j*m]); + Pivot[j] = j; + } + + /** Loop over columns of A. **/ +- assert(n <= m); ++ ++ assert( n <= m ); + for (j = 0; j < n; j++) { + + /** Bring the column of largest norm into the pivot position. **/ ++ + kmax = j; + for (k = j+1; k < n; k++) + if (Rdiag[k] > Rdiag[kmax]) +@@ -834,6 +915,7 @@ + + /** Compute the Householder reflection vector w_j to reduce the + j-th column of A to a multiple of the j-th unit vector. **/ ++ + ajnorm = lm_enorm(m-j, &A[j*m+j]); + if (ajnorm == 0) { + Rdiag[j] = 0; +@@ -850,6 +932,7 @@ + + /** Apply the Householder transformation U_w := 1 - 2*w_j.w_j/|w_j|^2 + to the remaining columns, and update the norms. **/ ++ + for (k = j+1; k < n; k++) { + /* Compute scalar product w_j * a_j. */ + sum = 0; +@@ -866,12 +949,12 @@ + /* No idea what happens here. */ + if (Rdiag[k] != 0) { + temp = A[m*k+j] / Rdiag[k]; +- if (fabs(temp) < 1) { +- Rdiag[k] *= sqrt(1 - SQR(temp)); ++ if ( fabs(temp)<1 ) { ++ Rdiag[k] *= sqrt(1-SQR(temp)); + temp = Rdiag[k] / W[k]; + } else + temp = 0; +- if (temp == 0 || 0.05 * SQR(temp) <= LM_MACHEP) { ++ if ( temp == 0 || 0.05 * SQR(temp) <= LM_MACHEP ) { + Rdiag[k] = lm_enorm(m-j-1, &A[m*k+j+1]); + W[k] = Rdiag[k]; + } +@@ -882,13 +965,16 @@ + } + } /*** lm_qrfac. ***/ + +-/******************************************************************************/ +-/* lm_qrsolv (linear least-squares) */ +-/******************************************************************************/ + +-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) ++/*****************************************************************************/ ++/* lm_qrsolv (linear least-squares) */ ++/*****************************************************************************/ ++ ++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) ++{ + /* + * 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 +@@ -921,145 +1007,152 @@ + * + * 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. ++ * 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. ++ * ldr is a positive integer INPUT variable not less than n ++ * which specifies the leading dimension of the array R. + * +- * Pivot 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 +- * Pivot(j) of the identity matrix. ++ * 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. ++ * 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^T*b. + * +- * x is an OUTPUT array of length n which contains the least-squares +- * solution of the system A*x = b, D*x = 0. ++ * 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. ++ * sdiag is an OUTPUT array of length n which contains the ++ * diagonal elements of the upper triangular matrix S. + * +- * W is a work array of length n. ++ * wa is a work array of length n. + * + */ +-{ + int i, kk, j, k, nsing; + double qtbpj, sum, temp; + double _sin, _cos, _tan, _cot; /* local variables, not functions */ + +- /*** Copy R and Q^T*b to preserve input and initialize S. +- In particular, save the diagonal elements of R in x. ***/ ++/*** qrsolv: copy R and Q^T*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[j*ldr+i] = r[i*ldr+j]; +- x[j] = r[j*ldr+j]; +- W[j] = qtb[j]; ++ r[j * ldr + i] = r[i * ldr + j]; ++ x[j] = r[j * ldr + j]; ++ wa[j] = qtb[j]; + } + +- /*** Eliminate the diagonal matrix D using a Givens rotation. ***/ ++/*** qrsolv: 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. ***/ ++/*** qrsolv: prepare the row of D to be eliminated, locating the ++ diagonal element using P from the QR factorization. ***/ + +- if (diag[Pivot[j]] != 0) { +- for (k = j; k < n; k++) +- Sdiag[k] = 0; +- Sdiag[j] = diag[Pivot[j]]; ++ if (diag[ipvt[j]] == 0) ++ goto L90; ++ for (k = j; k < n; k++) ++ sdiag[k] = 0; ++ sdiag[j] = diag[ipvt[j]]; + +- /*** The transformations to eliminate the row of D modify only +- a single element of Q^T*b beyond the first n, which is +- initially 0. ***/ ++/*** qrsolv: the transformations to eliminate the row of D modify only ++ a single element of Q^T*b beyond the first n, which is initially 0. ***/ + +- qtbpj = 0; +- for (k = j; k < n; k++) { ++ 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) +- continue; +- kk = k + ldr * k; +- if (fabs(r[kk]) < fabs(Sdiag[k])) { +- _cot = r[kk] / Sdiag[k]; +- _sin = 1 / hypot(1, _cot); +- _cos = _sin * _cot; +- } else { +- _tan = Sdiag[k] / r[kk]; +- _cos = 1 / hypot(1, _tan); +- _sin = _cos * _tan; +- } ++ /** determine a Givens rotation which eliminates the ++ appropriate element in the current row of D. **/ + +- /** Compute the modified diagonal element of R and +- the modified element of (Q^T*b,0). **/ +- r[kk] = _cos * r[kk] + _sin * Sdiag[k]; +- temp = _cos * W[k] + _sin * qtbpj; +- qtbpj = -_sin * W[k] + _cos * qtbpj; +- W[k] = temp; ++ if (sdiag[k] == 0) ++ continue; ++ kk = k + ldr * k; ++ if (fabs(r[kk]) < fabs(sdiag[k])) { ++ _cot = r[kk] / sdiag[k]; ++ _sin = 1 / sqrt(1 + SQR(_cot)); ++ _cos = _sin * _cot; ++ } else { ++ _tan = sdiag[k] / r[kk]; ++ _cos = 1 / sqrt(1 + SQR(_tan)); ++ _sin = _cos * _tan; ++ } + +- /** Accumulate the tranformation in the row of S. **/ +- for (i = k+1; i < n; i++) { +- temp = _cos * r[k*ldr+i] + _sin * Sdiag[i]; +- Sdiag[i] = -_sin * r[k*ldr+i] + _cos * Sdiag[i]; +- r[k*ldr+i] = temp; +- } ++ /** compute the modified diagonal element of R and ++ the modified element of (Q^T*b,0). **/ ++ ++ r[kk] = _cos * r[kk] + _sin * sdiag[k]; ++ temp = _cos * wa[k] + _sin * qtbpj; ++ qtbpj = -_sin * wa[k] + _cos * qtbpj; ++ wa[k] = temp; ++ ++ /** accumulate the tranformation in the row of S. **/ ++ ++ for (i = k + 1; i < n; i++) { ++ temp = _cos * r[k * ldr + i] + _sin * sdiag[i]; ++ sdiag[i] = -_sin * r[k * ldr + i] + _cos * sdiag[i]; ++ r[k * ldr + i] = temp; + } + } + +- /** Store the diagonal element of S and restore ++ L90: ++ /** store the diagonal element of S and restore + the corresponding diagonal element of R. **/ +- Sdiag[j] = r[j*ldr+j]; +- r[j*ldr+j] = x[j]; ++ ++ sdiag[j] = r[j * ldr + j]; ++ r[j * ldr + j] = x[j]; + } + +- /*** Solve the triangular system for z. If the system is singular, then +- obtain a least-squares solution. ***/ ++/*** qrsolv: 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) ++ if (sdiag[j] == 0 && nsing == n) + nsing = j; + if (nsing < n) +- W[j] = 0; ++ wa[j] = 0; + } + +- for (j = nsing-1; j >= 0; j--) { ++ for (j = nsing - 1; j >= 0; j--) { + sum = 0; +- for (i = j+1; i < nsing; i++) +- sum += r[j*ldr+i] * W[i]; +- W[j] = (W[j] - sum) / Sdiag[j]; ++ for (i = j + 1; i < nsing; i++) ++ sum += r[j * ldr + i] * wa[i]; ++ wa[j] = (wa[j] - sum) / sdiag[j]; + } + +- /*** Permute the components of z back to components of x. ***/ ++/*** qrsolv: permute the components of z back to components of x. ***/ + + for (j = 0; j < n; j++) +- x[Pivot[j]] = W[j]; ++ x[ipvt[j]] = wa[j]; + + } /*** lm_qrsolv. ***/ + +-/******************************************************************************/ +-/* lm_enorm (Euclidean norm) */ +-/******************************************************************************/ + +-static double lm_enorm(int n, const double* x) ++/*****************************************************************************/ ++/* lm_enorm (Euclidean norm) */ ++/*****************************************************************************/ ++ ++double lm_enorm(const int n, const double *const x) ++{ + /* This function calculates the Euclidean norm of an n-vector 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 ++ * 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, LM_SQRT_DWARF and LM_SQRT_GIANT. The main +- * restrictions on these constants are that LM_SQRT_DWARF**2 not underflow +- * and LM_SQRT_GIANT**2 not overflow. ++ * restrictions on these constants are that LM_SQRT_DWARF**2 not ++ * underflow and LM_SQRT_GIANT**2 not overflow. + * + * Parameters: + * +@@ -1067,9 +1160,8 @@ + * + * x is an INPUT array of length n. + */ +-{ + int i; +- double agiant, s1, s2, s3, xabs, x1max, x3max; ++ double agiant, s1, s2, s3, xabs, x1max, x3max, temp; + + s1 = 0; + s2 = 0; +@@ -1078,27 +1170,33 @@ + x3max = 0; + agiant = LM_SQRT_GIANT / n; + +- /** Sum squares. **/ ++ /** sum squares. **/ ++ + for (i = 0; i < n; i++) { + xabs = fabs(x[i]); + if (xabs > LM_SQRT_DWARF) { +- if (xabs < agiant) { +- s2 += SQR(xabs); +- } else if (xabs > x1max) { +- s1 = 1 + s1 * SQR(x1max / xabs); ++ if ( xabs < agiant ) { ++ s2 += xabs * xabs; ++ } else if ( xabs > x1max ) { ++ temp = x1max / xabs; ++ s1 = 1 + s1 * SQR(temp); + x1max = xabs; + } else { +- s1 += SQR(xabs / x1max); ++ temp = xabs / x1max; ++ s1 += SQR(temp); + } +- } else if (xabs > x3max) { +- s3 = 1 + s3 * SQR(x3max / xabs); ++ } else if ( xabs > x3max ) { ++ temp = x3max / xabs; ++ s3 = 1 + s3 * SQR(temp); + x3max = xabs; + } else if (xabs != 0) { +- s3 += SQR(xabs / x3max); ++ temp = xabs / x3max; ++ s3 += SQR(temp); + } + } + +- /** Calculate the norm. **/ ++ /** calculation of norm. **/ ++ + if (s1 != 0) + return x1max * sqrt(s1 + (s2 / x1max) / x1max); + else if (s2 != 0) +@@ -1110,3 +1208,80 @@ + return x3max * sqrt(s3); + + } /*** lm_enorm. ***/ ++ ++ ++/*****************************************************************************/ ++/* lm_fnorm (Euclidean norm of difference) */ ++/*****************************************************************************/ ++ ++double lm_fnorm(const int n, const double *const x, const double *const y) ++{ ++/* This function calculates the Euclidean norm of an n-vector x-y. ++ * ++ * 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, LM_SQRT_DWARF and LM_SQRT_GIANT. The main ++ * restrictions on these constants are that LM_SQRT_DWARF**2 not ++ * underflow and LM_SQRT_GIANT**2 not overflow. ++ * ++ * Parameters: ++ * ++ * n is a positive integer INPUT variable. ++ * ++ * x, y are INPUT arrays of length n. ++ */ ++ if (!y) ++ return lm_enorm(n, x); ++ int i; ++ double agiant, s1, s2, s3, xabs, x1max, x3max, temp; ++ ++ s1 = 0; ++ s2 = 0; ++ s3 = 0; ++ x1max = 0; ++ x3max = 0; ++ agiant = LM_SQRT_GIANT / n; ++ ++ /** sum squares. **/ ++ ++ for (i = 0; i < n; i++) { ++ xabs = fabs(x[i]-y[i]); ++ if (xabs > LM_SQRT_DWARF) { ++ if ( xabs < agiant ) { ++ s2 += xabs * xabs; ++ } else if ( xabs > x1max ) { ++ temp = x1max / xabs; ++ s1 = 1 + s1 * SQR(temp); ++ x1max = xabs; ++ } else { ++ temp = xabs / x1max; ++ s1 += SQR(temp); ++ } ++ } else if ( xabs > x3max ) { ++ temp = x3max / xabs; ++ s3 = 1 + s3 * SQR(temp); ++ x3max = xabs; ++ } else if (xabs != 0) { ++ temp = xabs / x3max; ++ s3 += SQR(temp); ++ } ++ } ++ ++ /** calculation of norm. **/ ++ ++ if (s1 != 0) ++ return x1max * sqrt(s1 + (s2 / x1max) / x1max); ++ else if (s2 != 0) ++ if (s2 >= x3max) ++ return sqrt(s2 * (1 + (x3max / s2) * (x3max * s3))); ++ else ++ return sqrt(x3max * ((s2 / x3max) + (x3max * s3))); ++ else ++ return x3max * sqrt(s3); ++ ++} /*** lm_fnorm. ***/ +diff --git a/src/external/lmfit/lmmin.h b/src/external/lmfit/lmmin.h +index 17b4880..d2ccfae 100644 +--- a/src/external/lmfit/lmmin.h ++++ b/src/external/lmfit/lmmin.h +@@ -17,9 +17,15 @@ + + #include "lmstruct.h" + +-/*! \brief +- * Levenberg-Marquardt minimization. +- * ++/* Levenberg-Marquardt minimization. */ ++void lmmin( ++ const int n_par, double* par, const int m_dat, const double* y, ++ const void* data, ++ void (*evaluate)( ++ const double* par, const int m_dat, const void* data, ++ double* fvec, int* userbreak), ++ const lm_control_struct* control, lm_status_struct* status); ++/* + * This routine contains the core algorithm of our library. + * + * It minimizes the sum of the squares of m nonlinear functions +@@ -29,14 +35,16 @@ + * + * Parameters: + * +- * \param[in] n_par The number of variables (INPUT, positive integer). +- * \param par The parameters to be fitted +- * x is the solution vector (INPUT/OUTPUT, array of length n). ++ * n_par is the number of variables (INPUT, positive integer). ++ * ++ * par is the solution vector (INPUT/OUTPUT, array of length n). + * On input it must be set to an estimated solution. + * On output it yields the final estimate of the solution. + * +- * m is the number of functions to be minimized (INPUT, positive integer). ++ * m_dat is the number of functions to be minimized (INPUT, positive integer). + * It must fulfill m>=n. ++ * ++ * y contains data to be fitted. Use a null pointer if there are no data. + * + * data is a pointer that is ignored by lmmin; it is however forwarded + * to the user-supplied functions evaluate and printout. +@@ -56,9 +64,9 @@ + * status contains OUTPUT variables that inform about the fit result, + * as declared and explained in lmstruct.h + */ +-void lmmin( const int n_par, double *par, const int m_dat, const void *data, +- void (*evaluate) (const double *par, const int m_dat, const void *data, +- double *fvec, int *userbreak), +- const lm_control_struct *control, lm_status_struct *status ); ++ ++/* Refined calculation of Eucledian norm. */ ++double lm_enorm(const int, const double*); ++double lm_fnorm(const int, const double*, const double*); + + #endif /* LMMIN_H */ +diff --git a/src/external/lmfit/lmstruct.h b/src/external/lmfit/lmstruct.h +index d545ff1..1f513b3 100644 +--- a/src/external/lmfit/lmstruct.h ++++ b/src/external/lmfit/lmstruct.h +@@ -39,23 +39,23 @@ + stepbound itself. In most cases stepbound should lie + in the interval (0.1,100.0). Generally, the value + 100.0 is recommended. */ +- int patience; /* Used to set the maximum number of function evaluations ++ int patience; /* Used to set the maximum number of function evaluations + to patience*(number_of_parameters+1). */ +- int scale_diag; /* If 1, the variables will be rescaled internally. ++ int scale_diag; /* If 1, the variables will be rescaled internally. + Recommended value is 1. */ + FILE* msgfile; /* Progress messages will be written to this file. */ +- int verbosity; /* OR'ed: 1: print some messages; 2: print Jacobian. */ +- int n_maxpri; /* -1, or max number of parameters to print. */ +- int m_maxpri; /* -1, or max number of residuals to print. */ ++ int verbosity; /* OR'ed: 1: print some messages; 2: print Jacobian. */ ++ int n_maxpri; /* -1, or max number of parameters to print. */ ++ int m_maxpri; /* -1, or max number of residuals to print. */ + } lm_control_struct; + + /* Collection of output parameters for status info. */ + typedef struct { +- double fnorm; /* norm of the residue vector fvec. */ +- int nfev; /* actual number of iterations. */ +- int outcome; /* Status indicator. Nonnegative values are used as index +- for the message text lm_infmsg, set in lmmin.c. */ +- int userbreak; /* Set when function evaluation requests termination. */ ++ double fnorm; /* norm of the residue vector fvec. */ ++ int nfev; /* actual number of iterations. */ ++ int outcome; /* Status indicator. Nonnegative values are used as index ++ for the message text lm_infmsg, set in lmmin.c. */ ++ int userbreak; /* Set when function evaluation requests termination. */ + } lm_status_struct; + + /* Preset (and recommended) control parameter settings. */ +diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt +index 6f1b3e9..7aab572 100644 +--- a/src/gromacs/CMakeLists.txt ++++ b/src/gromacs/CMakeLists.txt +@@ -138,10 +138,6 @@ + tmpi_get_source_list(THREAD_MPI_SOURCES ${CMAKE_SOURCE_DIR}/src/external/thread_mpi/src) + list(APPEND LIBGROMACS_SOURCES ${THREAD_MPI_SOURCES}) + +-get_lmfit_properties(LMFIT_SOURCES LMFIT_LIBRARIES_TO_LINK LMFIT_INCLUDE_DIRECTORY LMFIT_INCLUDE_DIR_ORDER) +-include_directories(${LMFIT_INCLUDE_DIR_ORDER} SYSTEM "${LMFIT_INCLUDE_DIRECTORY}") +-list(APPEND LIBGROMACS_SOURCES ${LMFIT_SOURCES}) +- + configure_file(version.h.cmakein version.h) + gmx_install_headers( + analysisdata.h +@@ -260,6 +256,7 @@ + VERSION ${LIBRARY_VERSION} + COMPILE_FLAGS "${OpenMP_C_FLAGS}") + ++manage_lmfit() + gmx_write_installed_header_list() + + # Only install the library in mdrun-only mode if it is actually necessary +diff --git a/src/gromacs/correlationfunctions/expfit.cpp b/src/gromacs/correlationfunctions/expfit.cpp +index 4fc9fd4..d2cbd4b 100644 +--- a/src/gromacs/correlationfunctions/expfit.cpp ++++ b/src/gromacs/correlationfunctions/expfit.cpp +@@ -3,7 +3,7 @@ + * + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2004, The GROMACS development team. +- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by ++ * Copyright (c) 2013,2014,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,15 +48,10 @@ + + #include + +-#include +- + #include +- +-#include + + #include "gromacs/correlationfunctions/integrate.h" + #include "gromacs/fileio/xvgr.h" +-#include "gromacs/math/functions.h" + #include "gromacs/math/vec.h" + #include "gromacs/utility/fatalerror.h" + #include "gromacs/utility/futil.h" +@@ -394,9 +389,6 @@ + } + } + +-/*! \brief function type for passing to fitting routine */ +-typedef double (*t_lmcurve)(double x, const double *a); +- + /*! \brief array of fitting functions corresponding to the pre-defined types */ + t_lmcurve lmcurves[effnNR+1] = { + lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm, +@@ -414,104 +406,6 @@ + return 0.0; + } + return lmcurves[eFitFn](x, parm); +-} +- +-/*! \brief lmfit_exp supports fitting of different functions +- * +- * This routine calls the Levenberg-Marquardt non-linear fitting +- * routine for fitting a data set with errors to a target function. +- * Fitting routines included in gromacs in src/external/lmfit. +- */ +-static gmx_bool lmfit_exp(int nfit, +- const double x[], +- const double y[], +- const double dy[], +- double parm[], +- gmx_bool bVerbose, +- int eFitFn, +- int nfix) +-{ +- double chisq, ochisq; +- gmx_bool bCont; +- int j; +- int maxiter = 100; +- lm_control_struct control; +- lm_status_struct *status; +- int nparam = effnNparams(eFitFn); +- int p2; +- gmx_bool bSkipLast; +- +- if ((eFitFn < 0) || (eFitFn >= effnNR)) +- { +- fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n", +- eFitFn, effnNR-1); +- return FALSE; +- } +- /* Using default control structure for double precision fitting that +- * comes with the lmfit package (i.e. from the include file). +- */ +- control = lm_control_double; +- control.verbosity = (bVerbose ? 1 : 0); +- control.n_maxpri = 0; +- control.m_maxpri = 0; +- +- snew(status, 1); +- /* Initial params */ +- chisq = 1e12; +- j = 0; +- if (bVerbose) +- { +- printf("%4s %10s Parameters\n", "Step", "chi^2"); +- } +- /* Check whether we have to skip some params */ +- if (nfix > 0) +- { +- do +- { +- p2 = 1 << (nparam-1); +- bSkipLast = ((p2 & nfix) == p2); +- if (bSkipLast) +- { +- nparam--; +- nfix -= p2; +- } +- } +- while ((nparam > 0) && (bSkipLast)); +- if (bVerbose) +- { +- printf("Using %d out of %d parameters\n", nparam, effnNparams(eFitFn)); +- } +- } +- do +- { +- ochisq = chisq; +- gmx_lmcurve(nparam, parm, nfit, x, y, dy, +- lmcurves[eFitFn], &control, status); +- chisq = gmx::square(status->fnorm); +- if (bVerbose) +- { +- printf("status: fnorm = %g, nfev = %d, userbreak = %d\noutcome = %s\n", +- status->fnorm, status->nfev, status->userbreak, +- lm_infmsg[status->outcome]); +- } +- if (bVerbose) +- { +- int mmm; +- printf("%4d %8g", j, chisq); +- for (mmm = 0; (mmm < effnNparams(eFitFn)); mmm++) +- { +- printf(" %8g", parm[mmm]); +- } +- printf("\n"); +- } +- j++; +- bCont = (fabs(ochisq - chisq) > fabs(control.ftol*chisq)); +- } +- while (bCont && (j < maxiter)); +- +- sfree(status); +- +- return TRUE; + } + + /*! \brief Ensure the fitting parameters are well-behaved. +diff --git a/src/gromacs/correlationfunctions/gmx_lmcurve.cpp b/src/gromacs/correlationfunctions/gmx_lmcurve.cpp +index f3a103c..8b3870e 100644 +--- a/src/gromacs/correlationfunctions/gmx_lmcurve.cpp ++++ b/src/gromacs/correlationfunctions/gmx_lmcurve.cpp +@@ -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. +@@ -38,14 +38,30 @@ + * Defines a driver routine for lmfit, and a callback for it to use. + * + * \author David van der Spoel ++ * \author Mark Abraham + * \ingroup module_correlationfunctions + */ + #include "gmxpre.h" + + #include "gmx_lmcurve.h" + ++#include "config.h" ++ ++#include ++ ++#if HAVE_LMFIT + #include + #include ++#endif ++ ++#include "gromacs/correlationfunctions/expfit.h" ++#include "gromacs/math/functions.h" ++#include "gromacs/utility/fatalerror.h" ++#include "gromacs/utility/smalloc.h" ++ ++#if HAVE_LMFIT ++ ++extern t_lmcurve lmcurves[effnNR+1]; + + typedef struct { + const double* t; +@@ -73,7 +89,8 @@ + *info = 0; + } + +-void gmx_lmcurve( ++//! Calls lmmin with the given data, with callback function \c f. ++static void gmx_lmcurve( + const int n_par, double* par, const int m_dat, + const double* t, const double* y, const double *dy, + double (*f)(double t, const double* par), +@@ -81,6 +98,114 @@ + { + lmcurve_data_struct data = { t, y, dy, f }; + +- lmmin(n_par, par, m_dat, (const void*)&data, lmcurve_evaluate, ++ lmmin(n_par, par, m_dat, nullptr, (const void*)&data, lmcurve_evaluate, + control, status); + } ++ ++#endif ++ ++bool lmfit_exp(int nfit, ++ const double x[], ++ const double y[], ++ const double dy[], ++ double parm[], ++ bool bVerbose, ++ int eFitFn, ++ int nfix) ++{ ++ if ((eFitFn < 0) || (eFitFn >= effnNR)) ++ { ++ fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n", ++ eFitFn, effnNR-1); ++ return false; ++ } ++#if HAVE_LMFIT ++ double chisq, ochisq; ++ gmx_bool bCont; ++ int j; ++ int maxiter = 100; ++ lm_control_struct control; ++ lm_status_struct *status; ++ int nparam = effnNparams(eFitFn); ++ int p2; ++ gmx_bool bSkipLast; ++ ++ /* Using default control structure for double precision fitting that ++ * comes with the lmfit package (i.e. from the include file). ++ */ ++ control = lm_control_double; ++ control.verbosity = (bVerbose ? 1 : 0); ++ control.n_maxpri = 0; ++ control.m_maxpri = 0; ++ ++ snew(status, 1); ++ /* Initial params */ ++ chisq = 1e12; ++ j = 0; ++ if (bVerbose) ++ { ++ printf("%4s %10s Parameters\n", "Step", "chi^2"); ++ } ++ /* Check whether we have to skip some params */ ++ if (nfix > 0) ++ { ++ do ++ { ++ p2 = 1 << (nparam-1); ++ bSkipLast = ((p2 & nfix) == p2); ++ if (bSkipLast) ++ { ++ nparam--; ++ nfix -= p2; ++ } ++ } ++ while ((nparam > 0) && (bSkipLast)); ++ if (bVerbose) ++ { ++ printf("Using %d out of %d parameters\n", nparam, effnNparams(eFitFn)); ++ } ++ } ++ do ++ { ++ ochisq = chisq; ++ gmx_lmcurve(nparam, parm, nfit, x, y, dy, ++ lmcurves[eFitFn], &control, status); ++ chisq = gmx::square(status->fnorm); ++ if (bVerbose) ++ { ++ printf("status: fnorm = %g, nfev = %d, userbreak = %d\noutcome = %s\n", ++ status->fnorm, status->nfev, status->userbreak, ++ lm_infmsg[status->outcome]); ++ } ++ if (bVerbose) ++ { ++ int mmm; ++ printf("%4d %8g", j, chisq); ++ for (mmm = 0; (mmm < effnNparams(eFitFn)); mmm++) ++ { ++ printf(" %8g", parm[mmm]); ++ } ++ printf("\n"); ++ } ++ j++; ++ bCont = (fabs(ochisq - chisq) > fabs(control.ftol*chisq)); ++ } ++ while (bCont && (j < maxiter)); ++ ++ sfree(status); ++#else ++ gmx_fatal(FARGS, "This build of GROMACS was not configured with support " ++ "for lmfit, so the requested fitting cannot be performed. " ++ "See the install guide for instructions on how to build " ++ "GROMACS with lmfit supported."); ++ GMX_UNUSED_VALUE(nfit); ++ GMX_UNUSED_VALUE(x); ++ GMX_UNUSED_VALUE(y); ++ GMX_UNUSED_VALUE(dy); ++ GMX_UNUSED_VALUE(parm); ++ GMX_UNUSED_VALUE(bVerbose); ++ GMX_UNUSED_VALUE(eFitFn); ++ GMX_UNUSED_VALUE(nfix); ++#endif ++ return true; ++} +diff --git a/src/gromacs/correlationfunctions/gmx_lmcurve.h b/src/gromacs/correlationfunctions/gmx_lmcurve.h +index 0f6c957..fa77671 100644 +--- a/src/gromacs/correlationfunctions/gmx_lmcurve.h ++++ b/src/gromacs/correlationfunctions/gmx_lmcurve.h +@@ -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. +@@ -43,13 +43,22 @@ + #ifndef GMX_CORRELATION_FUNCTIONS_GMX_LMCURVE_H + #define GMX_CORRELATION_FUNCTIONS_GMX_LMCURVE_H + +-#include ++/*! \brief function type for passing to fitting routine */ ++typedef double (*t_lmcurve)(double x, const double *a); + +-//! Calls lmmin with the given data, with callback function \c f. +-void gmx_lmcurve( const int n_par, double *par, const int m_dat, +- const double *t, const double *y, const double *dy, +- double (*f)(double t, const double *par ), +- const lm_control_struct *control, +- lm_status_struct *status ); ++/*! \brief lmfit_exp supports fitting of different functions ++ * ++ * This routine calls the Levenberg-Marquardt non-linear fitting ++ * routine for fitting a data set with errors to a target function. ++ * Fitting routines included in gromacs in src/external/lmfit. ++ */ ++bool lmfit_exp(int nfit, ++ const double x[], ++ const double y[], ++ const double dy[], ++ double parm[], ++ bool bVerbose, ++ int eFitFn, ++ int nfix); + + #endif +diff --git a/src/gromacs/correlationfunctions/tests/expfit.cpp b/src/gromacs/correlationfunctions/tests/expfit.cpp +index eb93463..5d5bd96 100644 +--- a/src/gromacs/correlationfunctions/tests/expfit.cpp ++++ b/src/gromacs/correlationfunctions/tests/expfit.cpp +@@ -1,7 +1,7 @@ + /* + * This file is part of the GROMACS molecular simulation package. + * +- * Copyright (c) 2014,2015,2017, by the GROMACS development team, led by ++ * Copyright (c) 2014,2015,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. +@@ -43,6 +43,8 @@ + #include "gmxpre.h" + + #include "gromacs/correlationfunctions/expfit.h" ++ ++#include "config.h" + + #include + +@@ -152,6 +154,10 @@ + //static var + std::vector ExpfitTest::data_; + ++// TODO calling test() leads to a fatal error, which we could in ++// principle test for. ++#if HAVE_LMFIT ++ + TEST_F (ExpfitTest, EffnEXP1) { + double param[] = {25}; + test(effnEXP1, param, 1e-5, 0); +@@ -202,6 +208,8 @@ + test(effnPRES, param, 1e-4, 1); + } + ++#endif ++ + } + + } diff --git a/gromacs.spec b/gromacs.spec index 1d70cd9..eff14cf 100644 --- a/gromacs.spec +++ b/gromacs.spec @@ -36,8 +36,8 @@ %endif Name: gromacs -Version: 2018.1 -Release: 1%{?_rcname}%{?dist}.1 +Version: 2018.2 +Release: 1%{?_rcname}%{?dist} Summary: Fast, Free and Flexible Molecular Dynamics License: GPLv2+ URL: http://www.gromacs.org @@ -55,7 +55,8 @@ BuildRequires: python2-sphinx %else Source0: ftp://ftp.gromacs.org/pub/gromacs/gromacs-%{version}%{?_rc}.tar.gz Source1: ftp://ftp.gromacs.org/pub/manual/manual-%{version}%{?_rc}.pdf -Source2: http://gerrit.gromacs.org/download/regressiontests-%{version}%{?_rc}.tar.gz +# Too britle sind 2018.2 +# Source2: http://gerrit.gromacs.org/download/regressiontests-{version}{?_rc}.tar.gz %endif Source6: gromacs-README.fedora # fix path to packaged dssp @@ -66,9 +67,17 @@ Patch0: gromacs-dssp-path.patch Patch2: gromacs-issue-2366.patch # fix building documentation Patch3: gromacs-sphinx-no-man.patch +# add support for lmfit-7.0 +# https://redmine.gromacs.org/issues/2533 +Patch4: facb927.diff BuildRequires: gcc-c++ -BuildRequires: cmake -BuildRequires: atlas-devel >= 3.10.1 +%if 0%{?rhel} +BuildRequires: cmake3 >= 3.4.3 +%else +BuildRequires: cmake >= 3.4.3 +%global cmake3 %{cmake} +%endif +BuildRequires: openblas-devel BuildRequires: fftw-devel BuildRequires: gsl-devel BuildRequires: hwloc @@ -165,7 +174,6 @@ This package the manual in PDF format. %package devel Summary: GROMACS header files and development libraries Requires: gromacs-libs = %{version}-%{release} -Requires: cmake Obsoletes: gromacs-mpich-devel < 2016-0.1.20160318gitbec9c87 Obsoletes: gromacs-openmpi-devel < 2016-0.1.20160318gitbec9c87 @@ -241,8 +249,9 @@ This package single and double precision binaries and libraries. %setup -q -n gromacs-%{commit} %patch3 -p1 -b .sphinx-no-man %else -%setup -q -a 2 -n gromacs-%{version}%{?_rc} +%setup -q %{?SOURCE2:-a 2} -n gromacs-%{version}%{?_rc} %patch2 -p1 +%patch4 -p1 install -Dpm644 %{SOURCE1} ./serial/docs/manual/gromacs.pdf %endif %patch0 -p1 -b .dssp @@ -261,12 +270,13 @@ export LDFLAGS="-L%{_libdir}/atlas" -DBUILD_TESTING:BOOL=ON \\\ -DCMAKE_SKIP_RPATH:BOOL=ON \\\ -DCMAKE_SKIP_BUILD_RPATH:BOOL=ON \\\ - -DGMX_BLAS_USER=satlas \\\ + -DGMX_BLAS_USER=openblas \\\ -DGMX_BUILD_UNITTESTS:BOOL=ON \\\ -DGMX_EXTERNAL_LMFIT:BOOL=ON \\\ + -DGMX_USE_LMFIT=external \\\ -DGMX_EXTERNAL_TNG:BOOL=ON \\\ -DGMX_EXTERNAL_TINYXML2:BOOL=OFF \\\ - -DGMX_LAPACK_USER=satlas \\\ + -DGMX_LAPACK_USER=openblas \\\ -DGMX_USE_RDTSCP=OFF \\\ -DGMX_SIMD=%{simd} \\\ @@ -280,7 +290,7 @@ export LDFLAGS="-L%{_libdir}/atlas" %{_openmpi_load} for p in '' _d ; do cd openmpi${p} -%cmake \ +%{cmake3} \ %{defopts} \ %{mpi} \ -DGMX_BINARY_SUFFIX=${MPI_SUFFIX}${p} -DGMX_LIBS_SUFFIX=${MPI_SUFFIX}${p} \ @@ -294,7 +304,7 @@ done %{_mpich_load} for p in '' _d ; do cd mpich${p} -%cmake \ +%{cmake3} \ %{defopts} \ %{mpi} \ -DGMX_BINARY_SUFFIX=${MPI_SUFFIX}${p} -DGMX_LIBS_SUFFIX=${MPI_SUFFIX}${p} \ @@ -307,10 +317,9 @@ done for p in '' _d ; do cd serial${p} -cp -al ../regressiontests* tests/ -%cmake \ +# cp -al ../regressiontests* tests/ # use with -DREGRESSIONTEST_PATH=${PWD}/tests below +%{cmake3} \ %{defopts} \ - -DREGRESSIONTEST_PATH=${PWD}/tests \ -DGMX_X11=ON \ $(test -n "$p" && echo %{double} || echo %{?single}) \ .. @@ -320,7 +329,7 @@ done %if %{git} cd serial -%cmake \ +%{cmake3} \ %{defopts} \ -DGMX_X11=ON \ %{?single} \ @@ -373,13 +382,8 @@ rm ./%{_bindir}/gmx-completion-${bin}.bash done rm ./%{_bindir}/gmx-completion.bash +%ldconfig_scriptlets libs -# Post install for libs. MPI packages don't need this. -%post libs -p /sbin/ldconfig - -%postun libs -p /sbin/ldconfig - -%if 1 %check %{_openmpi_load} for p in '' _d ; do @@ -405,7 +409,6 @@ for p in '' ; do LD_LIBRARY_PATH=%{buildroot}%{_libdir} make VERBOSE=1 %{?_smp_mflags} check cd .. done -%endif %files @@ -447,6 +450,12 @@ done %{_libdir}/mpich/bin/mdrun_mpich* %changelog +* Wed Jul 18 2018 Christoph Junghans - 2018.2-1 +- Version bump to 2018.2 (bug #1591052) +- Add support for lmfit-7 (patch will be part of v2019) +- Switch to OpenBlas (bug #1602822) +- Disable brittle regressiontests + * Fri Jul 13 2018 Fedora Release Engineering - 2018.1-1.1 - Rebuilt for https://fedoraproject.org/wiki/Fedora_29_Mass_Rebuild diff --git a/sources b/sources index 6c3ea67..73baa45 100644 --- a/sources +++ b/sources @@ -1,3 +1,2 @@ -SHA512 (gromacs-2018.1.tar.gz) = d29f152e9f115c7de07881c6af4cc05481e0a5520bd33142a09507e8c4df9f8b6c9d6d96efcb7adaaf4e0127b76ab247421d43918bd7e7779d6e71f5984db715 -SHA512 (manual-2018.1.pdf) = 486e2b14dc35040df355d8cb7c2f76290214ca0e93021d4cb6ad30969d94b2105b3cd7c545ba13656f4950e222e86d8f69307af0c43ce63858fc8fdf8acfa3d4 -SHA512 (regressiontests-2018.1.tar.gz) = 81c9a62fca859923c1e27214b32b0cff1dd48224dd4ad9301554036b842ccc400a2729752ba71b284e0c5b6c1769ce7de5fe2c9ba4fc7cf0917fd4ced9883112 +SHA512 (gromacs-2018.2.tar.gz) = d444b503e24a9875b0ab7622772946ef73ab2c897da6ff45ac908f147ea398ba2404b064a8784996fd34b25e188e36f12a492e0070427e0929f422d934205d28 +SHA512 (manual-2018.2.pdf) = 0afd1fbb00be7d887015843c90ea63617c3796a1a035e30b2c10d0b45e43ded12648c531ff7a20185d890496aafa56b83d5e723033a48ad16fe34eee97f6dd25