diff --git a/Bindings/LevMar.hsc b/Bindings/LevMar.hsc --- a/Bindings/LevMar.hsc +++ b/Bindings/LevMar.hsc @@ -23,6 +23,7 @@ , _LM_INFO_SZ -- * Errors. + , _LM_ERROR , _LM_ERROR_LAPACK_ERROR , _LM_ERROR_NO_JACOBIAN , _LM_ERROR_NO_BOX_CONSTRAINTS @@ -39,6 +40,24 @@ , _LM_STOP_THRESH , _LM_DIFF_DELTA + -- * Handy type synonyms + , Parameters + , Measurements + , Options + , LowerBounds + , UpperBounds + , ConstraintsMatrix + , ConstraintsVector + , Weights + , Info + , Work + , Covar + , AData + , NrOfParameters + , NrOfMeasurements + , NrOfConstraints + , MaxIterations + -- * Model & Jacobian. , Model , Jacobian @@ -73,12 +92,34 @@ , slevmar_blec_der , dlevmar_blec_dif , slevmar_blec_dif + + -- * Jacobian verification + , Errors + , LevMarChkJac + , dlevmar_chkjac + , slevmar_chkjac + + -- * Utils + , BestFitParameterIx + + , LevMarStddev + , LevMarCorCoef + , LevMarR2 + + , Result + + , dlevmar_stddev + , slevmar_stddev + , dlevmar_corcoef + , slevmar_corcoef + , dlevmar_R2 + , slevmar_R2 ) where -import Foreign.C.Types (CInt, CFloat, CDouble) -import Foreign.Ptr (Ptr, FunPtr, freeHaskellFunPtr) -import Control.Exception (bracket) +import Foreign.C.Types ( CInt, CFloat, CDouble ) +import Foreign.Ptr ( Ptr, FunPtr, freeHaskellFunPtr ) +import Control.Exception ( bracket ) #include @@ -106,6 +147,7 @@ -------------------------------------------------------------------------------- #{enum CInt, + , _LM_ERROR = LM_ERROR , _LM_ERROR_LAPACK_ERROR = LM_ERROR_LAPACK_ERROR , _LM_ERROR_NO_JACOBIAN = LM_ERROR_NO_JACOBIAN , _LM_ERROR_NO_BOX_CONSTRAINTS = LM_ERROR_NO_BOX_CONSTRAINTS @@ -133,15 +175,37 @@ -------------------------------------------------------------------------------- +-- Handy type synonyms +-------------------------------------------------------------------------------- + +type Parameters = Ptr +type Measurements = Ptr +type Options = Ptr +type LowerBounds = Ptr +type UpperBounds = Ptr +type ConstraintsMatrix = Ptr +type ConstraintsVector = Ptr +type Weights = Ptr +type Info = Ptr +type Work = Ptr +type Covar = Ptr +type AData = Ptr () +type NrOfParameters = CInt +type NrOfMeasurements = CInt +type NrOfConstraints = CInt +type MaxIterations = CInt + + +-------------------------------------------------------------------------------- -- Model & Jacobian. -------------------------------------------------------------------------------- -- | Functional relation describing measurements. -type Model r = Ptr r -- p - -> Ptr r -- hx - -> CInt -- m - -> CInt -- n - -> Ptr () -- adata +type Model r = Parameters r + -> Measurements r + -> NrOfParameters + -> NrOfMeasurements + -> AData -> IO () type Jacobian a = Model a @@ -162,136 +226,137 @@ -- Types of the Levenberg-Marquardt algorithms. -------------------------------------------------------------------------------- -type LevMarDer cr = FunPtr (Model cr) -- func - -> FunPtr (Jacobian cr) -- jacf - -> Ptr cr -- p - -> Ptr cr -- x - -> CInt -- m - -> CInt -- n - -> CInt -- itmax - -> Ptr cr -- opts - -> Ptr cr -- info - -> Ptr cr -- work - -> Ptr cr -- covar - -> Ptr () -- adata +type LevMarDer cr = FunPtr (Model cr) + -> FunPtr (Jacobian cr) + -> Parameters cr + -> Measurements cr + -> NrOfParameters + -> NrOfMeasurements + -> MaxIterations + -> Options cr + -> Info cr + -> Work cr + -> Covar cr + -> AData -> IO CInt -type LevMarDif cr = FunPtr (Model cr) -- func - -> Ptr cr -- p - -> Ptr cr -- x - -> CInt -- m - -> CInt -- n - -> CInt -- itmax - -> Ptr cr -- opts - -> Ptr cr -- info - -> Ptr cr -- work - -> Ptr cr -- covar - -> Ptr () -- adata +type LevMarDif cr = FunPtr (Model cr) + -> Parameters cr + -> Measurements cr + -> NrOfParameters + -> NrOfMeasurements + -> MaxIterations + -> Options cr + -> Info cr + -> Work cr + -> Covar cr + -> AData -> IO CInt -type LevMarBCDer cr = FunPtr (Model cr) -- func - -> FunPtr (Jacobian cr) -- jacf - -> Ptr cr -- p - -> Ptr cr -- x - -> CInt -- m - -> CInt -- n - -> Ptr cr -- lb - -> Ptr cr -- ub - -> CInt -- itmax - -> Ptr cr -- opts - -> Ptr cr -- info - -> Ptr cr -- work - -> Ptr cr -- covar - -> Ptr () -- adata +type LevMarBCDer cr = FunPtr (Model cr) + -> FunPtr (Jacobian cr) + -> Parameters cr + -> Measurements cr + -> NrOfParameters + -> NrOfMeasurements + -> LowerBounds cr + -> UpperBounds cr + -> MaxIterations + -> Options cr + -> Info cr + -> Work cr + -> Covar cr + -> AData -> IO CInt -type LevMarBCDif cr = FunPtr (Model cr) -- func - -> Ptr cr -- p - -> Ptr cr -- x - -> CInt -- m - -> CInt -- n - -> Ptr cr -- lb - -> Ptr cr -- ub - -> CInt -- itmax - -> Ptr cr -- opts - -> Ptr cr -- info - -> Ptr cr -- work - -> Ptr cr -- covar - -> Ptr () -- adata +type LevMarBCDif cr = FunPtr (Model cr) + -> Parameters cr + -> Measurements cr + -> NrOfParameters + -> NrOfMeasurements + -> LowerBounds cr + -> UpperBounds cr + -> MaxIterations + -> Options cr + -> Info cr + -> Work cr + -> Covar cr + -> AData -> IO CInt -type LevMarLecDer cr = FunPtr (Model cr) -- func - -> FunPtr (Jacobian cr) -- jacf - -> Ptr cr -- p - -> Ptr cr -- x - -> CInt -- m - -> CInt -- n - -> Ptr cr -- A - -> Ptr cr -- B - -> CInt -- k - -> CInt -- itmax - -> Ptr cr -- opts - -> Ptr cr -- info - -> Ptr cr -- work - -> Ptr cr -- covar - -> Ptr () -- adata +type LevMarLecDer cr = FunPtr (Model cr) + -> FunPtr (Jacobian cr) + -> Parameters cr + -> Measurements cr + -> NrOfParameters + -> NrOfMeasurements + -> ConstraintsMatrix cr + -> ConstraintsVector cr + -> NrOfConstraints + -> MaxIterations + -> Options cr + -> Info cr + -> Work cr + -> Covar cr + -> AData -> IO CInt -type LevMarLecDif cr = FunPtr (Model cr) -- func - -> Ptr cr -- p - -> Ptr cr -- x - -> CInt -- m - -> CInt -- n - -> Ptr cr -- A - -> Ptr cr -- B - -> CInt -- k - -> CInt -- itmax - -> Ptr cr -- opts - -> Ptr cr -- info - -> Ptr cr -- work - -> Ptr cr -- covar - -> Ptr () -- adata +type LevMarLecDif cr = FunPtr (Model cr) + -> Parameters cr + -> Measurements cr + -> NrOfParameters + -> NrOfMeasurements + -> ConstraintsMatrix cr + -> ConstraintsVector cr + -> NrOfConstraints + -> MaxIterations + -> Options cr + -> Info cr + -> Work cr + -> Covar cr + -> AData -> IO CInt -type LevMarBLecDer cr = FunPtr (Model cr) -- func - -> FunPtr (Jacobian cr) -- jacf - -> Ptr cr -- p - -> Ptr cr -- x - -> CInt -- m - -> CInt -- n - -> Ptr cr -- lb - -> Ptr cr -- ub - -> Ptr cr -- A - -> Ptr cr -- B - -> CInt -- k - -> Ptr cr -- wghts - -> CInt -- itmax - -> Ptr cr -- opts - -> Ptr cr -- info - -> Ptr cr -- work - -> Ptr cr -- covar - -> Ptr () -- adata +type LevMarBLecDer cr = FunPtr (Model cr) + -> FunPtr (Jacobian cr) + -> Parameters cr + -> Measurements cr + -> NrOfParameters + -> NrOfMeasurements + -> LowerBounds cr + -> UpperBounds cr + -> ConstraintsMatrix cr + -> ConstraintsVector cr + -> NrOfConstraints + -> Weights cr + -> MaxIterations + -> Options cr + -> Info cr + -> Work cr + -> Covar cr + -> AData -> IO CInt -type LevMarBLecDif cr = FunPtr (Model cr) -- func - -> Ptr cr -- p - -> Ptr cr -- x - -> CInt -- m - -> CInt -- n - -> Ptr cr -- lb - -> Ptr cr -- ub - -> Ptr cr -- A - -> Ptr cr -- B - -> CInt -- k - -> Ptr cr -- wghts - -> CInt -- itmax - -> Ptr cr -- opts - -> Ptr cr -- info - -> Ptr cr -- work - -> Ptr cr -- covar - -> Ptr () -- adata +type LevMarBLecDif cr = FunPtr (Model cr) + -> Parameters cr + -> Measurements cr + -> NrOfParameters + -> NrOfMeasurements + -> LowerBounds cr + -> UpperBounds cr + -> ConstraintsMatrix cr + -> ConstraintsVector cr + -> NrOfConstraints + -> Weights cr + -> MaxIterations + -> Options cr + -> Info cr + -> Work cr + -> Covar cr + -> AData -> IO CInt + -------------------------------------------------------------------------------- -- Levenberg-Marquardt algorithms. -------------------------------------------------------------------------------- @@ -312,6 +377,64 @@ foreign import ccall "dlevmar_blec_der" dlevmar_blec_der :: LevMarBLecDer CDouble foreign import ccall "slevmar_blec_dif" slevmar_blec_dif :: LevMarBLecDif CFloat foreign import ccall "dlevmar_blec_dif" dlevmar_blec_dif :: LevMarBLecDif CDouble + + +-------------------------------------------------------------------------------- +-- Jacobian verification +-------------------------------------------------------------------------------- + +type Errors = Ptr + +type LevMarChkJac cr = FunPtr (Model cr) + -> FunPtr (Jacobian cr) + -> Parameters cr + -> NrOfParameters + -> NrOfMeasurements + -> AData + -> Errors cr + -> IO () + +foreign import ccall "dlevmar_chkjac" dlevmar_chkjac :: LevMarChkJac CDouble +foreign import ccall "slevmar_chkjac" slevmar_chkjac :: LevMarChkJac CFloat + + +-------------------------------------------------------------------------------- +-- Utils +-------------------------------------------------------------------------------- + +type BestFitParameterIx = CInt + +-- | Standard deviation. +type LevMarStddev cr = Covar cr + -> NrOfParameters + -> BestFitParameterIx + -> IO cr + +-- | Pearson's correlation coefficient for best-fit parameters. +type LevMarCorCoef cr = Covar cr + -> NrOfParameters + -> BestFitParameterIx + -> BestFitParameterIx + -> IO cr + +-- | Coefficient of determination (R2). +type LevMarR2 cr = FunPtr (Model cr) + -> Parameters cr + -> Measurements cr + -> NrOfParameters + -> NrOfMeasurements + -> AData + -> Result cr + -> IO CInt + +type Result = Ptr + +foreign import ccall "dlevmar_stddev" dlevmar_stddev :: LevMarStddev CDouble +foreign import ccall "slevmar_stddev" slevmar_stddev :: LevMarStddev CFloat +foreign import ccall "dlevmar_corcoef" dlevmar_corcoef :: LevMarCorCoef CDouble +foreign import ccall "slevmar_corcoef" slevmar_corcoef :: LevMarCorCoef CFloat +foreign import ccall "dlevmar_R2" dlevmar_R2 :: LevMarR2 CDouble +foreign import ccall "slevmar_R2" slevmar_R2 :: LevMarR2 CFloat -- The End --------------------------------------------------------------------- diff --git a/Bindings/LevMar/CurryFriendly.hs b/Bindings/LevMar/CurryFriendly.hs --- a/Bindings/LevMar/CurryFriendly.hs +++ b/Bindings/LevMar/CurryFriendly.hs @@ -23,6 +23,7 @@ , LMA_C._LM_INFO_SZ -- * Errors + , LMA_C._LM_ERROR , LMA_C._LM_ERROR_LAPACK_ERROR , LMA_C._LM_ERROR_NO_JACOBIAN , LMA_C._LM_ERROR_NO_BOX_CONSTRAINTS @@ -39,6 +40,24 @@ , LMA_C._LM_STOP_THRESH , LMA_C._LM_DIFF_DELTA + -- * Handy type synonyms + , LMA_C.Parameters + , LMA_C.Measurements + , LMA_C.Options + , LMA_C.LowerBounds + , LMA_C.UpperBounds + , LMA_C.ConstraintsMatrix + , LMA_C.ConstraintsVector + , LMA_C.Weights + , LMA_C.Info + , LMA_C.Work + , LMA_C.Covar + , LMA_C.AData + , LMA_C.NrOfParameters + , LMA_C.NrOfMeasurements + , LMA_C.NrOfConstraints + , LMA_C.MaxIterations + -- * Model & Jacobian , LMA_C.Model , LMA_C.Jacobian @@ -49,7 +68,6 @@ -- * Handy type synonyms used in the curry friendly types. , BoxConstraints , LinearConstraints - , Weights -- * Curry friendly types of the Levenberg-Marquardt algorithms. , LevMarDer @@ -78,11 +96,33 @@ , slevmar_blec_der , dlevmar_blec_dif , slevmar_blec_dif + + -- * Jacobian verification + , LMA_C.Errors + , LMA_C.LevMarChkJac + , LMA_C.dlevmar_chkjac + , LMA_C.slevmar_chkjac + + -- * Utils + , LMA_C.BestFitParameterIx + + , LMA_C.LevMarStddev + , LMA_C.LevMarCorCoef + , LMA_C.LevMarR2 + + , LMA_C.Result + + , LMA_C.dlevmar_stddev + , LMA_C.slevmar_stddev + , LMA_C.dlevmar_corcoef + , LMA_C.slevmar_corcoef + , LMA_C.dlevmar_R2 + , LMA_C.slevmar_R2 ) where -import Foreign.C.Types (CInt, CFloat, CDouble) -import Foreign.Ptr (Ptr, FunPtr) +import Foreign.C.Types ( CFloat, CDouble ) +import Foreign.Ptr ( FunPtr ) import qualified Bindings.LevMar as LMA_C @@ -91,16 +131,13 @@ -- Handy type synonyms used in the curry friendly types. -------------------------------------------------------------------------------- -type BoxConstraints cr a = Ptr cr -- Lower bounds - -> Ptr cr -- Upper bounds - -> a - -type LinearConstraints cr a = Ptr cr -- Constraints matrix - -> Ptr cr -- Right hand constraints vector - -> CInt -- Number of constraints +type BoxConstraints cr a = LMA_C.LowerBounds cr + -> LMA_C.UpperBounds cr -> a -type Weights cr a = Ptr cr -- Weights +type LinearConstraints cr a = LMA_C.ConstraintsMatrix cr + -> LMA_C.ConstraintsVector cr + -> LMA_C.NrOfConstraints -> a @@ -114,8 +151,8 @@ type LevMarBCDer cr = BoxConstraints cr (LevMarDer cr) type LevMarLecDif cr = LinearConstraints cr (LevMarDif cr) type LevMarLecDer cr = LinearConstraints cr (LevMarDer cr) -type LevMarBLecDif cr = BoxConstraints cr (LinearConstraints cr (Weights cr (LevMarDif cr))) -type LevMarBLecDer cr = BoxConstraints cr (LinearConstraints cr (Weights cr (LevMarDer cr))) +type LevMarBLecDif cr = BoxConstraints cr (LinearConstraints cr (LMA_C.Weights cr -> (LevMarDif cr))) +type LevMarBLecDer cr = BoxConstraints cr (LinearConstraints cr (LMA_C.Weights cr -> (LevMarDer cr))) -------------------------------------------------------------------------------- diff --git a/Setup.hs b/Setup.hs --- a/Setup.hs +++ b/Setup.hs @@ -1,3 +1,6 @@ +#! /usr/bin/env runhaskell + import Distribution.Simple -main = defaultMain +main :: IO () +main = defaultMainWithHooks autoconfUserHooks diff --git a/bindings-levmar.cabal b/bindings-levmar.cabal --- a/bindings-levmar.cabal +++ b/bindings-levmar.cabal @@ -1,14 +1,15 @@ name: bindings-levmar -version: 0.1.0.1 +version: 0.1.1 cabal-version: >= 1.6 -build-type: Simple +build-type: Custom stability: experimental +tested-with: GHC ==6.10.4 author: Roel van Dijk & Bas van Dijk maintainer: vandijk.roel@gmail.com, v.dijk.bas@gmail.com copyright: (c) 2009 Roel van Dijk & Bas van Dijk license: OtherLicense license-file: LICENSE -category: numerical +category: numerical, FFI synopsis: A binding to the C levmar (Levenberg-Marquardt) library description: The Levenberg-Marquardt algorithm is an iterative technique that finds a local minimum of a function that @@ -43,8 +44,11 @@ program can only by distributed under the terms of the GPL. +extra-source-files: configure configure.hs +extra-tmp-files: bindings-levmar.buildinfo -extra-source-files: levmar-2.4/LICENSE +extra-source-files: LICENSE + , levmar-2.4/LICENSE , levmar-2.4/*.h , levmar-2.4/*.c , levmar-2.4/*.txt @@ -58,6 +62,14 @@ , levmar-2.4/matlab/Makefile , levmar-2.4/matlab/Makefile.w32 +flag mkl + description: Link with Intel's MKL optimized libraries. + default: False + +flag accelerate + description: Use the accelerate framework for LAPACK/BLAS on OS X + default: False + source-repository head type: darcs location: http://code.haskell.org/bindings-levmar @@ -68,7 +80,6 @@ , Bindings.LevMar.CurryFriendly ghc-options: -Wall -O2 cc-options: -D_OPENMP - extra-libraries: lapack include-dirs: levmar-2.4 c-sources: levmar-2.4/Axb.c @@ -77,3 +88,12 @@ levmar-2.4/lmblec.c levmar-2.4/lmlec.c levmar-2.4/misc.c + + if flag(mkl) + if arch(x86_64) + extra-libraries: mkl_lapack mkl_intel_lp64 mkl_sequential mkl_core + else + extra-libraries: mkl_lapack mkl_intel mkl_sequential mkl_core + + if flag(accelerate) + frameworks: Accelerate diff --git a/configure b/configure new file mode 100644 --- /dev/null +++ b/configure @@ -0,0 +1,3 @@ +#! /bin/sh + +runhaskell configure.hs $* diff --git a/configure.hs b/configure.hs new file mode 100644 --- /dev/null +++ b/configure.hs @@ -0,0 +1,90 @@ +#! /usr/bin/env runhaskell + +{- +configure.hs for bindings-levmar + +This file was completely copied from hmatrix-0.5.2.2 and then adjusted +for bindings-levmar. +-} + +import System +import Data.List (isPrefixOf) +import Distribution.Simple.LocalBuildInfo +import Distribution.Simple.Configure +import Distribution.PackageDescription + +-- possible additional dependencies for the desired libs (by default lapack) +opts :: [String] +opts = [ "" -- Ubuntu/Debian + , "blas" + , "blas cblas" + , "cblas" + , "f77blas cblas atlas gcc_s" -- Arch Linux (older version of atlas-lapack) + , "blas gfortran" -- Arch Linux with normal blas and lapack + ] + +-- Compile a simple program with symbols from LAPACK with the given libs. +testProg :: String -> String -> String +testProg libs fmks = "gcc levmar_lapack_test.c -o /tmp/dummy " + ++ f1 libs ++ " " ++ f2 fmks + ++ " > /dev/null 2> /dev/null" + +f1, f2 :: String -> String +f1 = unwords . map ("-l" ++) . words +f2 = unwords . map ("-framework "++) . words + +check :: String -> String -> IO Bool +check libs fmks = checkCommand (testProg libs fmks) + +checkCommand :: String -> IO Bool +checkCommand c = (ExitSuccess ==) `fmap` system c + +-- test different configurations until the first one works +try :: String -> String -> [String] -> IO (Maybe String) +try _ _ [] = return Nothing +try b f (opt:rest) = do + ok <- check (b ++ " " ++ opt) f + if ok then return (Just opt) + else try b f rest + +-- read --configure-option=link:lib1,lib2,lib3,etc +linkOp :: String +linkOp = "link:" + +getUserLink :: [String] -> String +getUserLink = concatMap (g . drop (length linkOp)) . filter (isPrefixOf linkOp) + where g = map cs + cs ',' = ' ' + cs x = x + +main :: IO () +main = do + putStr "Checking foreign libraries..." + + args <- getArgs + Just bInfo <- maybeGetPersistBuildConfig "dist" + + let Just lib = library . localPkgDescr $ bInfo + base = unwords . extraLibs . libBuildInfo $ lib + fwks = unwords . frameworks . libBuildInfo $ lib + auxpref = getUserLink args + + -- We extract the desired libs from bindings-levmar.cabal (using cabal flags) + -- and from a possible --configure-option=link:lib1,lib2,lib3 + -- by default the desired lib is lapack. + + let pref = if null (words (base ++ " " ++ auxpref)) then "lapack" else auxpref + fullOpts = map ((pref ++ " ") ++ ) opts + + r <- try base fwks fullOpts + case r of + Nothing -> do + putStrLn " FAIL" + putStrLn " *** Sorry, I can't link LAPACK." + putStrLn " *** Please make sure that the appropriate -dev packages are installed." + putStrLn " *** You can also specify the required libraries using" + putStrLn " *** cabal install bindings-levmar --configure-option=link:lib1,lib2,lib3,etc." + writeFile "bindings-levmar.buildinfo" "buildable: False\n" + Just ops -> do + putStrLn " OK" + writeFile "bindings-levmar.buildinfo" $ "extra-libraries: " ++ ops ++ "\n" diff --git a/levmar-2.4/Axb_core.c b/levmar-2.4/Axb_core.c --- a/levmar-2.4/Axb_core.c +++ b/levmar-2.4/Axb_core.c @@ -29,8 +29,10 @@ #ifdef LINSOLVERS_RETAIN_MEMORY #define __STATIC__ static +#define FREE_LINSOLVER_MEM(B) // empty #else #define __STATIC__ // empty +#define FREE_LINSOLVER_MEM(B) free(B) #endif /* LINSOLVERS_RETAIN_MEMORY */ #ifdef HAVE_LAPACK @@ -147,7 +149,7 @@ buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL)); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_QR) "() failed!\n"); - exit(1); + return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } } #else @@ -155,7 +157,7 @@ buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL)); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_QR) "() failed!\n"); - exit(1); + return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } #endif /* LINSOLVERS_RETAIN_MEMORY */ @@ -174,16 +176,13 @@ GEQRF((int *)&m, (int *)&m, a, (int *)&m, tau, work, (int *)&worksz, (int *)&info); /* error treatment */ if(info!=0){ + FREE_LINSOLVER_MEM(buf); if(info<0){ PRINT_ERROR(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", GEQRF) " in ", AX_EQ_B_QR) "()\n", -info); - exit(1); + return LM_ERROR_LAPACK_ERROR; } else{ PRINT_ERROR(RCAT(RCAT("Unknown LAPACK error %d for ", GEQRF) " in ", AX_EQ_B_QR) "()\n", info); -#ifndef LINSOLVERS_RETAIN_MEMORY - free(buf); -#endif - return 0; } } @@ -195,16 +194,13 @@ /* compute Q using the elementary reflectors computed by the above decomposition */ ORGQR((int *)&m, (int *)&m, (int *)&m, a, (int *)&m, tau, work, (int *)&worksz, (int *)&info); if(info!=0){ + FREE_LINSOLVER_MEM(buf); if(info<0){ PRINT_ERROR(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", ORGQR) " in ", AX_EQ_B_QR) "()\n", -info); - exit(1); + return LM_ERROR_LAPACK_ERROR; } else{ PRINT_ERROR(RCAT("Unknown LAPACK error (%d) in ", AX_EQ_B_QR) "()\n", info); -#ifndef LINSOLVERS_RETAIN_MEMORY - free(buf); -#endif - return 0; } } @@ -220,16 +216,13 @@ TRTRS("U", "N", "N", (int *)&m, (int *)&nrhs, r, (int *)&m, qtb, (int *)&m, &info); /* error treatment */ if(info!=0){ + FREE_LINSOLVER_MEM(buf); if(info<0){ PRINT_ERROR(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) " in ", AX_EQ_B_QR) "()\n", -info); - exit(1); + return LM_ERROR_LAPACK_ERROR; } else{ PRINT_ERROR(RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_QR) "()\n", info); -#ifndef LINSOLVERS_RETAIN_MEMORY - free(buf); -#endif - return 0; } } @@ -238,10 +231,7 @@ for(i=0; i=0.0)? (x) : -(x)) @@ -248,14 +247,14 @@ /* Jacobian verification, double & single precision */ #ifdef LM_DBL_PREC -extern void dlevmar_chkjac( +extern int 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( +extern int 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); @@ -265,14 +264,14 @@ #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); +extern int dlevmar_R2(void (*func)(double *p, double *hx, int m, int n, void *adata), double *p, double *x, int m, int n, void *adata, double *result); #endif /* LM_DBL_PREC */ #ifdef LM_SNGL_PREC extern float slevmar_stddev( float *covar, int m, int i); extern float slevmar_corcoef(float *covar, int m, int i, int j); -extern float slevmar_R2(void (*func)(float *p, float *hx, int m, int n, void *adata), float *p, float *x, int m, int n, void *adata); +extern int slevmar_R2(void (*func)(float *p, float *hx, int m, int n, void *adata), float *p, float *x, int m, int n, void *adata, float *result); #endif /* LM_SNGL_PREC */ #ifdef __cplusplus diff --git a/levmar-2.4/lm_core.c b/levmar-2.4/lm_core.c --- a/levmar-2.4/lm_core.c +++ b/levmar-2.4/lm_core.c @@ -309,6 +309,13 @@ issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU; #endif /* HAVE_LAPACK */ + if (issolved < 0) + { + k = issolved; + stop = 0; + goto levmar_der_end; + } + if(issolved){ /* compute p's new estimate and ||Dp||^2 */ for(i=0, Dp_L2=0.0; i