Skip to content
GitLab
Commit aa52c8bf authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Added support for the INTEL MKL library in place of ATLAS and FFTW (configure --enable-mkl option). 

Updated LevMar library to V2.6.
Pushed version number to 2.18.0.
parent 1b190b55
Hide whitespace changes
Inline Side-by-side
src/levmar/Makefile.in
View file @ aa52c8bf
# Makefile.in generated by automake 1.10.1 from Makefile.am.
# Makefile.in generated by automake 1.11.1 from Makefile.am.
# @configure_input@
# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
# 2003, 2004, 2005, 2006, 2007, 2008 Free Software Foundation, Inc.
# 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation,
# Inc.
# This Makefile.in is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.
......@@ -45,8 +46,9 @@
VPATH = @srcdir@
pkgdatadir = $(datadir)/@PACKAGE@
pkglibdir = $(libdir)/@PACKAGE@
pkgincludedir = $(includedir)/@PACKAGE@
pkglibdir = $(libdir)/@PACKAGE@
pkglibexecdir = $(libexecdir)/@PACKAGE@
am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd
install_sh_DATA = $(install_sh) -c -m 644
install_sh_PROGRAM = $(install_sh) -c
......@@ -65,7 +67,8 @@ subdir = src/levmar
DIST_COMMON = README $(srcdir)/Makefile.am $(srcdir)/Makefile.in
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/acx_atlas.m4 \
$(top_srcdir)/acx_fftw.m4 $(top_srcdir)/acx_prog_cc_optim.m4 \
$(top_srcdir)/acx_fftw.m4 $(top_srcdir)/acx_mkl.m4 \
$(top_srcdir)/acx_prog_cc_optim.m4 \
$(top_srcdir)/acx_pthread.m4 \
$(top_srcdir)/acx_urbi_resolve_dir.m4 \
$(top_srcdir)/configure.ac
......@@ -74,6 +77,7 @@ am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \
mkinstalldirs = $(SHELL) $(top_srcdir)/autoconf/mkinstalldirs
CONFIG_HEADER = $(top_builddir)/config.h
CONFIG_CLEAN_FILES =
CONFIG_CLEAN_VPATH_FILES =
LIBRARIES = $(noinst_LIBRARIES)
ARFLAGS = cru
liblevmar_a_AR = $(AR) $(ARFLAGS)
......@@ -85,6 +89,7 @@ liblevmar_a_OBJECTS = $(am_liblevmar_a_OBJECTS)
DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir)
depcomp = $(SHELL) $(top_srcdir)/autoconf/depcomp
am__depfiles_maybe = depfiles
am__mv = mv -f
COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS)
LTCOMPILE = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
......@@ -107,8 +112,8 @@ AM_LDFLAGS = @AM_LDFLAGS@
AR = @AR@
ATLAS_CFLAGS = @ATLAS_CFLAGS@
ATLAS_ERROR = @ATLAS_ERROR@
ATLAS_LIB = @ATLAS_LIB@
ATLAS_LIBPATH = @ATLAS_LIBPATH@
ATLAS_LIBS = @ATLAS_LIBS@
ATLAS_WARN = @ATLAS_WARN@
AUTOCONF = @AUTOCONF@
AUTOHEADER = @AUTOHEADER@
......@@ -119,47 +124,56 @@ CCDEPMODE = @CCDEPMODE@
CFLAGS = @CFLAGS@
CPP = @CPP@
CPPFLAGS = @CPPFLAGS@
CXX = @CXX@
CXXCPP = @CXXCPP@
CXXDEPMODE = @CXXDEPMODE@
CXXFLAGS = @CXXFLAGS@
CYGPATH_W = @CYGPATH_W@
DATE2 = @DATE2@
DATE3 = @DATE3@
DEFS = @DEFS@
DEPDIR = @DEPDIR@
ECHO = @ECHO@
DLLTOOL = @DLLTOOL@
DSYMUTIL = @DSYMUTIL@
DUMPBIN = @DUMPBIN@
ECHO_C = @ECHO_C@
ECHO_N = @ECHO_N@
ECHO_T = @ECHO_T@
EGREP = @EGREP@
EXEEXT = @EXEEXT@
F77 = @F77@
FFLAGS = @FFLAGS@
FFTW_ERROR = @FFTW_ERROR@
FFTW_LIBS = @FFTW_LIBS@
FFTW_WARN = @FFTW_WARN@
FGREP = @FGREP@
GREP = @GREP@
INSTALL = @INSTALL@
INSTALL_DATA = @INSTALL_DATA@
INSTALL_PROGRAM = @INSTALL_PROGRAM@
INSTALL_SCRIPT = @INSTALL_SCRIPT@
INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@
LD = @LD@
LDFLAGS = @LDFLAGS@
LIBOBJS = @LIBOBJS@
LIBS = @LIBS@
LIBTOOL = @LIBTOOL@
LIPO = @LIPO@
LN_S = @LN_S@
LTLIBOBJS = @LTLIBOBJS@
MAKEINFO = @MAKEINFO@
MANIFEST_TOOL = @MANIFEST_TOOL@
MKDIR_P = @MKDIR_P@
MKL_CFLAGS = @MKL_CFLAGS@
MKL_LIBS = @MKL_LIBS@
MKL_WARN = @MKL_WARN@
NM = @NM@
NMEDIT = @NMEDIT@
OBJDUMP = @OBJDUMP@
OBJEXT = @OBJEXT@
OTOOL = @OTOOL@
OTOOL64 = @OTOOL64@
PACKAGE = @PACKAGE@
PACKAGER = @PACKAGER@
PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@
PACKAGE_NAME = @PACKAGE_NAME@
PACKAGE_STRING = @PACKAGE_STRING@
PACKAGE_TARNAME = @PACKAGE_TARNAME@
PACKAGE_URL = @PACKAGE_URL@
PACKAGE_VERSION = @PACKAGE_VERSION@
PATH_SEPARATOR = @PATH_SEPARATOR@
PTHREAD_CC = @PTHREAD_CC@
......@@ -175,9 +189,9 @@ abs_builddir = @abs_builddir@
abs_srcdir = @abs_srcdir@
abs_top_builddir = @abs_top_builddir@
abs_top_srcdir = @abs_top_srcdir@
ac_ct_AR = @ac_ct_AR@
ac_ct_CC = @ac_ct_CC@
ac_ct_CXX = @ac_ct_CXX@
ac_ct_F77 = @ac_ct_F77@
ac_ct_DUMPBIN = @ac_ct_DUMPBIN@
am__include = @am__include@
am__leading_dot = @am__leading_dot@
am__quote = @am__quote@
......@@ -231,8 +245,7 @@ EXTRA_liblevmar_a_SOURCES = Axb_core.c lmbc_core.c lm_core.c \
lmblec_core.c lmbleic_core.c lmlec_core.c \
misc_core.c \
LICENSE README README.txt \
Makefile.icc Makefile.vc levmar.vcproj lmdemo.vcproj \
expfit.c lmdemo.c lm.h matlab
Makefile.icc Makefile.vc expfit.c lmdemo.c lm.h matlab
all: all-am
......@@ -242,14 +255,14 @@ $(srcdir)/Makefile.in: $(srcdir)/Makefile.am $(am__configure_deps)
@for dep in $?; do \
case '$(am__configure_deps)' in \
*$$dep*) \
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh \
&& exit 0; \
( cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh ) \
&& { if test -f $@; then exit 0; else break; fi; }; \
exit 1;; \
esac; \
done; \
echo ' cd $(top_srcdir) && $(AUTOMAKE) --foreign src/levmar/Makefile'; \
cd $(top_srcdir) && \
$(AUTOMAKE) --foreign src/levmar/Makefile
echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu src/levmar/Makefile'; \
$(am__cd) $(top_srcdir) && \
$(AUTOMAKE) --gnu src/levmar/Makefile
.PRECIOUS: Makefile
Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status
@case '$?' in \
......@@ -267,6 +280,7 @@ $(top_srcdir)/configure: $(am__configure_deps)
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh
$(ACLOCAL_M4): $(am__aclocal_m4_deps)
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh
$(am__aclocal_m4_deps):
clean-noinstLIBRARIES:
-test -z "$(noinst_LIBRARIES)" || rm -f $(noinst_LIBRARIES)
......@@ -300,21 +314,21 @@ distclean-compile:
.c.o:
@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@am__fastdepCC_TRUE@ $(am__mv) $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(COMPILE) -c $<
.c.obj:
@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ `$(CYGPATH_W) '$<'`
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@am__fastdepCC_TRUE@ $(am__mv) $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(COMPILE) -c `$(CYGPATH_W) '$<'`
.c.lo:
@am__fastdepCC_TRUE@ $(LTCOMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Plo
@am__fastdepCC_TRUE@ $(am__mv) $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Plo
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(LTCOMPILE) -c -o $@ $<
......@@ -330,14 +344,14 @@ ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES)
unique=`for i in $$list; do \
if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \
done | \
$(AWK) '{ files[$$0] = 1; nonemtpy = 1; } \
$(AWK) '{ files[$$0] = 1; nonempty = 1; } \
END { if (nonempty) { for (i in files) print i; }; }'`; \
mkid -fID $$unique
tags: TAGS
TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \
$(TAGS_FILES) $(LISP)
tags=; \
set x; \
here=`pwd`; \
list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
unique=`for i in $$list; do \
......@@ -345,29 +359,34 @@ TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \
done | \
$(AWK) '{ files[$$0] = 1; nonempty = 1; } \
END { if (nonempty) { for (i in files) print i; }; }'`; \
if test -z "$(ETAGS_ARGS)$$tags$$unique"; then :; else \
shift; \
if test -z "$(ETAGS_ARGS)$$*$$unique"; then :; else \
test -n "$$unique" || unique=$$empty_fix; \
$(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \
$$tags $$unique; \
if test $$# -gt 0; then \
$(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \
"$$@" $$unique; \
else \
$(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \
$$unique; \
fi; \
fi
ctags: CTAGS
CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \
$(TAGS_FILES) $(LISP)
tags=; \
list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
unique=`for i in $$list; do \
if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \
done | \
$(AWK) '{ files[$$0] = 1; nonempty = 1; } \
END { if (nonempty) { for (i in files) print i; }; }'`; \
test -z "$(CTAGS_ARGS)$$tags$$unique" \
test -z "$(CTAGS_ARGS)$$unique" \
|| $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \
$$tags $$unique
$$unique
GTAGS:
here=`$(am__cd) $(top_builddir) && pwd` \
&& cd $(top_srcdir) \
&& gtags -i $(GTAGS_ARGS) $$here
&& $(am__cd) $(top_srcdir) \
&& gtags -i $(GTAGS_ARGS) "$$here"
distclean-tags:
-rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags
......@@ -388,13 +407,17 @@ distdir: $(DISTFILES)
if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \
if test -d $$d/$$file; then \
dir=`echo "/$$file" | sed -e 's,/[^/]*$$,,'`; \
if test -d "$(distdir)/$$file"; then \
find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \
fi; \
if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \
cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \
cp -fpR $(srcdir)/$$file "$(distdir)$$dir" || exit 1; \
find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \
fi; \
cp -pR $$d/$$file $(distdir)$$dir || exit 1; \
cp -fpR $$d/$$file "$(distdir)$$dir" || exit 1; \
else \
test -f $(distdir)/$$file \
|| cp -p $$d/$$file $(distdir)/$$file \
test -f "$(distdir)/$$file" \
|| cp -p $$d/$$file "$(distdir)/$$file" \
|| exit 1; \
fi; \
done
......@@ -422,6 +445,7 @@ clean-generic:
distclean-generic:
-test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES)
-test . = "$(srcdir)" || test -z "$(CONFIG_CLEAN_VPATH_FILES)" || rm -f $(CONFIG_CLEAN_VPATH_FILES)
maintainer-clean-generic:
@echo "This command is intended for maintainers to use"
......@@ -443,6 +467,8 @@ dvi-am:
html: html-am
html-am:
info: info-am
info-am:
......@@ -451,18 +477,28 @@ install-data-am:
install-dvi: install-dvi-am
install-dvi-am:
install-exec-am:
install-html: install-html-am
install-html-am:
install-info: install-info-am
install-info-am:
install-man:
install-pdf: install-pdf-am
install-pdf-am:
install-ps: install-ps-am
install-ps-am:
installcheck-am:
maintainer-clean: maintainer-clean-am
......@@ -500,6 +536,7 @@ uninstall-am:
mostlyclean-compile mostlyclean-generic mostlyclean-libtool \
pdf pdf-am ps ps-am tags uninstall uninstall-am
# Tell versions [3.59,3.63) of GNU make to not export all variables.
# Otherwise a system limit (for SysV at least) may be exceeded.
.NOEXPORT:
src/levmar/Makefile.vc
View file @ aa52c8bf
......@@ -16,9 +16,9 @@ CC=cl /nologo
CONFIGFLAGS=#/ULINSOLVERS_RETAIN_MEMORY
# YOU MIGHT WANT TO UNCOMMENT THE FOLLOWING LINE
#SPOPTFLAGS=/GL /G7 /arch:SSE2 # special optimization: resp. whole program opt., Athlon/Pentium4 opt., SSE2 extensions
# /MD COMPILES WITH MULTIPLE THREADS SUPPORT. TO DISABLE IT, SUBSTITUTE WITH /ML
# /MD COMPILES WITH DYNAMIC LINK SUPPORT. FOR STATIC LINK, SUBSTITUTE WITH /MT
# FLAG /EHsc SUPERSEDED /GX IN MSVC'05. IF YOU HAVE AN EARLIER VERSION THAT COMPLAINS ABOUT IT, CHANGE /EHsc TO /GX
CFLAGS=$(CONFIGFLAGS) /I. /MD /W3 /EHsc /O2 $(SPOPTFLAGS) # /Wall
CFLAGS=$(CONFIGFLAGS) /I. /MD /TC /W3 /EHsc /D_CRT_SECURE_NO_WARNINGS /O2 $(SPOPTFLAGS) # /Wall
LAPACKLIBS_PATH=C:\src\lib # WHEN USING LAPACK, CHANGE THIS TO WHERE YOUR COMPILED LIBS ARE!
LDFLAGS=/link /subsystem:console /opt:ref /libpath:$(LAPACKLIBS_PATH) /libpath:.
LIBOBJS=lm.obj Axb.obj misc.obj lmlec.obj lmbc.obj lmblec.obj lmbleic.obj
......@@ -38,7 +38,7 @@ levmar.lib: $(LIBOBJS)
$(AR) /out:levmar.lib $(LIBOBJS)
lmdemo.exe: $(DEMOBJS) levmar.lib
$(CC) $(DEMOBJS) $(LDFLAGS) /out:lmdemo.exe $(LIBS)
$(CC) $(DEMOBJS) $(LDFLAGS) /out:lmdemo.exe /incremental:no $(LIBS)
lm.obj: lm.c lm_core.c levmar.h misc.h compiler.h
Axb.obj: Axb.c Axb_core.c levmar.h misc.h
......
src/levmar/README
View file @ aa52c8bf
The levmar v2.5 library has been included in this package untouched, except for
The levmar v2.6 library has been included in this package untouched, except for
three warnings removed, LU decomposition replaced with calls to ATLAS-Lapack
routines, Hessian matrix inversion done with SVD, and an AutoMakefile added.
Emmanuel Bertin <bertin@iap.fr>
src/levmar/README.txt
View file @ aa52c8bf
**************************************************************
LEVMAR
version 2.5
version 2.6
By Manolis Lourakis
Institute of Computer Science
......@@ -42,28 +42,29 @@ in using levmar in a proprietary commercial application, a commercial license fo
can be obtained by contacting the author using the email address at the end of this file.
COMPILATION
- You might first consider setting a few configuration options at the top of
levmar.h. See the accompanying comments for more details.
- On a Linux/Unix system, typing "make" will build both levmar and the demo
program using gcc. Alternatively, if Intel's C++ compiler is installed, it
can be used by typing "make -f Makefile.icc".
- Under Windows and if Visual C is installed & configured for command line
use, type "nmake /f Makefile.vc" in a cmd window to build levmar and the
demo program. In case of trouble, read the comments on top of Makefile.vc
Visual C++ project files (levmar.vcproj and lmdemo.vcproj) are also included,
however they are not supported and are only meant to serve as a starting point
for creating your own. Check http://www.arstdesign.com/articles/prjconverter.html
if you need to convert to .dsw/.dsp (i.e., Visual C++ 6.0) project files.
- levmar can also be built under various platforms using the CMake cross-platform
build system. The included CMakeLists.txt file can be used to generate makefiles
for Unix systems or project files for Windows systems. See http://www.cmake.org
for details.
- The preferred way to build levmar is through the CMake cross-platform build
system. The included CMakeLists.txt file can be used to generate makefiles
for Unix systems or project files for Windows systems. CMakeLists.txt defines
some configuration variables that control certain aspects of levmar and can
be modified from CMake's user interface. The values of these variables are
automatically propagated to levmar.h after CMake runs.
More information on how to use CMake can be found at http://www.cmake.org
- levmar can also be built using the supplied makefiles. Platform-specific
instructions are given below. Before compiling, you might consider setting
a few configuration options found at the top of levmar.h. See the
accompanying comments for more details.
-- On a Linux/Unix system, typing "make" will build both levmar and the demo
program using gcc. Alternatively, if Intel's C++ compiler is installed, it
can be used by typing "make -f Makefile.icc".
-- Under Windows and if Visual C is installed & configured for command line
use, type "nmake /f Makefile.vc" in a cmd window to build levmar and the
demo program. In case of trouble, read the comments on top of Makefile.vc
MATLAB INTERFACE
Since version 2.2, the levmar distribution includes a matlab interface.
Since version 2.2, the levmar distribution includes a matlab mex interface.
See the 'matlab' subdirectory for more information and examples of use.
Notice that *_core.c files are not to be compiled directly; For example,
......@@ -71,4 +72,4 @@ Axb_core.c is included by Axb.c, to provide single and double precision
routine versions.
Send your comments/bug reports to lourakis at ics forth gr
Send your comments/bug reports to lourakis (at) ics (dot) forth (dot) gr
src/levmar/compiler.h
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......@@ -63,10 +63,14 @@
#define LM_FINITE finite // ICC, GCC
#else
#define LM_FINITE finite // other than MSVC, ICC, GCC, let's hope this will work
#endif
#endif
#ifdef _MSC_VER // avoid deprecation warnings in VS2005
#define _CRT_SECURE_NO_WARNINGS
#ifdef _MSC_VER
#define LM_ISINF(x) (!_finite(x) && !_isnan(x)) // MSVC
#elif defined(__ICC) || defined(__INTEL_COMPILER) || defined(__GNUC__)
#define LM_ISINF(x) isinf(x) // ICC, GCC
#else
#define LM_ISINF(x) isinf(x) // other than MSVC, ICC, GCC, let's hope this will work
#endif
#endif /* _COMPILER_H_ */
src/levmar/levmar.h
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/*
......@@ -49,27 +49,48 @@
#ifndef _LEVMAR_H_
#define _LEVMAR_H_
/* Added by EB */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
/* End added by EB */
/************************************* Start of configuration options *************************************/
/* Note that when compiling with CMake, this configuration section is automatically generated
* based on the user's input, see levmar.h.in
*/
/* specify whether to use LAPACK or not. The first option is strongly recommended */
/* specifies whether to use LAPACK or not. Using LAPACK is strongly recommended */
/*#define HAVE_LAPACK*/ /* use LAPACK */
#undef HAVE_LAPACK /* uncomment this to force not using LAPACK */
/* specifies whether the PLASMA parallel library for multicore CPUs is available */
/* #undef HAVE_PLASMA */
/* to avoid the overhead of repeated mallocs(), routines in Axb.c can be instructed to
* retain working memory between calls. Such a choice, however, renders these routines
* non-reentrant and is not safe in a shared memory multiprocessing environment.
* Bellow, this option is turned on only when not compiling with OpenMP.
*/
#if !defined(_OPENMP)
#define LINSOLVERS_RETAIN_MEMORY /* comment this if you don't want routines in Axb.c retain working memory between calls */
#endif
/* determine the precision variants to be build. Default settings build
* both the single and double precision routines
* Bellow, an attempt is made to issue a warning if this option is turned on and OpenMP
* is being used (note that this will work only if omp.h is included before levmar.h)
* NOTE FROM EB: TURNED OFF IN ALL CASES
*/
#define LM_DBL_PREC /* comment this if you don't want the double precision routines to be compiled */
#define LM_SNGL_PREC /* comment this if you don't want the single precision routines to be compiled */
/*
#define LINSOLVERS_RETAIN_MEMORY
*/
#if (defined(_OPENMP))
# ifdef LINSOLVERS_RETAIN_MEMORY
# ifdef _MSC_VER
# pragma message("LINSOLVERS_RETAIN_MEMORY is not safe in a multithreaded environment and should be turned off!")
# else
# warning LINSOLVERS_RETAIN_MEMORY is not safe in a multithreaded environment and should be turned off!
# endif /* _MSC_VER */
# endif /* LINSOLVERS_RETAIN_MEMORY */
#endif /* _OPENMP */
/* specifies whether double precision routines will be compiled or not */
#define LM_DBL_PREC
/* specifies whether single precision routines will be compiled or not */
#define LM_SNGL_PREC
/****************** End of configuration options, no changes necessary beyond this point ******************/
......@@ -78,9 +99,6 @@
extern "C" {
#endif
#define FABS(x) (((x)>=0.0)? (x) : -(x))
/* work arrays size for ?levmar_der and ?levmar_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
......@@ -117,7 +135,7 @@ extern "C" {
#define LM_INIT_MU 1E-03
#define LM_STOP_THRESH 1E-17
#define LM_DIFF_DELTA 1E-06
#define LM_VERSION "2.5 (December 2009)"
#define LM_VERSION "2.6 (November 2011)"
#ifdef LM_DBL_PREC
/* double precision LM, with & without Jacobian */
......@@ -137,12 +155,12 @@ extern int dlevmar_dif(
extern int dlevmar_bc_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub,
double *p, double *x, int m, int n, double *lb, double *ub, double *dscl,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
extern int dlevmar_bc_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub,
double *p, double *x, int m, int n, double *lb, double *ub, double *dscl,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
#ifdef HAVE_LAPACK
......@@ -242,12 +260,12 @@ extern int slevmar_dif(
extern int slevmar_bc_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub,
float *p, float *x, int m, int n, float *lb, float *ub, float *dscl,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
extern int slevmar_bc_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub,
float *p, float *x, int m, int n, float *lb, float *ub, float *dscl,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
#ifdef HAVE_LAPACK
......@@ -361,6 +379,16 @@ extern int sAx_eq_b_LU_noLapack(float *A, float *B, float *x, int n);
#endif /* HAVE_LAPACK */
#ifdef HAVE_PLASMA
#ifdef LM_DBL_PREC
extern int dAx_eq_b_PLASMA_Chol(double *A, double *B, double *x, int m);
#endif
#ifdef LM_SNGL_PREC
extern int sAx_eq_b_PLASMA_Chol(float *A, float *B, float *x, int m);
#endif
extern void levmar_PLASMA_setnbcores(int cores);
#endif /* HAVE_PLASMA */
/* Jacobian verification, double & single precision */
#ifdef LM_DBL_PREC
extern void dlevmar_chkjac(
......@@ -376,7 +404,9 @@ extern void slevmar_chkjac(
float *p, int m, int n, void *adata, float *err);
#endif /* LM_SNGL_PREC */
/* standard deviation, coefficient of determination (R2) & Pearson's correlation coefficient for best-fit parameters */
/* miscellaneous: standard deviation, coefficient of determination (R2),
* Pearson's correlation coefficient for best-fit parameters
*/
#ifdef LM_DBL_PREC
extern double dlevmar_stddev( double *covar, int m, int i);
extern double dlevmar_corcoef(double *covar, int m, int i, int j);
......@@ -388,6 +418,14 @@ extern double dlevmar_R2(void (*func)(double *p, double *hx, int m, int n, void
extern float slevmar_stddev( float *covar, int m, int i);
extern float slevmar_corcoef(float *covar, int m, int i, int j);
extern float slevmar_R2(void (*func)(float *p, float *hx, int m, int n, void *adata), float *p, float *x, int m, int n, void *adata);
extern void slevmar_locscale(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, void *adata,
int howto, float locscl[2], float **residptr);
extern int slevmar_outlid(float *r, int n, float thresh, float ls[2], char *outlmap);
#endif /* LM_SNGL_PREC */
#ifdef __cplusplus
......
src/levmar/levmar.h.in 0 → 100644
View file @ aa52c8bf
/*
////////////////////////////////////////////////////////////////////////////////////
//
// Prototypes and definitions for the Levenberg - Marquardt minimization algorithm
// Copyright (C) 2004 Manolis Lourakis (lourakis at ics forth gr)
// Institute of Computer Science, Foundation for Research & Technology - Hellas
// Heraklion, Crete, Greece.
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
////////////////////////////////////////////////////////////////////////////////////
*/
#ifndef _LEVMAR_H_
#define _LEVMAR_H_
/************************************* Start of configuration options *************************************/
/* Note that when compiling with CMake, this configuration section is automatically generated
* based on the user's input, see levmar.h.in
*/
/* specifies whether to use LAPACK or not. Using LAPACK is strongly recommended */
#cmakedefine HAVE_LAPACK
/* specifies whether the PLASMA parallel library for multicore CPUs is available */
#cmakedefine HAVE_PLASMA
/* to avoid the overhead of repeated mallocs(), routines in Axb.c can be instructed to
* retain working memory between calls. Such a choice, however, renders these routines
* non-reentrant and is not safe in a shared memory multiprocessing environment.
* Bellow, an attempt is made to issue a warning if this option is turned on and OpenMP
* is being used (note that this will work only if omp.h is included before levmar.h)
*/
#cmakedefine LINSOLVERS_RETAIN_MEMORY
#if (defined(_OPENMP))
# ifdef LINSOLVERS_RETAIN_MEMORY
# ifdef _MSC_VER
# pragma message("LINSOLVERS_RETAIN_MEMORY is not safe in a multithreaded environment and should be turned off!")
# else
# warning LINSOLVERS_RETAIN_MEMORY is not safe in a multithreaded environment and should be turned off!
# endif /* _MSC_VER */
# endif /* LINSOLVERS_RETAIN_MEMORY */
#endif /* _OPENMP */
/* specifies whether double precision routines will be compiled or not */
#cmakedefine LM_DBL_PREC
/* specifies whether single precision routines will be compiled or not */
#cmakedefine LM_SNGL_PREC
/****************** End of configuration options, no changes necessary beyond this point ******************/
#ifdef __cplusplus
extern "C" {
#endif
/* work arrays size for ?levmar_der and ?levmar_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_DER_WORKSZ(npar, nmeas) (2*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar))
#define LM_DIF_WORKSZ(npar, nmeas) (4*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar))
/* work arrays size for ?levmar_bc_der and ?levmar_bc_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_BC_DER_WORKSZ(npar, nmeas) (2*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar))
#define LM_BC_DIF_WORKSZ(npar, nmeas) LM_BC_DER_WORKSZ((npar), (nmeas)) /* LEVMAR_BC_DIF currently implemented using LEVMAR_BC_DER()! */
/* work arrays size for ?levmar_lec_der and ?levmar_lec_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_LEC_DER_WORKSZ(npar, nmeas, nconstr) LM_DER_WORKSZ((npar)-(nconstr), (nmeas))
#define LM_LEC_DIF_WORKSZ(npar, nmeas, nconstr) LM_DIF_WORKSZ((npar)-(nconstr), (nmeas))
/* work arrays size for ?levmar_blec_der and ?levmar_blec_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_BLEC_DER_WORKSZ(npar, nmeas, nconstr) LM_LEC_DER_WORKSZ((npar), (nmeas)+(npar), (nconstr))
#define LM_BLEC_DIF_WORKSZ(npar, nmeas, nconstr) LM_LEC_DIF_WORKSZ((npar), (nmeas)+(npar), (nconstr))
/* work arrays size for ?levmar_bleic_der and ?levmar_bleic_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_BLEIC_DER_WORKSZ(npar, nmeas, nconstr1, nconstr2) LM_BLEC_DER_WORKSZ((npar)+(nconstr2), (nmeas)+(nconstr2), (nconstr1)+(nconstr2))
#define LM_BLEIC_DIF_WORKSZ(npar, nmeas, nconstr1, nconstr2) LM_BLEC_DIF_WORKSZ((npar)+(nconstr2), (nmeas)+(nconstr2), (nconstr1)+(nconstr2))
#define LM_OPTS_SZ 5 /* max(4, 5) */
#define LM_INFO_SZ 10
#define LM_ERROR -1
#define LM_INIT_MU 1E-03
#define LM_STOP_THRESH 1E-17
#define LM_DIFF_DELTA 1E-06
#define LM_VERSION "2.6 (November 2011)"
#ifdef LM_DBL_PREC
/* double precision LM, with & without Jacobian */
/* unconstrained minimization */
extern int dlevmar_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, int itmax, double *opts,
double *info, double *work, double *covar, void *adata);
extern int dlevmar_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, int itmax, double *opts,
double *info, double *work, double *covar, void *adata);
/* box-constrained minimization */
extern int dlevmar_bc_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *dscl,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
extern int dlevmar_bc_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *dscl,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
#ifdef HAVE_LAPACK
/* linear equation constrained minimization */
extern int dlevmar_lec_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *A, double *b, int k,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
extern int dlevmar_lec_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *A, double *b, int k,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
/* box & linear equation constrained minimization */
extern int dlevmar_blec_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *A, double *b, int k, double *wghts,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
extern int dlevmar_blec_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *A, double *b, int k, double *wghts,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
/* box, linear equations & inequalities constrained minimization */
extern int dlevmar_bleic_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub,
double *A, double *b, int k1, double *C, double *d, int k2,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
extern int dlevmar_bleic_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub,
double *A, double *b, int k1, double *C, double *d, int k2,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
/* box & linear inequality constraints */
extern int dlevmar_blic_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *C, double *d, int k2,
int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
extern int dlevmar_blic_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *C, double *d, int k2,
int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
/* linear equation & inequality constraints */
extern int dlevmar_leic_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *A, double *b, int k1, double *C, double *d, int k2,
int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
extern int dlevmar_leic_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *A, double *b, int k1, double *C, double *d, int k2,
int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
/* linear inequality constraints */
extern int dlevmar_lic_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *C, double *d, int k2,
int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
extern int dlevmar_lic_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *C, double *d, int k2,
int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
#endif /* HAVE_LAPACK */
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
/* single precision LM, with & without Jacobian */
/* unconstrained minimization */
extern int slevmar_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, int itmax, float *opts,
float *info, float *work, float *covar, void *adata);
extern int slevmar_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, int itmax, float *opts,
float *info, float *work, float *covar, void *adata);
/* box-constrained minimization */
extern int slevmar_bc_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *dscl,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
extern int slevmar_bc_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *dscl,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
#ifdef HAVE_LAPACK
/* linear equation constrained minimization */
extern int slevmar_lec_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *A, float *b, int k,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
extern int slevmar_lec_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *A, float *b, int k,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
/* box & linear equation constrained minimization */
extern int slevmar_blec_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *A, float *b, int k, float *wghts,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
extern int slevmar_blec_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *A, float *b, int k, float *wghts,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
/* box, linear equations & inequalities constrained minimization */
extern int slevmar_bleic_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub,
float *A, float *b, int k1, float *C, float *d, int k2,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
extern int slevmar_bleic_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub,
float *A, float *b, int k1, float *C, float *d, int k2,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
/* box & linear inequality constraints */
extern int slevmar_blic_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *C, float *d, int k2,
int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
extern int slevmar_blic_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *C, float *d, int k2,
int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
/* linear equality & inequality constraints */
extern int slevmar_leic_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *A, float *b, int k1, float *C, float *d, int k2,
int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
extern int slevmar_leic_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *A, float *b, int k1, float *C, float *d, int k2,
int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
/* linear inequality constraints */
extern int slevmar_lic_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *C, float *d, int k2,
int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
extern int slevmar_lic_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *C, float *d, int k2,
int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
#endif /* HAVE_LAPACK */
#endif /* LM_SNGL_PREC */
/* linear system solvers */
#ifdef HAVE_LAPACK
#ifdef LM_DBL_PREC
extern int dAx_eq_b_QR(double *A, double *B, double *x, int m);
extern int dAx_eq_b_QRLS(double *A, double *B, double *x, int m, int n);
extern int dAx_eq_b_Chol(double *A, double *B, double *x, int m);
extern int dAx_eq_b_LU(double *A, double *B, double *x, int m);
extern int dAx_eq_b_SVD(double *A, double *B, double *x, int m);
extern int dAx_eq_b_BK(double *A, double *B, double *x, int m);
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
extern int sAx_eq_b_QR(float *A, float *B, float *x, int m);
extern int sAx_eq_b_QRLS(float *A, float *B, float *x, int m, int n);
extern int sAx_eq_b_Chol(float *A, float *B, float *x, int m);
extern int sAx_eq_b_LU(float *A, float *B, float *x, int m);
extern int sAx_eq_b_SVD(float *A, float *B, float *x, int m);
extern int sAx_eq_b_BK(float *A, float *B, float *x, int m);
#endif /* LM_SNGL_PREC */
#else /* no LAPACK */
#ifdef LM_DBL_PREC
extern int dAx_eq_b_LU_noLapack(double *A, double *B, double *x, int n);
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
extern int sAx_eq_b_LU_noLapack(float *A, float *B, float *x, int n);
#endif /* LM_SNGL_PREC */
#endif /* HAVE_LAPACK */
#ifdef HAVE_PLASMA
#ifdef LM_DBL_PREC
extern int dAx_eq_b_PLASMA_Chol(double *A, double *B, double *x, int m);
#endif
#ifdef LM_SNGL_PREC
extern int sAx_eq_b_PLASMA_Chol(float *A, float *B, float *x, int m);
#endif
extern void levmar_PLASMA_setnbcores(int cores);
#endif /* HAVE_PLASMA */
/* Jacobian verification, double & single precision */
#ifdef LM_DBL_PREC
extern void dlevmar_chkjac(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, int m, int n, void *adata, double *err);
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
extern void slevmar_chkjac(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, int m, int n, void *adata, float *err);
#endif /* LM_SNGL_PREC */
/* miscellaneous: standard deviation, coefficient of determination (R2),
* Pearson's correlation coefficient for best-fit parameters
*/
#ifdef LM_DBL_PREC
extern double dlevmar_stddev( double *covar, int m, int i);
extern double dlevmar_corcoef(double *covar, int m, int i, int j);
extern double dlevmar_R2(void (*func)(double *p, double *hx, int m, int n, void *adata), double *p, double *x, int m, int n, void *adata);
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
extern float slevmar_stddev( float *covar, int m, int i);
extern float slevmar_corcoef(float *covar, int m, int i, int j);
extern float slevmar_R2(void (*func)(float *p, float *hx, int m, int n, void *adata), float *p, float *x, int m, int n, void *adata);
extern void slevmar_locscale(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, void *adata,
int howto, float locscl[2], float **residptr);
extern int slevmar_outlid(float *r, int n, float thresh, float ls[2], char *outlmap);
#endif /* LM_SNGL_PREC */
#ifdef __cplusplus
}
#endif
#endif /* _LEVMAR_H_ */
src/levmar/levmar.vcproj deleted 100644 → 0
View file @ 1b190b55
<?xml version="1.0" encoding="UTF-8"?>
<VisualStudioProject
ProjectType="Visual C++"
Version="8,00"
Name="levmar"
ProjectGUID="{F329E490-DB04-453A-A0BF-FEB90BD949D8}"
Keyword="Win32Proj"
>
<Platforms>
<Platform
Name="Win32"
/>
</Platforms>
<ToolFiles>
</ToolFiles>
<Configurations>
<Configuration
Name="Debug|Win32"
OutputDirectory="Debug"
IntermediateDirectory="Debug"
ConfigurationType="4"
>
<Tool
Name="VCPreBuildEventTool"
/>
<Tool
Name="VCCustomBuildTool"
/>
<Tool
Name="VCXMLDataGeneratorTool"
/>
<Tool
Name="VCWebServiceProxyGeneratorTool"
/>
<Tool
Name="VCMIDLTool"
/>
<Tool
Name="VCCLCompilerTool"
Optimization="0"
PreprocessorDefinitions="WIN32;_DEBUG;_CONSOLE"
MinimalRebuild="true"
BasicRuntimeChecks="3"
RuntimeLibrary="3"
UsePrecompiledHeader="0"
WarningLevel="3"
Detect64BitPortabilityProblems="true"
DebugInformationFormat="4"
CompileAs="1"
/>
<Tool
Name="VCManagedResourceCompilerTool"
/>
<Tool
Name="VCResourceCompilerTool"
/>
<Tool
Name="VCPreLinkEventTool"
/>
<Tool
Name="VCLibrarianTool"
/>
<Tool
Name="VCALinkTool"
/>
<Tool
Name="VCXDCMakeTool"
/>
<Tool
Name="VCBscMakeTool"
/>
<Tool
Name="VCFxCopTool"
/>
<Tool
Name="VCPostBuildEventTool"
/>
</Configuration>
<Configuration
Name="Release|Win32"
OutputDirectory="Release"
IntermediateDirectory="Release"
ConfigurationType="4"
>
<Tool
Name="VCPreBuildEventTool"
/>
<Tool
Name="VCCustomBuildTool"
/>
<Tool
Name="VCXMLDataGeneratorTool"
/>
<Tool
Name="VCWebServiceProxyGeneratorTool"
/>
<Tool
Name="VCMIDLTool"
/>
<Tool
Name="VCCLCompilerTool"
PreprocessorDefinitions="WIN32;NDEBUG;_CONSOLE"
RuntimeLibrary="2"
UsePrecompiledHeader="0"
WarningLevel="3"
Detect64BitPortabilityProblems="true"
DebugInformationFormat="3"
CompileAs="1"
/>
<Tool
Name="VCManagedResourceCompilerTool"
/>
<Tool
Name="VCResourceCompilerTool"
/>
<Tool
Name="VCPreLinkEventTool"
/>
<Tool
Name="VCLibrarianTool"
/>
<Tool
Name="VCALinkTool"
/>
<Tool
Name="VCXDCMakeTool"
/>
<Tool
Name="VCBscMakeTool"
/>
<Tool
Name="VCFxCopTool"
/>
<Tool
Name="VCPostBuildEventTool"
/>
</Configuration>
</Configurations>
<References>
</References>
<Files>
<Filter
Name="Header Files"
Filter="h;hpp;hxx;hm;inl;inc;xsd"
>
<File
RelativePath=".\compiler.h"
>
</File>
<File
RelativePath=".\levmar.h"
>
</File>
<File
RelativePath=".\misc.h"
>
</File>
</Filter>
<Filter
Name="Resource Files"
Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx"
>
</Filter>
<Filter
Name="Source Files"
Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"
>
<File
RelativePath=".\Axb.c"
>
</File>
<File
RelativePath=".\lm.c"
>
</File>
<File
RelativePath=".\lmbc.c"
>
</File>
<File
RelativePath=".\lmblec.c"
>
</File>
<File
RelativePath=".\lmlec.c"
>
</File>
<File
RelativePath=".\misc.c"
>
</File>
<File
RelativePath=".\lmbleic.c"
>
</File>
</Filter>
</Files>
<Globals>
</Globals>
</VisualStudioProject>
src/levmar/lm.c
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......
src/levmar/lm.h
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,376 +22,17 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/*
////////////////////////////////////////////////////////////////////////////////////
//
// Prototypes and definitions for the Levenberg - Marquardt minimization algorithm
// Copyright (C) 2004 Manolis Lourakis (lourakis at ics forth gr)
// Institute of Computer Science, Foundation for Research & Technology - Hellas
// Heraklion, Crete, Greece.
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
////////////////////////////////////////////////////////////////////////////////////
*/
#ifndef _DEPR_LM_H_
#define _DEPR_LM_H_
#ifndef _LEVMAR_H_
#define _LEVMAR_H_
#ifdef _MSC_VER
#pragma message("lm.h is deprecated, please use levmar.h instead!")
#else
#error lm.h is deprecated, please use levmar.h instead!
#endif /* _MSC_VER */
#endif /* _DEPR_LM_H_ */
/************************************* Start of configuration options *************************************/
/* specify whether to use LAPACK or not. The first option is strongly recommended */
/*#define HAVE_LAPACK*/ /* use LAPACK */
#undef HAVE_LAPACK /* uncomment this to force not using LAPACK */
/* to avoid the overhead of repeated mallocs(), routines in Axb.c can be instructed to
* retain working memory between calls. Such a choice, however, renders these routines
* non-reentrant and is not safe in a shared memory multiprocessing environment.
* Bellow, this option is turned on only when not compiling with OpenMP.
*/
#if !defined(_OPENMP)
#define LINSOLVERS_RETAIN_MEMORY /* comment this if you don't want routines in Axb.c retain working memory between calls */
#endif
/* determine the precision variants to be build. Default settings build
* both the single and double precision routines
*/
#define LM_DBL_PREC /* comment this if you don't want the double precision routines to be compiled */
#define LM_SNGL_PREC /* comment this if you don't want the single precision routines to be compiled */
/****************** End of configuration options, no changes necessary beyond this point ******************/
#ifdef __cplusplus
extern "C" {
#endif
#define FABS(x) (((x)>=0.0)? (x) : -(x))
/* work arrays size for ?levmar_der and ?levmar_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_DER_WORKSZ(npar, nmeas) (2*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar))
#define LM_DIF_WORKSZ(npar, nmeas) (4*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar))
/* work arrays size for ?levmar_bc_der and ?levmar_bc_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_BC_DER_WORKSZ(npar, nmeas) (2*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar))
#define LM_BC_DIF_WORKSZ(npar, nmeas) LM_BC_DER_WORKSZ((npar), (nmeas)) /* LEVMAR_BC_DIF currently implemented using LEVMAR_BC_DER()! */
/* work arrays size for ?levmar_lec_der and ?levmar_lec_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_LEC_DER_WORKSZ(npar, nmeas, nconstr) LM_DER_WORKSZ((npar)-(nconstr), (nmeas))
#define LM_LEC_DIF_WORKSZ(npar, nmeas, nconstr) LM_DIF_WORKSZ((npar)-(nconstr), (nmeas))
/* work arrays size for ?levmar_blec_der and ?levmar_blec_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_BLEC_DER_WORKSZ(npar, nmeas, nconstr) LM_LEC_DER_WORKSZ((npar), (nmeas)+(npar), (nconstr))
#define LM_BLEC_DIF_WORKSZ(npar, nmeas, nconstr) LM_LEC_DIF_WORKSZ((npar), (nmeas)+(npar), (nconstr))
/* work arrays size for ?levmar_bleic_der and ?levmar_bleic_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_BLEIC_DER_WORKSZ(npar, nmeas, nconstr1, nconstr2) LM_BLEC_DER_WORKSZ((npar)+(nconstr2), (nmeas)+(nconstr2), (nconstr1)+(nconstr2))
#define LM_BLEIC_DIF_WORKSZ(npar, nmeas, nconstr1, nconstr2) LM_BLEC_DIF_WORKSZ((npar)+(nconstr2), (nmeas)+(nconstr2), (nconstr1)+(nconstr2))
#define LM_OPTS_SZ 5 /* max(4, 5) */
#define LM_INFO_SZ 10
#define LM_ERROR -1
#define LM_INIT_MU 1E-03
#define LM_STOP_THRESH 1E-17
#define LM_DIFF_DELTA 1E-06
#define LM_VERSION "2.5 (December 2009)"
#ifdef LM_DBL_PREC
/* double precision LM, with & without Jacobian */
/* unconstrained minimization */
extern int dlevmar_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, int itmax, double *opts,
double *info, double *work, double *covar, void *adata);
extern int dlevmar_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, int itmax, double *opts,
double *info, double *work, double *covar, void *adata);
/* box-constrained minimization */
extern int dlevmar_bc_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
extern int dlevmar_bc_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
#ifdef HAVE_LAPACK
/* linear equation constrained minimization */
extern int dlevmar_lec_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *A, double *b, int k,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
extern int dlevmar_lec_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *A, double *b, int k,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
/* box & linear equation constrained minimization */
extern int dlevmar_blec_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *A, double *b, int k, double *wghts,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
extern int dlevmar_blec_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *A, double *b, int k, double *wghts,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
/* box, linear equations & inequalities constrained minimization */
extern int dlevmar_bleic_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub,
double *A, double *b, int k1, double *C, double *d, int k2,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
extern int dlevmar_bleic_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub,
double *A, double *b, int k1, double *C, double *d, int k2,
int itmax, double *opts, double *info, double *work, double *covar, void *adata);
/* box & linear inequality constraints */
extern int dlevmar_blic_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *C, double *d, int k2,
int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
extern int dlevmar_blic_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *lb, double *ub, double *C, double *d, int k2,
int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
/* linear equation & inequality constraints */
extern int dlevmar_leic_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *A, double *b, int k1, double *C, double *d, int k2,
int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
extern int dlevmar_leic_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *A, double *b, int k1, double *C, double *d, int k2,
int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
/* linear inequality constraints */
extern int dlevmar_lic_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, double *C, double *d, int k2,
int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
extern int dlevmar_lic_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, double *C, double *d, int k2,
int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata);
#endif /* HAVE_LAPACK */
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
/* single precision LM, with & without Jacobian */
/* unconstrained minimization */
extern int slevmar_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, int itmax, float *opts,
float *info, float *work, float *covar, void *adata);
extern int slevmar_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, int itmax, float *opts,
float *info, float *work, float *covar, void *adata);
/* box-constrained minimization */
extern int slevmar_bc_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
extern int slevmar_bc_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
#ifdef HAVE_LAPACK
/* linear equation constrained minimization */
extern int slevmar_lec_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *A, float *b, int k,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
extern int slevmar_lec_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *A, float *b, int k,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
/* box & linear equation constrained minimization */
extern int slevmar_blec_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *A, float *b, int k, float *wghts,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
extern int slevmar_blec_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *A, float *b, int k, float *wghts,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
/* box, linear equations & inequalities constrained minimization */
extern int slevmar_bleic_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub,
float *A, float *b, int k1, float *C, float *d, int k2,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
extern int slevmar_bleic_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub,
float *A, float *b, int k1, float *C, float *d, int k2,
int itmax, float *opts, float *info, float *work, float *covar, void *adata);
/* box & linear inequality constraints */
extern int slevmar_blic_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *C, float *d, int k2,
int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
extern int slevmar_blic_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *lb, float *ub, float *C, float *d, int k2,
int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
/* linear equality & inequality constraints */
extern int slevmar_leic_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *A, float *b, int k1, float *C, float *d, int k2,
int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
extern int slevmar_leic_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *A, float *b, int k1, float *C, float *d, int k2,
int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
/* linear inequality constraints */
extern int slevmar_lic_der(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, float *x, int m, int n, float *C, float *d, int k2,
int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
extern int slevmar_lic_dif(
void (*func)(float *p, float *hx, int m, int n, void *adata),
float *p, float *x, int m, int n, float *C, float *d, int k2,
int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata);
#endif /* HAVE_LAPACK */
#endif /* LM_SNGL_PREC */
/* linear system solvers */
#ifdef HAVE_LAPACK
#ifdef LM_DBL_PREC
extern int dAx_eq_b_QR(double *A, double *B, double *x, int m);
extern int dAx_eq_b_QRLS(double *A, double *B, double *x, int m, int n);
extern int dAx_eq_b_Chol(double *A, double *B, double *x, int m);
extern int dAx_eq_b_LU(double *A, double *B, double *x, int m);
extern int dAx_eq_b_SVD(double *A, double *B, double *x, int m);
extern int dAx_eq_b_BK(double *A, double *B, double *x, int m);
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
extern int sAx_eq_b_QR(float *A, float *B, float *x, int m);
extern int sAx_eq_b_QRLS(float *A, float *B, float *x, int m, int n);
extern int sAx_eq_b_Chol(float *A, float *B, float *x, int m);
extern int sAx_eq_b_LU(float *A, float *B, float *x, int m);
extern int sAx_eq_b_SVD(float *A, float *B, float *x, int m);
extern int sAx_eq_b_BK(float *A, float *B, float *x, int m);
#endif /* LM_SNGL_PREC */
#else /* no LAPACK */
#ifdef LM_DBL_PREC
extern int dAx_eq_b_LU_noLapack(double *A, double *B, double *x, int n);
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
extern int sAx_eq_b_LU_noLapack(float *A, float *B, float *x, int n);
#endif /* LM_SNGL_PREC */
#endif /* HAVE_LAPACK */
/* Jacobian verification, double & single precision */
#ifdef LM_DBL_PREC
extern void dlevmar_chkjac(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, int m, int n, void *adata, double *err);
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
extern void slevmar_chkjac(
void (*func)(float *p, float *hx, int m, int n, void *adata),
void (*jacf)(float *p, float *j, int m, int n, void *adata),
float *p, int m, int n, void *adata, float *err);
#endif /* LM_SNGL_PREC */
/* standard deviation, coefficient of determination (R2) & Pearson's correlation coefficient for best-fit parameters */
#ifdef LM_DBL_PREC
extern double dlevmar_stddev( double *covar, int m, int i);
extern double dlevmar_corcoef(double *covar, int m, int i, int j);
extern double dlevmar_R2(void (*func)(double *p, double *hx, int m, int n, void *adata), double *p, double *x, int m, int n, void *adata);
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
extern float slevmar_stddev( float *covar, int m, int i);
extern float slevmar_corcoef(float *covar, int m, int i, int j);
extern float slevmar_R2(void (*func)(float *p, float *hx, int m, int n, void *adata), float *p, float *x, int m, int n, void *adata);
#endif /* LM_SNGL_PREC */
#ifdef __cplusplus
}
#endif
#endif /* _LEVMAR_H_ */
src/levmar/lm_core.c
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......@@ -69,6 +69,10 @@
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack)
#endif /* HAVE_LAPACK */
#ifdef HAVE_PLASMA
#define AX_EQ_B_PLASMA_CHOL LM_ADD_PREFIX(Ax_eq_b_PLASMA_Chol)
#endif
/*
* This function seeks the parameter vector p that best describes the measurements vector x.
* More precisely, given a vector function func : R^m --> R^n with n>=m,
......@@ -238,8 +242,8 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
*/
/* looping downwards saves a few computations */
register int l, im;
register LM_REAL alpha, *jaclm;
register int l;
register LM_REAL alpha, *jaclm, *jacTjacim;
for(i=m*m; i-->0; )
jacTjac[i]=0.0;
......@@ -249,10 +253,10 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
for(l=n; l-->0; ){
jaclm=jac+l*m;
for(i=m; i-->0; ){
im=i*m;
jacTjacim=jacTjac+i*m;
alpha=jaclm[i]; //jac[l*m+i];
for(j=i+1; j-->0; ) /* j<=i computes lower triangular part only */
jacTjac[im+j]+=jaclm[j]*alpha; //jac[l*m+j]
jacTjacim[j]+=jaclm[j]*alpha; //jacTjac[i*m+j]+=jac[l*m+j]*alpha
/* J^T e */
jacTe[i]+=alpha*e[l];
......@@ -321,14 +325,19 @@ if(!(k%100)){
/* solve augmented equations */
#ifdef HAVE_LAPACK
/* 6 alternatives are available: LU, Cholesky, 2 variants of QR decomposition, SVD and LDLt.
* Cholesky is the fastest but might be inaccurate; QR is slower but more accurate;
* SVD is the slowest but most accurate; LU offers a tradeoff between accuracy and speed
/* 7 alternatives are available: LU, Cholesky + Cholesky with PLASMA, LDLt, 2 variants of QR decomposition and SVD.
* For matrices with dimensions of at least a few hundreds, the PLASMA implementation of Cholesky is the fastest.
* From the serial solvers, Cholesky is the fastest but might occasionally be inapplicable due to numerical round-off;
* QR is slower but more robust; SVD is the slowest but most robust; LU is quite robust but
* slower than LDLt; LDLt offers a good tradeoff between robustness and speed
*/
//issolved=AX_EQ_B_BK(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_BK;
issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU;
issolved=AX_EQ_B_BK(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_BK;
//issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU;
//issolved=AX_EQ_B_CHOL(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_CHOL;
#ifdef HAVE_PLASMA
//issolved=AX_EQ_B_PLASMA_CHOL(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_PLASMA_CHOL;
#endif
//issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR;
//issolved=AX_EQ_B_QRLS(jacTjac, jacTe, Dp, m, m); ++nlss; linsolver=(int (*)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m))AX_EQ_B_QRLS;
//issolved=AX_EQ_B_SVD(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_SVD;
......@@ -545,7 +554,7 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
if(!work){
worksz=LM_DIF_WORKSZ(m, n); //4*n+4*m + n*m + m*m;
work=(LM_REAL *)calloc(1,worksz*sizeof(LM_REAL)); /* allocate a big chunk in one step */
work=(LM_REAL *)malloc(worksz*sizeof(LM_REAL)); /* allocate a big chunk in one step */
if(!work){
fprintf(stderr, LCAT(LEVMAR_DIF, "(): memory allocation request failed\n"));
return LM_ERROR;
......@@ -628,8 +637,8 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
* Note that the non-blocking algorithm is faster on small
* problems since in this case it avoids the overheads of blocking.
*/
register int l, im;
register LM_REAL alpha, *jaclm;
register int l;
register LM_REAL alpha, *jaclm, *jacTjacim;
/* looping downwards saves a few computations */
for(i=m*m; i-->0; )
......@@ -640,10 +649,10 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
for(l=n; l-->0; ){
jaclm=jac+l*m;
for(i=m; i-->0; ){
im=i*m;
jacTjacim=jacTjac+i*m;
alpha=jaclm[i]; //jac[l*m+i];
for(j=i+1; j-->0; ) /* j<=i computes lower triangular part only */
jacTjac[im+j]+=jaclm[j]*alpha; //jac[l*m+j]
jacTjacim[j]+=jaclm[j]*alpha; //jacTjac[i*m+j]+=jac[l*m+j]*alpha
/* J^T e */
jacTe[i]+=alpha*e[l];
......@@ -712,14 +721,19 @@ if(!(k%100)){
/* solve augmented equations */
#ifdef HAVE_LAPACK
/* 6 alternatives are available: LU, Cholesky, 2 variants of QR decomposition, SVD and LDLt.
* Cholesky is the fastest but might be inaccurate; QR is slower but more accurate;
* SVD is the slowest but most accurate; LU offers a tradeoff between accuracy and speed
/* 7 alternatives are available: LU, Cholesky + Cholesky with PLASMA, LDLt, 2 variants of QR decomposition and SVD.
* For matrices with dimensions of at least a few hundreds, the PLASMA implementation of Cholesky is the fastest.
* From the serial solvers, Cholesky is the fastest but might occasionally be inapplicable due to numerical round-off;
* QR is slower but more robust; SVD is the slowest but most robust; LU is quite robust but
* slower than LDLt; LDLt offers a good tradeoff between robustness and speed
*/
//issolved=AX_EQ_B_BK(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_BK;
issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU;
issolved=AX_EQ_B_BK(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_BK;
//issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU;
//issolved=AX_EQ_B_CHOL(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_CHOL;
#ifdef HAVE_PLASMA
//issolved=AX_EQ_B_PLASMA_CHOL(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_PLASMA_CHOL;
#endif
//issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR;
//issolved=AX_EQ_B_QRLS(jacTjac, jacTe, Dp, m, m); ++nlss; linsolver=(int (*)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m))AX_EQ_B_QRLS;
//issolved=AX_EQ_B_SVD(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_SVD;
......@@ -864,6 +878,7 @@ if(!(k%100)){
#undef LEVMAR_L2NRMXMY
#undef AX_EQ_B_LU
#undef AX_EQ_B_CHOL
#undef AX_EQ_B_PLASMA_CHOL
#undef AX_EQ_B_QR
#undef AX_EQ_B_QRLS
#undef AX_EQ_B_SVD
......
src/levmar/lmbc.c
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......@@ -61,6 +61,8 @@
#define EPSILON 1E-12
#define ONE_THIRD 0.3333333334 /* 1.0/3.0 */
#define _LSITMAX_ 150 /* max #iterations for line search */
#define _POW_ 2.1
#if !defined(LM_DBL_PREC) && !defined(LM_SNGL_PREC)
#error At least one of LM_DBL_PREC, LM_SNGL_PREC should be defined!
......
src/levmar/lmbc_core.c
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......@@ -53,7 +53,9 @@
#define FUNC_STATE LM_ADD_PREFIX(func_state)
#define LNSRCH LM_ADD_PREFIX(lnsrch)
#define BOXPROJECT LM_ADD_PREFIX(boxProject)
#define BOXSCALE LM_ADD_PREFIX(boxScale)
#define LEVMAR_BOX_CHECK LM_ADD_PREFIX(levmar_box_check)
#define VECNORM LM_ADD_PREFIX(vecnorm)
#define LEVMAR_BC_DER LM_ADD_PREFIX(levmar_bc_der)
#define LEVMAR_BC_DIF LM_ADD_PREFIX(levmar_bc_dif)
#define LEVMAR_FDIF_FORW_JAC_APPROX LM_ADD_PREFIX(levmar_fdif_forw_jac_approx)
......@@ -76,31 +78,143 @@
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack)
#endif /* HAVE_LAPACK */
#ifdef HAVE_PLASMA
#define AX_EQ_B_PLASMA_CHOL LM_ADD_PREFIX(Ax_eq_b_PLASMA_Chol)
#endif
/* find the median of 3 numbers */
#define __MEDIAN3(a, b, c) ( ((a) >= (b))?\
( ((c) >= (a))? (a) : ( ((c) <= (b))? (b) : (c) ) ) : \
( ((c) >= (b))? (b) : ( ((c) <= (a))? (a) : (c) ) ) )
#define _POW_ LM_CNST(2.1)
/* Projections to feasible set \Omega: P_{\Omega}(y) := arg min { ||x - y|| : x \in \Omega}, y \in R^m */
/* project vector p to a box shaped feasible set. p is a mx1 vector.
* Either lb, ub can be NULL. If not NULL, they are mx1 vectors
*/
static void BOXPROJECT(LM_REAL *p, LM_REAL *lb, LM_REAL *ub, int m)
{
register int i;
if(!lb){ /* no lower bounds */
if(!ub) /* no upper bounds */
return;
else{ /* upper bounds only */
for(i=m; i-->0; )
if(p[i]>ub[i]) p[i]=ub[i];
}
}
else
if(!ub){ /* lower bounds only */
for(i=m; i-->0; )
if(p[i]<lb[i]) p[i]=lb[i];
}
else /* box bounds */
for(i=m; i-->0; )
p[i]=__MEDIAN3(lb[i], p[i], ub[i]);
}
#undef __MEDIAN3
/* pointwise scaling of bounds with the mx1 vector scl. If div=1 scaling is by 1./scl.
* Either lb, ub can be NULL. If not NULL, they are mx1 vectors
*/
static void BOXSCALE(LM_REAL *lb, LM_REAL *ub, LM_REAL *scl, int m, int div)
{
register int i;
if(!lb){ /* no lower bounds */
if(!ub) /* no upper bounds */
return;
else{ /* upper bounds only */
if(div){
for(i=m; i-->0; )
if(ub[i]!=LM_REAL_MAX)
ub[i]=ub[i]/scl[i];
}else{
for(i=m; i-->0; )
if(ub[i]!=LM_REAL_MAX)
ub[i]=ub[i]*scl[i];
}
}
}
else
if(!ub){ /* lower bounds only */
if(div){
for(i=m; i-->0; )
if(lb[i]!=LM_REAL_MIN)
lb[i]=lb[i]/scl[i];
}else{
for(i=m; i-->0; )
if(lb[i]!=LM_REAL_MIN)
lb[i]=lb[i]*scl[i];
}
}
else{ /* box bounds */
if(div){
for(i=m; i-->0; ){
if(ub[i]!=LM_REAL_MAX)
ub[i]=ub[i]/scl[i];
if(lb[i]!=LM_REAL_MIN)
lb[i]=lb[i]/scl[i];
}
}else{
for(i=m; i-->0; ){
if(ub[i]!=LM_REAL_MAX)
ub[i]=ub[i]*scl[i];
if(lb[i]!=LM_REAL_MIN)
lb[i]=lb[i]*scl[i];
}
}
}
}
/* compute the norm of a vector in a manner that avoids overflows
*/
static LM_REAL VECNORM(LM_REAL *x, int n)
{
#ifdef HAVE_LAPACK
#define NRM2 LM_MK_BLAS_NAME(nrm2)
extern LM_REAL NRM2(int *n, LM_REAL *dx, int *incx);
int one=1;
return NRM2(&n, x, &one);
#undef NRM2
#else // no LAPACK, use the simple method described by Blue in TOMS78
register int i;
LM_REAL max, sum, tmp;
#define __LSITMAX 150 // max #iterations for line search
for(i=n, max=0.0; i-->0; )
if(x[i]>max) max=x[i];
else if(x[i]<-max) max=-x[i];
for(i=n, sum=0.0; i-->0; ){
tmp=x[i]/max;
sum+=tmp*tmp;
}
return max*(LM_REAL)sqrt(sum);
#endif /* HAVE_LAPACK */
}
struct FUNC_STATE{
int n, *nfev;
LM_REAL *hx, *x;
LM_REAL *lb, *ub;
void *adata;
};
static void
LNSRCH(int m, LM_REAL *x, LM_REAL f, LM_REAL *g, LM_REAL *p, LM_REAL alpha, LM_REAL *xpls,
LM_REAL *ffpls, void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata), struct FUNC_STATE state,
LM_REAL *ffpls, void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata), struct FUNC_STATE *state,
int *mxtake, int *iretcd, LM_REAL stepmx, LM_REAL steptl, LM_REAL *sx)
{
/* Find a next newton iterate by backtracking line search.
* Specifically, finds a \lambda such that for a fixed alpha<0.5 (usually 1e-4),
* f(x + \lambda*p) <= f(x) + alpha * \lambda * g^T*p
*
* Translated (with minor changes) from Schnabel, Koontz & Weiss uncmin.f, v1.3
* Translated (with a few changes) from Schnabel, Koontz & Weiss uncmin.f, v1.3
* Main changes include the addition of box projection and modification of the scaling
* logic since uncmin.f operates in the original (unscaled) variable space.
* PARAMETERS :
......@@ -141,52 +255,48 @@ LNSRCH(int m, LM_REAL *x, LM_REAL f, LM_REAL *g, LM_REAL *p, LM_REAL alpha, LM_R
*mxtake = 0;
*iretcd = 2;
tmp1 = 0.;
if(!sx) /* no scaling */
for (i = 0; i < m; ++i)
tmp1 += p[i] * p[i];
else
for (i = 0; i < m; ++i)
tmp1 += sx[i] * sx[i] * p[i] * p[i];
for (i = m; i-- > 0; )
tmp1 += p[i] * p[i];
sln = (LM_REAL)sqrt(tmp1);
if (sln > stepmx) {
/* newton step longer than maximum allowed */
scl = stepmx / sln;
for(i=0; i<m; ++i) /* p * scl */
for (i = m; i-- > 0; ) /* p * scl */
p[i]*=scl;
sln = stepmx;
}
for(i=0, slp=0.; i<m; ++i) /* g^T * p */
slp+=g[i]*p[i];
rln = 0.;
if(!sx) /* no scaling */
for (i = 0; i < m; ++i) {
tmp1 = (FABS(x[i])>=LM_CNST(1.))? FABS(x[i]) : LM_CNST(1.);
tmp2 = FABS(p[i])/tmp1;
if(rln < tmp2) rln = tmp2;
}
else
for (i = 0; i < m; ++i) {
tmp1 = (FABS(x[i])>=LM_CNST(1.)/sx[i])? FABS(x[i]) : LM_CNST(1.)/sx[i];
tmp2 = FABS(p[i])/tmp1;
if(rln < tmp2) rln = tmp2;
}
for (i = m, slp = rln = 0.; i-- > 0; ){
slp+=g[i]*p[i]; /* g^T * p */
tmp1 = (FABS(x[i])>=LM_CNST(1.))? FABS(x[i]) : LM_CNST(1.);
tmp2 = FABS(p[i])/tmp1;
if(rln < tmp2) rln = tmp2;
}
rmnlmb = steptl / rln;
lambda = LM_CNST(1.0);
/* check if new iterate satisfactory. generate new lambda if necessary. */
for(j=__LSITMAX; j>=0; --j) {
for (i = 0; i < m; ++i)
xpls[i] = x[i] + lambda * p[i];
for(j = _LSITMAX_; j-- > 0; ) {
for (i = m; i-- > 0; )
xpls[i] = x[i] + lambda * p[i];
BOXPROJECT(xpls, state->lb, state->ub, m); /* project to feasible set */
/* evaluate function at new point */
(*func)(xpls, state.hx, m, state.n, state.adata); ++(*(state.nfev));
/* ### state.hx=state.x-state.hx, tmp1=||state.hx|| */
if(!sx){
(*func)(xpls, state->hx, m, state->n, state->adata); ++(*(state->nfev));
}
else{
for (i = m; i-- > 0; ) xpls[i] *= sx[i];
(*func)(xpls, state->hx, m, state->n, state->adata); ++(*(state->nfev));
for (i = m; i-- > 0; ) xpls[i] /= sx[i];
}
/* ### state->hx=state->x-state->hx, tmp1=||state->hx|| */
#if 1
tmp1=LEVMAR_L2NRMXMY(state.hx, state.x, state.hx, state.n);
tmp1=LEVMAR_L2NRMXMY(state->hx, state->x, state->hx, state->n);
#else
for(i=0, tmp1=0.0; i<state.n; ++i){
state.hx[i]=tmp2=state.x[i]-state.hx[i];
for(i=0, tmp1=0.0; i<state->n; ++i){
state->hx[i]=tmp2=state->x[i]-state->hx[i];
tmp1+=tmp2*tmp2;
}
#endif
......@@ -253,33 +363,6 @@ LNSRCH(int m, LM_REAL *x, LM_REAL f, LM_REAL *g, LM_REAL *p, LM_REAL alpha, LM_R
return;
} /* LNSRCH */
/* Projections to feasible set \Omega: P_{\Omega}(y) := arg min { ||x - y|| : x \in \Omega}, y \in R^m */
/* project vector p to a box shaped feasible set. p is a mx1 vector.
* Either lb, ub can be NULL. If not NULL, they are mx1 vectors
*/
static void BOXPROJECT(LM_REAL *p, LM_REAL *lb, LM_REAL *ub, int m)
{
register int i;
if(!lb){ /* no lower bounds */
if(!ub) /* no upper bounds */
return;
else{ /* upper bounds only */
for(i=0; i<m; ++i)
if(p[i]>ub[i]) p[i]=ub[i];
}
}
else
if(!ub){ /* lower bounds only */
for(i=0; i<m; ++i)
if(p[i]<lb[i]) p[i]=lb[i];
}
else /* box bounds */
for(i=0; i<m; ++i)
p[i]=__MEDIAN3(lb[i], p[i], ub[i]);
}
/*
* This function seeks the parameter vector p that best describes the measurements
* vector x under box constraints.
......@@ -299,6 +382,15 @@ register int i;
* Journal of Computational and Applied Mathematics 172, 2004, pp. 375-397.
* Also, see K. Madsen, H.B. Nielsen and O. Tingleff's lecture notes on
* unconstrained Levenberg-Marquardt at http://www.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
*
* The algorithm implemented by this function employs projected gradient steps. Since steepest descent
* is very sensitive to poor scaling, diagonal scaling has been implemented through the dscl argument:
* Instead of minimizing f(p) for p, f(D*q) is minimized for q=D^-1*p, D being a diagonal scaling
* matrix whose diagonal equals dscl (see Nocedal-Wright p.27). dscl should contain "typical" magnitudes
* for the parameters p. A NULL value for dscl implies no scaling. i.e. D=I.
* To account for scaling, the code divides the starting point and box bounds pointwise by dscl. Moreover,
* before calling func and jacf the scaling has to be undone (by multiplying), as should be done with
* the final point. Note also that jac_q=jac_p*D, where jac_q, jac_p are the jacobians w.r.t. q & p, resp.
*/
int LEVMAR_BC_DER(
......@@ -310,6 +402,7 @@ int LEVMAR_BC_DER(
int n, /* I: measurement vector dimension */
LM_REAL *lb, /* I: vector of lower bounds. If NULL, no lower bounds apply */
LM_REAL *ub, /* I: vector of upper bounds. If NULL, no upper bounds apply */
LM_REAL *dscl, /* I: diagonal scaling constants. NULL implies no scaling */
int itmax, /* I: maximum number of iterations */
LM_REAL opts[4], /* I: minim. options [\mu, \epsilon1, \epsilon2, \epsilon3]. Respectively the scale factor for initial \mu,
* stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2. Set to NULL for defaults to be used.
......@@ -347,7 +440,8 @@ LM_REAL *e, /* nx1 */
*jacTjac, /* mxm */
*Dp, /* mx1 */
*diag_jacTjac, /* diagonal of J^T J, mx1 */
*pDp; /* p + Dp, mx1 */
*pDp, /* p + Dp, mx1 */
*sp_pDp=NULL; /* dscl*p or dscl*pDp, mx1 */
register LM_REAL mu, /* damping constant */
tmp; /* mainly used in matrix & vector multiplications */
......@@ -360,9 +454,8 @@ const int nm=n*m;
/* variables for constrained LM */
struct FUNC_STATE fstate;
LM_REAL alpha=LM_CNST(1e-4), beta=LM_CNST(0.9), gamma=LM_CNST(0.99995), gamma_sq=gamma*gamma, rho=LM_CNST(1e-8);
LM_REAL t, t0;
LM_REAL steptl=LM_CNST(1e3)*(LM_REAL)sqrt(LM_REAL_EPSILON), jacTeDp;
LM_REAL alpha=LM_CNST(1e-4), beta=LM_CNST(0.9), gamma=LM_CNST(0.99995), rho=LM_CNST(1e-8);
LM_REAL t, t0, jacTeDp;
LM_REAL tmin=LM_CNST(1e-12), tming=LM_CNST(1e-18); /* minimum step length for LS and PG steps */
const LM_REAL tini=LM_CNST(1.0); /* initial step length for LS and PG steps */
int nLMsteps=0, nLSsteps=0, nPGsteps=0, gprevtaken=0;
......@@ -387,6 +480,20 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
return LM_ERROR;
}
if(dscl){ /* check that scaling consts are valid */
for(i=m; i-->0; )
if(dscl[i]<=0.0){
fprintf(stderr, LCAT(LEVMAR_BC_DER, "(): scaling constants should be positive (scale %d: %g <= 0)\n"), i, dscl[i]);
return LM_ERROR;
}
sp_pDp=(LM_REAL *)malloc(m*sizeof(LM_REAL));
if(!sp_pDp){
fprintf(stderr, LCAT(LEVMAR_BC_DER, "(): memory allocation request failed\n"));
return LM_ERROR;
}
}
if(opts){
tau=opts[0];
eps1=opts[1];
......@@ -425,10 +532,12 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
fstate.n=n;
fstate.hx=hx;
fstate.x=x;
fstate.lb=lb;
fstate.ub=ub;
fstate.adata=adata;
fstate.nfev=&nfev;
/* see if starting point is within the feasile set */
/* see if starting point is within the feasible set */
for(i=0; i<m; ++i)
pDp[i]=p[i];
BOXPROJECT(p, lb, ub, m); /* project to feasible set */
......@@ -451,6 +560,12 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
init_p_eL2=p_eL2;
if(!LM_FINITE(p_eL2)) stop=7;
if(dscl){
/* scale starting point and constraints */
for(i=m; i-->0; ) p[i]/=dscl[i];
BOXSCALE(lb, ub, dscl, m, 1);
}
for(k=0; k<itmax && !stop; ++k){
/* Note that p and e have been updated at a previous iteration */
......@@ -464,7 +579,22 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
* only its upper triangular part and copying it to the lower part
*/
(*jacf)(p, jac, m, n, adata); ++njev;
if(!dscl){
(*jacf)(p, jac, m, n, adata); ++njev;
}
else{
for(i=m; i-->0; ) sp_pDp[i]=p[i]*dscl[i];
(*jacf)(sp_pDp, jac, m, n, adata); ++njev;
/* compute jac*D */
for(i=n; i-->0; ){
register LM_REAL *jacim;
jacim=jac+i*m;
for(j=m; j-->0; )
jacim[j]*=dscl[j]; // jac[i*m+j]*=dscl[j];
}
}
/* J^T J, J^T e */
if(nm<__BLOCKSZ__SQ){ // this is a small problem
......@@ -486,8 +616,7 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
* Note that the non-blocking algorithm is faster on small
* problems since in this case it avoids the overheads of blocking.
*/
register int l, im;
register LM_REAL alpha, *jaclm;
register LM_REAL alpha, *jaclm, *jacTjacim;
/* looping downwards saves a few computations */
for(i=m*m; i-->0; )
......@@ -498,10 +627,10 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
for(l=n; l-->0; ){
jaclm=jac+l*m;
for(i=m; i-->0; ){
im=i*m;
jacTjacim=jacTjac+i*m;
alpha=jaclm[i]; //jac[l*m+i];
for(j=i+1; j-->0; ) /* j<=i computes lower triangular part only */
jacTjac[im+j]+=jaclm[j]*alpha; //jac[l*m+j]
jacTjacim[j]+=jaclm[j]*alpha; //jacTjac[i*m+j]+=jac[l*m+j]*alpha
/* J^T e */
jacTe[i]+=alpha*e[l];
......@@ -579,14 +708,19 @@ if(!(k%100)){
/* solve augmented equations */
#ifdef HAVE_LAPACK
/* 6 alternatives are available: LU, Cholesky, 2 variants of QR decomposition, SVD and LDLt.
* Cholesky is the fastest but might be inaccurate; QR is slower but more accurate;
* SVD is the slowest but most accurate; LU offers a tradeoff between accuracy and speed
/* 7 alternatives are available: LU, Cholesky + Cholesky with PLASMA, LDLt, 2 variants of QR decomposition and SVD.
* For matrices with dimensions of at least a few hundreds, the PLASMA implementation of Cholesky is the fastest.
* From the serial solvers, Cholesky is the fastest but might occasionally be inapplicable due to numerical round-off;
* QR is slower but more robust; SVD is the slowest but most robust; LU is quite robust but
* slower than LDLt; LDLt offers a good tradeoff between robustness and speed
*/
//issolved=AX_EQ_B_BK(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_BK;
issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU;
issolved=AX_EQ_B_BK(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_BK;
//issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU;
//issolved=AX_EQ_B_CHOL(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_CHOL;
#ifdef HAVE_PLASMA
//issolved=AX_EQ_B_PLASMA_CHOL(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_PLASMA_CHOL;
#endif
//issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR;
//issolved=AX_EQ_B_QRLS(jacTjac, jacTe, Dp, m, m); ++nlss; linsolver=(int (*)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m))AX_EQ_B_QRLS;
//issolved=AX_EQ_B_SVD(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_SVD;
......@@ -618,7 +752,14 @@ if(!(k%100)){
break;
}
(*func)(pDp, hx, m, n, adata); ++nfev; /* evaluate function at p + Dp */
if(!dscl){
(*func)(pDp, hx, m, n, adata); ++nfev; /* evaluate function at p + Dp */
}
else{
for(i=m; i-->0; ) sp_pDp[i]=pDp[i]*dscl[i];
(*func)(sp_pDp, hx, m, n, adata); ++nfev; /* evaluate function at p + Dp */
}
/* ### hx=x-hx, pDp_eL2=||hx|| */
#if 1
pDp_eL2=LEVMAR_L2NRMXMY(hx, x, hx, n);
......@@ -628,12 +769,15 @@ if(!(k%100)){
pDp_eL2+=tmp*tmp;
}
#endif
if(!LM_FINITE(pDp_eL2)){
/* the following test ensures that the computation of pDp_eL2 has not overflowed.
* Such an overflow does no harm here, thus it is not signalled as an error
*/
if(!LM_FINITE(pDp_eL2) && !LM_FINITE(VECNORM(hx, n))){
stop=7;
break;
}
if(pDp_eL2<=gamma_sq*p_eL2){
if(pDp_eL2<=gamma*p_eL2){
for(i=0, dL=0.0; i<m; ++i)
dL+=Dp[i]*(mu*Dp[i]+jacTe[i]);
......@@ -644,11 +788,14 @@ if(!(k%100)){
tmp=LM_CNST(1.0)-tmp*tmp*tmp;
mu=mu*( (tmp>=LM_CNST(ONE_THIRD))? tmp : LM_CNST(ONE_THIRD) );
}
else
mu=(mu>=pDp_eL2)? pDp_eL2 : mu; /* pDp_eL2 is the new pDp_eL2 */
else{
tmp=LM_CNST(0.1)*pDp_eL2; /* pDp_eL2 is the new p_eL2 */
mu=(mu>=tmp)? tmp : mu;
}
#else
mu=(mu>=pDp_eL2)? pDp_eL2 : mu; /* pDp_eL2 is the new pDp_eL2 */
tmp=LM_CNST(0.1)*pDp_eL2; /* pDp_eL2 is the new p_eL2 */
mu=(mu>=tmp)? tmp : mu;
#endif
nu=2;
......@@ -663,6 +810,7 @@ if(!(k%100)){
gprevtaken=0;
break;
}
/* note that if the LM step is not taken, code falls through to the LM line search below */
}
else{
......@@ -692,37 +840,52 @@ if(!(k%100)){
jacTeDp+=jacTe[i]*Dp[i];
}
if(jacTeDp<=-rho*pow(Dp_L2, _POW_/LM_CNST(2.0))){
if(jacTeDp<=-rho*pow(Dp_L2, LM_CNST(_POW_)/LM_CNST(2.0))){
/* Dp is a descent direction; do a line search along it */
#if 1
/* use Schnabel's backtracking line search; it requires fewer "func" evaluations */
{
int mxtake, iretcd;
LM_REAL stepmx;
LM_REAL stepmx, steptl=LM_CNST(1e3)*(LM_REAL)sqrt(LM_REAL_EPSILON);
tmp=(LM_REAL)sqrt(p_L2); stepmx=LM_CNST(1e3)*( (tmp>=LM_CNST(1.0))? tmp : LM_CNST(1.0) );
#if 1
/* use Schnabel's backtracking line search; it requires fewer "func" evaluations */
LNSRCH(m, p, p_eL2, jacTe, Dp, alpha, pDp, &pDp_eL2, func, fstate,
&mxtake, &iretcd, stepmx, steptl, NULL); /* NOTE: LNSRCH() updates hx */
if(iretcd!=0) goto gradproj; /* rather inelegant but effective way to handle LNSRCH() failures... */
LNSRCH(m, p, p_eL2, jacTe, Dp, alpha, pDp, &pDp_eL2, func, &fstate,
&mxtake, &iretcd, stepmx, steptl, dscl); /* NOTE: LNSRCH() updates hx */
if(iretcd!=0 || !LM_FINITE(pDp_eL2)) goto gradproj; /* rather inelegant but effective way to handle LNSRCH() failures... */
}
#else
/* use the simpler (but slower!) line search described by Kanzow et al */
for(t=tini; t>tmin; t*=beta){
for(i=0; i<m; ++i){
for(i=0; i<m; ++i)
pDp[i]=p[i] + t*Dp[i];
//pDp[i]=__MEDIAN3(lb[i], pDp[i], ub[i]); /* project to feasible set */
BOXPROJECT(pDp, lb, ub, m); /* project to feasible set */
if(!dscl){
(*func)(pDp, hx, m, n, adata); ++nfev; /* evaluate function at p + t*Dp */
}
else{
for(i=m; i-->0; ) sp_pDp[i]=pDp[i]*dscl[i];
(*func)(sp_pDp, hx, m, n, adata); ++nfev; /* evaluate function at p + t*Dp */
}
(*func)(pDp, hx, m, n, adata); ++nfev; /* evaluate function at p + t*Dp */
for(i=0, pDp_eL2=0.0; i<n; ++i){ /* compute ||e(pDp)||_2 */
/* compute ||e(pDp)||_2 */
/* ### hx=x-hx, pDp_eL2=||hx|| */
#if 1
pDp_eL2=LEVMAR_L2NRMXMY(hx, x, hx, n);
#else
for(i=0, pDp_eL2=0.0; i<n; ++i){
hx[i]=tmp=x[i]-hx[i];
pDp_eL2+=tmp*tmp;
}
#endif /* ||e(pDp)||_2 */
if(!LM_FINITE(pDp_eL2)) goto gradproj; /* treat as line search failure */
//if(LM_CNST(0.5)*pDp_eL2<=LM_CNST(0.5)*p_eL2 + t*alpha*jacTeDp) break;
if(pDp_eL2<=p_eL2 + LM_CNST(2.0)*t*alpha*jacTeDp) break;
}
#endif
#endif /* line search alternatives */
++nLSsteps;
gprevtaken=0;
......@@ -731,11 +894,13 @@ if(!(k%100)){
*/
}
else{
gradproj: /* Note that this point can also be reached via a goto when LNSRCH() fails */
/* Note that this point can also be reached via a goto when LNSRCH() fails. */
gradproj:
/* jacTe is a descent direction; make a projected gradient step */
/* jacTe has been negated above. Being a descent direction, it is next used
* to make a projected gradient step
*/
/* if the previous step was along the gradient descent, try to use the t employed in that step */
/* compute ||g|| */
for(i=0, tmp=0.0; i<m; ++i)
tmp+=jacTe[i]*jacTe[i];
......@@ -743,14 +908,24 @@ gradproj: /* Note that this point can also be reached via a goto when LNSRCH() f
tmp=LM_CNST(100.0)/(LM_CNST(1.0)+tmp);
t0=(tmp<=tini)? tmp : tini; /* guard against poor scaling & large steps; see (3.50) in C.T. Kelley's book */
/* if the previous step was along the gradient descent, try to use the t employed in that step */
for(t=(gprevtaken)? t : t0; t>tming; t*=beta){
for(i=0; i<m; ++i)
pDp[i]=p[i] - t*jacTe[i];
BOXPROJECT(pDp, lb, ub, m); /* project to feasible set */
for(i=0; i<m; ++i)
Dp[i]=pDp[i]-p[i];
for(i=0, Dp_L2=0.0; i<m; ++i){
Dp[i]=tmp=pDp[i]-p[i];
Dp_L2+=tmp*tmp;
}
if(!dscl){
(*func)(pDp, hx, m, n, adata); ++nfev; /* evaluate function at p - t*g */
}
else{
for(i=m; i-->0; ) sp_pDp[i]=pDp[i]*dscl[i];
(*func)(sp_pDp, hx, m, n, adata); ++nfev; /* evaluate function at p - t*g */
}
(*func)(pDp, hx, m, n, adata); ++nfev; /* evaluate function at p - t*g */
/* compute ||e(pDp)||_2 */
/* ### hx=x-hx, pDp_eL2=||hx|| */
#if 1
......@@ -761,22 +936,36 @@ gradproj: /* Note that this point can also be reached via a goto when LNSRCH() f
pDp_eL2+=tmp*tmp;
}
#endif
if(!LM_FINITE(pDp_eL2)){
/* the following test ensures that the computation of pDp_eL2 has not overflowed.
* Such an overflow does no harm here, thus it is not signalled as an error
*/
if(!LM_FINITE(pDp_eL2) && !LM_FINITE(VECNORM(hx, n))){
stop=7;
goto breaknested;
}
for(i=0, tmp=0.0; i<m; ++i) /* compute ||g^T * Dp|| */
tmp+=jacTe[i]*Dp[i];
/* compute ||g^T * Dp||. Note that if pDp has not been altered by projection
* (i.e. BOXPROJECT), jacTeDp=-t*||g||^2
*/
for(i=0, jacTeDp=0.0; i<m; ++i)
jacTeDp+=jacTe[i]*Dp[i];
if(gprevtaken && pDp_eL2<=p_eL2 + LM_CNST(2.0)*LM_CNST(0.99999)*tmp){ /* starting t too small */
if(gprevtaken && pDp_eL2<=p_eL2 + LM_CNST(2.0)*LM_CNST(0.99999)*jacTeDp){ /* starting t too small */
t=t0;
gprevtaken=0;
continue;
}
//if(LM_CNST(0.5)*pDp_eL2<=LM_CNST(0.5)*p_eL2 + alpha*tmp) break;
if(pDp_eL2<=p_eL2 + LM_CNST(2.0)*alpha*tmp) break;
//if(LM_CNST(0.5)*pDp_eL2<=LM_CNST(0.5)*p_eL2 + alpha*jacTeDp) terminatePGLS;
if(pDp_eL2<=p_eL2 + LM_CNST(2.0)*alpha*jacTeDp) goto terminatePGLS;
//if(pDp_eL2<=p_eL2 - LM_CNST(2.0)*alpha/t*Dp_L2) goto terminatePGLS; // sufficient decrease condition proposed by Kelley in (5.13)
}
/* if this point is reached then the gradient line search has failed */
gprevtaken=0;
break;
terminatePGLS:
++nPGsteps;
gprevtaken=1;
......@@ -831,6 +1020,12 @@ breaknested: /* NOTE: this point is also reached via an explicit goto! */
/* covariance matrix */
if(covar){
LEVMAR_COVAR(jacTjac, covar, p_eL2, m, n);
if(dscl){ /* correct for the scaling */
for(i=m; i-->0; )
for(j=m; j-->0; )
covar[i*m+j]*=(dscl[i]*dscl[j]);
}
}
if(freework) free(work);
......@@ -843,6 +1038,13 @@ breaknested: /* NOTE: this point is also reached via an explicit goto! */
printf("%d LM steps, %d line search, %d projected gradient\n", nLMsteps, nLSsteps, nPGsteps);
#endif
if(dscl){
/* scale final point and constraints */
for(i=0; i<m; ++i) p[i]*=dscl[i];
BOXSCALE(lb, ub, dscl, m, 0);
free(sp_pDp);
}
return (stop!=4 && stop!=7)? k : LM_ERROR;
}
......@@ -892,6 +1094,7 @@ int LEVMAR_BC_DIF(
int n, /* I: measurement vector dimension */
LM_REAL *lb, /* I: vector of lower bounds. If NULL, no lower bounds apply */
LM_REAL *ub, /* I: vector of upper bounds. If NULL, no upper bounds apply */
LM_REAL *dscl, /* I: diagonal scaling constants. NULL implies no scaling */
int itmax, /* I: maximum number of iterations */
LM_REAL opts[5], /* I: opts[0-4] = minim. options [\mu, \epsilon1, \epsilon2, \epsilon3, \delta]. Respectively the
* scale factor for initial \mu, stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2 and
......@@ -938,7 +1141,7 @@ int ret;
data.adata=adata;
data.delta=(opts)? FABS(opts[4]) : (LM_REAL)LM_DIFF_DELTA;
ret=LEVMAR_BC_DER(LMBC_DIF_FUNC, LMBC_DIF_JACF, p, x, m, n, lb, ub, itmax, opts, info, work, covar, (void *)&data);
ret=LEVMAR_BC_DER(LMBC_DIF_FUNC, LMBC_DIF_JACF, p, x, m, n, lb, ub, dscl, itmax, opts, info, work, covar, (void *)&data);
if(info){ /* correct the number of function calls */
if(data.ffdif)
......@@ -956,7 +1159,9 @@ int ret;
#undef FUNC_STATE
#undef LNSRCH
#undef BOXPROJECT
#undef BOXSCALE
#undef LEVMAR_BOX_CHECK
#undef VECNORM
#undef LEVMAR_BC_DER
#undef LMBC_DIF_DATA
#undef LMBC_DIF_FUNC
......@@ -969,6 +1174,7 @@ int ret;
#undef LEVMAR_L2NRMXMY
#undef AX_EQ_B_LU
#undef AX_EQ_B_CHOL
#undef AX_EQ_B_PLASMA_CHOL
#undef AX_EQ_B_QR
#undef AX_EQ_B_QRLS
#undef AX_EQ_B_SVD
......
src/levmar/lmblec.c
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......
src/levmar/lmblec_core.c
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......
src/levmar/lmbleic.c
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......
src/levmar/lmbleic_core.c
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......
src/levmar/lmdemo.c
View file @ aa52c8bf
......@@ -1038,7 +1038,7 @@ char *probname[]={
lb[0]=-DBL_MAX; lb[1]=-1.5;
ub[0]=ub[1]=DBL_MAX;
ret=dlevmar_bc_der(hs01, jachs01, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
ret=dlevmar_bc_der(hs01, jachs01, p, x, m, n, lb, ub, NULL, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
......@@ -1054,7 +1054,7 @@ char *probname[]={
lb[0]=2.0; lb[1]=-50.0;
ub[0]=50.0; ub[1]=50.0;
ret=dlevmar_bc_der(hs21, jachs21, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
ret=dlevmar_bc_der(hs21, jachs21, p, x, m, n, lb, ub, NULL, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
......@@ -1072,7 +1072,7 @@ char *probname[]={
ub[0]=ub[2]=ub[3]=DBL_MAX;
ub[1]=0.8;
ret=dlevmar_bc_der(hatfldb, jachatfldb, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
ret=dlevmar_bc_der(hatfldb, jachatfldb, p, x, m, n, lb, ub, NULL, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
......@@ -1089,7 +1089,7 @@ char *probname[]={
ub[0]=ub[1]=ub[2]=ub[3]=10.0;
ret=dlevmar_bc_der(hatfldc, jachatfldc, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
ret=dlevmar_bc_der(hatfldc, jachatfldc, p, x, m, n, lb, ub, NULL, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
......@@ -1106,7 +1106,7 @@ char *probname[]={
ub[0]=ub[1]=ub[2]=ub[3]=ub[4]=100.0;
ret=dlevmar_bc_der(combust, jaccombust, p, x, m, n, lb, ub, 5000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
ret=dlevmar_bc_der(combust, jaccombust, p, x, m, n, lb, ub, NULL, 5000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
......
src/levmar/lmdemo.vcproj deleted 100644 → 0
View file @ 1b190b55
<?xml version="1.0" encoding="UTF-8"?>
<VisualStudioProject
ProjectType="Visual C++"
Version="8,00"
Name="lmdemo"
ProjectGUID="{4A085982-23EF-4B63-8D72-B901D0521326}"
Keyword="Win32Proj"
>
<Platforms>
<Platform
Name="Win32"
/>
</Platforms>
<ToolFiles>
</ToolFiles>
<Configurations>
<Configuration
Name="Debug|Win32"
OutputDirectory="Debug"
IntermediateDirectory="Debug"
ConfigurationType="1"
>
<Tool
Name="VCPreBuildEventTool"
/>
<Tool
Name="VCCustomBuildTool"
/>
<Tool
Name="VCXMLDataGeneratorTool"
/>
<Tool
Name="VCWebServiceProxyGeneratorTool"
/>
<Tool
Name="VCMIDLTool"
/>
<Tool
Name="VCCLCompilerTool"
Optimization="0"
AdditionalIncludeDirectories="."
PreprocessorDefinitions="WIN32;_DEBUG;_CONSOLE"
MinimalRebuild="true"
BasicRuntimeChecks="3"
RuntimeLibrary="3"
UsePrecompiledHeader="0"
WarningLevel="3"
Detect64BitPortabilityProblems="true"
DebugInformationFormat="4"
CompileAs="1"
/>
<Tool
Name="VCManagedResourceCompilerTool"
/>
<Tool
Name="VCResourceCompilerTool"
/>
<Tool
Name="VCPreLinkEventTool"
/>
<Tool
Name="VCLinkerTool"
AdditionalDependencies="levmar.lib clapack.lib blas.lib libF77.lib libI77.lib"
AdditionalLibraryDirectories=".\Debug;c:\RECOVER\src\lib"
GenerateManifest="true"
/>
<Tool
Name="VCALinkTool"
/>
<Tool
Name="VCManifestTool"
/>
<Tool
Name="VCXDCMakeTool"
/>
<Tool
Name="VCBscMakeTool"
/>
<Tool
Name="VCFxCopTool"
/>
<Tool
Name="VCAppVerifierTool"
/>
<Tool
Name="VCWebDeploymentTool"
/>
<Tool
Name="VCPostBuildEventTool"
/>
</Configuration>
<Configuration
Name="Release|Win32"
OutputDirectory="Release"
IntermediateDirectory="Release"
ConfigurationType="1"
>
<Tool
Name="VCPreBuildEventTool"
/>
<Tool
Name="VCCustomBuildTool"
/>
<Tool
Name="VCXMLDataGeneratorTool"
/>
<Tool
Name="VCWebServiceProxyGeneratorTool"
/>
<Tool
Name="VCMIDLTool"
/>
<Tool
Name="VCCLCompilerTool"
AdditionalIncludeDirectories="."
PreprocessorDefinitions="WIN32;NDEBUG;_CONSOLE"
RuntimeLibrary="2"
UsePrecompiledHeader="0"
WarningLevel="3"
Detect64BitPortabilityProblems="true"
DebugInformationFormat="3"
CompileAs="1"
/>
<Tool
Name="VCManagedResourceCompilerTool"
/>
<Tool
Name="VCResourceCompilerTool"
/>
<Tool
Name="VCPreLinkEventTool"
/>
<Tool
Name="VCLinkerTool"
AdditionalDependencies="levmar.lib clapack.lib blas.lib libF77.lib libI77.lib"
AdditionalLibraryDirectories=".\Release;c:\RECOVER\src\lib"
GenerateManifest="true"
/>
<Tool
Name="VCALinkTool"
/>
<Tool
Name="VCManifestTool"
/>
<Tool
Name="VCXDCMakeTool"
/>
<Tool
Name="VCBscMakeTool"
/>
<Tool
Name="VCFxCopTool"
/>
<Tool
Name="VCAppVerifierTool"
/>
<Tool
Name="VCWebDeploymentTool"
/>
<Tool
Name="VCPostBuildEventTool"
/>
</Configuration>
</Configurations>
<References>
</References>
<Files>
<Filter
Name="Header Files"
Filter="h;hpp;hxx;hm;inl;inc;xsd"
>
<File
RelativePath=".\levmar.h"
>
</File>
</Filter>
<Filter
Name="Resource Files"
Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx"
>
</Filter>
<Filter
Name="Source Files"
Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"
>
<File
RelativePath=".\lmdemo.c"
>
</File>
</Filter>
</Files>
<Globals>
</Globals>
</VisualStudioProject>
src/levmar/lmlec.c
View file @ aa52c8bf
......@@ -5,8 +5,8 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004 Manolis Lourakis (original version)
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* (C) 2004-2011 Manolis Lourakis (orig. version)
*
* Licenses: GNU General Public License
*
......@@ -22,7 +22,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/10/2010
* Last modified: 09/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/////////////////////////////////////////////////////////////////////////////////
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment