bindings-levmar 0.1.0.1 → 0.1.1
raw patch · 11 files changed
+513/−263 lines, 11 filesbuild-type:Customsetup-changed
Files
- Bindings/LevMar.hsc +245/−122
- Bindings/LevMar/CurryFriendly.hs +50/−13
- Setup.hs +4/−1
- bindings-levmar.cabal +25/−5
- configure +3/−0
- configure.hs +90/−0
- levmar-2.4/Axb_core.c +46/−98
- levmar-2.4/lm.h +15/−16
- levmar-2.4/lm_core.c +17/−0
- levmar-2.4/lmbc_core.c +10/−1
- levmar-2.4/misc_core.c +8/−7
Bindings/LevMar.hsc view
@@ -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 <lm.h> @@ -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 ---------------------------------------------------------------------
Bindings/LevMar/CurryFriendly.hs view
@@ -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))) --------------------------------------------------------------------------------
Setup.hs view
@@ -1,3 +1,6 @@+#! /usr/bin/env runhaskell+ import Distribution.Simple -main = defaultMain+main :: IO ()+main = defaultMainWithHooks autoconfUserHooks
bindings-levmar.cabal view
@@ -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
+ configure view
@@ -0,0 +1,3 @@+#! /bin/sh++runhaskell configure.hs $*
+ configure.hs view
@@ -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"
levmar-2.4/Axb_core.c view
@@ -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<m; i++) x[i]=qtb[i]; -#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif-+ FREE_LINSOLVER_MEM(buf); return 1; } @@ -292,7 +282,7 @@ if(m<n){ PRINT_ERROR(RCAT("Normal equations require that the number of rows is greater than number of columns in ", AX_EQ_B_QRLS) "() [%d x %d]! -- try transposing\n", m, n);- exit(1);+ return LM_ERROR; } /* calculate required memory size */@@ -318,7 +308,7 @@ buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL)); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_QRLS) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } } #else@@ -326,7 +316,7 @@ buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL)); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_QRLS) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } #endif /* LINSOLVERS_RETAIN_MEMORY */ @@ -352,16 +342,13 @@ GEQRF((int *)&m, (int *)&n, 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_QRLS) "()\n", -info);- exit(1);+ return LM_ERROR_LAPACK_ERROR; } else{ PRINT_ERROR(RCAT(RCAT("Unknown LAPACK error %d for ", GEQRF) " in ", AX_EQ_B_QRLS) "()\n", info);-#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif- return 0; } }@@ -380,16 +367,13 @@ TRTRS("U", "T", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, atb, (int *)&n, &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_QRLS) "()\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_QRLS) "()\n", info);-#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif- return 0; } }@@ -398,16 +382,13 @@ TRTRS("U", "N", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, atb, (int *)&n, &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_QRLS) "()\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_QRLS) "()\n", info);-#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif- return 0; } }@@ -416,10 +397,7 @@ for(i=0; i<n; i++) x[i]=atb[i]; -#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif-+ FREE_LINSOLVER_MEM(buf); return 1; } @@ -477,7 +455,7 @@ buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL)); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_CHOL) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } } #else@@ -485,7 +463,7 @@ buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL)); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_CHOL) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } #endif /* LINSOLVERS_RETAIN_MEMORY */ @@ -507,17 +485,14 @@ POTRF("U", (int *)&m, a, (int *)&m, (int *)&info); /* error treatment */ if(info!=0){+ FREE_LINSOLVER_MEM(buf); if(info<0){ PRINT_ERROR(RCAT(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", POTF2) "/", POTRF) " in ", AX_EQ_B_CHOL) "()\n", -info);- exit(1);+ return LM_ERROR; } else{ PRINT_ERROR(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- return 0; } }@@ -526,7 +501,7 @@ POTRS("U", (int *)&m, (int *)&nrhs, a, (int *)&m, b, (int *)&m, &info); if(info<0){ PRINT_ERROR(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", POTRS) " in ", AX_EQ_B_CHOL) "()\n", -info);- exit(1);+ return LM_ERROR_LAPACK_ERROR; } #if 0@@ -534,16 +509,13 @@ TRTRS("U", "T", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, b, (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_CHOL) "()\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_CHOL) "()\n", info);-#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif- return 0; } }@@ -552,16 +524,13 @@ TRTRS("U", "N", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, b, (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_CHOL) "()\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_CHOL) "()\n", info);-#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif- return 0; } }@@ -571,10 +540,7 @@ for(i=0; i<m; i++) x[i]=b[i]; -#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif-+ FREE_LINSOLVER_MEM(buf); return 1; } @@ -632,7 +598,7 @@ buf=(LM_REAL *)malloc(buf_sz); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } } #else@@ -640,7 +606,7 @@ buf=(LM_REAL *)malloc(buf_sz); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } #endif /* LINSOLVERS_RETAIN_MEMORY */ @@ -659,16 +625,13 @@ /* LU decomposition for A */ GETRF((int *)&m, (int *)&m, a, (int *)&m, ipiv, (int *)&info); if(info!=0){+ FREE_LINSOLVER_MEM(buf); if(info<0){ PRINT_ERROR(RCAT(RCAT("argument %d of ", GETRF) " illegal in ", AX_EQ_B_LU) "()\n", -info);- exit(1);+ return LM_ERROR; } else{ PRINT_ERROR(RCAT(RCAT("singular matrix A for ", GETRF) " in ", AX_EQ_B_LU) "()\n");-#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif- return 0; } }@@ -676,16 +639,13 @@ /* solve the system with the computed LU */ GETRS("N", (int *)&m, (int *)&nrhs, a, (int *)&m, ipiv, b, (int *)&m, (int *)&info); if(info!=0){+ FREE_LINSOLVER_MEM(buf); if(info<0){ PRINT_ERROR(RCAT(RCAT("argument %d of ", GETRS) " illegal in ", AX_EQ_B_LU) "()\n", -info);- exit(1);+ return LM_ERROR; } else{ PRINT_ERROR(RCAT(RCAT("unknown error for ", GETRS) " in ", AX_EQ_B_LU) "()\n");-#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif- return 0; } }@@ -695,10 +655,7 @@ x[i]=b[i]; } -#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif-+ FREE_LINSOLVER_MEM(buf); return 1; } @@ -770,7 +727,7 @@ buf=(LM_REAL *)malloc(buf_sz); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_SVD) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } } #else@@ -778,7 +735,7 @@ buf=(LM_REAL *)malloc(buf_sz); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_SVD) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } #endif /* LINSOLVERS_RETAIN_MEMORY */ @@ -800,16 +757,13 @@ /* error treatment */ if(info!=0){+ FREE_LINSOLVER_MEM(buf); if(info<0){ PRINT_ERROR(RCAT(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", GESVD), "/" GESDD) " in ", AX_EQ_B_SVD) "()\n", -info);- exit(1);+ return LM_ERROR_LAPACK_ERROR; } else{ PRINT_ERROR(RCAT("LAPACK error: dgesdd (dbdsdc)/dgesvd (dbdsqr) failed to converge in ", AX_EQ_B_SVD) "() [info=%d]\n", info);-#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif- return 0; } }@@ -840,10 +794,7 @@ x[i]=sum; } -#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif-+ FREE_LINSOLVER_MEM(buf); return 1; } @@ -921,7 +872,7 @@ buf=(void *)malloc(tot_sz); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } } #else@@ -929,7 +880,7 @@ buf=(void *)malloc(tot_sz); if(!buf){ PRINT_ERROR(RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } #endif /* LINSOLVERS_RETAIN_MEMORY */ @@ -1027,10 +978,7 @@ x[i]=sum/a[i*m+i]; } -#ifndef LINSOLVERS_RETAIN_MEMORY- free(buf);-#endif-+ FREE_LINSOLVER_MEM(buf); return 1; }
levmar-2.4/lm.h view
@@ -22,7 +22,6 @@ #ifndef _LM_H_ #define _LM_H_ - /************************************* Start of configuration options *************************************/ /* specify whether to use LAPACK or not. The first option is strongly recommended */@@ -49,7 +48,6 @@ /****************** End of configuration options, no changes necessary beyond this point ******************/ - #ifdef __cplusplus extern "C" { #endif@@ -61,16 +59,17 @@ #endif enum lmerror-{ LM_ERROR_LAPACK_ERROR = -1-, LM_ERROR_NO_JACOBIAN = -2-, LM_ERROR_NO_BOX_CONSTRAINTS = -3-, LM_ERROR_FAILED_BOX_CHECK = -4-, LM_ERROR_MEMORY_ALLOCATION_FAILURE = -5-, LM_ERROR_CONSTRAINT_MATRIX_ROWS_GT_COLS = -6-, LM_ERROR_CONSTRAINT_MATRIX_NOT_FULL_ROW_RANK = -7-, LM_ERROR_TOO_FEW_MEASUREMENTS = -8-, LM_ERROR_SINGULAR_MATRIX = -9-, LM_ERROR_SUM_OF_SQUARES_NOT_FINITE = -10+{ LM_ERROR = -1+, LM_ERROR_LAPACK_ERROR = -2+, LM_ERROR_NO_JACOBIAN = -3+, LM_ERROR_NO_BOX_CONSTRAINTS = -4+, LM_ERROR_FAILED_BOX_CHECK = -5+, LM_ERROR_MEMORY_ALLOCATION_FAILURE = -6+, LM_ERROR_CONSTRAINT_MATRIX_ROWS_GT_COLS = -7+, LM_ERROR_CONSTRAINT_MATRIX_NOT_FULL_ROW_RANK = -8+, LM_ERROR_TOO_FEW_MEASUREMENTS = -9+, LM_ERROR_SINGULAR_MATRIX = -10+, LM_ERROR_SUM_OF_SQUARES_NOT_FINITE = -11 }; #define FABS(x) (((x)>=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
levmar-2.4/lm_core.c view
@@ -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<m; ++i){@@ -411,6 +418,8 @@ LEVMAR_COVAR(jacTjac, covar, p_eL2, m, n); } +levmar_der_end: /* NOTE: this point is also reached via an explicit goto! */+ if(freework) free(work); #ifdef LINSOLVERS_RETAIN_MEMORY@@ -702,6 +711,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_dif_end;+ }+ if(issolved){ /* compute p's new estimate and ||Dp||^2 */ for(i=0, Dp_L2=0.0; i<m; ++i){@@ -818,6 +834,7 @@ LEVMAR_COVAR(jacTjac, covar, p_eL2, m, n); } +levmar_dif_end: /* NOTE: this point is also reached via an explicit goto! */ if(freework) free(work);
levmar-2.4/lmbc_core.c view
@@ -567,7 +567,14 @@ issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU; #endif /* HAVE_LAPACK */ - if(issolved){+ if (issolved < 0)+ {+ k = issolved;+ stop = 0;+ goto levmar_bcder_end;+ }++ if(issolved) { for(i=0; i<m; ++i) pDp[i]=p[i] + Dp[i]; @@ -803,6 +810,8 @@ if(covar){ LEVMAR_COVAR(jacTjac, covar, p_eL2, m, n); }++levmar_bcder_end: /* NOTE: this pointis also reached via an explicit goto! */ if(freework) free(work);
levmar-2.4/misc_core.c view
@@ -245,7 +245,7 @@ * other value which may cause loss of significance. */ -void LEVMAR_CHKJAC(+int LEVMAR_CHKJAC( 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, int m, int n, void *adata, LM_REAL *err)@@ -266,7 +266,7 @@ buf=(LM_REAL *)malloc((fvec_sz + fjac_sz + pp_sz + fvecp_sz)*sizeof(LM_REAL)); if(!buf){ PRINT_ERROR(LCAT(LEVMAR_CHKJAC, "(): memory allocation request failed\n"));- exit(1);+ return 0; } fvec=buf; fjac=fvec+fvec_sz;@@ -315,7 +315,7 @@ free(buf); - return;+ return 1; } #ifdef HAVE_LAPACK@@ -611,8 +611,8 @@ /* coefficient of determination. * see http://en.wikipedia.org/wiki/Coefficient_of_determination */-LM_REAL LEVMAR_R2(void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),- LM_REAL *p, LM_REAL *x, int m, int n, void *adata)+int LEVMAR_R2(void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),+ LM_REAL *p, LM_REAL *x, int m, int n, void *adata, LM_REAL *result) { register int i; register LM_REAL tmp;@@ -623,7 +623,7 @@ if((hx=(LM_REAL *)malloc(n*sizeof(LM_REAL)))==NULL){ PRINT_ERROR(RCAT("memory allocation request failed in ", LEVMAR_R2) "()\n");- exit(1);+ return LM_ERROR_MEMORY_ALLOCATION_FAILURE; } /* hx=f(p) */@@ -643,7 +643,8 @@ free(hx); - return LM_CNST(1.0) - SSerr/SStot;+ *result = LM_CNST(1.0) - SSerr/SStot;+ return 0; } /* check box constraints for consistency */