public software.sextractor

Compare Revisions

Rev 223 → Rev 232

/gpu/xsl/sextractor.xsl
16,8 → 16,8
<xsl:variable name="time" select="/VOTABLE/RESOURCE/RESOURCE[@name='MetaData']/PARAM[@name='Time']/@value"/>
<HTML>
<HEAD>
<link rel="shortcut icon" type="image/x-icon" href="http://astromatic.iap.fr/xsl/favicon.ico" />
<script type="text/javascript" src="http://astromatic.iap.fr/xsl/sorttable.js"/>
<link rel="shortcut icon" type="image/x-icon" href="http://astromatic.net/xsl/favicon.ico" />
<script type="text/javascript" src="http://astromatic.net/xsl/sorttable.js"/>
 
<style type="text/css">
p {
27,7 → 27,7
body {
margin: 10px;
background-color: #e0e0e0;
background-image: url("http://astromatic.iap.fr/xsl/body_bg.jpg");
background-image: url("http://astromatic.net/xsl/body_bg.jpg");
background-repeat: repeat-x;
background-position: top;
min-width:662px;
65,7 → 65,7
#header {
padding: 5px;
min-width: 662px;
background-image: url("http://astromatic.iap.fr/xsl/astromaticleft.png");
background-image: url("http://astromatic.net/xsl/astromaticleft.png");
background-repeat: repeat-x;
background-position: left top;
text-align: left;
132,7 → 132,7
</HEAD>
<BODY>
<div id="header">
<a href="/"><img style="vertical-align: middle; border:0px" src="http://astromatic.iap.fr/xsl/astromatic.png" title="Astromatic home" alt="Astromatic.net" /></a> Processing summary
<a href="/"><img style="vertical-align: middle; border:0px" src="http://astromatic.net/xsl/astromatic.png" title="Astromatic home" alt="Astromatic.net" /></a> Processing summary
</div>
<xsl:call-template name="VOTable"/>
</BODY>
gpu/xsl Property changes : Modified: svn:mergeinfo Merged /trunk/xsl:r222-231
/gpu/man/sex.1.in
18,7 → 18,7
image. Although it is particularly oriented towards reduction of large scale
galaxy-survey data, it performs rather well on moderately crowded star fields.
.RE
See http://terapix.iap.fr/soft/sextractor for more details.
See http://astromatic.net/software/sextractor for more details.
.SS "Operation modes:"
.TP
\fB\-h\fR, \fB\-\-help\fR
43,5 → 43,5
.PP
The full documentation for
.B SExtractor
is maintained as a Postscript manual available at
.B http://terapix.iap.fr/soft/sextractor
is maintained as a PDF manual available at
.B http://astromatic.net/software/sextractor
/gpu/man/sex.1
1,4 → 1,4
.TH SEXTRACTOR "1" "July 2010" "SWarp 2.11.7" "User Commands"
.TH SEXTRACTOR "1" "September 2010" "SWarp 2.12.3" "User Commands"
.SH NAME
sex \- extract a source catalog from an astronomical FITS image
.SH SYNOPSIS
18,7 → 18,7
image. Although it is particularly oriented towards reduction of large scale
galaxy-survey data, it performs rather well on moderately crowded star fields.
.RE
See http://terapix.iap.fr/soft/sextractor for more details.
See http://astromatic.net/software/sextractor for more details.
.SS "Operation modes:"
.TP
\fB\-h\fR, \fB\-\-help\fR
43,5 → 43,5
.PP
The full documentation for
.B SExtractor
is maintained as a Postscript manual available at
.B http://terapix.iap.fr/soft/sextractor
is maintained as a PDF manual available at
.B http://astromatic.net/software/sextractor
/gpu/configure
1,6 → 1,6
#! /bin/sh
# Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.63 for sextractor 2.11.7.
# Generated by GNU Autoconf 2.63 for sextractor 2.12.3.
#
# Report bugs to <bertin@iap.fr>.
#
750,8 → 750,8
# Identity of this package.
PACKAGE_NAME='sextractor'
PACKAGE_TARNAME='sextractor'
PACKAGE_VERSION='2.11.7'
PACKAGE_STRING='sextractor 2.11.7'
PACKAGE_VERSION='2.12.3'
PACKAGE_STRING='sextractor 2.12.3'
PACKAGE_BUGREPORT='bertin@iap.fr'
 
ac_unique_file="src/makeit.c"
808,6 → 808,8
PTHREAD_CFLAGS
PTHREAD_LIBS
PTHREAD_CC
USE_MODEL_FALSE
USE_MODEL_TRUE
LIBOBJS
LIBTOOL
ac_ct_F77
934,6 → 936,7
with_fftw
with_fftw_incdir
with_xsl_url
enable_model_fitting
enable_threads
enable_gprof
enable_best_link
1505,7 → 1508,7
# Omit some internal or obsolete options to make the list less imposing.
# This message is too long to be a string in the A/UX 3.1 sh.
cat <<_ACEOF
\`configure' configures sextractor 2.11.7 to adapt to many kinds of systems.
\`configure' configures sextractor 2.12.3 to adapt to many kinds of systems.
 
Usage: $0 [OPTION]... [VAR=VALUE]...
 
1575,7 → 1578,7
 
if test -n "$ac_init_help"; then
case $ac_init_help in
short | recursive ) echo "Configuration of sextractor 2.11.7:";;
short | recursive ) echo "Configuration of sextractor 2.12.3:";;
esac
cat <<\_ACEOF
 
1595,6 → 1598,8
optimize for fast installation [default=yes]
--disable-libtool-lock avoid locking (might break parallel builds)
--disable-largefile omit support for large files
--disable-model-fitting Disable model-fitting and ATLAS/FFTW dependency
(enabled by default)
--enable-threads[=<max_number_of_threads>]
Enable multhreading (on with up to 16 threads by
default)
1706,7 → 1711,7
test -n "$ac_init_help" && exit $ac_status
if $ac_init_version; then
cat <<\_ACEOF
sextractor configure 2.11.7
sextractor configure 2.12.3
generated by GNU Autoconf 2.63
 
Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
1720,7 → 1725,7
This file contains any messages produced by compilers while
running configure, to aid debugging if configure makes a mistake.
 
It was created by sextractor $as_me 2.11.7, which was
It was created by sextractor $as_me 2.12.3, which was
generated by GNU Autoconf 2.63. Invocation command line was
 
$ $0 $@
2423,7 → 2428,7
 
# Define the identity of the package.
PACKAGE='sextractor'
VERSION='2.11.7'
VERSION='2.12.3'
 
 
cat >>confdefs.h <<_ACEOF
5316,7 → 5321,7
;;
*-*-irix6*)
# Find out which ABI we are using.
echo '#line 5319 "configure"' > conftest.$ac_ext
echo '#line 5324 "configure"' > conftest.$ac_ext
if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5
(eval $ac_compile) 2>&5
ac_status=$?
8091,11 → 8096,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:8094: $lt_compile\"" >&5)
(eval echo "\"\$as_me:8099: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
echo "$as_me:8098: \$? = $ac_status" >&5
echo "$as_me:8103: \$? = $ac_status" >&5
if (exit $ac_status) && test -s "$ac_outfile"; then
# The compiler can only warn and ignore the option if not recognized
# So say no if there are warnings other than the usual output.
8381,11 → 8386,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:8384: $lt_compile\"" >&5)
(eval echo "\"\$as_me:8389: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
echo "$as_me:8388: \$? = $ac_status" >&5
echo "$as_me:8393: \$? = $ac_status" >&5
if (exit $ac_status) && test -s "$ac_outfile"; then
# The compiler can only warn and ignore the option if not recognized
# So say no if there are warnings other than the usual output.
8485,11 → 8490,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:8488: $lt_compile\"" >&5)
(eval echo "\"\$as_me:8493: $lt_compile\"" >&5)
(eval "$lt_compile" 2>out/conftest.err)
ac_status=$?
cat out/conftest.err >&5
echo "$as_me:8492: \$? = $ac_status" >&5
echo "$as_me:8497: \$? = $ac_status" >&5
if (exit $ac_status) && test -s out/conftest2.$ac_objext
then
# The compiler can only warn and ignore the option if not recognized
10872,7 → 10877,7
lt_dlunknown=0; lt_dlno_uscore=1; lt_dlneed_uscore=2
lt_status=$lt_dlunknown
cat > conftest.$ac_ext <<EOF
#line 10875 "configure"
#line 10880 "configure"
#include "confdefs.h"
 
#if HAVE_DLFCN_H
10972,7 → 10977,7
lt_dlunknown=0; lt_dlno_uscore=1; lt_dlneed_uscore=2
lt_status=$lt_dlunknown
cat > conftest.$ac_ext <<EOF
#line 10975 "configure"
#line 10980 "configure"
#include "confdefs.h"
 
#if HAVE_DLFCN_H
13400,11 → 13405,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:13403: $lt_compile\"" >&5)
(eval echo "\"\$as_me:13408: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
echo "$as_me:13407: \$? = $ac_status" >&5
echo "$as_me:13412: \$? = $ac_status" >&5
if (exit $ac_status) && test -s "$ac_outfile"; then
# The compiler can only warn and ignore the option if not recognized
# So say no if there are warnings other than the usual output.
13504,11 → 13509,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:13507: $lt_compile\"" >&5)
(eval echo "\"\$as_me:13512: $lt_compile\"" >&5)
(eval "$lt_compile" 2>out/conftest.err)
ac_status=$?
cat out/conftest.err >&5
echo "$as_me:13511: \$? = $ac_status" >&5
echo "$as_me:13516: \$? = $ac_status" >&5
if (exit $ac_status) && test -s out/conftest2.$ac_objext
then
# The compiler can only warn and ignore the option if not recognized
15068,11 → 15073,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:15071: $lt_compile\"" >&5)
(eval echo "\"\$as_me:15076: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
echo "$as_me:15075: \$? = $ac_status" >&5
echo "$as_me:15080: \$? = $ac_status" >&5
if (exit $ac_status) && test -s "$ac_outfile"; then
# The compiler can only warn and ignore the option if not recognized
# So say no if there are warnings other than the usual output.
15172,11 → 15177,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:15175: $lt_compile\"" >&5)
(eval echo "\"\$as_me:15180: $lt_compile\"" >&5)
(eval "$lt_compile" 2>out/conftest.err)
ac_status=$?
cat out/conftest.err >&5
echo "$as_me:15179: \$? = $ac_status" >&5
echo "$as_me:15184: \$? = $ac_status" >&5
if (exit $ac_status) && test -s out/conftest2.$ac_objext
then
# The compiler can only warn and ignore the option if not recognized
17369,11 → 17374,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:17372: $lt_compile\"" >&5)
(eval echo "\"\$as_me:17377: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
echo "$as_me:17376: \$? = $ac_status" >&5
echo "$as_me:17381: \$? = $ac_status" >&5
if (exit $ac_status) && test -s "$ac_outfile"; then
# The compiler can only warn and ignore the option if not recognized
# So say no if there are warnings other than the usual output.
17659,11 → 17664,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:17662: $lt_compile\"" >&5)
(eval echo "\"\$as_me:17667: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
echo "$as_me:17666: \$? = $ac_status" >&5
echo "$as_me:17671: \$? = $ac_status" >&5
if (exit $ac_status) && test -s "$ac_outfile"; then
# The compiler can only warn and ignore the option if not recognized
# So say no if there are warnings other than the usual output.
17763,11 → 17768,11
-e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \
-e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
-e 's:$: $lt_compiler_flag:'`
(eval echo "\"\$as_me:17766: $lt_compile\"" >&5)
(eval echo "\"\$as_me:17771: $lt_compile\"" >&5)
(eval "$lt_compile" 2>out/conftest.err)
ac_status=$?
cat out/conftest.err >&5
echo "$as_me:17770: \$? = $ac_status" >&5
echo "$as_me:17775: \$? = $ac_status" >&5
if (exit $ac_status) && test -s out/conftest2.$ac_objext
then
# The compiler can only warn and ignore the option if not recognized
23353,6 → 23358,24
_ACEOF
 
 
# Provide special option to disable model-fitting (enabled by default)
{ $as_echo "$as_me:$LINENO: checking if model-fitting should be disabled" >&5
$as_echo_n "checking if model-fitting should be disabled... " >&6; }
# Check whether --enable-model-fitting was given.
if test "${enable_model_fitting+set}" = set; then
enableval=$enable_model_fitting; use_model="$enableval"
else
use_model="yes"
fi
 
if test "$use_model" = "no"; then
{ $as_echo "$as_me:$LINENO: result: yes" >&5
$as_echo "yes" >&6; }
else
{ $as_echo "$as_me:$LINENO: result: no" >&5
$as_echo "no" >&6; }
fi
 
# Set flags for multithreading
n_pthreads=16
# Check whether --enable-threads was given.
23399,6 → 23422,30
fi
 
 
############# Actions to complete if model-fitting is activated ##############
{ $as_echo "$as_me:$LINENO: checking for model-fitting configure option" >&5
$as_echo_n "checking for model-fitting configure option... " >&6; }
if test "$use_model" = "yes"; then
{ $as_echo "$as_me:$LINENO: result: enabled" >&5
$as_echo "enabled" >&6; }
 
cat >>confdefs.h <<\_ACEOF
#define USE_MODEL 1
_ACEOF
 
else
{ $as_echo "$as_me:$LINENO: result: disabled" >&5
$as_echo "disabled" >&6; }
fi
if test $use_model = "yes"; then
USE_MODEL_TRUE=
USE_MODEL_FALSE='#'
else
USE_MODEL_TRUE='#'
USE_MODEL_FALSE=
fi
 
 
################# Actions to complete in case of multhreading ################
 
cat >>confdefs.h <<_ACEOF
23993,6 → 24040,7
 
 
################ handle the FFTW library (Fourier transforms) ################
if test "$use_model" = "yes"; then
 
 
 
25676,19 → 25724,21
fi
 
 
if test "$use_fftw" = "yes"; then
LIBS="$FFTW_LIBS $LIBS"
if test "$FFTW_WARN" != ""; then
{ $as_echo "$as_me:$LINENO: WARNING: $FFTW_WARN" >&5
if test "$use_fftw" = "yes"; then
LIBS="$FFTW_LIBS $LIBS"
if test "$FFTW_WARN" != ""; then
{ $as_echo "$as_me:$LINENO: WARNING: $FFTW_WARN" >&5
$as_echo "$as_me: WARNING: $FFTW_WARN" >&2;}
fi
else
{ { $as_echo "$as_me:$LINENO: error: $FFTW_ERROR Exiting." >&5
fi
else
{ { $as_echo "$as_me:$LINENO: error: $FFTW_ERROR Exiting." >&5
$as_echo "$as_me: error: $FFTW_ERROR Exiting." >&2;}
{ (exit 1); exit 1; }; }
fi
fi
 
################## handle the ATLAS library(linear algebra) ##################
if test "$use_model" = "yes"; then
 
 
 
28109,16 → 28159,17
fi
 
 
if test "$use_atlas" = "yes"; then
LIBS="$ATLAS_LIB $LIBS"
if test "$ATLAS_WARN" != ""; then
{ $as_echo "$as_me:$LINENO: WARNING: $ATLAS_WARN" >&5
if test "$use_atlas" = "yes"; then
LIBS="$ATLAS_LIB $LIBS"
if test "$ATLAS_WARN" != ""; then
{ $as_echo "$as_me:$LINENO: WARNING: $ATLAS_WARN" >&5
$as_echo "$as_me: WARNING: $ATLAS_WARN" >&2;}
fi
else
{ { $as_echo "$as_me:$LINENO: error: $ATLAS_ERROR Exiting." >&5
fi
else
{ { $as_echo "$as_me:$LINENO: error: $ATLAS_ERROR Exiting." >&5
$as_echo "$as_me: error: $ATLAS_ERROR Exiting." >&2;}
{ (exit 1); exit 1; }; }
fi
fi
 
# Link with gprof option
28206,6 → 28257,13
Usually this means the macro was only invoked conditionally." >&2;}
{ (exit 1); exit 1; }; }
fi
if test -z "${USE_MODEL_TRUE}" && test -z "${USE_MODEL_FALSE}"; then
{ { $as_echo "$as_me:$LINENO: error: conditional \"USE_MODEL\" was never defined.
Usually this means the macro was only invoked conditionally." >&5
$as_echo "$as_me: error: conditional \"USE_MODEL\" was never defined.
Usually this means the macro was only invoked conditionally." >&2;}
{ (exit 1); exit 1; }; }
fi
if test -z "${USE_THREADS_TRUE}" && test -z "${USE_THREADS_FALSE}"; then
{ { $as_echo "$as_me:$LINENO: error: conditional \"USE_THREADS\" was never defined.
Usually this means the macro was only invoked conditionally." >&5
28535,7 → 28593,7
# report actual input values of CONFIG_FILES etc. instead of their
# values after options handling.
ac_log="
This file was extended by sextractor $as_me 2.11.7, which was
This file was extended by sextractor $as_me 2.12.3, which was
generated by GNU Autoconf 2.63. Invocation command line was
 
CONFIG_FILES = $CONFIG_FILES
28598,7 → 28656,7
_ACEOF
cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
ac_cs_version="\\
sextractor config.status 2.11.7
sextractor config.status 2.12.3
configured by $0, generated by GNU Autoconf 2.63,
with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"
 
/gpu/src/define.h
27,7 → 27,7
#define MYVERSION VERSION
#define EXECUTABLE "sex"
#define COPYRIGHT "Emmanuel BERTIN <bertin@iap.fr>"
#define WEBSITE "http://astromatic.iap.fr/software/sextractor"
#define WEBSITE "http://astromatic.net/software/sextractor"
#define INSTITUTE "IAP http://www.iap.fr"
 
/*--------------------------- Internal constants ----------------------------*/
/gpu/src/param.h
9,7 → 9,7
*
* Contents: parameter list for catalog data.
*
* Last modify: 18/05/2010
* Last modify: 23/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
886,7 → 886,9
"src.morph.param", "", 1, &prefs.pc_vectorsize},
*/
 
#ifdef USE_MODEL
#include "paramprofit.h"
#endif
 
{""}
};
/gpu/src/prefs.c
9,7 → 9,7
*
* Contents: Functions to handle the configuration file.
*
* Last modify: 08/03/2010
* Last modify: 22/09/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
600,7 → 600,7
prefs.prof_flag = 1;
 
/*--------------------------- Adaptive class-star ---------------------------*/
if (prefs.seeing_fwhm == 0)
if (prefs.seeing_fwhm == 0 && FLAG(obj2.sprob))
prefs.psf_flag = 1;
 
/*-------------------------- Pattern-fitting -------------------------------*/
/gpu/src/fitswcs.c
9,7 → 9,7
*
* Contents: Read and write WCS header info.
*
* Last modify: 20/02/2009
* Last modify: 30/07/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
315,7 → 315,7
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 20/02/2009
VERSION 30/07/2010
***/
wcsstruct *read_wcs(tabstruct *tab)
 
433,7 → 433,7
else
{
/*---- Search for an observation date expressed in "civilian" format */
FITSREADS(buf, "DATE-OBS ", str, "");
FITSREADS(buf, "DATE-OBS", str, "");
if (*str)
{
/*------ Decode DATE-OBS format: DD/MM/YY or YYYY-MM-DD */
/gpu/src/globals.h
9,7 → 9,7
*
* Contents: global declarations.
*
* Last modify: 01/10/2009
* Last modify: 20/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
54,6 → 54,8
neurclose(void),
neurresp(double *, double *),
preanalyse(int, objliststruct *, int),
propagate_covar(double *vi, double *d, double *vo,
int ni, int no, double *temp),
readcatparams(char *),
readdata(picstruct *, PIXTYPE *, int),
readidata(picstruct *, FLAGTYPE *, int),
/gpu/src/xml.c
9,7 → 9,7
*
* Contents: XML logging.
*
* Last modify: 25/05/2010
* Last modify: 03/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
219,7 → 219,7
OUTPUT RETURN_OK if everything went fine, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 05/02/2010
VERSION 03/08/2010
***/
int write_xml_meta(FILE *file, char *error)
{
388,8 → 388,8
xmlstack[n].ident[0], xmlstack[n].ident[1],
xmlstack[n].backmean[0], xmlstack[n].backmean[1],
xmlstack[n].backsig[0], xmlstack[n].backsig[1],
xmlstack[n].thresh[0], xmlstack[n].thresh[1],
xmlstack[n].sigfac[0], xmlstack[n].sigfac[1],
xmlstack[n].thresh[0], xmlstack[n].thresh[1],
xmlstack[n].pixscale[0], xmlstack[n].pixscale[1],
xmlstack[n].epoch[0], xmlstack[n].epoch[1],
xmlstack[n].gain[0], xmlstack[n].gain[1],
412,8 → 412,8
xmlstack[n].ident[0],
xmlstack[n].backmean[0],
xmlstack[n].backsig[0],
xmlstack[n].thresh[0],
xmlstack[n].sigfac[0],
xmlstack[n].thresh[0],
xmlstack[n].pixscale[0],
xmlstack[n].epoch[0],
xmlstack[n].gain[0],
/gpu/src/levmar/README
1,4 → 1,4
The levmar v2.4 library has been included in this package untouched, except for
The levmar v2.5 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>
/gpu/src/levmar/lm.h
19,15 → 19,15
////////////////////////////////////////////////////////////////////////////////////
*/
 
#ifndef _LM_H_
#define _LM_H_
#ifndef _LEVMAR_H_
#define _LEVMAR_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 */
/*#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
78,13 → 78,19
#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.4 (April 2009)"
#define LM_VERSION "2.5 (December 2009)"
 
#ifdef LM_DBL_PREC
/* double precision LM, with & without Jacobian */
136,6 → 142,56
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 */
191,6 → 247,56
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 */
204,6 → 310,7
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
212,6 → 319,7
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 */
259,4 → 367,4
}
#endif
 
#endif /* _LM_H_ */
#endif /* _LEVMAR_H_ */
/gpu/src/levmar/lmbc.c
28,7 → 28,7
#include <math.h>
#include <float.h>
 
#include "lm.h"
#include "levmar.h"
#include "compiler.h"
#include "misc.h"
 
/gpu/src/levmar/lmdemo.c
28,7 → 28,7
#include <math.h>
#include <float.h>
 
#include "lm.h"
#include "levmar.h"
 
#ifndef LM_DBL_PREC
#error Demo program assumes that levmar has been compiled with double precision, see LM_DBL_PREC!
162,6 → 162,36
}
}
 
/* Osborne's problem, minimum at (0.3754, 1.9358, -1.4647, 0.0129, 0.0221) */
void osborne(double *p, double *x, int m, int n, void *data)
{
register int i;
double t;
 
for(i=0; i<n; ++i){
t=10*i;
x[i]=p[0] + p[1]*exp(-p[3]*t) + p[2]*exp(-p[4]*t);
}
}
 
void jacosborne(double *p, double *jac, int m, int n, void *data)
{
register int i, j;
double t, tmp1, tmp2;
 
for(i=j=0; i<n; ++i){
t=10*i;
tmp1=exp(-p[3]*t);
tmp2=exp(-p[4]*t);
 
jac[j++]=1.0;
jac[j++]=tmp1;
jac[j++]=tmp2;
jac[j++]=-p[1]*t*tmp1;
jac[j++]=-p[2]*t*tmp2;
}
}
 
/* helical valley function, minimum at (1.0, 0.0, 0.0) */
#ifndef M_PI
#define M_PI 3.14159265358979323846 /* pi */
471,7 → 501,7
jac[j++]=1.0;
}
 
/* Hock - Schittkowski (modified) problem 52 (box/linearly constrained), minimum at (-0.09, 0.03, 0.25, -0.19, 0.03)
/* Hock - Schittkowski (modified #1) problem 52 (box/linearly constrained), minimum at (-0.09, 0.03, 0.25, -0.19, 0.03)
* constr1: p[0] + 3*p[1] = 0;
* constr2: p[2] + p[3] - 2*p[4] = 0;
* constr3: p[1] - p[4] = 0;
484,7 → 514,7
* constr8: 0.0 <= p[4] <= 0.3;
*
*/
void modhs52(double *p, double *x, int m, int n, void *data)
void mod1hs52(double *p, double *x, int m, int n, void *data)
{
x[0]=4.0*p[0]-p[1];
x[1]=p[1]+p[2]-2.0;
492,7 → 522,7
x[3]=p[4]-1.0;
}
 
void jacmodhs52(double *p, double *jac, int m, int n, void *data)
void jacmod1hs52(double *p, double *jac, int m, int n, void *data)
{
register int j=0;
 
521,6 → 551,60
jac[j++]=1.0;
}
 
 
/* Hock - Schittkowski (modified #2) problem 52 (linear inequality constrained), minimum at (0.5, 2.0, 0.0, 1.0, 1.0)
* A fifth term [(p[0]-0.5)^2] is added to the objective function and
* the equality contraints are replaced by the following inequalities:
* constr1: p[0] + 3*p[1] >= -1.0;
* constr2: p[2] + p[3] - 2*p[4] >= -2.0;
* constr3: p[1] - p[4] <= 7.0;
*
*
*/
void mod2hs52(double *p, double *x, int m, int n, void *data)
{
x[0]=4.0*p[0]-p[1];
x[1]=p[1]+p[2]-2.0;
x[2]=p[3]-1.0;
x[3]=p[4]-1.0;
x[4]=p[0]-0.5;
}
 
void jacmod2hs52(double *p, double *jac, int m, int n, void *data)
{
register int j=0;
 
jac[j++]=4.0;
jac[j++]=-1.0;
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=0.0;
 
jac[j++]=0.0;
jac[j++]=1.0;
jac[j++]=1.0;
jac[j++]=0.0;
jac[j++]=0.0;
 
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=1.0;
jac[j++]=0.0;
 
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=1.0;
 
jac[j++]=1.0;
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=0.0;
}
 
/* Schittkowski (modified) problem 235 (box/linearly constrained), minimum at (-1.725, 2.9, 0.725)
* constr1: p[0] + p[2] = -1.0;
*
653,13 → 737,55
jac[j+3]=R9*p[1]+2*p[3];
}
 
/* Hock - Schittkowski (modified) problem 76 (linear inequalities & equations constrained), minimum at (0.0, 0.00909091, 0.372727, 0.354545)
* The non-squared terms in the objective function have been removed, the rhs of constr2 has been changed to 0.4 (from 4)
* and constr3 has been changed to an equation.
*
* constr1: p[0] + 2*p[1] + p[2] + p[3] <= 5;
* constr2: 3*p[0] + p[1] + 2*p[2] - p[3] <= 0.4;
* constr3: p[1] + 4*p[2] = 1.5;
*
*/
void modhs76(double *p, double *x, int m, int n, void *data)
{
x[0]=p[0];
x[1]=sqrt(0.5)*p[1];
x[2]=p[2];
x[3]=sqrt(0.5)*p[3];
}
 
void jacmodhs76(double *p, double *jac, int m, int n, void *data)
{
register int j=0;
 
jac[j++]=1.0;
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=0.0;
 
jac[j++]=0.0;
jac[j++]=sqrt(0.5);
jac[j++]=0.0;
jac[j++]=0.0;
 
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=1.0;
jac[j++]=0.0;
 
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=0.0;
jac[j++]=sqrt(0.5);
}
 
 
 
int main()
{
register int i, j;
int problem, ret;
double p[5], // 6 is max(2, 3, 5)
double p[5], // 5 is max(2, 3, 5)
x[16]; // 16 is max(2, 3, 5, 6, 16)
int m, n;
double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
669,6 → 795,7
"Powell's function",
"Wood's function",
"Meyer's (reformulated) problem",
"Osborne's problem",
"helical valley function",
"Boggs & Tolle's problem #3",
"Hock - Schittkowski problem #28",
679,13 → 806,15
"hatfldb problem",
"hatfldc problem",
"equilibrium combustion problem",
"Hock - Schittkowski modified problem #52",
"Hock - Schittkowski modified #1 problem #52",
"Schittkowski modified problem #235",
"Boggs & Tolle modified problem #7",
"Hock - Schittkowski modified #2 problem #52",
"Hock - Schittkowski modified problem #76",
};
 
opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
opts[4]=LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing
opts[4]= LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing
//opts[4]=-LM_DIFF_DELTA; // specifies central differencing to approximate Jacobian; more accurate but more expensive to compute!
 
/* uncomment the appropriate line below to select a minimization problem */
695,12 → 824,13
//2; // Powell's function
//3; // Wood's function
4; // Meyer's (reformulated) problem
//5; // helical valley function
//5; // Osborne's problem
//6; // helical valley function
#ifdef HAVE_LAPACK
//6; // Boggs & Tolle's problem 3
//7; // Hock - Schittkowski problem 28
//8; // Hock - Schittkowski problem 48
//9; // Hock - Schittkowski problem 51
//7; // Boggs & Tolle's problem 3
//8; // Hock - Schittkowski problem 28
//9; // Hock - Schittkowski problem 48
//10; // Hock - Schittkowski problem 51
#else // no LAPACK
#ifdef _MSC_VER
#pragma message("LAPACK not available, some test problems cannot be used")
709,21 → 839,24
#endif // _MSC_VER
 
#endif /* HAVE_LAPACK */
//10; // Hock - Schittkowski problem 01
//11; // Hock - Schittkowski modified problem 21
//12; // hatfldb problem
//13; // hatfldc problem
//14; // equilibrium combustion problem
//11; // Hock - Schittkowski problem 01
//12; // Hock - Schittkowski modified problem 21
//13; // hatfldb problem
//14; // hatfldc problem
//15; // equilibrium combustion problem
#ifdef HAVE_LAPACK
//15; // Hock - Schittkowski modified problem 52
//16; // Schittkowski modified problem 235
//17; // Boggs & Tolle modified problem #7
//16; // Hock - Schittkowski modified #1 problem 52
//17; // Schittkowski modified problem 235
//18; // Boggs & Tolle modified problem #7
//19; // Hock - Schittkowski modified #2 problem 52
//20; // Hock - Schittkowski modified problem #76"
#endif /* HAVE_LAPACK */
switch(problem){
default: fprintf(stderr, "unknown problem specified (#%d)! Note that some minimization problems require LAPACK.\n", problem);
exit(1);
break;
 
case 0:
/* Rosenbrock function */
m=2; n=2;
798,10 → 931,27
for(i=0; i<n; ++i) printf("gradient %d, err %g\n", i, err[i]);
}
*/
break;
 
case 5:
/* Osborne's data fitting problem */
{
double x33[]={
8.44E-1, 9.08E-1, 9.32E-1, 9.36E-1, 9.25E-1, 9.08E-1, 8.81E-1,
8.5E-1, 8.18E-1, 7.84E-1, 7.51E-1, 7.18E-1, 6.85E-1, 6.58E-1,
6.28E-1, 6.03E-1, 5.8E-1, 5.58E-1, 5.38E-1, 5.22E-1, 5.06E-1,
4.9E-1, 4.78E-1, 4.67E-1, 4.57E-1, 4.48E-1, 4.38E-1, 4.31E-1,
4.24E-1, 4.2E-1, 4.14E-1, 4.11E-1, 4.06E-1};
 
m=5; n=33;
p[0]=0.5; p[1]=1.5; p[2]=-1.0; p[3]=1.0E-2; p[4]=2.0E-2;
 
ret=dlevmar_der(osborne, jacosborne, p, x33, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
//ret=dlevmar_dif(osborne, p, x33, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
}
break;
 
case 5:
case 6:
/* helical valley function */
m=3; n=3;
p[0]=-1.0; p[1]=0.0; p[2]=0.0;
811,7 → 961,7
break;
 
#ifdef HAVE_LAPACK
case 6:
case 7:
/* Boggs-Tolle problem 3 */
m=5; n=5;
p[0]=2.0; p[1]=2.0; p[2]=2.0;
826,7 → 976,8
//ret=dlevmar_lec_dif(bt3, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
}
break;
case 7:
 
case 8:
/* Hock - Schittkowski problem 28 */
m=3; n=3;
p[0]=-4.0; p[1]=1.0; p[2]=1.0;
840,7 → 991,8
//ret=dlevmar_lec_dif(hs28, p, x, m, n, A, b, 1, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
}
break;
case 8:
 
case 9:
/* Hock - Schittkowski problem 48 */
m=5; n=5;
p[0]=3.0; p[1]=5.0; p[2]=-3.0;
855,7 → 1007,8
//ret=dlevmar_lec_dif(hs48, p, x, m, n, A, b, 2, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
}
break;
case 9:
 
case 10:
/* Hock - Schittkowski problem 51 */
m=5; n=5;
p[0]=2.5; p[1]=0.5; p[2]=2.0;
870,9 → 1023,10
//ret=dlevmar_lec_dif(hs51, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
}
break;
 
#endif /* HAVE_LAPACK */
 
case 10:
case 11:
/* Hock - Schittkowski problem 01 */
m=2; n=2;
p[0]=-2.0; p[1]=1.0;
887,7 → 1041,8
ret=dlevmar_bc_der(hs01, jachs01, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
case 11:
 
case 12:
/* Hock - Schittkowski (modified) problem 21 */
m=2; n=2;
p[0]=-1.0; p[1]=-1.0;
902,7 → 1057,8
ret=dlevmar_bc_der(hs21, jachs21, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
case 12:
 
case 13:
/* hatfldb problem */
m=4; n=4;
p[0]=p[1]=p[2]=p[3]=0.1;
919,7 → 1075,8
ret=dlevmar_bc_der(hatfldb, jachatfldb, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
case 13:
 
case 14:
/* hatfldc problem */
m=4; n=4;
p[0]=p[1]=p[2]=p[3]=0.9;
935,7 → 1092,8
ret=dlevmar_bc_der(hatfldc, jachatfldc, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
case 14:
 
case 15:
/* equilibrium combustion problem */
m=5; n=5;
p[0]=p[1]=p[2]=p[3]=p[4]=0.0001;
951,9 → 1109,10
ret=dlevmar_bc_der(combust, jaccombust, p, x, m, n, lb, ub, 5000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
}
break;
 
#ifdef HAVE_LAPACK
case 15:
/* Hock - Schittkowski modified problem 52 */
case 16:
/* Hock - Schittkowski modified #1 problem 52 */
m=5; n=4;
p[0]=2.0; p[1]=2.0; p[2]=2.0;
p[3]=2.0; p[4]=2.0;
970,11 → 1129,12
lb[0]=-0.09; lb[1]=0.0; lb[2]=-DBL_MAX; lb[3]=-0.2; lb[4]=0.0;
ub[0]=DBL_MAX; ub[1]=0.3; ub[2]=0.25; ub[3]=0.3; ub[4]=0.3;
 
ret=dlevmar_blec_der(modhs52, jacmodhs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, analytic Jacobian
//ret=dlevmar_blec_dif(modhs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
ret=dlevmar_blec_der(mod1hs52, jacmod1hs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, analytic Jacobian
//ret=dlevmar_blec_dif(mod1hs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
}
break;
case 16:
 
case 17:
/* Schittkowski modified problem 235 */
m=3; n=2;
p[0]=-2.0; p[1]=3.0; p[2]=1.0;
993,7 → 1153,8
//ret=dlevmar_blec_dif(mods235, p, x, m, n, lb, ub, A, b, 2, NULL, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
}
break;
case 17:
 
case 18:
/* Boggs & Tolle modified problem 7 */
m=5; n=5;
p[0]=-2.0; p[1]=1.0; p[2]=1.0; p[3]=1.0; p[4]=1.0;
1012,6 → 1173,47
//ret=dlevmar_blec_dif(modbt7, p, x, m, n, lb, ub, A, b, 3, NULL, 10000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
}
break;
 
case 19:
/* Hock - Schittkowski modified #2 problem 52 */
m=5; n=5;
p[0]=2.0; p[1]=2.0; p[2]=2.0;
p[3]=2.0; p[4]=2.0;
for(i=0; i<n; i++) x[i]=0.0;
 
{
double C[3*5]={1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, -2.0, 0.0, -1.0, 0.0, 0.0, 1.0},
d[3]={-1.0, -2.0, -7.0};
 
ret=dlevmar_bleic_der(mod2hs52, jacmod2hs52, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, 3, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, analytic Jacobian
//ret=dlevmar_bleic_dif(mod2hs52, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, 3, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, no Jacobian
}
break;
 
case 20:
/* Hock - Schittkowski modified problem 76 */
m=4; n=4;
p[0]=0.5; p[1]=0.5; p[2]=0.5; p[3]=0.5;
for(i=0; i<n; i++) x[i]=0.0;
 
{
double A[1*4]={0.0, 1.0, 4.0, 0.0},
b[1]={1.5};
 
double C[2*4]={-1.0, -2.0, -1.0, -1.0, -3.0, -1.0, -2.0, 1.0},
d[2]={-5.0, -0.4};
 
double lb[4]={0.0, 0.0, 0.0, 0.0};
 
ret=dlevmar_bleic_der(modhs76, jacmodhs76, p, x, m, n, lb, NULL, A, b, 1, C, d, 2, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, analytic Jacobian
//ret=dlevmar_bleic_dif(modhs76, p, x, m, n, lb, NULL, A, b, 1, C, d, 2, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, no Jacobian
/* variations:
* if no lb is used, the minimizer is (-0.1135922 0.1330097 0.3417476 0.07572816)
* if the rhs of constr2 is 4.0, the minimizer is (0.0, 0.166667, 0.333333, 0.0)
*/
}
break;
 
#endif /* HAVE_LAPACK */
} /* switch */
/gpu/src/levmar/compiler.h
38,4 → 38,8
#define LM_FINITE finite // other than MSVC, ICC, GCC, let's hope this will work
#endif
 
#ifdef _MSC_VER // avoid deprecation warnings in VS2005
#define _CRT_SECURE_NO_WARNINGS
#endif
 
#endif /* _COMPILER_H_ */
/gpu/src/levmar/lmbleic_core.c New file
0,0 → 1,506
/////////////////////////////////////////////////////////////////////////////////
//
// Levenberg - Marquardt non-linear minimization algorithm
// Copyright (C) 2009 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 LM_REAL // not included by lmbleic.c
#error This file should not be compiled directly!
#endif
 
 
/* precision-specific definitions */
#define LMBLEIC_DATA LM_ADD_PREFIX(lmbleic_data)
#define LMBLEIC_ELIM LM_ADD_PREFIX(lmbleic_elim)
#define LMBLEIC_FUNC LM_ADD_PREFIX(lmbleic_func)
#define LMBLEIC_JACF LM_ADD_PREFIX(lmbleic_jacf)
#define LEVMAR_BLEIC_DER LM_ADD_PREFIX(levmar_bleic_der)
#define LEVMAR_BLEIC_DIF LM_ADD_PREFIX(levmar_bleic_dif)
#define LEVMAR_BLIC_DER LM_ADD_PREFIX(levmar_blic_der)
#define LEVMAR_BLIC_DIF LM_ADD_PREFIX(levmar_blic_dif)
#define LEVMAR_LEIC_DER LM_ADD_PREFIX(levmar_leic_der)
#define LEVMAR_LEIC_DIF LM_ADD_PREFIX(levmar_leic_dif)
#define LEVMAR_LIC_DER LM_ADD_PREFIX(levmar_lic_der)
#define LEVMAR_LIC_DIF LM_ADD_PREFIX(levmar_lic_dif)
#define LEVMAR_BLEC_DER LM_ADD_PREFIX(levmar_blec_der)
#define LEVMAR_BLEC_DIF LM_ADD_PREFIX(levmar_blec_dif)
#define LEVMAR_TRANS_MAT_MAT_MULT LM_ADD_PREFIX(levmar_trans_mat_mat_mult)
#define LEVMAR_COVAR LM_ADD_PREFIX(levmar_covar)
#define LEVMAR_FDIF_FORW_JAC_APPROX LM_ADD_PREFIX(levmar_fdif_forw_jac_approx)
 
struct LMBLEIC_DATA{
LM_REAL *jac;
int nineqcnstr; // #inequality constraints
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata);
void (*jacf)(LM_REAL *p, LM_REAL *jac, int m, int n, void *adata);
void *adata;
};
 
 
/* wrapper ensuring that the user-supplied function is called with the right number of variables (i.e. m) */
static void LMBLEIC_FUNC(LM_REAL *pext, LM_REAL *hx, int mm, int n, void *adata)
{
struct LMBLEIC_DATA *data=(struct LMBLEIC_DATA *)adata;
int m;
 
m=mm-data->nineqcnstr;
(*(data->func))(pext, hx, m, n, data->adata);
}
 
 
/* wrapper for computing the Jacobian at pext. The Jacobian is nxmm */
static void LMBLEIC_JACF(LM_REAL *pext, LM_REAL *jacext, int mm, int n, void *adata)
{
struct LMBLEIC_DATA *data=(struct LMBLEIC_DATA *)adata;
int m;
register int i, j;
LM_REAL *jac, *jacim, *jacextimm;
 
m=mm-data->nineqcnstr;
jac=data->jac;
 
(*(data->jacf))(pext, jac, m, n, data->adata);
 
for(i=0; i<n; ++i){
jacextimm=jacext+i*mm;
jacim=jac+i*m;
for(j=0; j<m; ++j)
jacextimm[j]=jacim[j]; //jacext[i*mm+j]=jac[i*m+j];
 
for(j=m; j<mm; ++j)
jacextimm[j]=0.0; //jacext[i*mm+j]=0.0;
}
}
 
 
/*
* This function is similar to LEVMAR_DER except that the minimization is
* performed subject to the box constraints lb[i]<=p[i]<=ub[i], the linear
* equation constraints A*p=b, A being k1xm, b k1x1, and the linear inequality
* constraints C*p>=d, C being k2xm, d k2x1.
*
* The inequalities are converted to equations by introducing surplus variables,
* i.e. c^T*p >= d becomes c^T*p - y = d, with y>=0. To transform all inequalities
* to equations, a total of k2 surplus variables are introduced; a problem with only
* box and linear constraints results then and is solved with LEVMAR_BLEC_DER()
* Note that opposite direction inequalities should be converted to the desired
* direction by negating, i.e. c^T*p <= d becomes -c^T*p >= -d
*
* This function requires an analytic Jacobian. In case the latter is unavailable,
* use LEVMAR_BLEIC_DIF() bellow
*
*/
int LEVMAR_BLEIC_DER(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata), /* functional relation describing measurements. A p \in R^m yields a \hat{x} \in R^n */
void (*jacf)(LM_REAL *p, LM_REAL *j, int m, int n, void *adata), /* function to evaluate the Jacobian \part x / \part p */
LM_REAL *p, /* I/O: initial parameter estimates. On output has the estimated solution */
LM_REAL *x, /* I: measurement vector. NULL implies a zero vector */
int m, /* I: parameter vector dimension (i.e. #unknowns) */
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 *A, /* I: equality constraints matrix, k1xm. If NULL, no linear equation constraints apply */
LM_REAL *b, /* I: right hand constraints vector, k1x1 */
int k1, /* I: number of constraints (i.e. A's #rows) */
LM_REAL *C, /* I: inequality constraints matrix, k2xm */
LM_REAL *d, /* I: right hand constraints vector, k2x1 */
int k2, /* I: number of inequality constraints (i.e. C's #rows) */
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
*/
LM_REAL info[LM_INFO_SZ],
/* O: information regarding the minimization. Set to NULL if don't care
* info[0]= ||e||_2 at initial p.
* info[1-4]=[ ||e||_2, ||J^T e||_inf, ||Dp||_2, mu/max[J^T J]_ii ], all computed at estimated p.
* info[5]= # iterations,
* info[6]=reason for terminating: 1 - stopped by small gradient J^T e
* 2 - stopped by small Dp
* 3 - stopped by itmax
* 4 - singular matrix. Restart from current p with increased mu
* 5 - no further error reduction is possible. Restart with increased mu
* 6 - stopped by small ||e||_2
* 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error
* info[7]= # function evaluations
* info[8]= # Jacobian evaluations
* info[9]= # linear systems solved, i.e. # attempts for reducing error
*/
LM_REAL *work, /* working memory at least LM_BLEIC_DER_WORKSZ() reals large, allocated if NULL */
LM_REAL *covar, /* O: Covariance matrix corresponding to LS solution; mxm. Set to NULL if not needed. */
void *adata) /* pointer to possibly additional data, passed uninterpreted to func & jacf.
* Set to NULL if not needed
*/
{
struct LMBLEIC_DATA data;
LM_REAL *ptr, *pext, *Aext, *bext, *covext; /* corresponding to p, A, b, covar for the full set of variables;
pext=[p, surplus], pext is mm, Aext is (k1+k2)xmm, bext (k1+k2), covext is mmxmm
*/
LM_REAL *lbext, *ubext; // corresponding to lb, ub for the full set of variables
int mm, ret, k12;
register int i, j, ii;
register LM_REAL tmp;
LM_REAL locinfo[LM_INFO_SZ];
 
if(!jacf){
fprintf(stderr, RCAT("No function specified for computing the Jacobian in ", LEVMAR_BLEIC_DER)
RCAT("().\nIf no such function is available, use ", LEVMAR_BLEIC_DIF) RCAT("() rather than ", LEVMAR_BLEIC_DER) "()\n");
return LM_ERROR;
}
 
if(!C || !d){
fprintf(stderr, RCAT(LCAT(LEVMAR_BLEIC_DER, "(): missing inequality constraints, use "), LEVMAR_BLEC_DER) "() in this case!\n");
return LM_ERROR;
}
 
if(!A || !b) k1=0; // sanity check
 
mm=m+k2;
 
if(n<m-k1){
fprintf(stderr, LCAT(LEVMAR_BLEIC_DER, "(): cannot solve a problem with fewer measurements + equality constraints [%d + %d] than unknowns [%d]\n"), n, k1, m);
return LM_ERROR;
}
 
k12=k1+k2;
ptr=(LM_REAL *)malloc((3*mm + k12*mm + k12 + n*m + (covar? mm*mm : 0))*sizeof(LM_REAL));
if(!ptr){
fprintf(stderr, LCAT(LEVMAR_BLEIC_DER, "(): memory allocation request failed\n"));
return LM_ERROR;
}
pext=ptr;
lbext=pext+mm;
ubext=lbext+mm;
Aext=ubext+mm;
bext=Aext+k12*mm;
data.jac=bext+k12;
covext=covar? data.jac+n*m : NULL;
data.nineqcnstr=k2;
data.func=func;
data.jacf=jacf;
data.adata=adata;
 
/* compute y s.t. C*p - y=d, i.e. y=C*p-d.
* y is stored in the last k2 elements of pext
*/
for(i=0; i<k2; ++i){
for(j=0, tmp=0.0; j<m; ++j)
tmp+=C[i*m+j]*p[j];
pext[j=i+m]=tmp-d[i];
 
/* surplus variables must be >=0 */
lbext[j]=0.0;
ubext[j]=LM_REAL_MAX;
}
/* set the first m elements of pext equal to p */
for(i=0; i<m; ++i){
pext[i]=p[i];
lbext[i]=lb? lb[i] : LM_REAL_MIN;
ubext[i]=ub? ub[i] : LM_REAL_MAX;
}
 
/* setup the constraints matrix */
/* original linear equation constraints */
for(i=0; i<k1; ++i){
for(j=0; j<m; ++j)
Aext[i*mm+j]=A[i*m+j];
 
for(j=m; j<mm; ++j)
Aext[i*mm+j]=0.0;
 
bext[i]=b[i];
}
/* linear equation constraints resulting from surplus variables */
for(i=0, ii=k1; i<k2; ++i, ++ii){
for(j=0; j<m; ++j)
Aext[ii*mm+j]=C[i*m+j];
 
for(j=m; j<mm; ++j)
Aext[ii*mm+j]=0.0;
 
Aext[ii*mm+m+i]=-1.0;
 
bext[ii]=d[i];
}
 
if(!info) info=locinfo; /* make sure that LEVMAR_BLEC_DER() is called with non-null info */
/* note that the default weights for the penalty terms are being used below */
ret=LEVMAR_BLEC_DER(LMBLEIC_FUNC, LMBLEIC_JACF, pext, x, mm, n, lbext, ubext, Aext, bext, k12, NULL, itmax, opts, info, work, covext, (void *)&data);
 
/* copy back the minimizer */
for(i=0; i<m; ++i)
p[i]=pext[i];
 
#if 0
printf("Surplus variables for the minimizer:\n");
for(i=m; i<mm; ++i)
printf("%g ", pext[i]);
printf("\n\n");
#endif
 
if(covar){
for(i=0; i<m; ++i){
for(j=0; j<m; ++j)
covar[i*m+j]=covext[i*mm+j];
}
}
 
free(ptr);
 
return ret;
}
 
/* Similar to the LEVMAR_BLEIC_DER() function above, except that the Jacobian is approximated
* with the aid of finite differences (forward or central, see the comment for the opts argument)
*/
int LEVMAR_BLEIC_DIF(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata), /* functional relation describing measurements. A p \in R^m yields a \hat{x} \in R^n */
LM_REAL *p, /* I/O: initial parameter estimates. On output has the estimated solution */
LM_REAL *x, /* I: measurement vector. NULL implies a zero vector */
int m, /* I: parameter vector dimension (i.e. #unknowns) */
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 *A, /* I: equality constraints matrix, k1xm. If NULL, no linear equation constraints apply */
LM_REAL *b, /* I: right hand constraints vector, k1x1 */
int k1, /* I: number of constraints (i.e. A's #rows) */
LM_REAL *C, /* I: inequality constraints matrix, k2xm */
LM_REAL *d, /* I: right hand constraints vector, k2x1 */
int k2, /* I: number of inequality constraints (i.e. C's #rows) */
int itmax, /* I: maximum number of iterations */
LM_REAL opts[5], /* I: opts[0-3] = 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
* the step used in difference approximation to the Jacobian. Set to NULL for defaults to be used.
* If \delta<0, the Jacobian is approximated with central differences which are more accurate
* (but slower!) compared to the forward differences employed by default.
*/
LM_REAL info[LM_INFO_SZ],
/* O: information regarding the minimization. Set to NULL if don't care
* info[0]= ||e||_2 at initial p.
* info[1-4]=[ ||e||_2, ||J^T e||_inf, ||Dp||_2, mu/max[J^T J]_ii ], all computed at estimated p.
* info[5]= # iterations,
* info[6]=reason for terminating: 1 - stopped by small gradient J^T e
* 2 - stopped by small Dp
* 3 - stopped by itmax
* 4 - singular matrix. Restart from current p with increased mu
* 5 - no further error reduction is possible. Restart with increased mu
* 6 - stopped by small ||e||_2
* 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error
* info[7]= # function evaluations
* info[8]= # Jacobian evaluations
* info[9]= # linear systems solved, i.e. # attempts for reducing error
*/
LM_REAL *work, /* working memory at least LM_BLEIC_DIF_WORKSZ() reals large, allocated if NULL */
LM_REAL *covar, /* O: Covariance matrix corresponding to LS solution; mxm. Set to NULL if not needed. */
void *adata) /* pointer to possibly additional data, passed uninterpreted to func.
* Set to NULL if not needed
*/
{
struct LMBLEIC_DATA data;
LM_REAL *ptr, *pext, *Aext, *bext, *covext; /* corresponding to p, A, b, covar for the full set of variables;
pext=[p, surplus], pext is mm, Aext is (k1+k2)xmm, bext (k1+k2), covext is mmxmm
*/
LM_REAL *lbext, *ubext; // corresponding to lb, ub for the full set of variables
int mm, ret, k12;
register int i, j, ii;
register LM_REAL tmp;
LM_REAL locinfo[LM_INFO_SZ];
 
if(!C || !d){
fprintf(stderr, RCAT(LCAT(LEVMAR_BLEIC_DIF, "(): missing inequality constraints, use "), LEVMAR_BLEC_DIF) "() in this case!\n");
return LM_ERROR;
}
if(!A || !b) k1=0; // sanity check
 
mm=m+k2;
 
if(n<m-k1){
fprintf(stderr, LCAT(LEVMAR_BLEIC_DIF, "(): cannot solve a problem with fewer measurements + equality constraints [%d + %d] than unknowns [%d]\n"), n, k1, m);
return LM_ERROR;
}
 
k12=k1+k2;
ptr=(LM_REAL *)malloc((3*mm + k12*mm + k12 + (covar? mm*mm : 0))*sizeof(LM_REAL));
if(!ptr){
fprintf(stderr, LCAT(LEVMAR_BLEIC_DIF, "(): memory allocation request failed\n"));
return LM_ERROR;
}
pext=ptr;
lbext=pext+mm;
ubext=lbext+mm;
Aext=ubext+mm;
bext=Aext+k12*mm;
data.jac=NULL;
covext=covar? bext+k12 : NULL;
data.nineqcnstr=k2;
data.func=func;
data.jacf=NULL;
data.adata=adata;
 
/* compute y s.t. C*p - y=d, i.e. y=C*p-d.
* y is stored in the last k2 elements of pext
*/
for(i=0; i<k2; ++i){
for(j=0, tmp=0.0; j<m; ++j)
tmp+=C[i*m+j]*p[j];
pext[j=i+m]=tmp-d[i];
 
/* surplus variables must be >=0 */
lbext[j]=0.0;
ubext[j]=LM_REAL_MAX;
}
/* set the first m elements of pext equal to p */
for(i=0; i<m; ++i){
pext[i]=p[i];
lbext[i]=lb? lb[i] : LM_REAL_MIN;
ubext[i]=ub? ub[i] : LM_REAL_MAX;
}
 
/* setup the constraints matrix */
/* original linear equation constraints */
for(i=0; i<k1; ++i){
for(j=0; j<m; ++j)
Aext[i*mm+j]=A[i*m+j];
 
for(j=m; j<mm; ++j)
Aext[i*mm+j]=0.0;
 
bext[i]=b[i];
}
/* linear equation constraints resulting from surplus variables */
for(i=0, ii=k1; i<k2; ++i, ++ii){
for(j=0; j<m; ++j)
Aext[ii*mm+j]=C[i*m+j];
 
for(j=m; j<mm; ++j)
Aext[ii*mm+j]=0.0;
 
Aext[ii*mm+m+i]=-1.0;
 
bext[ii]=d[i];
}
 
if(!info) info=locinfo; /* make sure that LEVMAR_BLEC_DIF() is called with non-null info */
/* note that the default weights for the penalty terms are being used below */
ret=LEVMAR_BLEC_DIF(LMBLEIC_FUNC, pext, x, mm, n, lbext, ubext, Aext, bext, k12, NULL, itmax, opts, info, work, covext, (void *)&data);
 
/* copy back the minimizer */
for(i=0; i<m; ++i)
p[i]=pext[i];
 
#if 0
printf("Surplus variables for the minimizer:\n");
for(i=m; i<mm; ++i)
printf("%g ", pext[i]);
printf("\n\n");
#endif
 
if(covar){
for(i=0; i<m; ++i){
for(j=0; j<m; ++j)
covar[i*m+j]=covext[i*mm+j];
}
}
 
free(ptr);
 
return ret;
}
 
 
/* convenience wrappers to LEVMAR_BLEIC_DER/LEVMAR_BLEIC_DIF */
 
/* box & linear inequality constraints */
int LEVMAR_BLIC_DER(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
void (*jacf)(LM_REAL *p, LM_REAL *j, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *lb, LM_REAL *ub,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[4], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DER(func, jacf, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, k2, itmax, opts, info, work, covar, adata);
}
 
int LEVMAR_BLIC_DIF(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *lb, LM_REAL *ub,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[5], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DIF(func, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, k2, itmax, opts, info, work, covar, adata);
}
 
/* linear equation & inequality constraints */
int LEVMAR_LEIC_DER(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
void (*jacf)(LM_REAL *p, LM_REAL *j, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *A, LM_REAL *b, int k1,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[4], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DER(func, jacf, p, x, m, n, NULL, NULL, A, b, k1, C, d, k2, itmax, opts, info, work, covar, adata);
}
 
int LEVMAR_LEIC_DIF(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *A, LM_REAL *b, int k1,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[5], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DIF(func, p, x, m, n, NULL, NULL, A, b, k1, C, d, k2, itmax, opts, info, work, covar, adata);
}
 
/* linear inequality constraints */
int LEVMAR_LIC_DER(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
void (*jacf)(LM_REAL *p, LM_REAL *j, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[4], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DER(func, jacf, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, k2, itmax, opts, info, work, covar, adata);
}
 
int LEVMAR_LIC_DIF(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[5], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DIF(func, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, k2, itmax, opts, info, work, covar, adata);
}
 
/* undefine all. THIS MUST REMAIN AT THE END OF THE FILE */
#undef LMBLEIC_DATA
#undef LMBLEIC_ELIM
#undef LMBLEIC_FUNC
#undef LMBLEIC_JACF
#undef LEVMAR_FDIF_FORW_JAC_APPROX
#undef LEVMAR_COVAR
#undef LEVMAR_TRANS_MAT_MAT_MULT
#undef LEVMAR_BLEIC_DER
#undef LEVMAR_BLEIC_DIF
#undef LEVMAR_BLIC_DER
#undef LEVMAR_BLIC_DIF
#undef LEVMAR_LEIC_DER
#undef LEVMAR_LEIC_DIF
#undef LEVMAR_LIC_DER
#undef LEVMAR_LIC_DIF
#undef LEVMAR_BLEC_DER
#undef LEVMAR_BLEC_DIF
/gpu/src/levmar/matlab/jacmodhs52.m
1,4 → 1,4
function jac = modhs52_jac(p)
function jac = jacmodhs52(p)
m=5;
 
jac(1, 1:m)=[4.0, -1.0, 0.0, 0.0, 0.0];
/gpu/src/levmar/matlab/jacmeyer.m
1,4 → 1,4
function jac = meyer_jac(p, data1, data2)
function jac = jacmeyer(p, data1, data2)
n=16;
m=3;
 
/gpu/src/levmar/matlab/jacmodhs76.m New file
0,0 → 1,7
function jac = jacmodhs76(p)
m=4;
 
jac(1, 1:m)=[1.0, 0.0, 0.0, 0.0];
jac(2, 1:m)=[0.0, sqrt(0.5), 0.0, 0.0];
jac(3, 1:m)=[0.0, 0.0, 1.0, 0.0];
jac(4, 1:m)=[0.0, 0.0, 0.0, sqrt(0.5)];
/gpu/src/levmar/matlab/jacbt3.m
1,4 → 1,4
function jac = bt3_jac(p, adata)
function jac = jacbt3(p, adata)
n=5;
m=5;
 
/gpu/src/levmar/matlab/osborne.m New file
0,0 → 1,7
function x = osborne(p)
n=33;
 
for i=1:n
t=10*(i-1);
x(i)=p(1) + p(2)*exp(-p(4)*t) + p(3)*exp(-p(5)*t);
end
/gpu/src/levmar/matlab/levmar.c
24,7 → 24,7
#include <string.h>
#include <ctype.h>
 
#include <lm.h>
#include <levmar.h>
 
#include <mex.h>
 
46,6 → 46,10
#define MIN_CONSTRAINED_BC 1
#define MIN_CONSTRAINED_LEC 2
#define MIN_CONSTRAINED_BLEC 3
#define MIN_CONSTRAINED_BLEIC 4
#define MIN_CONSTRAINED_BLIC 5
#define MIN_CONSTRAINED_LEIC 6
#define MIN_CONSTRAINED_LIC 7
 
struct mexdata {
/* matlab names of the fitting function & its Jacobian */
68,9 → 72,9
char buf[256];
va_list args;
 
va_start(args, fmt);
vsprintf(buf, fmt, args);
va_end(args);
va_start(args, fmt);
vsprintf(buf, fmt, args);
va_end(args);
 
mexErrMsgTxt(buf);
}
81,9 → 85,9
char buf[256];
va_list args;
 
va_start(args, fmt);
vsprintf(buf, fmt, args);
va_end(args);
va_start(args, fmt);
vsprintf(buf, fmt, args);
va_end(args);
 
mexWarnMsgTxt(buf);
}
183,7 → 187,7
else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || !(mxGetM(lhs[0])==1 || mxGetN(lhs[0])==1) ||
__MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))!=n){
fprintf(stderr, "levmar: '%s' should produce a real vector with %d elements (got %d).\n",
dat->fname, m, __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0])));
dat->fname, n, __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0])));
ret=1;
}
/* delete the matrix created by matlab */
220,20 → 224,26
[ret, p, info, covar]=levmar_bc (f, j, p0, x, itmax, opts, 'bc', lb, ub, ...)
[ret, p, info, covar]=levmar_lec (f, j, p0, x, itmax, opts, 'lec', A, b, ...)
[ret, p, info, covar]=levmar_blec(f, j, p0, x, itmax, opts, 'blec', lb, ub, A, b, wghts, ...)
 
[ret, p, info, covar]=levmar_bleic(f, j, p0, x, itmax, opts, 'bleic', lb, ub, A, b, C, d, ...)
[ret, p, info, covar]=levmar_blic (f, j, p0, x, itmax, opts, 'blic', lb, ub, C, d, ...)
[ret, p, info, covar]=levmar_leic (f, j, p0, x, itmax, opts, 'leic', A, b, C, d, ...)
[ret, p, info, covar]=levmar_lic (f, j, p0, x, itmax, opts, 'lic', C, d, ...)
 
*/
 
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *Prhs[])
{
register int i;
register double *pdbl;
mxArray **prhs=(mxArray **)&Prhs[0], *At;
mxArray **prhs=(mxArray **)&Prhs[0], *At, *Ct;
struct mexdata mdata;
int len, status;
double *p, *p0, *ret, *x;
int m, n, havejac, Arows, itmax, nopts, mintype, nextra;
int m, n, havejac, Arows, Crows, itmax, nopts, mintype, nextra;
double opts[LM_OPTS_SZ]={LM_INIT_MU, LM_STOP_THRESH, LM_STOP_THRESH, LM_STOP_THRESH, LM_DIFF_DELTA};
double info[LM_INFO_SZ];
double *lb=NULL, *ub=NULL, *A=NULL, *b=NULL, *wghts=NULL, *covar=NULL;
double *lb=NULL, *ub=NULL, *A=NULL, *b=NULL, *wghts=NULL, *C=NULL, *d=NULL, *covar=NULL;
 
/* parse input args; start by checking their number */
if((nrhs<5))
362,6 → 372,10
else if(!strncmp(minhowto, "bc", 2)) mintype=MIN_CONSTRAINED_BC;
else if(!strncmp(minhowto, "lec", 3)) mintype=MIN_CONSTRAINED_LEC;
else if(!strncmp(minhowto, "blec", 4)) mintype=MIN_CONSTRAINED_BLEC;
else if(!strncmp(minhowto, "bleic", 5)) mintype=MIN_CONSTRAINED_BLEIC;
else if(!strncmp(minhowto, "blic", 4)) mintype=MIN_CONSTRAINED_BLIC;
else if(!strncmp(minhowto, "leic", 4)) mintype=MIN_CONSTRAINED_LEIC;
else if(!strncmp(minhowto, "lic", 3)) mintype=MIN_CONSTRAINED_BLIC;
else matlabFmtdErrMsgTxt("levmar: unknown minimization type '%s'.", minhowto);
 
mxFree(minhowto);
378,7 → 392,7
* upon the minimization type determined above
*/
/** lb, ub **/
if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BC || mintype==MIN_CONSTRAINED_BLEC)){
if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BC || mintype==MIN_CONSTRAINED_BLEC || mintype==MIN_CONSTRAINED_BLIC || mintype==MIN_CONSTRAINED_BLEIC)){
/* check if the next two arguments are real row or column vectors */
if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){
if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
397,7 → 411,7
}
 
/** A, b **/
if(nrhs>=7 && (mintype==MIN_CONSTRAINED_LEC || mintype==MIN_CONSTRAINED_BLEC)){
if(nrhs>=7 && (mintype==MIN_CONSTRAINED_LEC || mintype==MIN_CONSTRAINED_BLEC || mintype==MIN_CONSTRAINED_LEIC || mintype==MIN_CONSTRAINED_BLEIC)){
/* check if the next two arguments are a real matrix and a real row or column vector */
if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && mxGetM(prhs[5])>=1 && mxGetN(prhs[5])>=1){
if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
428,6 → 442,27
}
}
}
 
/** C, d **/
if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BLEIC || mintype==MIN_CONSTRAINED_BLIC || mintype==MIN_CONSTRAINED_LEIC || mintype==MIN_CONSTRAINED_LIC)){
/* check if the next two arguments are a real matrix and a real row or column vector */
if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && mxGetM(prhs[5])>=1 && mxGetN(prhs[5])>=1){
if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
if((i=mxGetN(prhs[5]))!=m)
matlabFmtdErrMsgTxt("levmar: C must have %d columns, got %d.", m, i);
if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=(Crows=mxGetM(prhs[5])))
matlabFmtdErrMsgTxt("levmar: d must have %d elements, got %d.", Crows, i);
 
Ct=prhs[5];
d=mxGetPr(prhs[6]);
C=getTranspose(Ct);
 
prhs+=2;
nrhs-=2;
}
}
}
 
/* arguments below this point are assumed to be extra arguments passed
* to every invocation of the fitting function and its Jacobian
*/
471,8 → 506,8
covar=mxMalloc(m*m*sizeof(double));
 
/* invoke levmar */
if(!lb && !ub){
if(!A && !b){ /* no constraints */
switch(mintype){
case MIN_UNCONSTRAINED: /* no constraints */
if(havejac)
status=dlevmar_der(func, jacfunc, p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata);
else
481,8 → 516,18
fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_der()/dlevmar_dif()\n");
#endif /* DEBUG */
}
else{ /* linear constraints */
break;
case MIN_CONSTRAINED_BC: /* box constraints */
if(havejac)
status=dlevmar_bc_der(func, jacfunc, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata);
else
status=dlevmar_bc_dif(func, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata);
#ifdef DEBUG
fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_bc_der()/dlevmar_bc_dif()\n");
#endif /* DEBUG */
break;
case MIN_CONSTRAINED_LEC: /* linear equation constraints */
#ifdef HAVE_LAPACK
if(havejac)
status=dlevmar_lec_der(func, jacfunc, p, x, m, n, A, b, Arows, itmax, opts, info, NULL, covar, (void *)&mdata);
496,35 → 541,86
fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_lec_der()/dlevmar_lec_dif()\n");
#endif /* DEBUG */
}
}
else{
if(!A && !b){ /* box constraints */
break;
case MIN_CONSTRAINED_BLEC: /* box & linear equation constraints */
#ifdef HAVE_LAPACK
if(havejac)
status=dlevmar_bc_der(func, jacfunc, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata);
status=dlevmar_blec_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
else
status=dlevmar_bc_dif(func, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata);
status=dlevmar_blec_dif(func, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
#else
mexErrMsgTxt("levmar: no box & linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
#endif /* HAVE_LAPACK */
 
#ifdef DEBUG
fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_bc_der()/dlevmar_bc_dif()\n");
fprintf(stderr, "LEVMAR: calling dlevmar_blec_der()/dlevmar_blec_dif()\n");
#endif /* DEBUG */
}
else{ /* box & linear constraints */
break;
case MIN_CONSTRAINED_BLEIC: /* box, linear equation & inequalities constraints */
#ifdef HAVE_LAPACK
if(havejac)
status=dlevmar_blec_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
else
status=dlevmar_blec_dif(func, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
status=dlevmar_bleic_dif(func, p, x, m, n, lb, ub, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
#else
mexErrMsgTxt("levmar: no box & linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
mexErrMsgTxt("levmar: no box, linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
#endif /* HAVE_LAPACK */
 
#ifdef DEBUG
fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_blec_der()/dlevmar_blec_dif()\n");
fprintf(stderr, "LEVMAR: calling dlevmar_bleic_der()/dlevmar_bleic_dif()\n");
#endif /* DEBUG */
}
break;
case MIN_CONSTRAINED_BLIC: /* box, linear inequalities constraints */
#ifdef HAVE_LAPACK
if(havejac)
status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
else
status=dlevmar_bleic_dif(func, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
#else
mexErrMsgTxt("levmar: no box & linear inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
#endif /* HAVE_LAPACK */
 
#ifdef DEBUG
fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_blic_der()/dlevmar_blic_dif()\n");
#endif /* DEBUG */
break;
case MIN_CONSTRAINED_LEIC: /* linear equation & inequalities constraints */
#ifdef HAVE_LAPACK
if(havejac)
status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, NULL, NULL, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
else
status=dlevmar_bleic_dif(func, p, x, m, n, NULL, NULL, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
#else
mexErrMsgTxt("levmar: no linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
#endif /* HAVE_LAPACK */
 
#ifdef DEBUG
fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_leic_der()/dlevmar_leic_dif()\n");
#endif /* DEBUG */
break;
case MIN_CONSTRAINED_LIC: /* linear inequalities constraints */
#ifdef HAVE_LAPACK
if(havejac)
status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
else
status=dlevmar_bleic_dif(func, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
#else
mexErrMsgTxt("levmar: no linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
#endif /* HAVE_LAPACK */
 
#ifdef DEBUG
fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_lic_der()/dlevmar_lic_dif()\n");
#endif /* DEBUG */
break;
default:
mexErrMsgTxt("levmar: unexpected internal error.");
}
 
#ifdef DEBUG
fflush(stderr);
printf("LEVMAR: minimization returned %d in %g iter, reason %g\n\tSolution: ", status, info[5], info[6]);
568,6 → 664,7
/* cleanup */
mxDestroyArray(mdata.rhs[0]);
if(A) mxFree(A);
if(C) mxFree(C);
 
mxFree(mdata.fname);
if(havejac) mxFree(mdata.jacname);
/gpu/src/levmar/matlab/lmdemo.m
1,3 → 1,8
% Demo program for levmar's MEX-file interface
% Performs minimization of several test problems
 
format long;
 
% Unconstrained minimization
 
% fitting the exponential model x_i=p(1)*exp(-p(2)*i)+p(3) of expfit.c to noisy measurements obtained with (5.0 0.1 1.0)
38,6 → 43,24
popt
 
 
% Osborne's problem
p0=[0.5, 1.5, -1.0, 1.0E-2, 2.0E-2];
 
x=[8.44E-1, 9.08E-1, 9.32E-1, 9.36E-1, 9.25E-1, 9.08E-1, 8.81E-1,...
8.5E-1, 8.18E-1, 7.84E-1, 7.51E-1, 7.18E-1, 6.85E-1, 6.58E-1,...
6.28E-1, 6.03E-1, 5.8E-1, 5.58E-1, 5.38E-1, 5.22E-1, 5.06E-1,...
4.9E-1, 4.78E-1, 4.67E-1, 4.57E-1, 4.48E-1, 4.38E-1, 4.31E-1,...
4.24E-1, 4.2E-1, 4.14E-1, 4.11E-1, 4.06E-1];
 
 
options=[1E-03, 1E-15, 1E-15, 1E-20, 1E-06];
 
[ret, popt, info, covar]=levmar('osborne', 'jacosborne', p0, x, 200, options);
%[ret, popt, info, covar]=levmar('osborne', p0, x, 200, options, 'unc');
disp('Osborne''s problem');
popt
 
 
% Linear constraints
 
% Boggs-Tolle problem 3
72,7 → 95,7
 
% Box and linear constraints
 
% Hock-Schittkowski modified problem 52
% Hock-Schittkowski modified problem 52 (#1)
p0=[2.0, 2.0, 2.0, 2.0, 2.0];
x=[0.0, 0.0, 0.0, 0.0];
lb=[-0.09, 0.0, -realmax, -0.2, 0.0];
84,10 → 107,10
options=[1E-03, 1E-15, 1E-15, 1E-20];
 
[ret, popt, info, covar]=levmar('modhs52', 'jacmodhs52', p0, x, 200, options, 'blec', lb, ub, A, b);
disp('Hock-Schittkowski modified problem 52');
disp('Hock-Schittkowski modified problem 52 (#1)');
popt
 
% Hock-Schittkowski modified problem 235
% Schittkowski modified problem 235
p0=[-2.0, 3.0, 1.0];
x=[0.0, 0.0];
lb=[-realmax, 0.1, 0.7];
101,3 → 124,17
disp('Hock-Schittkowski modified problem 235');
popt
 
% Box, linear equation & inequality constraints
p0=[0.5, 0.5, 0.5, 0.5];
x=[0.0, 0.0, 0.0, 0.0];
lb=[0.0, 0.0, 0.0, 0.0];
ub=[realmax, realmax, realmax, realmax];
A=[0.0, 1.0, 4.0, 0.0];
b=[1.5]';
C=[-1.0, -2.0, -1.0, -1.0;
-3.0, -1.0, -2.0, 1.0];
d=[-5.0, -0.4]';
 
[ret, popt, info, covar]=levmar('modhs76', 'jacmodhs76', p0, x, 200, options, 'bleic', lb, ub, A, b, C, d);
disp('Hock-Schittkowski modified problem 76');
popt
/gpu/src/levmar/matlab/jacexpfit.m
1,4 → 1,4
function jac = expfit_jac(p, data)
function jac = jacexpfit(p, data)
n=data;
m=max(size(p));
 
/gpu/src/levmar/matlab/README.txt
1,5 → 1,6
This directory contains a matlab MEX interface to levmar. This interface
has been tested with Matlab v. 6.5 R13 under linux and v. 7.4 R2007 under Windows.
Users have also reported success with GNU Octave.
 
FILES
The following files are included:
/gpu/src/levmar/matlab/CMakeLists.txt New file
0,0 → 1,58
# CMake file for levmar's MEX-file; see http://www.cmake.org
# Requires FindMatlab.cmake included with cmake
 
PROJECT(LEVMARMEX)
#CMAKE_MINIMUM_REQUIRED(VERSION 1.4)
 
INCLUDE("C:/Program Files/CMake 2.4/share/cmake-2.4/Modules/FindMatlab.cmake")
 
# f2c is sometimes equivalent to libF77 & libI77; in that case, set HAVE_F2C to 0
SET(HAVE_F2C 1 CACHE BOOL "Do we have f2c or F77/I77?" )
 
# the directory where the lapack/blas/f2c libraries reside
SET(LAPACKBLAS_DIR /usr/lib CACHE PATH "Path to lapack/blas libraries")
 
# the directory where levmar.h resides
SET(LM_H_DIR .. CACHE PATH "Path to levmar.h")
# the directory where the levmar library resides
SET(LEVMAR_DIR .. CACHE PATH "Path to levmar library")
 
# actual names for the lapack/blas/f2c libraries
SET(LAPACK_LIB lapack CACHE STRING "The name of the lapack library")
SET(BLAS_LIB blas CACHE STRING "The name of the blas library")
IF(HAVE_F2C)
SET(F2C_LIB f2c CACHE STRING "The name of the f2c library")
ELSE(HAVE_F2C)
SET(F77_LIB libF77 CACHE STRING "The name of the F77 library")
SET(I77_LIB libI77 CACHE STRING "The name of the I77 library")
ENDIF(HAVE_F2C)
 
########################## NO CHANGES BEYOND THIS POINT ##########################
 
INCLUDE_DIRECTORIES(${LM_H_DIR})
LINK_DIRECTORIES(${LAPACKBLAS_DIR} ${LEVMAR_DIR})
 
SET(SRC levmar.c)
 
# naming conventions for the generated file's suffix
IF(WIN32)
SET(SUFFIX ".mexw32")
ELSE(WIN32)
SET(SUFFIX ".mexglx")
ENDIF(WIN32)
 
SET(OUTNAME "levmar${SUFFIX}")
 
ADD_LIBRARY(${OUTNAME} MODULE ${SRC})
 
IF(HAVE_F2C)
ADD_CUSTOM_COMMAND(OUTPUT ${OUTNAME}
DEPENDS ${SRC}
COMMAND mex
ARGS -O -I${LM_H_DIR} ${SRC} -I${MATLAB_INCLUDE_DIR} -L${LAPACKBLAS_DIR} -L${LEVMAR_DIR} -L${MATLAB_MEX_LIBRARY} -llevmar -l${LAPACK_LIB} -l${BLAS_LIB} -l${F2C_LIB} -output ${MATLAB_LIBRARIES} ${OUTNAME})
ELSE(HAVE_F2C)
ADD_CUSTOM_COMMAND(OUTPUT ${OUTNAME}
DEPENDS ${SRC}
COMMAND mex
ARGS -O -I${LM_H_DIR} ${SRC} -I${MATLAB_INCLUDE_DIR} -L${LAPACKBLAS_DIR} -L${LEVMAR_DIR} -L${MATLAB_MEX_LIBRARY} -llevmar -l${LAPACK_LIB} -l${BLAS_LIB} -l${F77_LIB} -l${I77_LIB} ${MATLAB_LIBRARIES} -output ${OUTNAME})
ENDIF(HAVE_F2C)
/gpu/src/levmar/matlab/jacosborne.m New file
0,0 → 1,11
function jac = jacosborne(p)
n=33;
m=5;
 
for i=1:n
t=10*(i-1);
tmp1=exp(-p(4)*t);
tmp2=exp(-p(5)*t);
 
jac(i, 1:m)=[1.0, tmp1, tmp2, -p(2)*t*tmp1, -p(3)*t*tmp2];
end
/gpu/src/levmar/matlab/levmar.m
2,11 → 2,17
% LEVMAR matlab MEX interface to the levmar non-linear least squares minimization
% library available from http://www.ics.forth.gr/~lourakis/levmar/
%
% levmar can be used in any of the 4 following ways:
% levmar can be used in any of the 8 following ways:
% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'unc', ...)
% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'bc', lb, ub, ...)
% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'lec', A, b, ...)
% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'blec', lb, ub, A, b, wghts, ...)
%
% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'bleic', lb, ub, A, b, C, d, ...)
% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'blic', lb, ub, C, d, ...)
% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'leic', A, b, C, d, ...)
% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'lic', C, d, ...)
%
%
% The dots at the end denote any additional, problem specific data that are passed uniterpreted to
% all invocations of fname and jacname, see below for details.
44,14 → 50,21
% 'bc' specifies minimization subject to box constraints.
% 'lec' specifies minimization subject to linear equation constraints.
% 'blec' specifies minimization subject to box and linear equation constraints.
% 'bleic' specifies minimization subject to box, linear equation and inequality constraints.
% 'blic' specifies minimization subject to box and linear inequality constraints.
% 'leic' specifies minimization subject to linear equation and inequality constraints.
% 'lic' specifies minimization subject to linear inequality constraints.
% If omitted, a default of 'unc' is assumed. Depending on the minimization type, the MEX
% interface will invoke one of dlevmar_XXX, dlevmar_bc_XXX, dlevmar_lec_XXX or dlevmar_blec_XXX
% interface will invoke one of dlevmar_XXX, dlevmar_bc_XXX, dlevmar_lec_XXX, dlevmar_blec_XXX or dlevmar_bleic_XXX
%
% - lb, ub: vectors of doubles specifying lower and upper bounds for p, respectively
%
% - A, b: k x m matrix and k vector specifying linear equation constraints for p, i.e. A*p=b
% A should have full rank.
%
% - C, d: k x m matrix and k vector specifying linear inequality constraints for p, i.e. C*p>=d
% A should have full rank.
%
% - wghts: vector of doubles specifying the weights for the penalty terms corresponding to
% the box constraints, see lmblec_core.c for more details. If omitted and a 'blec' type
% minimization is to be carried out, default weights are used.
/gpu/src/levmar/matlab/modhs76.m New file
0,0 → 1,7
function x = modhs76(p)
 
x(1)=p(1);
x(2)=sqrt(0.5)*p(2);
x(3)=p(3);
x(4)=sqrt(0.5)*p(4);
 
/gpu/src/levmar/matlab/jachs01.m
1,4 → 1,4
function jac = hs01_jac(p)
function jac = jachs01(p)
m=2;
 
jac(1, 1:m)=[-20.0*p(1), 10.0];
/gpu/src/levmar/lm_core.c
37,6 → 37,7
#define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR)
#define AX_EQ_B_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS)
#define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD)
#define AX_EQ_B_BK LM_ADD_PREFIX(Ax_eq_b_BK)
#else
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack)
#endif /* HAVE_LAPACK */
293,11 → 294,12
 
/* solve augmented equations */
#ifdef HAVE_LAPACK
/* 5 alternatives are available: LU, Cholesky, 2 variants of QR decomposition and SVD.
/* 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
*/
 
//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;
//issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR;
516,7 → 518,7
 
if(!work){
worksz=LM_DIF_WORKSZ(m, n); //4*n+4*m + n*m + m*m;
work=(LM_REAL *)malloc(worksz*sizeof(LM_REAL)); /* allocate a big chunk in one step */
work=(LM_REAL *)calloc(1,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;
683,11 → 685,12
 
/* solve augmented equations */
#ifdef HAVE_LAPACK
/* 5 alternatives are available: LU, Cholesky, 2 variants of QR decomposition and SVD.
/* 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
*/
 
//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;
//issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR;
837,3 → 840,4
#undef AX_EQ_B_QR
#undef AX_EQ_B_QRLS
#undef AX_EQ_B_SVD
#undef AX_EQ_B_BK
/gpu/src/levmar/Axb_core.c
26,6 → 26,7
#error This file should not be compiled directly!
#endif
 
 
#ifdef LINSOLVERS_RETAIN_MEMORY
#define __STATIC__ static
#else
36,7 → 37,6
 
/* prototypes of LAPACK routines */
 
 
#define GEQRF LM_MK_LAPACK_NAME(geqrf)
#define ORGQR LM_MK_LAPACK_NAME(orgqr)
#define TRTRS LM_MK_LAPACK_NAME(trtrs)
47,6 → 47,8
#define GETRS LM_MK_LAPACK_NAME(getrs)
#define GESVD LM_MK_LAPACK_NAME(gesvd)
#define GESDD LM_MK_LAPACK_NAME(gesdd)
#define SYTRF LM_MK_LAPACK_NAME(sytrf)
#define SYTRS LM_MK_LAPACK_NAME(sytrs)
 
/* QR decomposition */
extern int GEQRF(int *m, int *n, LM_REAL *a, int *lda, LM_REAL *tau, LM_REAL *work, int *lwork, int *info);
74,12 → 76,17
extern int GESDD(char *jobz, int *m, int *n, LM_REAL *a, int *lda, LM_REAL *s, LM_REAL *u, int *ldu, LM_REAL *vt, int *ldvt,
LM_REAL *work, int *lwork, int *iwork, int *info);
 
/* LDLt/UDUt factorization and systems solution */
extern int SYTRF(char *uplo, int *n, LM_REAL *a, int *lda, int *ipiv, LM_REAL *work, int *lwork, int *info);
extern int SYTRS(char *uplo, int *n, int *nrhs, LM_REAL *a, int *lda, int *ipiv, LM_REAL *b, int *ldb, int *info);
 
/* precision-specific definitions */
#define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR)
#define AX_EQ_B_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS)
#define AX_EQ_B_CHOL LM_ADD_PREFIX(Ax_eq_b_Chol)
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU)
#define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD)
#define AX_EQ_B_BK LM_ADD_PREFIX(Ax_eq_b_BK)
 
/*
* This function returns the solution of Ax = b
105,8 → 112,8
 
static int nb=0; /* no __STATIC__ decl. here! */
 
LM_REAL *a, *qtb, *tau, *r, *work;
int a_sz, qtb_sz, tau_sz, r_sz, tot_sz;
LM_REAL *a, *tau, *r, *work;
int a_sz, tau_sz, r_sz, tot_sz;
register int i, j;
int info, worksz, nrhs=1;
register LM_REAL sum;
115,8 → 122,9
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
buf=NULL;
 
return 1;
}
#else
125,10 → 133,8
/* calculate required memory size */
a_sz=m*m;
qtb_sz=m;
tau_sz=m;
r_sz=m*m; /* only the upper triangular part really needed */
 
if(!nb){
LM_REAL tmp;
 
137,7 → 143,7
nb=((int)tmp)/m; // optimal worksize is m*nb
}
worksz=nb*m;
tot_sz=a_sz + qtb_sz + tau_sz + r_sz + worksz;
tot_sz=a_sz + tau_sz + r_sz + worksz;
 
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
160,8 → 166,7
#endif /* LINSOLVERS_RETAIN_MEMORY */
 
a=buf;
qtb=a+a_sz;
tau=qtb+qtb_sz;
tau=a+a_sz;
r=tau+tau_sz;
work=r+r_sz;
 
209,15 → 214,15
}
}
 
/* Q is now in a; compute Q^T b in qtb */
/* Q is now in a; compute Q^T b in x */
for(i=0; i<m; i++){
for(j=0, sum=0.0; j<m; j++)
sum+=a[i*m+j]*B[j];
qtb[i]=sum;
x[i]=sum;
}
 
/* solve the linear system R x = Q^t b */
TRTRS("U", "N", "N", (int *)&m, (int *)&nrhs, r, (int *)&m, qtb, (int *)&m, &info);
TRTRS("U", "N", "N", (int *)&m, (int *)&nrhs, r, (int *)&m, x, (int *)&m, &info);
/* error treatment */
if(info!=0){
if(info<0){
234,10 → 239,6
}
}
 
/* copy the result in x */
for(i=0; i<m; i++)
x[i]=qtb[i];
 
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
271,18 → 272,19
 
static int nb=0; /* no __STATIC__ decl. here! */
 
LM_REAL *a, *atb, *tau, *r, *work;
int a_sz, atb_sz, tau_sz, r_sz, tot_sz;
LM_REAL *a, *tau, *r, *work;
int a_sz, tau_sz, r_sz, tot_sz;
register int i, j;
int info, worksz, nrhs=1;
register LM_REAL sum;
 
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
 
return 1;
}
#else
296,7 → 298,6
/* calculate required memory size */
a_sz=m*n;
atb_sz=n;
tau_sz=n;
r_sz=n*n;
if(!nb){
307,7 → 308,7
nb=((int)tmp)/m; // optimal worksize is m*nb
}
worksz=nb*m;
tot_sz=a_sz + atb_sz + tau_sz + r_sz + worksz;
tot_sz=a_sz + tau_sz + r_sz + worksz;
 
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
330,8 → 331,7
#endif /* LINSOLVERS_RETAIN_MEMORY */
 
a=buf;
atb=a+a_sz;
tau=atb+atb_sz;
tau=a+a_sz;
r=tau+tau_sz;
work=r+r_sz;
 
340,11 → 340,11
for(j=0; j<n; j++)
a[i+j*m]=A[i*n+j];
 
/* compute A^T b in atb */
/* compute A^T b in x */
for(i=0; i<n; i++){
for(j=0, sum=0.0; j<m; j++)
sum+=A[j*n+i]*B[j];
atb[i]=sum;
x[i]=sum;
}
 
/* QR decomposition of A */
376,7 → 376,7
}
 
/* solve the linear system R^T y = A^t b */
TRTRS("U", "T", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, atb, (int *)&n, &info);
TRTRS("U", "T", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, x, (int *)&n, &info);
/* error treatment */
if(info!=0){
if(info<0){
394,7 → 394,7
}
 
/* solve the linear system R x = y */
TRTRS("U", "N", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, atb, (int *)&n, &info);
TRTRS("U", "N", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, x, (int *)&n, &info);
/* error treatment */
if(info!=0){
if(info<0){
411,10 → 411,6
}
}
 
/* copy the result in x */
for(i=0; i<n; i++)
x[i]=atb[i];
 
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
445,17 → 441,18
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
 
LM_REAL *a, *b;
int a_sz, b_sz, tot_sz;
LM_REAL *a;
int a_sz, tot_sz;
register int i;
int info, nrhs=1;
 
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
 
return 1;
}
#else
464,8 → 461,7
/* calculate required memory size */
a_sz=m*m;
b_sz=m;
tot_sz=a_sz + b_sz;
tot_sz=a_sz;
 
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
488,16 → 484,15
#endif /* LINSOLVERS_RETAIN_MEMORY */
 
a=buf;
b=a+a_sz;
 
/* store A into a anb B into b. A is assumed symmetric,
/* store A into a and B into x. A is assumed symmetric,
* hence no transposition is needed
*/
for(i=0; i<m; i++){
a[i]=A[i];
b[i]=B[i];
x[i]=B[i];
}
for(i=m; i<m*m; i++)
for(i=m; i<m*m; i++)
a[i]=A[i];
 
/* Cholesky decomposition of A */
511,8 → 506,7
exit(1);
}
else{
fprintf(stderr, RCAT(RCAT(RCAT("LAPACK error: the leading minor of order %d is not positive definite,\nthe factorization could not be completed for ", POTF2) "/", POTRF), " in ",
AX_EQ_B_CHOL) "()\n", info);
fprintf(stderr, RCAT(RCAT(RCAT("LAPACK error: the leading minor of order %d is not positive definite,\nthe factorization could not be completed for ", POTF2) "/", POTRF) " in ", AX_EQ_B_CHOL) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
522,7 → 516,7
}
 
/* solve using the computed Cholesky in one lapack call */
POTRS("U", (int *)&m, (int *)&nrhs, a, (int *)&m, b, (int *)&m, &info);
POTRS("U", (int *)&m, (int *)&nrhs, a, (int *)&m, x, (int *)&m, &info);
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", POTRS) " in ", AX_EQ_B_CHOL) "()\n", -info);
exit(1);
530,7 → 524,7
 
#if 0
/* alternative: solve the linear system U^T y = b ... */
TRTRS("U", "T", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, b, (int *)&m, &info);
TRTRS("U", "T", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, x, (int *)&m, &info);
/* error treatment */
if(info!=0){
if(info<0){
548,7 → 542,7
}
 
/* ... solve the linear system U x = y */
TRTRS("U", "N", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, b, (int *)&m, &info);
TRTRS("U", "N", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, x, (int *)&m, &info);
/* error treatment */
if(info!=0){
if(info<0){
566,10 → 560,6
}
#endif /* 0 */
 
/* copy the result in x */
for(i=0; i<m; i++)
x[i]=b[i];
 
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
599,17 → 589,18
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
 
int a_sz, ipiv_sz, b_sz, tot_sz;
int a_sz, ipiv_sz, tot_sz;
register int i, j;
int info, *ipiv, nrhs=1;
LM_REAL *a, *b;
 
LM_REAL *a;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
 
return 1;
}
#else
619,8 → 610,7
/* calculate required memory size */
ipiv_sz=m;
a_sz=m*m;
b_sz=m;
tot_sz=(a_sz + b_sz)*sizeof(LM_REAL) + ipiv_sz*sizeof(int); /* should be arranged in that order for proper doubles alignment */
tot_sz=a_sz*sizeof(LM_REAL) + ipiv_sz*sizeof(int); /* should be arranged in that order for proper doubles alignment */
 
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
643,15 → 633,14
#endif /* LINSOLVERS_RETAIN_MEMORY */
 
a=buf;
b=a+a_sz;
ipiv=(int *)(b+b_sz);
ipiv=(int *)(a+a_sz);
 
/* store A (column major!) into a and B into b */
/* store A (column major!) into a and B into x */
for(i=0; i<m; i++){
for(j=0; j<m; j++)
a[i+j*m]=A[i*m+j];
 
b[i]=B[i];
x[i]=B[i];
}
 
/* LU decomposition for A */
672,7 → 661,7
}
 
/* solve the system with the computed LU */
GETRS("N", (int *)&m, (int *)&nrhs, a, (int *)&m, ipiv, b, (int *)&m, (int *)&info);
GETRS("N", (int *)&m, (int *)&nrhs, a, (int *)&m, ipiv, x, (int *)&m, (int *)&info);
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("argument %d of ", GETRS) " illegal in ", AX_EQ_B_LU) "()\n", -info);
688,11 → 677,6
}
}
 
/* copy the result in x */
for(i=0; i<m; i++){
x[i]=b[i];
}
 
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
729,13 → 713,14
LM_REAL thresh, one_over_denom;
register LM_REAL sum;
int info, rank, worksz, *iwork, iworksz;
 
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
 
return 1;
}
#else
844,22 → 829,145
return 1;
}
 
/*
* This function returns the solution of Ax = b for a real symmetric matrix A
*
* The function is based on UDUT factorization with the pivoting
* strategy of Bunch and Kaufman:
* A is factored as U*D*U^T where U is upper triangular and
* D symmetric and block diagonal (aka spectral decomposition,
* Banachiewicz factorization, modified Cholesky factorization)
*
* A is mxm, b is mx1.
*
* The function returns 0 in case of error, 1 if successfull
*
* This function is often called repetitively to solve problems of identical
* dimensions. To avoid repetitive malloc's and free's, allocated memory is
* retained between calls and free'd-malloc'ed when not of the appropriate size.
* A call with NULL as the first argument forces this memory to be released.
*/
int AX_EQ_B_BK(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0, nb=0;
 
LM_REAL *a, *work;
int a_sz, ipiv_sz, work_sz, tot_sz;
register int i, j;
int info, *ipiv, nrhs=1;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
 
return 1;
}
#else
return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */
 
/* calculate required memory size */
ipiv_sz=m;
a_sz=m*m;
if(!nb){
LM_REAL tmp;
 
work_sz=-1; // workspace query; optimal size is returned in tmp
SYTRF("U", (int *)&m, NULL, (int *)&m, NULL, (LM_REAL *)&tmp, (int *)&work_sz, (int *)&info);
nb=((int)tmp)/m; // optimal worksize is m*nb
}
work_sz=(nb!=-1)? nb*m : 1;
tot_sz=(a_sz + work_sz)*sizeof(LM_REAL) + ipiv_sz*sizeof(int); /* should be arranged in that order for proper doubles alignment */
 
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
if(buf) free(buf); /* free previously allocated memory */
 
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz);
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_BK) "() failed!\n");
exit(1);
}
}
#else
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz);
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_BK) "() failed!\n");
exit(1);
}
#endif /* LINSOLVERS_RETAIN_MEMORY */
 
a=buf;
work=a+a_sz;
ipiv=(int *)(work+work_sz);
 
/* store A into a and B into x; A is assumed to be symmetric, hence
* the column and row major order representations are the same
*/
for(i=0; i<m; ++i){
a[i]=A[i];
x[i]=B[i];
}
for(j=m*m; i<j; ++i) // copy remaining rows; note that i is not re-initialized
a[i]=A[i];
 
/* UDUt factorization for A */
SYTRF("U", (int *)&m, a, (int *)&m, ipiv, work, (int *)&work_sz, (int *)&info);
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", SYTRF) " in ", AX_EQ_B_BK) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT(RCAT("LAPACK error: singular block diagonal matrix D for", SYTRF) " in ", AX_EQ_B_BK)"() [D(%d, %d) is zero]\n", info, info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
 
return 0;
}
}
 
/* solve the system with the computed factorization */
SYTRS("U", (int *)&m, (int *)&nrhs, a, (int *)&m, ipiv, x, (int *)&m, (int *)&info);
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", SYTRS) " in ", AX_EQ_B_BK) "()\n", -info);
exit(1);
}
 
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
 
return 1;
}
 
/* undefine all. IT MUST REMAIN IN THIS POSITION IN FILE */
#undef AX_EQ_B_QR
#undef AX_EQ_B_QRLS
#undef AX_EQ_B_CHOL
#undef AX_EQ_B_LU
#undef AX_EQ_B_SVD
#undef AX_EQ_B_BK
 
#undef GEQRF
#undef ORGQR
#undef TRTRS
#undef POTF2
#undef POTRF
#undef POTRS
#undef GETRF
#undef GETRS
#undef GESVD
#undef GESDD
#undef SYTRF
#undef SYTRS
 
#else // no LAPACK
 
876,14 → 984,15
extern int GETRF(int *m, int *n, LM_REAL *a, int *lda, int *ipiv, int *info);
extern int GETRS(char *trans, int *n, int *nrhs, LM_REAL *a, int *lda, int *ipiv, LM_REAL *b, int *ldb, int *info);
/* End added by EB */
 
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack)
 
/*
* This function returns the solution of Ax = b
*
* The function employs LU decomposition followed by forward/back substitution (see
* also the LAPACK-based LU solver above)
* The function employs LU decomposition:
* If A=L U with L lower and U upper triangular, then the original system
* amounts to solving
* L y = b, U x = y
*
* A is mxm, b is mx1
*
896,20 → 1005,21
*/
int AX_EQ_B_LU(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
 
__STATIC__ void *buf=NULL;
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
 
int a_sz, ipiv_sz, b_sz, tot_sz;
int a_sz, ipiv_sz, tot_sz;
register int i, j;
int info, *ipiv;
LM_REAL *a, *b;
int info, *ipiv, nrhs=1;
LM_REAL *a;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
 
return 1;
}
#else
919,15 → 1029,14
/* calculate required memory size */
ipiv_sz=m;
a_sz=m*m;
b_sz=m;
tot_sz=(ipiv_sz + a_sz + b_sz)*sizeof(LM_REAL);
tot_sz=a_sz*sizeof(LM_REAL) + ipiv_sz*sizeof(int); /* should be arranged in that order for proper doubles alignment */
 
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
if(buf) free(buf); /* free previously allocated memory */
 
buf_sz=tot_sz;
buf=(void *)malloc(buf_sz);
buf=(LM_REAL *)malloc(buf_sz);
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
exit(1);
935,26 → 1044,26
}
#else
buf_sz=tot_sz;
buf=(void *)malloc(buf_sz);
buf=(LM_REAL *)malloc(buf_sz);
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
exit(1);
}
#endif /* LINSOLVERS_RETAIN_MEMORY */
 
ipiv=(int *)buf;
a=(LM_REAL *)(ipiv + ipiv_sz);
b=a+a_sz;
a=buf;
ipiv=(int *)(a+a_sz);
 
/* store A (column major!) into a and B into b */
/* store A (column major!) into a and B into x */
for(i=0; i<m; i++){
for(j=0; j<m; j++)
a[i+j*m]=A[i*m+j];
 
b[i]=B[i];
x[i]=B[i];
}
 
/* LU decomposition for A */
//GETRF((int *)&m, (int *)&m, a, (int *)&m, ipiv, (int *)&info);
info = ATLAS_POTRF(CblasRowMajor, CblasUpper, m, a, m);
if(info!=0){
if(info<0){
972,7 → 1081,8
}
 
/* solve the system with the computed LU */
info = ATLAS_POTRS(CblasRowMajor, CblasUpper, m, 1, a, m, b, m);
//GETRS("N", (int *)&m, (int *)&nrhs, a, (int *)&m, ipiv, x, (int *)&m, (int *)&info);
info = ATLAS_POTRS(CblasRowMajor, CblasUpper, m, 1, a, m, x, m);
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("argument %d of ", GETRS) " illegal in ", AX_EQ_B_LU) "()\n", -info);
988,21 → 1098,177
}
}
 
/* copy the result in x */
for(i=0; i<m; i++){
x[i]=b[i];
}
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
 
return 1;
return 1;
}
 
///* This function returns the solution of Ax = b
// *
// * The function employs LU decomposition followed by forward/back substitution (see
// * also the LAPACK-based LU solver above)
// *
// * A is mxm, b is mx1
// *
// * The function returns 0 in case of error, 1 if successful
// *
// * This function is often called repetitively to solve problems of identical
// * dimensions. To avoid repetitive malloc's and free's, allocated memory is
// * retained between calls and free'd-malloc'ed when not of the appropriate size.
// * A call with NULL as the first argument forces this memory to be released.
// */
//int AX_EQ_B_LU(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
//{
//__STATIC__ void *buf=NULL;
//__STATIC__ int buf_sz=0;
//
//register int i, j, k;
//int *idx, maxi=-1, idx_sz, a_sz, work_sz, tot_sz;
//LM_REAL *a, *work, max, sum, tmp;
//
// if(!A)
//#ifdef LINSOLVERS_RETAIN_MEMORY
// {
// if(buf) free(buf);
// buf=NULL;
// buf_sz=0;
//
// return 1;
// }
//#else
// return 1; /* NOP */
//#endif /* LINSOLVERS_RETAIN_MEMORY */
//
// /* calculate required memory size */
// idx_sz=m;
// a_sz=m*m;
// work_sz=m;
// tot_sz=(a_sz+work_sz)*sizeof(LM_REAL) + idx_sz*sizeof(int); /* should be arranged in that order for proper doubles alignment */
//
//#ifdef LINSOLVERS_RETAIN_MEMORY
// if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
// if(buf) free(buf); /* free previously allocated memory */
//
// buf_sz=tot_sz;
// buf=(void *)malloc(tot_sz);
// if(!buf){
// fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
// exit(1);
// }
// }
//#else
// buf_sz=tot_sz;
// buf=(void *)malloc(tot_sz);
// if(!buf){
// fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
// exit(1);
// }
//#endif /* LINSOLVERS_RETAIN_MEMORY */
//
// a=buf;
// work=a+a_sz;
// idx=(int *)(work+work_sz);
//
// /* avoid destroying A, B by copying them to a, x resp. */
// for(i=0; i<m; ++i){ // B & 1st row of A
// a[i]=A[i];
// x[i]=B[i];
// }
// for( ; i<a_sz; ++i) a[i]=A[i]; // copy A's remaining rows
// /****
// for(i=0; i<m; ++i){
// for(j=0; j<m; ++j)
// a[i*m+j]=A[i*m+j];
// x[i]=B[i];
// }
// ****/
//
// /* compute the LU decomposition of a row permutation of matrix a; the permutation itself is saved in idx[] */
// for(i=0; i<m; ++i){
// max=0.0;
// for(j=0; j<m; ++j)
// if((tmp=FABS(a[i*m+j]))>max)
// max=tmp;
// if(max==0.0){
// fprintf(stderr, RCAT("Singular matrix A in ", AX_EQ_B_LU) "()!\n");
//#ifndef LINSOLVERS_RETAIN_MEMORY
// free(buf);
//#endif
//
// return 0;
// }
// work[i]=LM_CNST(1.0)/max;
// }
//
// for(j=0; j<m; ++j){
// for(i=0; i<j; ++i){
// sum=a[i*m+j];
// for(k=0; k<i; ++k)
// sum-=a[i*m+k]*a[k*m+j];
// a[i*m+j]=sum;
// }
// max=0.0;
// for(i=j; i<m; ++i){
// sum=a[i*m+j];
// for(k=0; k<j; ++k)
// sum-=a[i*m+k]*a[k*m+j];
// a[i*m+j]=sum;
// if((tmp=work[i]*FABS(sum))>=max){
// max=tmp;
// maxi=i;
// }
// }
// if(j!=maxi){
// for(k=0; k<m; ++k){
// tmp=a[maxi*m+k];
// a[maxi*m+k]=a[j*m+k];
// a[j*m+k]=tmp;
// }
// work[maxi]=work[j];
// }
// idx[j]=maxi;
// if(a[j*m+j]==0.0)
// a[j*m+j]=LM_REAL_EPSILON;
// if(j!=m-1){
// tmp=LM_CNST(1.0)/(a[j*m+j]);
// for(i=j+1; i<m; ++i)
// a[i*m+j]*=tmp;
// }
// }
//
// /* The decomposition has now replaced a. Solve the linear system using
// * forward and back substitution
// */
// for(i=k=0; i<m; ++i){
// j=idx[i];
// sum=x[j];
// x[j]=x[i];
// if(k!=0)
// for(j=k-1; j<i; ++j)
// sum-=a[i*m+j]*x[j];
// else
// if(sum!=0.0)
// k=i+1;
// x[i]=sum;
// }
//
// for(i=m-1; i>=0; --i){
// sum=x[i];
// for(j=i+1; j<m; ++j)
// sum-=a[i*m+j]*x[j];
// x[i]=sum/a[i*m+i];
// }
//
//#ifndef LINSOLVERS_RETAIN_MEMORY
// free(buf);
//#endif
//
// return 1;
//}
//
/* undefine all. IT MUST REMAIN IN THIS POSITION IN FILE */
#undef AX_EQ_B_LU
/* Added by EB */
#undef GETRF
#undef GETRS
/* End Added by EB */
 
#endif /* HAVE_LAPACK */
/gpu/src/levmar/lmbc_core.c
44,6 → 44,7
#define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR)
#define AX_EQ_B_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS)
#define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD)
#define AX_EQ_B_BK LM_ADD_PREFIX(Ax_eq_b_BK)
#else
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack)
#endif /* HAVE_LAPACK */
551,11 → 552,12
 
/* solve augmented equations */
#ifdef HAVE_LAPACK
/* 5 alternatives are available: LU, Cholesky, 2 variants of QR decomposition and SVD.
/* 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
*/
 
//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;
//issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR;
943,3 → 945,4
#undef AX_EQ_B_QR
#undef AX_EQ_B_QRLS
#undef AX_EQ_B_SVD
#undef AX_EQ_B_BK
/gpu/src/levmar/Makefile.am
1,10 → 1,13
# Makefile.am for the levmar Levenberg-Marquardt fitting library
# Copyright (C) 2008 Emmanuel Bertin.
# Copyright (C) 2007-2010 Emmanuel Bertin.
# levmar code by Manolis Lourakis <lourakis@ics.forth.gr>
# http://www.ics.forth.gr/~lourakis/levmar/
noinst_LIBRARIES = liblevmar.a
liblevmar_a_SOURCES = Axb.c lmbc.c lm.c lmblec.c lmlec.c misc.c \
compiler.h lm.h misc.h
liblevmar_a_SOURCES = Axb.c lmbc.c lm.c lmblec.c lmbleic.c lmlec.c misc.c \
compiler.h levmar.h misc.h
EXTRA_liblevmar_a_SOURCES = Axb_core.c lmbc_core.c lm_core.c \
lmblec_core.c lmlec_core.c misc_core.c \
LICENSE README README.txt lmdemo.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
/gpu/src/levmar/Makefile.icc
11,8 → 11,8
CFLAGS=$(CONFIGFLAGS) $(ARCHFLAGS) -O3 -tpp7 -xW -ip -ipo -unroll #-g
LAPACKLIBS_PATH=/usr/local/lib # WHEN USING LAPACK, CHANGE THIS TO WHERE YOUR COMPILED LIBS ARE!
LDFLAGS=-L$(LAPACKLIBS_PATH) -L.
LIBOBJS=lm.o Axb.o misc.o lmlec.o lmbc.o lmblec.o
LIBSRCS=lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c
LIBOBJS=lm.o Axb.o misc.o lmlec.o lmbc.o lmblec.o lmbleic.o
LIBSRCS=lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c lmbleic.c
DEMOBJS=lmdemo.o
DEMOSRCS=lmdemo.c
AR=xiar
24,6 → 24,9
# The following works with the ATLAS updated lapack and Linux_P4SSE2 from http://www.netlib.org/atlas/archives/linux/
#LAPACKLIBS=-L/usr/local/atlas/lib -llapack -lcblas -lf77blas -latlas -lf2c
 
#LAPACKLIBS=-llapack -lgoto2 -lpthread -lf2c # This works with GotoBLAS
# from http://www.tacc.utexas.edu/research-development/tacc-projects/
 
LIBS=$(LAPACKLIBS)
 
all: liblevmar.a lmdemo
35,14 → 38,15
lmdemo: $(DEMOBJS) liblevmar.a
$(CC) $(ARCHFLAGS) $(LDFLAGS) $(DEMOBJS) -o lmdemo -llevmar $(LIBS) -lm
 
lm.o: lm.c lm_core.c lm.h misc.h compiler.h
Axb.o: Axb.c Axb_core.c lm.h misc.h
misc.o: misc.c misc_core.c lm.h misc.h
lmlec.o: lmlec.c lmlec_core.c lm.h misc.h
lmbc.o: lmbc.c lmbc_core.c lm.h misc.h compiler.h
lmblec.o: lmblec.c lmblec_core.c lm.h misc.h
lm.o: lm.c lm_core.c levmar.h misc.h compiler.h
Axb.o: Axb.c Axb_core.c levmar.h misc.h
misc.o: misc.c misc_core.c levmar.h misc.h
lmlec.o: lmlec.c lmlec_core.c levmar.h misc.h
lmbc.o: lmbc.c lmbc_core.c levmar.h misc.h compiler.h
lmblec.o: lmblec.c lmblec_core.c levmar.h misc.h
lmbleic.o: lmbleic.c lmbleic_core.c levmar.h misc.h
 
lmdemo.o: lm.h
lmdemo.o: levmar.h
 
clean:
@rm -f $(LIBOBJS) $(DEMOBJS)
/gpu/src/levmar/README.txt
1,6 → 1,6
**************************************************************
LEVMAR
version 2.4
version 2.5
By Manolis Lourakis
 
Institute of Computer Science
12,15 → 12,15
GENERAL
This is levmar, a copylefted C/C++ implementation of the Levenberg-Marquardt non-linear
least squares algorithm. levmar includes double and single precision LM versions, both
with analytic and finite difference approximated jacobians. levmar also has some support
for constrained non-linear least squares, allowing linear equation and box constraints.
You have the following options regarding the solution of the underlying augmented normal
equations:
with analytic and finite difference approximated Jacobians. levmar also has some support
for constrained non-linear least squares, allowing linear equation, box and linear
inequality constraints. The following options regarding the solution of the underlying
augmented normal equations are offered:
 
1) Assuming that you have LAPACK (or an equivalent vendor library such as ESSL, MKL,
NAG, ...) installed, you can use the included LAPACK-based solvers (default).
 
2) If you don't have LAPACK or decide not to use it, undefine HAVE_LAPACK in lm.h
2) If you don't have LAPACK or decide not to use it, undefine HAVE_LAPACK in levmar.h
and a LAPACK-free, LU-based linear systems solver will by used. Also, the line
setting the variable LAPACKLIBS in the Makefile should be commented out.
 
38,12 → 38,12
levmar is released under the GNU Public License (GPL), which can be found in the included
LICENSE file. Note that under the terms of GPL, commercial use is allowed only if a software
employing levmar is also published in source under the GPL. However, if you are interested
in using levmar in a proprietary commercial apprlication, a commercial license for levmar
in using levmar in a proprietary commercial application, a commercial license for levmar
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
lm.h. See the accompanying comments for more details.
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
63,7 → 63,7
for details.
 
MATLAB INTERFACE
Since version 2.2, the levmar distrubution includes a matlab interface.
Since version 2.2, the levmar distribution includes a matlab 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,
/gpu/src/levmar/misc.c
28,7 → 28,7
#include <math.h>
#include <float.h>
 
#include "lm.h"
#include "levmar.h"
#include "misc.h"
 
#if !defined(LM_DBL_PREC) && !defined(LM_SNGL_PREC)
/gpu/src/levmar/levmar.h New file
0,0 → 1,370
/*
////////////////////////////////////////////////////////////////////////////////////
//
// 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 *************************************/
 
/* 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_ */
/gpu/src/levmar/lmlec.c
27,7 → 27,7
#include <stdlib.h>
#include <math.h>
 
#include "lm.h"
#include "levmar.h"
#include "misc.h"
 
 
35,10 → 35,8
 
#ifdef _MSC_VER
#pragma message("Linearly constrained optimization requires LAPACK and was not compiled!")
/*
#else
#warning Linearly constrained optimization requires LAPACK and was not compiled!
*/
//#warning Linearly constrained optimization requires LAPACK and was not compiled!
#endif // _MSC_VER
 
#else // LAPACK present
/gpu/src/levmar/CMakeLists.txt
5,8 → 5,8
#CMAKE_MINIMUM_REQUIRED(VERSION 1.4)
 
# compiler flags
ADD_DEFINITIONS(-DLINSOLVERS_RETAIN_MEMORY) # do not free memory between linear solvers calls
#REMOVE_DEFINITIONS(-DLINSOLVERS_RETAIN_MEMORY)
#ADD_DEFINITIONS(-DLINSOLVERS_RETAIN_MEMORY) # do not free memory between linear solvers calls
#REMOVE_DEFINITIONS(-DLINSOLVERS_RETAIN_MEMORY) # free memory between calls
 
# f2c is sometimes equivalent to libF77 & libI77; in that case, set HAVE_F2C to 0
SET(HAVE_F2C 1 CACHE BOOL "Do we have f2c or F77/I77?" )
26,17 → 26,19
 
########################## NO CHANGES BEYOND THIS POINT ##########################
 
INCLUDE_DIRECTORIES(.)
#INCLUDE_DIRECTORIES(/usr/include)
LINK_DIRECTORIES(${LAPACKBLAS_DIR})
 
# levmar library source files
ADD_LIBRARY(levmar STATIC
lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c
lm.h misc.h compiler.h
lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c lmbleic.c
levmar.h misc.h compiler.h
)
 
# demo program
ADD_EXECUTABLE(lmdemo lmdemo.c lm.h)
LINK_DIRECTORIES(${LAPACKBLAS_DIR})
LINK_DIRECTORIES(.)
ADD_EXECUTABLE(lmdemo lmdemo.c levmar.h)
# libraries the demo depends on
IF(HAVE_F2C)
TARGET_LINK_LIBRARIES(lmdemo levmar ${LAPACK_LIB} ${BLAS_LIB} ${F2C_LIB})
/gpu/src/levmar/lmblec.c
28,17 → 28,15
#include <math.h>
#include <float.h>
 
#include "lm.h"
#include "levmar.h"
#include "misc.h"
 
#ifndef HAVE_LAPACK
 
#ifdef _MSC_VER
#pragma message("Combined box and linearly constrained optimization requires LAPACK and was not compiled!")
/*
#else
#warning Combined box and linearly constrained optimization requires LAPACK and was not compiled!
*/
//#warning Combined box and linearly constrained optimization requires LAPACK and was not compiled!
#endif // _MSC_VER
 
#else // LAPACK present
/gpu/src/levmar/expfit.c
23,7 → 23,7
#include <stdlib.h>
#include <math.h>
 
#include <lm.h>
#include <levmar.h>
 
#ifndef LM_DBL_PREC
#error Example program assumes that levmar has been compiled with double precision, see LM_DBL_PREC!
/gpu/src/levmar/Makefile.in
50,7 → 50,8
liblevmar_a_AR = $(AR) $(ARFLAGS)
liblevmar_a_LIBADD =
am_liblevmar_a_OBJECTS = Axb.$(OBJEXT) lmbc.$(OBJEXT) lm.$(OBJEXT) \
lmblec.$(OBJEXT) lmlec.$(OBJEXT) misc.$(OBJEXT)
lmblec.$(OBJEXT) lmbleic.$(OBJEXT) lmlec.$(OBJEXT) \
misc.$(OBJEXT)
liblevmar_a_OBJECTS = $(am_liblevmar_a_OBJECTS)
DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir)
depcomp = $(SHELL) $(top_srcdir)/autoconf/depcomp
195,16 → 196,19
top_srcdir = @top_srcdir@
 
# Makefile.am for the levmar Levenberg-Marquardt fitting library
# Copyright (C) 2008 Emmanuel Bertin.
# Copyright (C) 2007-2010 Emmanuel Bertin.
# levmar code by Manolis Lourakis <lourakis@ics.forth.gr>
# http://www.ics.forth.gr/~lourakis/levmar/
noinst_LIBRARIES = liblevmar.a
liblevmar_a_SOURCES = Axb.c lmbc.c lm.c lmblec.c lmlec.c misc.c \
compiler.h lm.h misc.h
liblevmar_a_SOURCES = Axb.c lmbc.c lm.c lmblec.c lmbleic.c lmlec.c misc.c \
compiler.h levmar.h misc.h
 
EXTRA_liblevmar_a_SOURCES = Axb_core.c lmbc_core.c lm_core.c \
lmblec_core.c lmlec_core.c misc_core.c \
LICENSE README README.txt lmdemo.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
 
all: all-am
 
255,12 → 259,15
 
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Axb.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Axb_core.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/expfit.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lm.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lm_core.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmbc.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmbc_core.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmblec.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmblec_core.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmbleic.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmbleic_core.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmdemo.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmlec.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmlec_core.Po@am__quote@
/gpu/src/levmar/misc_core.c
543,6 → 543,16
}
#endif /* HAVE_LAPACK */
 
/* precision-specific definitions */
/* Added by EB */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include ATLAS_LAPACK_H
#define ATLAS_POTRF LM_CAT_(clapack_, LM_ADD_PREFIX(potrf))
#define ATLAS_POTRI LM_CAT_(clapack_, LM_ADD_PREFIX(potri))
/* End added by EB */
 
/*
* This function computes in C the covariance matrix corresponding to a least
* squares fit. JtJ is the approximate Hessian at the solution (i.e. J^T*J, where
564,7 → 574,7
*/
int LEVMAR_COVAR(LM_REAL *JtJ, LM_REAL *C, LM_REAL sumsq, int m, int n)
{
register int i;
register int i,j;
int rnk;
LM_REAL fact;
 
581,24 → 591,27
#endif // _MSC_VER
 
// rnk=LEVMAR_LUINVERSE(JtJ, C, m);
 
rnk = SVDINV(JtJ, C, m);
if(!rnk) return 0;
 
// rnk=m; /* assume full rank */
if (!rnk) return 0;
 
// rnk=m; /* assume full rank */
#endif /* HAVE_LAPACK */
fact=(LM_REAL)n/(LM_REAL)(n-rnk);
for(i=0; i<m*m; ++i)
C[i]*=fact;
 
// fact=sumsq/(LM_REAL)(n-rnk);
// for(i=0; i<m*m; ++i)
// C[i]*=fact;
for(i=0; i<m; ++i)
if (C[i*(m+1)]<0.0 || fabs(C[i*(m+1)])>1e30)
{
for(j=0; j<m; ++j)
C[i*m+j] = 0.0;
for(j=0; j<m; ++j)
C[j*m+i] = 0.0;
}
 
fact=(LM_REAL)n/(LM_REAL)(n-rnk);
for(i=0; i<m*m; ++i)
C[i]*=fact;
 
for(i=0; i<m; ++i)
if (C[i*(m+1)]<0.0)
C[i*(m+1)] = 0.0;
 
return rnk;
}
 
681,12 → 694,6
register int i, j;
int info;
 
/* compute the Cholesky decomposition of C in W, s.t. C=W^t W and W is upper triangular */
int LEVMAR_CHOLESKY(LM_REAL *C, LM_REAL *W, int m)
{
register int i, j;
int info;
 
/* copy weights array C to W so that LAPACK won't destroy it;
* C is assumed symmetric, hence no transposition is needed
*/
697,13 → 704,13
POTF2("U", (int *)&m, W, (int *)&m, (int *)&info);
/* error treatment */
if(info!=0){
if(info<0){
if(info<0){
fprintf(stderr, "LAPACK error: illegal value for argument %d of dpotf2 in %s\n", -info, LCAT(LEVMAR_CHOLESKY, "()"));
}
else{
fprintf(stderr, "LAPACK error: the leading minor of order %d is not positive definite,\n%s()\n", info,
RCAT("and the Cholesky factorization could not be completed in ", LEVMAR_CHOLESKY));
}
}
else{
fprintf(stderr, "LAPACK error: the leading minor of order %d is not positive definite,\n%s()\n", info,
RCAT("and the Cholesky factorization could not be completed in ", LEVMAR_CHOLESKY));
}
return LM_ERROR;
}
 
769,12 → 776,12
 
switch(n - i){
case 7 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i;
case 6 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i;
case 5 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i;
case 4 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i;
case 6 : e[i]=x[i]-y[i]; sum1+=e[i]*e[i]; ++i;
case 5 : e[i]=x[i]-y[i]; sum2+=e[i]*e[i]; ++i;
case 4 : e[i]=x[i]-y[i]; sum3+=e[i]*e[i]; ++i;
case 3 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i;
case 2 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i;
case 1 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i;
case 2 : e[i]=x[i]-y[i]; sum1+=e[i]*e[i]; ++i;
case 1 : e[i]=x[i]-y[i]; sum2+=e[i]*e[i]; //++i;
}
}
}
804,12 → 811,12
 
switch(n - i){
case 7 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i;
case 6 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i;
case 5 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i;
case 4 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i;
case 6 : e[i]=-y[i]; sum1+=e[i]*e[i]; ++i;
case 5 : e[i]=-y[i]; sum2+=e[i]*e[i]; ++i;
case 4 : e[i]=-y[i]; sum3+=e[i]*e[i]; ++i;
case 3 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i;
case 2 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i;
case 1 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i;
case 2 : e[i]=-y[i]; sum1+=e[i]*e[i]; ++i;
case 1 : e[i]=-y[i]; sum2+=e[i]*e[i]; //++i;
}
}
}
827,284 → 834,270
NOTES Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671). The a
and v matrices are transposed with respect to the N.R. convention.
AUTHOR E. Bertin (IAP)
VERSION 30/05/2008
VERSION 02/09/2010
***/
 
static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
{
#define MIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) < (maxarg2) ?\
(maxarg1) : (maxarg2))
#define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
(maxarg1) : (maxarg2))
#define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \
(ct=bt/at,at*sqrt(1.0+ct*ct)) \
: (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
#define TOL 1.0e-6
 
int flag,i,its,j,jj,k,l,mmi,nm, nml, rank;
int flag,i,its,j,jj,k,l,nm, rank;
LM_REAL *vmat,*wmat,
*w,*ap,*ap0,*ap1,*ap10,*rv1p,*vp,*vp0,*vp1,*vp10,
*rv1,*tmp,
*w,*rv1,*tmp,
c,f,h,s,x,y,z,
anorm, g, scale,
at,bt,ct,maxarg1,maxarg2,
thresh, wmax;
thresh, wmax, tol,tanorm;
anorm = g = scale = 0.0;
 
tol = sizeof(anorm)>4? 1.0e-9 : 1.0e-6;
 
rv1=(LM_REAL *)malloc(m*sizeof(LM_REAL));
tmp=(LM_REAL *)malloc(m*sizeof(LM_REAL));
vmat=(LM_REAL *)malloc(m*m*sizeof(LM_REAL));
wmat=(LM_REAL *)malloc(m*sizeof(LM_REAL));
 
l = nm = nml = 0; /* To avoid gcc -Wall warnings */
for (i=0;i<m;i++)
l = nm = 0; /* To avoid gcc -Wall warnings */
 
for (i=0; i<m; i++)
{
l = i+1;
nml = m-l;
rv1[i] = scale*g;
g = s = scale = 0.0;
mmi = m - i;
ap = ap0 = a+i*(m+1);
for (k=mmi;k--;)
scale += fabs(*(ap++));
for (k=i; k<m; k++)
scale += fabs(a[k*m+i]);
if (scale)
{
for (ap=ap0,k=mmi; k--; ap++)
for (k=i; k<m; k++)
{
*ap /= scale;
s += *ap**ap;
a[k*m+i] /= scale;
s += a[k*m+i]*a[k*m+i];
}
f = *ap0;
f = a[i*m+i];
g = -SIGN(sqrt(s),f);
h = f*g-s;
*ap0 = f-g;
ap10 = a+l*m+i;
for (j=nml;j--; ap10+=m)
h = f*g - s;
a[i*m+i] = f-g;
for (j=l; j<m; j++)
{
s = 0.0;
for (ap=ap0,ap1=ap10,k=mmi; k--;)
s += *(ap1++)**(ap++);
for (k=i; k<m; k++)
s += a[k*m+i] * a[k*m+j];
f = s/h;
for (ap=ap0,ap1=ap10,k=mmi; k--;)
*(ap1++) += f**(ap++);
for (k=i; k<m; k++)
a[k*m+j] += f * a[k*m+i];
}
for (ap=ap0,k=mmi; k--;)
*(ap++) *= scale;
for (k=i; k<m; k++)
a[k*m+i] *= scale;
}
wmat[i] = scale*g;
wmat[i]= scale * g;
g = s = scale = 0.0;
if (i+1 != m)
if (l != m)
{
ap = ap0 = a+i+m*l;
for (k=nml;k--; ap+=m)
scale += fabs(*ap);
for (k=l; k<m; k++)
scale += fabs(a[i*m+k]);
if (scale)
{
for (ap=ap0,k=nml;k--; ap+=m)
for (k=l; k<m; k++)
{
*ap /= scale;
s += *ap**ap;
a[i*m+k] /= scale;
s += a[i*m+k] * a[i*m+k];
}
f=*ap0;
f = a[i*m+l];
g = -SIGN(sqrt(s),f);
h=f*g-s;
*ap0=f-g;
rv1p = rv1+l;
for (ap=ap0,k=nml;k--; ap+=m)
*(rv1p++) = *ap/h;
ap10 = a+l+m*l;
for (j=m-l; j--; ap10++)
h = f*g - s;
a[i*m+l] = f - g;
for (k=l; k<m; k++)
rv1[k] = a[i*m+k] / h;
for (j=l; j<m; j++)
{
for (s=0.0,ap=ap0,ap1=ap10,k=nml; k--; ap+=m,ap1+=m)
s += *ap1**ap;
rv1p = rv1+l;
for (ap1=ap10,k=nml;k--; ap1+=m)
*ap1 += s**(rv1p++);
s = 0.0;
for (k=l; k<m; k++)
s += a[j*m+k] * a[i*m+k];
for (k=l; k<m; k++)
a[j*m+k] += s * rv1[k];
}
for (ap=ap0,k=nml;k--; ap+=m)
*ap *= scale;
for (k=l; k<m; k++)
a[i*m+k] *= scale;
}
}
anorm=MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
anorm = MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
tanorm = tol*anorm;
}
 
for (i=m-1;i>=0;i--)
for (i=m;i--;)
{
if (i < m-1)
{
if (g)
{
ap0 = a+l*m+i;
vp0 = vmat+i*m+l;
vp10 = vmat+l*m+l;
g *= *ap0;
for (ap=ap0,vp=vp0,j=nml; j--; ap+=m)
*(vp++) = *ap/g;
for (j=nml; j--; vp10+=m)
for (j=l; j<m; j++)
vmat[j*m+i]=(a[i*m+j]/a[i*m+l]) / g;
for (j=l; j<m; j++)
{
for (s=0.0,ap=ap0,vp1=vp10,k=nml; k--; ap+=m)
s += *ap**(vp1++);
for (vp=vp0,vp1=vp10,k=nml; k--;)
*(vp1++) += s**(vp++);
s = 0.0;
for (k=l; k<m; k++)
s += a[i*m+k] * vmat[k*m+j];
for (k=l; k<m; k++)
vmat[k*m+j] += s * vmat[k*m+i];
}
}
vp = vmat+l*m+i;
vp1 = vmat+i*m+l;
for (j=nml; j--; vp+=m)
*vp = *(vp1++) = 0.0;
for (j=l; j<m; j++)
vmat[i*m+j] = vmat[j*m+i] =0.0;
}
vmat[i*m+i]=1.0;
g=rv1[i];
l=i;
nml = m-l;
vmat[i*m+i] = 1.0;
g = rv1[i];
l = i;
}
 
for (i=m; --i>=0;)
for (i=m; i--;)
{
l=i+1;
nml = m-l;
mmi=m-i;
g=wmat[i];
ap0 = a+i*m+i;
ap10 = ap0 + m;
for (ap=ap10,j=nml;j--;ap+=m)
*ap=0.0;
l = i+1;
g = wmat[i];
for (j=l; j<m;j++)
a[i*m+j] = 0.0;
if (g)
{
g=1.0/g;
for (j=nml;j--; ap10+=m)
g = 1.0/g;
for (j=l; j<m; j++)
{
for (s=0.0,ap=ap0,ap1=ap10,k=mmi; --k;)
s += *(++ap)**(++ap1);
f = (s/(*ap0))*g;
for (ap=ap0,ap1=ap10,k=mmi;k--;)
*(ap1++) += f**(ap++);
s = 0.0;
for (k=l; k<m; k++)
s += a[k*m+i] * a[k*m+j];
f = (s / a[i*m+i]) * g;
for (k=i; k<m; k++)
a[k*m+j] += f * a[k*m+i];
}
for (ap=ap0,j=mmi;j--;)
*(ap++) *= g;
for (j=i; j<m; j++)
a[j*m+i] *= g;
}
else
for (ap=ap0,j=mmi;j--;)
*(ap++)=0.0;
++(*ap0);
for (j=i; j<m; j++)
a[j*m+i] = 0.0;
++a[i*m+i];
}
 
for (k=m; --k>=0;)
for (k=m; k--;)
{
for (its=0; its<30; its++)
{
for (its=0;its<100;its++)
flag = 1;
for (l=k+1; l--; )
{
flag=1;
for (l=k;l>=0;l--)
nm=l-1;
if (fabs(rv1[l]) <= tanorm || !l)
{
nm=l-1;
if (!l || fabs(rv1[l]) <= anorm*TOL)
{
flag=0;
break;
}
if (fabs(wmat[nm]) <= anorm*TOL)
break;
flag=0;
break;
}
if (flag)
if (fabs(wmat[nm]) <= tanorm)
break;
}
if (flag)
{
c = 0.0;
s = 1.0;
for (i=l; i<=k; i++)
{
c=0.0;
s=1.0;
ap0 = a+nm*m;
ap10 = a+l*m;
for (i=l; i<=k; i++,ap10+=m)
f = s * rv1[i];
rv1[i] = c * rv1[i];
if (fabs(f) <= tanorm)
break;
g = wmat[i];
h = PYTHAG(f,g);
wmat[i] = h;
h = 1.0 / h;
c = g * h;
s = -f * h;
for (j=0; j<m; j++)
{
f=s*rv1[i];
if (fabs(f) <= anorm*TOL)
break;
g=wmat[i];
h=PYTHAG(f,g);
wmat[i]=h;
h=1.0/h;
c=g*h;
s=(-f*h);
for (ap=ap0,ap1=ap10,j=m; j--;)
{
z = *ap1;
y = *ap;
*(ap++) = y*c+z*s;
*(ap1++) = z*c-y*s;
}
y = a[j*m+nm];
z = a[j*m+i];
a[j*m+nm] = y*c + z*s;
a[j*m+i] = z*c - y*s;
}
}
z=wmat[k];
if (l == k)
}
z = wmat[k];
if (l==k)
{
if (z < 0.0)
{
if (z < 0.0)
{
wmat[k] = -z;
vp = vmat+k*m;
for (j=m; j--; vp++)
*vp = (-*vp);
}
break;
wmat[k] = -z;
for (j=0; j<m; j++)
vmat[j*m+k] = -vmat[j*m+k];
}
x=wmat[l];
nm=k-1;
y=wmat[nm];
g=rv1[nm];
h=rv1[k];
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=PYTHAG(f,1.0);
f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
c=s=1.0;
ap10 = a+l*m;
vp10 = vmat+l*m;
for (j=l;j<=nm;j++,ap10+=m,vp10+=m)
break;
}
if (its == 29)
return 0;
x = wmat[l];
nm = k-1;
y = wmat[nm];
g = rv1[nm];
h = rv1[k];
f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y);
g = PYTHAG(f, 1.0);
f = ((x-z)*(x+z) + h * ((y / (f+SIGN(g,f))) - h)) / x;
c = s = 1.0;
for (j=l; j<=nm; j++)
{
i = j+1;
g = rv1[i];
y = wmat[i];
h = s * g;
g = c * g;
z = PYTHAG(f, h);
rv1[j] = z;
c = f / z;
s = h / z;
f = x*c + g*s;
g = g*c - x*s;
h = y*s;
y *= c;
for (jj=0; jj<m; jj++)
{
i=j+1;
g=rv1[i];
y=wmat[i];
h=s*g;
g=c*g;
z=PYTHAG(f,h);
rv1[j]=z;
c=f/z;
s=h/z;
f=x*c+g*s;
g=g*c-x*s;
h=y*s;
y=y*c;
for (vp=(vp1=vp10)+m,jj=m; jj--;)
{
z = *vp;
x = *vp1;
*(vp1++) = x*c+z*s;
*(vp++) = z*c-x*s;
}
z=PYTHAG(f,h);
wmat[j]=z;
if (z)
{
z=1.0/z;
c=f*z;
s=h*z;
}
f=c*g+s*y;
x=c*y-s*g;
for (ap=(ap1=ap10)+m,jj=m; jj--;)
{
z = *ap;
y = *ap1;
*(ap1++) = y*c+z*s;
*(ap++) = z*c-y*s;
}
x = vmat[jj*m+j];
z = vmat[jj*m+i];
vmat[jj*m+j] = x*c + z*s;
vmat[jj*m+i] = z*c - x*s;
}
rv1[l]=0.0;
rv1[k]=f;
wmat[k]=x;
z = PYTHAG(f, h);
wmat[j] = z;
if (z)
{
z = 1.0 / z;
c = f*z;
s = h*z;
}
f = c*g + s*y;
x = c*y - s*g;
for (jj=0; jj<m; jj++)
{
y = a[jj*m+j];
z = a[jj*m+i];
a[jj*m+j] = y*c + z*s;
a[jj*m+i] = z*c - y*s;
}
}
rv1[l] = 0.0;
rv1[k] = f;
wmat[k] = x;
}
}
 
wmax=0.0;
w = wmat;
for (j=m;j--; w++)
if (*w > wmax)
wmax=*w;
thresh=TOL*wmax;
thresh=tol*wmax;
rank = 0;
w = wmat;
for (j=m;j--; w++)
1121,7 → 1114,7
{
s = 0.0;
for (k=0; k<m; k++)
s += vmat[j+k*m]*a[i+k*m]*wmat[k];
s += vmat[k+j*m]*a[k+i*m]*wmat[k];
b[j+i*m] = s;
}
 
1137,7 → 1130,6
#undef SIGN
#undef MAX
#undef PYTHAG
#undef TOL
/* undefine everything. THIS MUST REMAIN AT THE END OF THE FILE */
#undef LEVMAR_PSEUDOINVERSE
#undef LEVMAR_LUINVERSE
/gpu/src/levmar/Makefile.vc
21,8 → 21,8
CFLAGS=$(CONFIGFLAGS) /I. /MD /W3 /EHsc /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
LIBSRCS=lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c
LIBOBJS=lm.obj Axb.obj misc.obj lmlec.obj lmbc.obj lmblec.obj lmbleic.obj
LIBSRCS=lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c lmbleic.c
DEMOBJS=lmdemo.obj
DEMOSRCS=lmdemo.c
AR=lib /nologo
40,14 → 40,15
lmdemo.exe: $(DEMOBJS) levmar.lib
$(CC) $(DEMOBJS) $(LDFLAGS) /out:lmdemo.exe $(LIBS)
 
lm.obj: lm.c lm_core.c lm.h misc.h compiler.h
Axb.obj: Axb.c Axb_core.c lm.h misc.h
misc.obj: misc.c misc_core.c lm.h misc.h
lmlec.obj: lmlec.c lmlec_core.c lm.h misc.h
lmbc.obj: lmbc.c lmbc_core.c lm.h misc.h compiler.h
lmblec.obj: lmblec.c lmblec_core.c lm.h misc.h
lm.obj: lm.c lm_core.c levmar.h misc.h compiler.h
Axb.obj: Axb.c Axb_core.c levmar.h misc.h
misc.obj: misc.c misc_core.c levmar.h misc.h
lmlec.obj: lmlec.c lmlec_core.c levmar.h misc.h
lmbc.obj: lmbc.c lmbc_core.c levmar.h misc.h compiler.h
lmblec.obj: lmblec.c lmblec_core.c levmar.h misc.h
lmbleic.obj: lmbleic.c lmbleic_core.c levmar.h misc.h
 
lmdemo.obj: lm.h
lmdemo.obj: levmar.h
 
clean:
-del $(LIBOBJS) $(DEMOBJS)
/gpu/src/levmar/lmbleic.c New file
0,0 → 1,89
/////////////////////////////////////////////////////////////////////////////////
//
// Levenberg - Marquardt non-linear minimization algorithm
// Copyright (C) 2009 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.
//
/////////////////////////////////////////////////////////////////////////////////
 
/*******************************************************************************
* Wrappers for linear inequality constrained Levenberg-Marquardt minimization.
* The same core code is used with appropriate #defines to derive single and
* double precision versions, see also lmbleic_core.c
*******************************************************************************/
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
 
#include "levmar.h"
#include "misc.h"
 
 
#ifndef HAVE_LAPACK
 
#ifdef _MSC_VER
#pragma message("Linear inequalities constrained optimization requires LAPACK and was not compiled!")
#else
//#warning Linear inequalities constrained optimization requires LAPACK and was not compiled!
#endif // _MSC_VER
 
#else // LAPACK present
 
#if !defined(LM_DBL_PREC) && !defined(LM_SNGL_PREC)
#error At least one of LM_DBL_PREC, LM_SNGL_PREC should be defined!
#endif
 
 
#ifdef LM_SNGL_PREC
/* single precision (float) definitions */
#define LM_REAL float
#define LM_PREFIX s
 
#define LM_REAL_MAX FLT_MAX
#define LM_REAL_MIN -FLT_MAX
#define __SUBCNST(x) x##F
#define LM_CNST(x) __SUBCNST(x) // force substitution
 
#include "lmbleic_core.c" // read in core code
 
#undef LM_REAL
#undef LM_PREFIX
#undef LM_REAL_MAX
#undef LM_REAL_MIN
#undef __SUBCNST
#undef LM_CNST
#endif /* LM_SNGL_PREC */
 
#ifdef LM_DBL_PREC
/* double precision definitions */
#define LM_REAL double
#define LM_PREFIX d
 
#define LM_REAL_MAX DBL_MAX
#define LM_REAL_MIN -DBL_MAX
#define LM_CNST(x) (x)
 
#include "lmbleic_core.c" // read in core code
 
#undef LM_REAL
#undef LM_PREFIX
#undef LM_REAL_MAX
#undef LM_REAL_MIN
#undef LM_CNST
#endif /* LM_DBL_PREC */
 
#endif /* HAVE_LAPACK */
 
/gpu/src/levmar/lm.c
28,7 → 28,7
#include <math.h>
#include <float.h>
 
#include "lm.h"
#include "levmar.h"
#include "compiler.h"
#include "misc.h"
 
/gpu/src/levmar/Axb.c
28,7 → 28,7
#include <stdlib.h>
#include <math.h>
 
#include "lm.h"
#include "levmar.h"
#include "misc.h"
 
#if !defined(LM_DBL_PREC) && !defined(LM_SNGL_PREC)
/gpu/src/astrom.c
9,7 → 9,7
*
* Contents: Astrometrical computations.
*
* Last modify: 15/12/2009
* Last modify: 20/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
854,7 → 854,9
void astrom_profshapeparam(picstruct *field, objstruct *obj)
{
wcsstruct *wcs;
double dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3, ct,st;
double mat[9], tempmat[9], mx2wcov[9], dpdmx2[6], cov[4],
dx2,dy2,dxy, xm2,ym2,xym, pm2, lm0,lm1,lm2,lm3, ct,st,
temp, invstemp, den, invden, dval;
int lng,lat, naxis;
 
wcs = field->wcs;
1047,9 → 1049,18
dx2 = obj2->prof_mx2;
dy2 = obj2->prof_my2;
dxy = obj2->prof_mxy;
obj2->prof_mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
obj2->prof_my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
obj2->prof_mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
mat[0] = lm0*lm0;
mat[3] = lm1*lm1;
mat[6] = lm0*lm1;
mat[1] = lm2*lm2;
mat[4] = lm3*lm3;
mat[7] = lm2*lm3;
mat[2] = lm0*lm2;
mat[5] = lm1*lm3;
mat[8] = lm0*lm3+lm1*lm2;
obj2->prof_mx2w = xm2 = mat[0]*dx2 + mat[3]*dy2 + mat[6]*dxy;
obj2->prof_my2w = ym2 = mat[1]*dx2 + mat[4]*dy2 + mat[7]*dxy;
obj2->prof_mxyw = xym = mat[2]*dx2 + mat[5]*dy2 + mat[8]*dxy;
temp=xm2-ym2;
if (FLAG(obj2.prof_thetaw))
{
1093,24 → 1104,70
obj2->prof_cxyw = (float)(-2*xym/temp);
}
if (FLAG(obj2.prof_e1w))
/*-- Use the Jacobians to compute the moment covariance matrix */
if (FLAG(obj2.prof_pol1errw) || FLAG(obj2.prof_e1errw))
propagate_covar(obj2->prof_mx2cov, mat, mx2wcov, 3, 3, tempmat);
 
if (FLAG(obj2.prof_pol1w))
{
if (xm2+ym2 > 1.0/BIG)
{
obj2->prof_pol1w = (xm2 - ym2) / (xm2+ym2);
obj2->prof_pol2w = 2.0*xym / (xm2 + ym2);
temp = xm2*ym2-xym*xym;
if (temp>=0.0)
temp = xm2+ym2+2.0*sqrt(temp);
else
temp = xm2+ym2;
obj2->prof_e1w = (xm2 - ym2) / temp;
obj2->prof_e2w = 2.0*xym / temp;
if (FLAG(obj2.prof_pol1errw))
{
/*-------- Compute the Jacobian of polarisation */
invden = 1.0/(xm2+ym2);
dpdmx2[0] = 2.0*ym2*invden*invden;
dpdmx2[1] = -2.0*xm2*invden*invden;
dpdmx2[2] = 0.0;
dpdmx2[3] = -2.0*xym*invden*invden;
dpdmx2[4] = -2.0*xym*invden*invden;
dpdmx2[5] = 2.0*invden;
propagate_covar(mx2wcov, dpdmx2, cov, 3, 2, tempmat);
obj2->prof_pol1errw = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
obj2->prof_pol2errw = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
obj2->prof_pol12corrw = (dval=cov[0]*cov[3]) > 0.0?
(float)(cov[1]/sqrt(dval)) : 0.0;
}
}
else
obj2->prof_pol1w = obj2->prof_pol2w
= obj2->prof_e1w = obj2->prof_e2w = 0.0;
obj2->prof_pol1w = obj2->prof_pol2w = obj2->prof_pol1errw
= obj2->prof_pol2errw = obj2->prof_pol12corrw = 0.0;
}
 
if (FLAG(obj2.prof_e1w))
{
if (xm2+ym2 > 1.0/BIG)
{
temp = xm2*ym2 - xym*xym;
den = (temp>=0.0) ? xm2+ym2+2.0*sqrt(temp) : xm2+ym2;
invden = 1.0/den;
obj2->prof_e1w = (float)(invden*(xm2 - ym2));
obj2->prof_e2w = (float)(2.0 * invden * xym);
if (FLAG(obj2.prof_e1errw))
{
/*------ Compute the Jacobian of ellipticity */
invstemp = (temp>=0.0) ? 1.0/sqrt(temp) : 0.0;
dpdmx2[0] = ( den - (1.0+ym2*invstemp)*(xm2-ym2))*invden*invden;
dpdmx2[1] = (-den - (1.0+xm2*invstemp)*(xm2-ym2))*invden*invden;
dpdmx2[2] = 2.0*xym*invstemp*(xm2-ym2)*invden*invden;
dpdmx2[3] = -2.0*xym*(1.0+ym2*invstemp)*invden*invden;
dpdmx2[4] = -2.0*xym*(1.0+xm2*invstemp)*invden*invden;
dpdmx2[5] = (2.0*den+4.0*xym*xym*invstemp)*invden*invden;
 
/*------ Use the Jacobian to compute the ellipticity covariance matrix */
propagate_covar(mx2wcov, dpdmx2, cov, 3, 2, tempmat);
obj2->prof_e1errw = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
obj2->prof_e2errw = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
obj2->prof_e12corrw = (dval=cov[0]*cov[3]) > 0.0?
(float)(cov[1]/sqrt(dval)) : 0.0;
}
}
else
obj2->prof_e1w = obj2->prof_e2w = obj2->prof_e1errw
= obj2->prof_e2errw = obj2->prof_e12corrw = 0.0;
}
}
 
return;
/gpu/src/types.h
9,7 → 9,7
*
* Contents: global type definitions.
*
* Last modify: 23/03/2010
* Last modify: 20/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
350,12 → 350,17
float poserrcxxw_prof, poserrcyyw_prof,
poserrcxyw_prof; /* WORLD error ellipse */
double prof_mx2, prof_my2, prof_mxy; /* Model-fitting moments */
double prof_mx2cov[9]; /* Mod-fit moment cov. matrix */
float prof_a, prof_b,
prof_theta; /* Model-fitting ellip. params*/
float prof_cxx, prof_cyy,
prof_cxy; /* Model-fitting ellip. params*/
float prof_pol1, prof_pol2; /* Model-fitting pol. vector*/
float prof_pol1err, prof_pol2err,
prof_pol12corr; /* Model-fitting pol. errors */
float prof_e1, prof_e2; /* Model-fitting ellip.vector*/
float prof_e1err, prof_e2err,
prof_e12corr; /* Model-fitting ellip. errors */
double prof_mx2w, prof_my2w,
prof_mxyw; /* WORLD model-fitting moments*/
float prof_aw, prof_bw,
366,7 → 371,11
float prof_cxxw, prof_cyyw,
prof_cxyw; /* WORLD ellipse parameters */
float prof_pol1w, prof_pol2w; /* WORLD polarisation vector*/
float prof_pol1errw, prof_pol2errw,
prof_pol12corrw; /* WORLD polarisation errors */
float prof_e1w, prof_e2w; /* WORLD ellipticity vector*/
float prof_e1errw, prof_e2errw,
prof_e12corrw; /* WORLD ellipticity errors */
float prof_class_star; /* Mod.-fitting star/gal class*/
float prof_concentration; /* Model-fitting concentration*/
float prof_concentrationerr; /* RMS error */
/gpu/src/winpos.c
9,7 → 9,7
*
* Contents: Compute windowed barycenter
*
* Last modify: 19/12/2007
* Last modify: 16/07/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
38,16 → 38,16
OUTPUT -.
NOTES obj->posx and obj->posy are taken as initial centroid guesses.
AUTHOR E. Bertin (IAP)
VERSION 19/12/2007
VERSION 16/07/2010
***/
void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
 
{
float r2, raper,raper2, rintlim,rintlim2,rextlim2,
dx,dx1,dy,dy2, sig, pdbkg,
offsetx,offsety,scalex,scaley,scale2, ngamma, locarea;
double tv, tv2, pix, var, backnoise2, gain, locpix,
dxpos,dypos, twosig2, err,err2, emx2,emy2,emxy,
float r2,invtwosig2, raper,raper2, rintlim,rintlim2,rextlim2,
dx,dx1,dy,dy2, sig, invngamma, pdbkg,
offsetx,offsety,scalex,scaley,scale2, locarea;
double tv, norm, pix, var, backnoise2, invgain, locpix,
dxpos,dypos, err,err2, emx2,emy2,emxy,
esum, temp,temp2, mx2, my2,mxy,pmx2, theta, mx,my,
mx2ph, my2ph;
int i,x,y, x2,y2, xmin,xmax,ymin,ymax, sx,sy, w,h,
70,9 → 70,9
errflag = FLAG(obj2.winposerr_mx2);
momentflag = FLAG(obj2.win_mx2) | FLAG(obj2.winposerr_mx2);
var = backnoise2 = field->backsig*field->backsig;
gain = field->gain;
invgain = field->gain>0.0? 1.0/field->gain : 0.0;
sig = obj2->hl_radius*2.0/2.35; /* From half-FWHM to sigma */
twosig2 = 2.0*sig*sig;
invtwosig2 = 1.0/(2.0*sig*sig);
 
/* Integration radius */
raper = WINPOS_NSIG*sig;
80,12 → 80,12
/* For photographic data */
if (pflag)
{
ngamma = field->ngamma;
pdbkg = exp(obj->dbkg/ngamma);
invngamma = 1.0/field->ngamma;
pdbkg = expf(obj->dbkg*invngamma);
}
else
{
ngamma = 0.0;
invngamma = 0.0;
pdbkg = 0.0;
}
raper2 = raper*raper;
165,7 → 165,7
}
else
locarea = 1.0;
locarea *= exp(-r2/twosig2);
locarea *= expf(-r2*invtwosig2);
/*-------- Here begin tests for pixel and/or weight overflows. Things are a */
/*-------- bit intricated to have it running as fast as possible in the most */
/*-------- common cases */
191,7 → 191,7
}
}
if (pflag)
pix=exp(pix/ngamma);
pix = expf(pix*invngamma);
dx = x - mx;
dy = y - my;
locpix = locarea*pix;
202,13 → 202,13
{
err = var;
if (pflag)
err *= locpix*pix/(ngamma*ngamma);
else if (gain>0.0 && pix>0.0)
err *= locpix*pix*invngamma*invngamma;
else if (invgain>0.0 && pix>0.0)
{
if (gainflag)
err += pix/gain*var/backnoise2;
err += pix*invgain*var/backnoise2;
else
err += pix/gain;
err += pix*invgain;
}
err2 = locarea*locarea*err;
esum += err2;
228,8 → 228,8
 
if (tv>0.0)
{
mx += (dxpos /= tv)*WINPOS_GRADFAC;
my += (dypos /= tv)*WINPOS_GRADFAC;
mx += (dxpos /= tv)*WINPOS_FAC;
my += (dypos /= tv)*WINPOS_FAC;
}
else
break;
252,8 → 252,7
obj2->fluxerr_win = sqrt(esum);
}
temp2=mx2*my2-mxy*mxy;
obj2->win_flag = (tv <= 0.0)*4 + (mx2 < 0.0 || my2 < 0.0)*2
+ (temp2<0.0);
obj2->win_flag = (tv <= 0.0)*4 + (mx2 < 0.0 || my2 < 0.0)*2 + (temp2<0.0);
if (obj2->win_flag)
{
/*--- Negative values: revert to isophotal estimates */
299,12 → 298,12
{
if (errflag)
{
tv2 = tv*tv;
emx2 /= tv2;
emy2 /= tv2;
emxy /= tv2;
norm = WINPOS_FAC*WINPOS_FAC/(tv*tv);
emx2 *= norm;
emy2 *= norm;
emxy *= norm;
/*-- Handle fully correlated profiles (which cause a singularity...) */
esum *= 0.08333/tv2;
esum *= 0.08333*norm;
if (obj->singuflag && (emx2*emy2-emxy*emxy) < esum*esum)
{
emx2 += esum;
/gpu/src/winpos.h
9,7 → 9,7
*
* Contents: Include file for winpos.c.
*
* Last modify: 25/08/2005
* Last modify: 16/07/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
18,9 → 18,9
 
#define WINPOS_NITERMAX 16 /* Maximum number of steps */
#define WINPOS_NSIG 4 /* Measurement radius */
#define WINPOS_OVERSAMP 3 /* oversampling in each dimension */
#define WINPOS_STEPMIN 0.001 /* Minimum change in position for continueing*/
#define WINPOS_GRADFAC 2.0 /* Gradient descent acceleration factor */
#define WINPOS_OVERSAMP 11 /* oversampling in each dimension */
#define WINPOS_STEPMIN 0.0001 /* Minimum change in position for continueing*/
#define WINPOS_FAC 2.0 /* Centroid offset factor (2 for a Gaussian) */
 
/* NOTES:
One must have:
/gpu/src/scan.c
10,7 → 10,7
* Contents: functions for extraction of connected pixels from
* a pixmap.
*
* Last modify: 21/01/2010
* Last modify: 28/09/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
48,7 → 48,7
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 01/12/2009
VERSION 28/09/2010
***/
void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
int nffield, picstruct *wfield, picstruct *dwfield)
101,7 → 101,7
w = cfield->width;
h = cfield->height;
objlist.dthresh = cfield->dthresh;
objlist.thresh = cfield->thresh;
objlist.thresh = field->thresh;
cfield->yblank = 1;
field->y = field->stripy = 0;
field->ymin = field->stripylim = 0;
/gpu/src/analyse.c
9,7 → 9,7
*
* Contents: analyse(), endobject()...: measurements on detections.
*
* Last modify: 15/12/2009
* Last modify: 28/09/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
93,9 → 93,9
pospeakflag, minarea, gainflag;
double tv,sigtv, ngamma,
esum, emx2,emy2,emxy, err,gain,backnoise2,dbacknoise2,
xm,ym, x,y,var,var2;
xm,ym, x,y,var,var2, threshfac;
float *heap,*heapt,*heapj,*heapk, swap;
PIXTYPE pix, cdpix, tpix, peak,cdpeak, thresh,dthresh;
PIXTYPE pix, cdpix, tpix, peak,cdpeak, thresh,dthresh,minthresh;
static PIXTYPE threshs[NISO];
 
 
149,6 → 149,8
peak = -BIG;
cdpeak = -BIG;
thresh = field->thresh;
minthresh = (PLISTEXIST(var))? BIG : thresh;
threshfac = field->backsig > 0.0 ? field->thresh / field->backsig : 1.0;
dthresh = dfield->dthresh;
area = 0;
for (pixt=pixel+obj->firstpix; pixt>=pixel; pixt=pixel+PLIST(pixt,nextpix))
165,7 → 167,13
obj->peaky = PLIST(pixt,y) + 1;
}
if (PLISTEXIST(var))
{
var = PLISTPIX(pixt, var);
thresh = threshfac*sqrt(var);
if (thresh < minthresh)
minthresh = thresh;
}
 
if (photoflag)
{
pix = exp(pix/ngamma);
178,7 → 186,6
var2 += pix/gain*var/backnoise2;
 
sigtv += var2;
 
if (pix>thresh)
area++;
tv += pix;
283,8 → 290,7
obj->flux = tv;
obj->fluxerr = sigtv;
obj->peak = peak;
 
obj->thresh -= obj->dbkg;
obj->thresh = minthresh - obj->dbkg;
obj->peak -= obj->dbkg;
 
/* Initialize isophotal thresholds so as to sample optimally the full profile*/
684,6 → 690,7
}
 
/*----------------------------- Profile fitting -----------------------------*/
#ifdef USE_MODEL
if (prefs.prof_flag)
{
profit_fit(theprofit, field, wfield, obj, obj2);
697,7 → 704,7
if (FLAG(obj2.poserrmx2w_prof))
astrom_proferrparam(field, obj);
}
 
#endif
/*--- Express everything in magnitude units */
computemags(field, obj);
 
/gpu/src/preflist.h
9,7 → 9,7
*
* Contents: Keywords for the configuration file.
*
* Last modify: 18/05/2010
* Last modify: 23/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
70,9 → 70,12
"MINIBACK_RMS", "-BACKGROUND",
"FILTERED", "OBJECTS", "APERTURES", "SEGMENTATION", "ASSOC",
"-OBJECTS", "PSFS", "-PSFS",
"PC_CONVPROTOS", "-PC_CONVPROTOS", "PC_PROTOS",
"MAP_SOM", "MODELS", "-MODELS", "SPHEROIDS", "-SPHEROIDS",
"DISKS", "-DISKS", "PATTERNS", "OTHER", ""},
"PC_CONVPROTOS", "-PC_CONVPROTOS", "PC_PROTOS", "MAP_SOM",
#ifdef USE_MODEL
"MODELS", "-MODELS", "SPHEROIDS", "-SPHEROIDS",
"DISKS", "-DISKS", "PATTERNS",
#endif
"OTHER", ""},
0, 17, &prefs.ncheck_type},
{"CLEAN", P_BOOL, &prefs.clean_flag},
{"CLEAN_PARAM", P_FLOAT, &prefs.clean_param, 0,0, 0.1,10.0},
/gpu/src/ldactoasc.h
30,7 → 30,7
#define MYVERSION VERSION
#endif
#define COPYRIGHT "Emmanuel BERTIN <bertin@iap.fr>"
#define WEBSITE "http://astromatic.iap.fr/software/sextractor/"
#define WEBSITE "http://astromatic.net/software/sextractor/"
#define INSTITUTE "IAP http://www.iap.fr"
 
/*----------------------------- Physical constants --------------------------*/
/gpu/src/Makefile.am
1,15 → 1,23
# Program Makefile for SEx
# Copyright (C) 2004-2010 Emmanuel Bertin.
SUBDIRS = fits levmar wcs
if USE_MODEL
FFTSOURCE = fft.c
PATTERNSOURCE = pattern.c
PROFITSOURCE = profit.c
LEVLIB = $(top_builddir)/src/levmar/liblevmar.a
LEVDIR = levmar
endif
SUBDIRS = fits $(LEVDIR) wcs
bin_PROGRAMS = sex ldactoasc
check_PROGRAMS = sex
sex_SOURCES = analyse.c assoc.c astrom.c back.c bpro.c catout.c \
check.c clean.c extract.c fft.c field.c filter.c \
fitswcs.c flag.c graph.c growth.c header.c image.c \
interpolate.c main.c makeit.c manobjlist.c misc.c \
neurro.c pattern.c pc.c photom.c plist.c prefs.c \
profit.c psf.c readimage.c refine.c retina.c scan.c \
som.c weight.c winpos.c xml.c \
check.c clean.c extract.c $(FFTSOURCE) field.c \
filter.c fitswcs.c flag.c graph.c growth.c header.c \
image.c interpolate.c main.c makeit.c manobjlist.c \
misc.c neurro.c $(PATTERNSOURCE) pc.c photom.c \
plist.c prefs.c $(PROFITSOURCE) psf.c readimage.c \
refine.c retina.c scan.c som.c weight.c winpos.c \
xml.c \
assoc.h astrom.h back.h bpro.h check.h clean.h \
define.h extract.h fft.h field.h filter.h fitswcs.h \
flag.h globals.h growth.h header.h image.h \
21,7 → 29,7
ldactoasc_SOURCES = ldactoasc.c ldactoasc.h
sex_LDADD = $(top_builddir)/src/fits/libfits.a \
$(top_builddir)/src/wcs/libwcs_c.a \
$(top_builddir)/src/levmar/liblevmar.a
$(LEVLIB)
ldactoasc_LDADD = $(top_builddir)/src/fits/libfits.a
DATE=`date +"%Y-%m-%d"`
 
/gpu/src/paramprofit.h
9,7 → 9,7
*
* Contents: Model-fitting parameter list for catalog data.
*
* Last modify: 25/05/2010
* Last modify: 20/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
197,6 → 197,26
&outobj2.prof_pol2, H_FLOAT, T_FLOAT, "%10.6f", "",
"src.ellipticity;stat.fit;instr.det", ""},
 
{"ELLIP1ERRMODEL_IMAGE", "Ellipticity component std.error from model-fitting",
&outobj2.prof_e1err, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit;instr.det", ""},
{"ELLIP2ERRMODEL_IMAGE", "Ellipticity component std.error from model-fitting",
&outobj2.prof_e2err, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit;instr.det", ""},
{"ELLIPCORRMODEL_IMAGE", "Corr.coeff between ellip.components from model-fitting",
&outobj2.prof_e12corr, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.correlation;src.ellipticity;stat.fit;instr.det", ""},
 
{"POLAR1ERRMODEL_IMAGE", "Polarisation component std.error from model-fitting",
&outobj2.prof_pol1err, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit;instr.det", ""},
{"POLAR2ERRMODEL_IMAGE", "Polarisation component std.error from model-fitting",
&outobj2.prof_pol2err, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit;instr.det", ""},
{"POLARCORRMODEL_IMAGE", "Corr.coeff between polar. components from fitting",
&outobj2.prof_pol12corr, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.correlation;src.ellipticity;stat.fit;instr.det", ""},
 
{"X2MODEL_WORLD", "Variance along X-WORLD (alpha) from model-fitting",
&outobj2.prof_mx2w, H_EXPO, T_DOUBLE, "%18.10e", "deg**2",
"src.impactParam;stat.fit", "deg2"},
212,13 → 232,33
{"ELLIP2MODEL_WORLD", "Ellipticity component from model-fitting",
&outobj2.prof_e2w, H_FLOAT, T_FLOAT, "%10.6f", "",
"src.ellipticity;stat.fit", ""},
{"POLAR1MODEL_WORLD", "Ellipticity component (quadratic) from model-fitting",
{"POLAR1MODEL_WORLD", "Polarisation component from model-fitting",
&outobj2.prof_pol1w, H_FLOAT, T_FLOAT, "%10.6f", "",
"src.ellipticity;stat.fit", ""},
{"POLAR2MODEL_WORLD", "Ellipticity component (quadratic) from model-fitting",
{"POLAR2MODEL_WORLD", "Polarisation component from model-fitting",
&outobj2.prof_pol2w, H_FLOAT, T_FLOAT, "%10.6f", "",
"src.ellipticity;stat.fit", ""},
 
{"ELLIP1ERRMODEL_WORLD", "Ellipticity component std.error from model-fitting",
&outobj2.prof_e1errw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit", ""},
{"ELLIP2ERRMODEL_WORLD", "Ellipticity component std.error from model-fitting",
&outobj2.prof_e2errw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit", ""},
{"ELLIPCORRMODEL_WORLD", "Corr.coeff between ellip.components from model-fitting",
&outobj2.prof_e12corrw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.correlation;src.ellipticity;stat.fit", ""},
 
{"POLAR1ERRMODEL_WORLD", "Polarisation component std.error from model-fitting",
&outobj2.prof_pol1errw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit", ""},
{"POLAR2ERRMODEL_WORLD", "Polarisation component std.error from model-fitting",
&outobj2.prof_pol2errw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit", ""},
{"POLARCORRMODEL_WORLD", "Corr.coeff between polar. components from fitting",
&outobj2.prof_pol12corrw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.correlation;src.ellipticity;stat.fit", ""},
 
{"CXXMODEL_IMAGE", "Cxx ellipse parameter from model-fitting",
&outobj2.prof_cxx, H_EXPO, T_FLOAT, "%15.7e", "pixel**(-2)",
"src.impactParam;stat.fit;instr.det", "pix-2"},
/gpu/src/catout.c
9,7 → 9,7
*
* Contents: functions for output of catalog data.
*
* Last modify: 08/07/2010
* Last modify: 20/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
233,8 → 233,11
FLAG(obj2.prof_bar_aspect) |= FLAG(obj2.prof_bar_aspecterr);
FLAG(obj2.prof_bar_theta) |= FLAG(obj2.prof_bar_thetaerr);
FLAG(obj2.prof_arms_mag) |= FLAG(obj2.prof_arms_magerr);
FLAG(obj2.prof_e1w) |= FLAG(obj2.prof_e2w)
| FLAG(obj2.prof_pol1w) | FLAG(obj2.prof_pol2w);
FLAG(obj2.prof_e1errw) |= FLAG(obj2.prof_e2errw) | FLAG(obj2.prof_e12corrw);
FLAG(obj2.prof_e1w) |= FLAG(obj2.prof_e2w) | FLAG(obj2.prof_e1errw);
FLAG(obj2.prof_pol1errw) |= FLAG(obj2.prof_pol2errw)
| FLAG(obj2.prof_pol12corrw);
FLAG(obj2.prof_pol1w) |= FLAG(obj2.prof_pol2w) | FLAG(obj2.prof_pol1errw);
FLAG(obj2.prof_aw) |= FLAG(obj2.prof_bw);
FLAG(obj2.prof_cxxw) |= FLAG(obj2.prof_cyyw) | FLAG(obj2.prof_cxyw);
FLAG(obj2.prof_thetas) |= FLAG(obj2.prof_theta1950)
245,7 → 248,7
| FLAG(obj2.prof_mxyw)
| FLAG(obj2.prof_thetaw) | FLAG(obj2.prof_aw)
| FLAG(obj2.prof_cxxw)
| FLAG(obj2.prof_e1w);
| FLAG(obj2.prof_e1w) | FLAG(obj2.prof_pol1w);
 
FLAG(obj2.dtheta1950) |= FLAG(obj2.prof_theta1950)
| FLAG(obj2.prof_spheroid_theta1950)
313,10 → 316,6
| FLAG(obj2.prof_spheroid_aspecterrw)
| FLAG(obj2.prof_spheroid_thetaerrw);
 
FLAG(obj2.prof_mx2w) |= FLAG(obj2.prof_my2w) | FLAG(obj2.prof_mxyw)
| FLAG(obj2.prof_e1w) |FLAG(obj2.prof_e2w)
| FLAG(obj2.prof_pol1w) |FLAG(obj2.prof_pol2w);
 
FLAG(obj2.prof_spheroid_fluxmean) |= FLAG(obj2.prof_spheroid_mumean);
FLAG(obj2.prof_spheroid_fluxeff) |= FLAG(obj2.prof_spheroid_mueff);
FLAG(obj2.prof_spheroid_peak) |= FLAG(obj2.prof_spheroid_mumax)
328,16 → 327,22
| FLAG(obj2.prof_bar_lengthw)
| FLAG(obj2.prof_arms_scalew)
| FLAG(obj2.prof_mx2w);
 
FLAG(obj2.prof_e1err) |= FLAG(obj2.prof_e2err) | FLAG(obj2.prof_e12corr)
| FLAG(obj2.prof_e1errw);
FLAG(obj2.prof_e1) |= FLAG(obj2.prof_e2)
| FLAG(obj2.prof_pol1) | FLAG(obj2.prof_pol2)
| FLAG(obj2.prof_e1w);
| FLAG(obj2.prof_e1err) | FLAG(obj2.prof_e1w);
FLAG(obj2.prof_pol1err) |= FLAG(obj2.prof_pol2err)
| FLAG(obj2.prof_pol12corr) | FLAG(obj2.prof_pol1errw);
FLAG(obj2.prof_pol1) |= FLAG(obj2.prof_pol2)
| FLAG(obj2.prof_pol1err) | FLAG(obj2.prof_pol1w);
FLAG(obj2.prof_a) |= FLAG(obj2.prof_b) | FLAG(obj2.prof_theta)
| FLAG(obj2.prof_aw);
FLAG(obj2.prof_cxx) |= FLAG(obj2.prof_cyy)
| FLAG(obj2.prof_cxy) | FLAG(obj2.prof_cxxw);
FLAG(obj2.prof_mx2) |= FLAG(obj2.prof_my2) | FLAG(obj2.prof_mxy)
| FLAG(obj2.prof_e1)
| FLAG(obj2.prof_e1) | FLAG(obj2.prof_pol1)
| FLAG(obj2.prof_a) | FLAG(obj2.prof_cxx)
| FLAG(obj2.prof_mx2w);
FLAG(obj2.fluxmean_prof) |= FLAG(obj2.mumean_prof);
439,7 → 444,8
FLAG(obj2.winposerr_mx2) |= FLAG(obj2.winposerr_my2)
| FLAG(obj2.winposerr_mxy)
| FLAG(obj2.winposerr_a) | FLAG(obj2.winposerr_cxx)
| FLAG(obj2.winposerr_mx2w);
| FLAG(obj2.winposerr_mx2w)
| FLAG(obj2.fluxerr_win) | FLAG(obj2.magerr_win);
 
FLAG(obj2.winpos_alpha1950) |= FLAG(obj2.winpos_delta1950)
| FLAG(obj2.win_theta1950)
/gpu/src/Makefile.in
53,22 → 53,36
am_ldactoasc_OBJECTS = ldactoasc.$(OBJEXT)
ldactoasc_OBJECTS = $(am_ldactoasc_OBJECTS)
ldactoasc_DEPENDENCIES = $(top_builddir)/src/fits/libfits.a
am__sex_SOURCES_DIST = analyse.c assoc.c astrom.c back.c bpro.c \
catout.c check.c clean.c extract.c fft.c field.c filter.c \
fitswcs.c flag.c graph.c growth.c header.c image.c \
interpolate.c main.c makeit.c manobjlist.c misc.c neurro.c \
pattern.c pc.c photom.c plist.c prefs.c profit.c psf.c \
readimage.c refine.c retina.c scan.c som.c weight.c winpos.c \
xml.c assoc.h astrom.h back.h bpro.h check.h clean.h define.h \
extract.h fft.h field.h filter.h fitswcs.h flag.h globals.h \
growth.h header.h image.h interpolate.h key.h neurro.h param.h \
paramprofit.h pattern.h photom.h plist.h prefs.h preflist.h \
profit.h psf.h retina.h sexhead1.h sexhead.h sexheadsc.h som.h \
threads.h types.h wcscelsys.h weight.h winpos.h xml.h
@USE_MODEL_TRUE@am__objects_1 = fft.$(OBJEXT)
@USE_MODEL_TRUE@am__objects_2 = pattern.$(OBJEXT)
@USE_MODEL_TRUE@am__objects_3 = profit.$(OBJEXT)
am_sex_OBJECTS = analyse.$(OBJEXT) assoc.$(OBJEXT) astrom.$(OBJEXT) \
back.$(OBJEXT) bpro.$(OBJEXT) catout.$(OBJEXT) check.$(OBJEXT) \
clean.$(OBJEXT) extract.$(OBJEXT) fft.$(OBJEXT) \
clean.$(OBJEXT) extract.$(OBJEXT) $(am__objects_1) \
field.$(OBJEXT) filter.$(OBJEXT) fitswcs.$(OBJEXT) \
flag.$(OBJEXT) graph.$(OBJEXT) growth.$(OBJEXT) \
header.$(OBJEXT) image.$(OBJEXT) interpolate.$(OBJEXT) \
main.$(OBJEXT) makeit.$(OBJEXT) manobjlist.$(OBJEXT) \
misc.$(OBJEXT) neurro.$(OBJEXT) pattern.$(OBJEXT) pc.$(OBJEXT) \
misc.$(OBJEXT) neurro.$(OBJEXT) $(am__objects_2) pc.$(OBJEXT) \
photom.$(OBJEXT) plist.$(OBJEXT) prefs.$(OBJEXT) \
profit.$(OBJEXT) psf.$(OBJEXT) readimage.$(OBJEXT) \
$(am__objects_3) psf.$(OBJEXT) readimage.$(OBJEXT) \
refine.$(OBJEXT) retina.$(OBJEXT) scan.$(OBJEXT) som.$(OBJEXT) \
weight.$(OBJEXT) winpos.$(OBJEXT) xml.$(OBJEXT)
sex_OBJECTS = $(am_sex_OBJECTS)
sex_DEPENDENCIES = $(top_builddir)/src/fits/libfits.a \
$(top_builddir)/src/wcs/libwcs_c.a \
$(top_builddir)/src/levmar/liblevmar.a
$(top_builddir)/src/wcs/libwcs_c.a $(LEVLIB)
DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir)
depcomp = $(SHELL) $(top_srcdir)/autoconf/depcomp
am__depfiles_maybe = depfiles
82,7 → 96,7
--mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
$(LDFLAGS) -o $@
SOURCES = $(ldactoasc_SOURCES) $(sex_SOURCES)
DIST_SOURCES = $(ldactoasc_SOURCES) $(sex_SOURCES)
DIST_SOURCES = $(ldactoasc_SOURCES) $(am__sex_SOURCES_DIST)
RECURSIVE_TARGETS = all-recursive check-recursive dvi-recursive \
html-recursive info-recursive install-data-recursive \
install-dvi-recursive install-exec-recursive \
94,7 → 108,7
distclean-recursive maintainer-clean-recursive
ETAGS = etags
CTAGS = ctags
DIST_SUBDIRS = $(SUBDIRS)
DIST_SUBDIRS = fits levmar wcs
DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
ACLOCAL = @ACLOCAL@
AMTAR = @AMTAR@
223,14 → 237,20
 
# Program Makefile for SEx
# Copyright (C) 2004-2010 Emmanuel Bertin.
SUBDIRS = fits levmar wcs
@USE_MODEL_TRUE@FFTSOURCE = fft.c
@USE_MODEL_TRUE@PATTERNSOURCE = pattern.c
@USE_MODEL_TRUE@PROFITSOURCE = profit.c
@USE_MODEL_TRUE@LEVLIB = $(top_builddir)/src/levmar/liblevmar.a
@USE_MODEL_TRUE@LEVDIR = levmar
SUBDIRS = fits $(LEVDIR) wcs
sex_SOURCES = analyse.c assoc.c astrom.c back.c bpro.c catout.c \
check.c clean.c extract.c fft.c field.c filter.c \
fitswcs.c flag.c graph.c growth.c header.c image.c \
interpolate.c main.c makeit.c manobjlist.c misc.c \
neurro.c pattern.c pc.c photom.c plist.c prefs.c \
profit.c psf.c readimage.c refine.c retina.c scan.c \
som.c weight.c winpos.c xml.c \
check.c clean.c extract.c $(FFTSOURCE) field.c \
filter.c fitswcs.c flag.c graph.c growth.c header.c \
image.c interpolate.c main.c makeit.c manobjlist.c \
misc.c neurro.c $(PATTERNSOURCE) pc.c photom.c \
plist.c prefs.c $(PROFITSOURCE) psf.c readimage.c \
refine.c retina.c scan.c som.c weight.c winpos.c \
xml.c \
assoc.h astrom.h back.h bpro.h check.h clean.h \
define.h extract.h fft.h field.h filter.h fitswcs.h \
flag.h globals.h growth.h header.h image.h \
243,7 → 263,7
ldactoasc_SOURCES = ldactoasc.c ldactoasc.h
sex_LDADD = $(top_builddir)/src/fits/libfits.a \
$(top_builddir)/src/wcs/libwcs_c.a \
$(top_builddir)/src/levmar/liblevmar.a
$(LEVLIB)
 
ldactoasc_LDADD = $(top_builddir)/src/fits/libfits.a
DATE = `date +"%Y-%m-%d"`
/gpu/src/profit.c
9,7 → 9,7
*
* Contents: Fit an arbitrary profile combination to a detection.
*
* Last modify: 07/07/2010
* Last modify: 26/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
31,20 → 31,15
#include "globals.h"
#include "prefs.h"
#include "fits/fitscat.h"
#include "levmar/lm.h"
#include "levmar/levmar.h"
#include "fft.h"
#include "fitswcs.h"
#include "check.h"
#include "image.h"
#include "pattern.h"
#include "psf.h"
#include "profit.h"
 
#define INTERPW 6 /* Interpolation function range (x) */
#define INTERPH 6 /* Interpolation function range (y) */
 
#define INTERPF(x) (x==0.0?1.0:sinf(PI*x)*sinf(PI*x/3.0)/(PI*PI/3.0*x*x))
/* Lanczos approximation */
 
static double prof_gammainc(double x, double a),
prof_gamma(double x);
static float prof_interpolate(profstruct *prof, float *posin);
199,14 → 194,13
fit).
NOTES It is a modified version of the lm_minimize() of lmfit.
AUTHOR E. Bertin (IAP)
VERSION 08/03/2010
VERSION 26/08/2010
***/
void profit_fit(profitstruct *profit,
picstruct *field, picstruct *wfield,
objstruct *obj, obj2struct *obj2)
{
profitstruct pprofit;
profitstruct hdprofit;
patternstruct *pattern;
psfstruct *psf;
checkstruct *check;
235,9 → 229,9
profit->ix = (int)(obj->mx + 0.49999);/* internal convention: 1st pix = 0 */
profit->iy = (int)(obj->my + 0.49999);/* internal convention: 1st pix = 0 */
psf_fwhm = psf->masksize[0]*psf->pixstep;
profit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499)
profit->objnaxisn[0] = (((int)(obj2->hl_radius*6.0 + psf_fwhm + 0.499)
*1.2)/2)*2 + 1;
profit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499)
profit->objnaxisn[1] = (((int)(obj2->hl_radius*6.0 + psf_fwhm + 0.499)
*1.2)/2)*2 + 1;
if (profit->objnaxisn[1]<profit->objnaxisn[0])
profit->objnaxisn[1] = profit->objnaxisn[0];
1335,7 → 1329,7
OUTPUT RETURN_ERROR if the rasters don't overlap, RETURN_OK otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 01/05/2010
VERSION 12/09/2010
***/
int profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
float factor)
1344,7 → 1338,7
float *pixin,*pixin0, *mask,*maskt, *pixinout, *dpixin,*dpixin0,
*dpixout,*dpixout0, *dx,*dy,
xcin,xcout,ycin,ycout, xsin,ysin, xin,yin, x,y, dxm,dym, val,
invpixstep;
invpixstep, norm;
int *start,*startt, *nmask,*nmaskt,
i,j,k,n,t,
ixsout,iysout, ixout,iyout, dixout,diyout, nxout,nyout,
1404,10 → 1398,10
 
/* Set the yrange for the x-resampling with some margin for interpolation */
iysina = (int)ysin; /* Int. part of Input start y-coord with margin */
hmh = INTERPH/2 - 1; /* Interpolant start */
hmh = INTERPW/2 - 1; /* Interpolant start */
if (iysina<0 || ((iysina -= hmh)< 0))
iysina = 0;
nyin = (int)(ysin+nyout*invpixstep)+INTERPH-hmh;/* Interpolated Input y size*/
nyin = (int)(ysin+nyout*invpixstep)+INTERPW-hmh;/* Interpolated Input y size*/
if (nyin>profit->modnaxisn[1]) /* with margin */
nyin = profit->modnaxisn[1];
/* Express everything relative to the effective Input start (with margin) */
1440,8 → 1434,13
n=t;
*(startt++) = ix;
*(nmaskt++) = n;
norm = 0.0;
for (x=dxm, i=n; i--; x+=1.0)
*(maskt++) = INTERPF(x);
norm += (*(maskt++) = INTERPF(x));
norm = norm>0.0? 1.0/norm : 1.0;
maskt -= n;
for (i=n; i--;)
*(maskt++) *= norm;
}
 
QCALLOC(pixinout, float, nxout*nyin); /* Intermediary frame-buffer */
1466,12 → 1465,12
}
 
/* Reallocate interpolant stuff for the y direction */
QREALLOC(mask, float, nyout*INTERPH); /* Interpolation masks */
QREALLOC(mask, float, nyout*INTERPW); /* Interpolation masks */
QREALLOC(nmask, int, nyout); /* Interpolation mask sizes */
QREALLOC(start, int, nyout); /* Int. part of Input conv starts */
 
/* Compute the local interpolant and data starting points in y */
hmh = INTERPH/2 - 1;
hmh = INTERPW/2 - 1;
yin = ysin;
maskt = mask;
nmaskt = nmask;
1482,18 → 1481,23
dym = iyin - yin - hmh; /* starting point in the interpolation func */
if (iy < 0)
{
n = INTERPH+iy;
n = INTERPW+iy;
dym -= (float)iy;
iy = 0;
}
else
n = INTERPH;
n = INTERPW;
if (n>(t=nyin-iy))
n=t;
*(startt++) = iy;
*(nmaskt++) = n;
norm = 0.0;
for (y=dym, i=n; i--; y+=1.0)
*(maskt++) = INTERPF(y);
norm += (*(maskt++) = INTERPF(y));
norm = norm>0.0? 1.0/norm : 1.0;
maskt -= n;
for (i=n; i--;)
*(maskt++) *= norm;
}
 
/* Initialize destination buffer to zero */
2027,13 → 2031,20
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 07/07/2010
VERSION 20/08/2010
***/
void profit_moments(profitstruct *profit, obj2struct *obj2)
{
profstruct *prof;
double m0, mx2,my2,mxy, den, temp,temp2,pmx2,theta;
int p;
double dpdmx2[6], cov[4],
*jac,*jact, *pjac,*pjact, *dcovar,*dcovart,
*dmx2,*dmy2,*dmxy,
m0,invm0, mx2,my2,mxy, den,invden,
temp, temp2,invtemp2,invstemp2,
pmx2,theta, flux, dval;
float *covart;
int findex[PROF_NPROF],
i,j,p, nparam;
 
/* hw = (float)(profit->modnaxisn[0]/2);*/
/* hh = (float)(profit->modnaxisn[1]/2);*/
2068,19 → 2079,60
/* obj2->prof_my2 = my2 = my2/sum - my*my;*/
/* obj2->prof_mxy = mxy = mxy/sum - mx*my;*/
 
nparam = profit->nparam;
if (FLAG(obj2.prof_e1err) || FLAG(obj2.prof_pol1err))
{
/*-- Set up Jacobian matrices */
QCALLOC(jac, double, nparam*3);
QMALLOC(pjac, double, (nparam<2? 6 : nparam*3));
QMALLOC(dcovar, double, nparam*nparam);
dcovart = dcovar;
covart = profit->covar;
for (i=nparam*nparam; i--;)
*(dcovart++) = (double)(*(covart++));
dmx2 = jac;
dmy2 = jac+nparam;
dmxy = jac+2*nparam;
}
else
jac = pjac = dcovar = dmx2 = dmy2 = dmxy = NULL;
 
m0 = mx2 = my2 = mxy = 0.0;
for (p=0; p<profit->nprof; p++)
{
prof = profit->prof[p];
prof_moments(profit, prof);
m0 += *prof->flux;
mx2 += prof->mx2**prof->flux;
my2 += prof->my2**prof->flux;
mxy += prof->mxy**prof->flux;
findex[p] = prof_moments(profit, prof, pjac);
flux = *prof->flux;
m0 += flux;
mx2 += prof->mx2*flux;
my2 += prof->my2*flux;
mxy += prof->mxy*flux;
if (jac)
{
jact = jac;
pjact = pjac;
for (j=nparam*3; j--;)
*(jact++) += flux * *(pjact++);
}
}
obj2->prof_mx2 = mx2 / m0;
obj2->prof_my2 = my2 / m0;
obj2->prof_mxy = mxy / m0;
invm0 = 1.0 / m0;
obj2->prof_mx2 = (mx2 *= invm0);
obj2->prof_my2 = (my2 *= invm0);
obj2->prof_mxy = (mxy *= invm0);
/* Complete the flux derivative of moments */
if (jac)
{
for (p=0; p<profit->nprof; p++)
{
prof = profit->prof[p];
dmx2[findex[p]] = prof->mx2 - mx2;
dmy2[findex[p]] = prof->my2 - my2;
dmxy[findex[p]] = prof->mxy - mxy;
}
jact = jac;
for (j=nparam*3; j--;)
*(jact++) *= invm0;
}
 
/* Handle fully correlated profiles (which cause a singularity...) */
if ((temp2=mx2*my2-mxy*mxy)<0.00694)
2090,28 → 2142,83
temp2 = mx2*my2-mxy*mxy;
}
 
if (FLAG(obj2.prof_e1))
/* Use the Jacobians to compute the moment covariance matrix */
if (jac)
propagate_covar(dcovar, jac, obj2->prof_mx2cov, nparam, 3,
pjac); /* We re-use pjac */
 
if (FLAG(obj2.prof_pol1))
{
/*--- "Polarisation", i.e. module = (a^2-b^2)/(a^2+b^2) */
if (mx2+my2 > 1.0/BIG)
{
obj2->prof_pol1 = (mx2 - my2) / (mx2+my2);
obj2->prof_pol2 = 2.0*mxy / (mx2 + my2);
if (temp2>=0.0)
den = mx2+my2+2.0*sqrt(temp2);
else
den = mx2+my2;
obj2->prof_e1 = (mx2 - my2) / den;
obj2->prof_e2 = 2.0*mxy / den;
if (FLAG(obj2.prof_pol1err))
{
/*------ Compute the Jacobian of polarisation */
invden = 1.0/(mx2+my2);
dpdmx2[0] = 2.0*my2*invden*invden;
dpdmx2[1] = -2.0*mx2*invden*invden;
dpdmx2[2] = 0.0;
dpdmx2[3] = -2.0*mxy*invden*invden;
dpdmx2[4] = -2.0*mxy*invden*invden;
dpdmx2[5] = 2.0*invden;
 
/*------ Use the Jacobian to compute the polarisation covariance matrix */
propagate_covar(obj2->prof_mx2cov, dpdmx2, cov, 3, 2,
pjac); /* We re-use pjac */
obj2->prof_pol1err = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
obj2->prof_pol2err = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
obj2->prof_pol12corr = (dval=cov[0]*cov[3]) > 0.0?
(float)(cov[1]/sqrt(dval)) : 0.0;
}
}
else
obj2->prof_pol1 = obj2->prof_pol2 = obj2->prof_e1 = obj2->prof_e2 = 0.0;
obj2->prof_pol1 = obj2->prof_pol2
= obj2->prof_pol1err = obj2->prof_pol2err = obj2->prof_pol12corr = 0.0;
}
 
if (FLAG(obj2.prof_e1))
{
/*--- "Ellipticity", i.e. module = (a-b)/(a+b) */
if (mx2+my2 > 1.0/BIG)
{
den = (temp2>=0.0) ? mx2+my2+2.0*sqrt(temp2) : mx2+my2;
invden = 1.0/den;
obj2->prof_e1 = (float)(invden * (mx2 - my2));
obj2->prof_e2 = (float)(2.0 * invden * mxy);
if (FLAG(obj2.prof_e1err))
{
/*------ Compute the Jacobian of ellipticity */
invstemp2 = (temp2>=0.0) ? 1.0/sqrt(temp2) : 0.0;
dpdmx2[0] = ( den - (1.0+my2*invstemp2)*(mx2-my2))*invden*invden;
dpdmx2[1] = (-den - (1.0+mx2*invstemp2)*(mx2-my2))*invden*invden;
dpdmx2[2] = 2.0*mxy*invstemp2*(mx2-my2)*invden*invden;
dpdmx2[3] = -2.0*mxy*(1.0+my2*invstemp2)*invden*invden;
dpdmx2[4] = -2.0*mxy*(1.0+mx2*invstemp2)*invden*invden;
dpdmx2[5] = (2.0*den+4.0*mxy*mxy*invstemp2)*invden*invden;
 
/*------ Use the Jacobian to compute the ellipticity covariance matrix */
propagate_covar(obj2->prof_mx2cov, dpdmx2, cov, 3, 2,
pjac); /* We re-use pjac */
obj2->prof_e1err = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
obj2->prof_e2err = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
obj2->prof_e12corr = (dval=cov[0]*cov[3]) > 0.0?
(float)(cov[1]/sqrt(dval)) : 0.0;
}
}
else
obj2->prof_e1 = obj2->prof_e2
= obj2->prof_e1err = obj2->prof_e2err = obj2->prof_e12corr = 0.0;
}
 
if (FLAG(obj2.prof_cxx))
{
obj2->prof_cxx = (float)(my2/temp2);
obj2->prof_cyy = (float)(mx2/temp2);
obj2->prof_cxy = (float)(-2*mxy/temp2);
invtemp2 = (temp2>=0.0) ? 1.0/temp2 : 0.0;
obj2->prof_cxx = (float)(my2*invtemp2);
obj2->prof_cyy = (float)(mx2*invtemp2);
obj2->prof_cxy = (float)(-2*mxy*invtemp2);
}
 
if (FLAG(obj2.prof_a))
2128,6 → 2235,11
obj2->prof_theta = theta*180.0/PI;
}
 
/* Free memory used by Jacobians */
free(jac);
free(pjac);
free(dcovar);
 
return;
}
 
2342,7 → 2454,7
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 02/07/2010
VERSION 26/08/2010
***/
void profit_resetparam(profitstruct *profit, paramenum paramtype)
{
2523,6 → 2635,8
 
if (parammin!=parammax && (param<=parammin || param>=parammax))
param = (parammin+parammax)/2.0;
else if (parammin==0.0 && parammax==0.0)
parammax = 1.0;
profit_setparam(profit, paramtype, param, parammin, parammax);
 
return;
2695,13 → 2809,12
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 08/07/2010
VERSION 27/07/2010
***/
void profit_covarunboundtobound(profitstruct *profit,
double *dcovar, float *covar)
{
double *dxdy,
dxmin,dxmax;
double *dxdy;
float *x,*xmin,*xmax;
int *fflag,
f,f1,f2, nfree, p,p1,p2, nparam;
2718,12 → 2831,7
if (fflag[p])
{
if (xmin[p]!=xmax[p])
{
dxmin = x[p] - xmin[p];
dxmax= xmax[p] - x[p];
dxdy[f++] = (fabs(dxmin) < 1.0/BIG && fabs(dxmax) < 1.0/BIG) ?
0.0 : dxmin*dxmax/(dxmin+dxmax);
}
dxdy[f++] = (x[p] - xmin[p]) * (xmax[p] - x[p]) / (xmax[p] - xmin[p]);
else
dxdy[f++] = xmax[p];
}
3546,61 → 3654,167
 
 
/****** prof_moments **********************************************************
PROTO void prof_moments(profitstruct *profit, profstruct *prof)
PURPOSE Computes (analytically or numerically) the 2nd moments of a profile
model.
PROTO int prof_moments(profitstruct *profit, profstruct *prof)
PURPOSE Computes (analytically or numerically) the 2nd moments of a profile.
INPUT Profile-fitting structure,
profile structure.
OUTPUT -.
NOTES -.
profile structure,
optional pointer to 3xnparam Jacobian matrix.
OUTPUT Index to the profile flux for further processing.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 08/07/2010
VERSION 20/08/2010
***/
void prof_moments(profitstruct *profit, profstruct *prof)
int prof_moments(profitstruct *profit, profstruct *prof, double *jac)
{
double m20,m02, ct,st, bn, n;
double *dmx2,*dmy2,*dmxy,
m20, a2, ct,st, mx2fac, my2fac,mxyfac, dc2,ds2,dcs,
bn,bn2, n,n2, nfac,nfac2, hscale2, dmdn;
int nparam, index;
 
if (jac)
/*-- Clear output Jacobian */
{
nparam = profit->nparam;
memset(jac, 0, nparam*3*sizeof(double));
dmx2 = jac;
dmy2 = jac + nparam;
dmxy = jac + 2*nparam;
}
else
dmx2 = dmy2 = dmxy = NULL; /* To avoid gcc -Wall warnings */
 
m20 = 0.0; /* to avoid gcc -Wall warnings */
index = 0; /* to avoid gcc -Wall warnings */
 
 
if (prof->posangle)
{
a2 = *prof->aspect**prof->aspect;
ct = cos(*prof->posangle*DEG);
st = sin(*prof->posangle*DEG);
mx2fac = ct*ct + st*st*a2;
my2fac = st*st + ct*ct*a2;
mxyfac = ct*st * (1.0 - a2);
if (jac)
{
dc2 = -2.0*ct*st*DEG;
ds2 = 2.0*ct*st*DEG;
dcs = (ct*ct - st*st)*DEG;
}
else
dc2 = ds2 = dcs = 0.0; /* To avoid gcc -Wall warnings */
switch(prof->code)
{
case PROF_SERSIC:
n = fabs(*prof->extra[0]);
bn = 2.0*n - 1.0/3.0 + 4.0/(405.0*n) + 46.0/(25515.0*n*n)
+ 131.0/(1148175*n*n*n); /* Ciotti & Bertin 1999 */
m20 = 0.5 * *prof->scale**prof->scale
* prof_gamma(4.0*n) / (prof_gamma(2.0*n) * pow(bn, 2.0*n));
nfac = prof_gamma(4.0*n) / (prof_gamma(2.0*n)*pow(bn, 2.0*n));
hscale2 = 0.5 * *prof->scale**prof->scale;
m20 = hscale2 * nfac;
if (jac)
{
dmx2[profit->paramindex[PARAM_SPHEROID_REFF]]
= *prof->scale * nfac * mx2fac;
dmy2[profit->paramindex[PARAM_SPHEROID_REFF]]
= *prof->scale * nfac * my2fac;
dmxy[profit->paramindex[PARAM_SPHEROID_REFF]]
= *prof->scale * nfac * mxyfac;
n2 = n+0.01;
bn2 = 2.0*n2 - 1.0/3.0 + 4.0/(405.0*n2) + 46.0/(25515.0*n2*n2)
+ 131.0/(1148175*n2*n2*n2); /* Ciotti & Bertin 1999 */
nfac2 = prof_gamma(4.0*n2) / (prof_gamma(2.0*n2)*pow(bn2, 2.0*n2));
dmdn = 100.0 * hscale2 * (nfac2-nfac);
dmx2[profit->paramindex[PARAM_SPHEROID_SERSICN]] = dmdn * mx2fac;
dmy2[profit->paramindex[PARAM_SPHEROID_SERSICN]] = dmdn * my2fac;
dmxy[profit->paramindex[PARAM_SPHEROID_SERSICN]] = dmdn * mxyfac;
dmx2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
= 2.0 * m20 * st*st * *prof->aspect;
dmy2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
= 2.0 * m20 * ct*ct * *prof->aspect;
dmxy[profit->paramindex[PARAM_SPHEROID_ASPECT]]
= -2.0 * m20 * ct*st * *prof->aspect;
dmx2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (dc2+ds2*a2);
dmy2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (ds2+dc2*a2);
dmxy[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (1.0-a2)*dcs;
}
index = profit->paramindex[PARAM_SPHEROID_FLUX];
break;
case PROF_DEVAUCOULEURS:
m20 = 10.83995 * *prof->scale**prof->scale;
if (jac)
{
dmx2[profit->paramindex[PARAM_SPHEROID_REFF]]
= 21.680 * *prof->scale * mx2fac;
dmy2[profit->paramindex[PARAM_SPHEROID_REFF]]
= 21.680 * *prof->scale * my2fac;
dmxy[profit->paramindex[PARAM_SPHEROID_REFF]]
= 21.680 * *prof->scale * mxyfac;
dmx2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
= 2.0 * m20 * st*st * *prof->aspect;
dmy2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
= 2.0 * m20 * ct*ct * *prof->aspect;
dmxy[profit->paramindex[PARAM_SPHEROID_ASPECT]]
= -2.0 * m20 * ct*st * *prof->aspect;
dmx2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (dc2+ds2*a2);
dmy2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (ds2+dc2*a2);
dmxy[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (1.0-a2)*dcs;
}
index = profit->paramindex[PARAM_SPHEROID_FLUX];
break;
case PROF_EXPONENTIAL:
m20 = 3.0 * *prof->scale**prof->scale;
if (jac)
{
dmx2[profit->paramindex[PARAM_DISK_SCALE]]
= 6.0 * *prof->scale * mx2fac;
dmy2[profit->paramindex[PARAM_DISK_SCALE]]
= 6.0 * *prof->scale * my2fac;
dmxy[profit->paramindex[PARAM_DISK_SCALE]]
= 6.0 * *prof->scale * mxyfac;
dmx2[profit->paramindex[PARAM_DISK_ASPECT]]
= 2.0 * m20 * st*st * *prof->aspect;
dmy2[profit->paramindex[PARAM_DISK_ASPECT]]
= 2.0 * m20 * ct*ct * *prof->aspect;
dmxy[profit->paramindex[PARAM_DISK_ASPECT]]
= -2.0 * m20 * ct*st * *prof->aspect;
dmx2[profit->paramindex[PARAM_DISK_POSANG]] = m20 * (dc2 + ds2*a2);
dmy2[profit->paramindex[PARAM_DISK_POSANG]] = m20 * (ds2 + dc2*a2);
dmxy[profit->paramindex[PARAM_DISK_POSANG]] = m20 * (1.0 - a2) * dcs;
}
index = profit->paramindex[PARAM_DISK_FLUX];
break;
case PROF_ARMS:
m20 = 1.0;
index = profit->paramindex[PARAM_ARMS_FLUX];
break;
case PROF_BAR:
m20 = 1.0;
index = profit->paramindex[PARAM_BAR_FLUX];
break;
case PROF_INRING:
m20 = 1.0;
index = profit->paramindex[PARAM_INRING_FLUX];
break;
case PROF_OUTRING:
m20 = 1.0 / 1.0;
m20 = 1.0;
index = profit->paramindex[PARAM_OUTRING_FLUX];
break;
default:
m20 = 0.0; /* to avoid gcc -Wall warnings */
error(EXIT_FAILURE, "*Internal Error*: Unknown oriented model in ",
"prof_moments()");
break;
}
 
m02 = m20**prof->aspect**prof->aspect;
ct = cos(*prof->posangle*DEG);
st = sin(*prof->posangle*DEG);
prof->mx2 = ct*ct*m20 + st*st*m02;
prof->my2 = st*st*m20 + ct*ct*m02;
prof->mxy = ct*st*(m20 - m02);
prof->mx2 = m20*mx2fac;
prof->my2 = m20*my2fac;
prof->mxy = m20*mxyfac;
}
else
prof->mx2 = prof->my2 = prof->mxy = 0.0;
 
return;
return index;
}
 
 
/gpu/src/profit.h
9,7 → 9,7
*
* Contents: Include file for profit.c.
*
* Last modify: 08/07/2010
* Last modify: 21/07/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
178,13 → 178,14
int profit_copyobjpix(profitstruct *profit, picstruct *field,
picstruct *wfield),
profit_minimize(profitstruct *profit, int niter),
prof_moments(profitstruct *profit, profstruct *prof,
double *jac),
profit_resample(profitstruct *profit, float *inpix,
PIXTYPE *outpix, float factor),
profit_setparam(profitstruct *profit, paramenum paramtype,
float param, float parammin, float parammax);
 
void prof_end(profstruct *prof),
prof_moments(profitstruct *profit, profstruct *prof),
profit_addparam(profitstruct *profit, paramenum paramindex,
float **param),
profit_boundtounbound(profitstruct *profit,
/gpu/src/makeit.c
9,7 → 9,7
*
* Contents: main program.
*
* Last modify: 21/01/2010
* Last modify: 23/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
108,6 → 108,7
 
if (prefs.prof_flag)
{
#ifdef USE_MODEL
fft_init(prefs.nthreads);
/* Create profiles at full resolution */
NFPRINTF(OUTPUT, "Preparing profile models");
151,6 → 152,11
QPRINTF(OUTPUT, "%s", profname[theprofit->prof[i]->code]);
}
QPRINTF(OUTPUT, "\n");
#else
error(EXIT_FAILURE,
"*Error*: model-fitting is not supported in this build.\n",
" Please check your configure options");
#endif
}
 
if (prefs.filter_flag)
512,11 → 518,13
if (prefs.growth_flag)
endgrowth();
 
#ifdef USE_MODEL
if (prefs.prof_flag)
{
profit_end(theprofit);
fft_end();
}
#endif
 
if (prefs.psf_flag || prefs.prof_flag)
psf_end(thepsf,thepsfit); /*?*/
/gpu/src/image.c
9,7 → 9,7
*
* Contents: Function related to image manipulations.
*
* Last modify: 13/09/2009
* Last modify: 12/09/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
27,7 → 27,7
#include "prefs.h"
#include "image.h"
 
static float interpm[INTERPW*INTERPH];
static float interpm[INTERPW*INTERPW];
 
/********************************* copyimage *********************************/
/*
158,10 → 158,10
 
/* Compute the interpolation mask */
ddx0 = -(idmx=(INTERPW-(dx>0.0?1:0))/2)-dx;
ddy = -(idmy=(INTERPH-(dy>0.0?1:0))/2)-dy;
ddy = -(idmy=(INTERPW-(dy>0.0?1:0))/2)-dy;
sum = 0.0;
m = interpm;
for (my=INTERPH; my--; ddy+=1.0)
for (my=INTERPW; my--; ddy+=1.0)
{
ddx = ddx0;
fy = INTERPF(ddy);
171,7 → 171,7
 
/* Normalize it */
m = interpm;
for (i=INTERPW*INTERPH; i--;)
for (i=INTERPW*INTERPW; i--;)
*(m++) /= sum;
 
/* Do the interpolation */
180,7 → 180,7
ymin = iy - h/2 - idmy;
sw = field->width;
sh = field->stripheight;
for (my=INTERPH; my--; ymin++)
for (my=INTERPW; my--; ymin++)
{
/*-- Set the image boundaries in y */
if ((idy = field->ymin-ymin) > 0)
269,10 → 269,10
 
/* Compute the interpolation mask */
ddx0 = -(idmx=(INTERPW-(dx>0.0?1:0))/2)-dx;
ddy = -(idmy=(INTERPH-(dy>0.0?1:0))/2)-dy;
ddy = -(idmy=(INTERPW-(dy>0.0?1:0))/2)-dy;
sum = 0.0;
m = interpm;
for (my=INTERPH; my--; ddy+=1.0)
for (my=INTERPW; my--; ddy+=1.0)
{
ddx = ddx0;
fy = INTERPF(ddy);
282,7 → 282,7
 
/* Normalize it */
m = interpm;
for (i=INTERPW*INTERPH; i--;)
for (i=INTERPW*INTERPW; i--;)
*(m++) /= sum;
 
/* Do the interpolation */
291,7 → 291,7
ymin = iy - h/2 - idmy;
sw = field->width;
sh = field->stripheight;
for (my=INTERPH; my--; ymin++)
for (my=INTERPW; my--; ymin++)
{
/*-- Set the image boundaries in y */
if ((idy = field->ymin-ymin) > 0)
466,7 → 466,7
float dx, float dy, float step2)
{
float *mask,*maskt, xc1,xc2,yc1,yc2, xs1,ys1, x1,y1, x,y, dxm,dym,
val,
val, norm,
*pix12, *pixin,*pixin0, *pixout,*pixout0;
int i,j,k,n,t, *start,*startt, *nmask,*nmaskt,
ixs2,iys2, ix2,iy2, dix2,diy2, nx2,ny2, iys1a, ny1, hmw,hmh,
520,7 → 520,7
 
/* Set the yrange for the x-resampling with some margin for interpolation */
iys1a = (int)ys1; /* Int part of Im1 start y-coord with margin */
hmh = INTERPH/2 - 1; /* Interpolant start */
hmh = INTERPW/2 - 1; /* Interpolant start */
if (iys1a<0 || ((iys1a -= hmh)< 0))
iys1a = 0;
ny1 = (int)(ys1+ny2*step2)+INTERPW-hmh; /* Interpolated Im1 y size */
556,8 → 556,13
n=t;
*(startt++) = ix;
*(nmaskt++) = n;
norm = 0.0;
for (x=dxm, i=n; i--; x+=1.0)
*(maskt++) = INTERPF(x);
norm += (*(maskt++) = INTERPF(x));
norm = norm>0.0? 1.0/norm : 1.0;
maskt -= n;
for (i=n; i--;)
*(maskt++) *= norm;
}
 
QCALLOC(pix12, float, nx2*ny1); /* Intermediary frame-buffer */
582,12 → 587,12
}
 
/* Reallocate interpolant stuff for the y direction */
QREALLOC(mask, float, ny2*INTERPH); /* Interpolation masks */
QREALLOC(mask, float, ny2*INTERPW); /* Interpolation masks */
QREALLOC(nmask, int, ny2); /* Interpolation mask sizes */
QREALLOC(start, int, ny2); /* Int part of Im1 conv starts */
 
/* Compute the local interpolant and data starting points in y */
hmh = INTERPH/2 - 1;
hmh = INTERPW/2 - 1;
y1 = ys1;
maskt = mask;
nmaskt = nmask;
598,18 → 603,23
dym = iy1 - y1 - hmh; /* starting point in the interpolation func */
if (iy < 0)
{
n = INTERPH+iy;
n = INTERPW+iy;
dym -= (float)iy;
iy = 0;
}
else
n = INTERPH;
n = INTERPW;
if (n>(t=ny1-iy))
n=t;
*(startt++) = iy;
*(nmaskt++) = n;
norm = 0.0;
for (y=dym, i=n; i--; y+=1.0)
*(maskt++) = INTERPF(y);
norm += (*(maskt++) = INTERPF(y));
norm = norm>0.0? 1.0/norm : 1.0;
maskt -= n;
for (i=n; i--;)
*(maskt++) *= norm;
}
 
/* Make the interpolation in y and transpose once again */
/gpu/src/fits/fitshead.c
9,7 → 9,7
*
* Contents: general functions for handling FITS file headers.
*
* Last modify: 28/10/2009
* Last modify: 20/07/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
187,7 → 187,7
otherwise.
NOTES -.
AUTHOR E. Bertin (IAP & Leiden observatory)
VERSION 28/10/2009
VERSION 20/07/2010
***/
int readbintabparam_head(tabstruct *tab)
 
271,7 → 271,7
key->htype = H_STRING;
break;
default:
error(EXIT_FAILURE, "*Internal Error*: Unkwown T_TYPE for ", str);
error(EXIT_FAILURE, "*Error*: Unknown TFORM in ", cat->filename);
}
 
/*--handle the special case of multimensional arrays*/
/gpu/src/image.h
9,17 → 9,19
*
* Contents: Include file for image.c.
*
* Last modify: 13/09/2009
* Last modify: 12/09/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
 
/*----------------------------- Internal constants --------------------------*/
 
#define INTERPW 6 /* Interpolation function range (x) */
#define INTERPH 6 /* Interpolation function range (y) */
#define INTERPW 8 /* Interpolation function range */
#define INTERPFAC 4.0 /* Interpolation envelope factor */
 
#define INTERPF(x) (x==0.0?1.0:sin(PI*x)*sin(PI*x/3.0)/(PI*PI/3.0*x*x))
#define INTERPF(x) (x<1e-5 && x>-1e-5? 1.0 \
:(x>INTERPFAC?0.0:(x<-INTERPFAC?0.0 \
:sinf(PI*x)*sinf(PI/INTERPFAC*x)/(PI*PI/INTERPFAC*x*x))))
/* Lanczos approximation */
 
/*--------------------------- structure definitions -------------------------*/
/gpu/src/misc.c
9,7 → 9,7
*
* Contents: miscellaneous functions.
*
* Last modify: 24/09/2009
* Last modify: 20/08/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
73,6 → 73,60
}
 
 
/****** propagate_covar ******************************************************
PROTO void propagate_covar(double *vi, double *d, double *vo,
int ni, int no, double *temp)
PURPOSE Compute Dt.V.D (propagate covariance matrix errors)
INPUT Pointer to the original covariance matrix,
pointer to the matrix of derivatives,
input number of parameters,
output number of parameters,
pointer to a ni*no work array.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 20/08/2010
***/
void propagate_covar(double *vi, double *d, double *vo,
int ni, int no, double *temp)
{
double *vit,*dt,*vot,*tempt,
dval;
int i,j,k;
 
tempt = temp;
vit = vi;
for (j=0; j<ni; j++)
{
dt = d;
for (i=no; i--;)
{
vit = vi + j*ni;
dval = 0.0;
for (k=ni; k--;)
dval += *(vit++)**(dt++);
*(tempt++) = dval;
}
}
 
vot = vo;
for (j=0; j<no; j++)
{
for (i=0; i<no; i++)
{
dt = d + j*ni;
tempt = temp + i;
dval = 0.0;
for (k=ni; k--; tempt+=no)
dval += *(dt++)**tempt;
*(vot++) = dval;
}
}
 
return;
}
 
 
/****** counter_seconds *******************************************************
PROTO double counter_seconds(void)
PURPOSE Count the number of seconds (with an arbitrary offset).
/gpu/README
9,9 → 9,10
configuration files like default.param are still needed, though.
 
The SExtractor homepage is
http://terapix.iap.fr/soft/sextractor
http://astromatic.net/software/sextractor
In case of problems, questions or suggestions related to the software, please
refer to the TERAPIX forum:
http://terapix.iap.fr/forum/
refer to the SExtractor forum:
http://astromatic.net/forum/forumdisplay.php?fid=4
 
Emmanuel Bertin.
 
/gpu/config.h.in
197,6 → 197,9
/* Define to 1 if your <sys/time.h> declares `struct tm'. */
#undef TM_IN_SYS_TIME
 
/* Triggers model-fitting and linking with ATLAS/FFTW */
#undef USE_MODEL
 
/* Triggers multhreading */
#undef USE_THREADS
 
/gpu/configure.ac
6,7 → 6,7
define([AC_CACHE_SAVE],)
 
# This is your standard Bertin source code...
AC_INIT(sextractor, 2.11.7, [bertin@iap.fr])
AC_INIT(sextractor, 2.12.3, [bertin@iap.fr])
AC_CONFIG_SRCDIR(src/makeit.c)
AC_CONFIG_AUX_DIR(autoconf)
AM_CONFIG_HEADER(config.h)
142,6 → 142,20
 
AC_DEFINE_UNQUOTED([XSL_URL], "$xsl_url",[Default URL of the XSLT filter])
 
# Provide special option to disable model-fitting (enabled by default)
AC_MSG_CHECKING([if model-fitting should be disabled])
AC_ARG_ENABLE([model-fitting],
[AC_HELP_STRING([--disable-model-fitting],
[Disable model-fitting and ATLAS/FFTW dependency \
(enabled by default)])],
[use_model="$enableval"],
[use_model="yes"])
if test "$use_model" = "no"; then
AC_MSG_RESULT([yes])
else
AC_MSG_RESULT([no])
fi
 
# Set flags for multithreading
n_pthreads=16
AC_ARG_ENABLE(threads,
178,6 → 192,16
AC_MSG_RESULT([yes]),
AC_MSG_RESULT([no]))
 
############# Actions to complete if model-fitting is activated ##############
AC_MSG_CHECKING([for model-fitting configure option])
if test "$use_model" = "yes"; then
AC_MSG_RESULT([enabled])
AC_DEFINE(USE_MODEL, 1, [Triggers model-fitting and linking with ATLAS/FFTW])
else
AC_MSG_RESULT([disabled])
fi
AM_CONDITIONAL(USE_MODEL, [test $use_model = "yes"])
 
################# Actions to complete in case of multhreading ################
AC_DEFINE_UNQUOTED(THREADS_NMAX, $n_pthreads,[Maximum number of POSIX threads])
if test "$use_pthreads" = "yes"; then
190,30 → 214,34
[AM_CFLAGS="$AM_CFLAGS $PTHREAD_CFLAGS -D_REENTRANT"]
LIBS="$PTHREAD_LIBS $LIBS"
fi
AM_CONDITIONAL(USE_THREADS, test $use_pthreads = "yes")
AM_CONDITIONAL(USE_THREADS, [test $use_pthreads = "yes"])
 
################ handle the FFTW library (Fourier transforms) ################
ACX_FFTW($fftw_libdir,$fftw_incdir,$use_pthreads,yes,
if test "$use_model" = "yes"; then
ACX_FFTW($fftw_libdir,$fftw_incdir,$use_pthreads,yes,
[use_fftw=yes],[use_fftw=no])
if test "$use_fftw" = "yes"; then
LIBS="$FFTW_LIBS $LIBS"
if test "$FFTW_WARN" != ""; then
AC_MSG_WARN([$FFTW_WARN])
if test "$use_fftw" = "yes"; then
LIBS="$FFTW_LIBS $LIBS"
if test "$FFTW_WARN" != ""; then
AC_MSG_WARN([$FFTW_WARN])
fi
else
AC_MSG_ERROR([$FFTW_ERROR Exiting.])
fi
else
AC_MSG_ERROR([$FFTW_ERROR Exiting.])
fi
 
################## handle the ATLAS library(linear algebra) ##################
ACX_ATLAS($atlas_libdir,$atlas_incdir,$use_pthreads,
if test "$use_model" = "yes"; then
ACX_ATLAS($atlas_libdir,$atlas_incdir,$use_pthreads,
[use_atlas=yes],[use_atlas=no])
if test "$use_atlas" = "yes"; then
LIBS="$ATLAS_LIB $LIBS"
if test "$ATLAS_WARN" != ""; then
AC_MSG_WARN([$ATLAS_WARN])
if test "$use_atlas" = "yes"; then
LIBS="$ATLAS_LIB $LIBS"
if test "$ATLAS_WARN" != ""; then
AC_MSG_WARN([$ATLAS_WARN])
fi
else
AC_MSG_ERROR([$ATLAS_ERROR Exiting.])
fi
else
AC_MSG_ERROR([$ATLAS_ERROR Exiting.])
fi
 
# Link with gprof option
gpu Property changes : Modified: svn:mergeinfo Merged /trunk:r222-231