# albertoruiz/hmatrix

merge

1 parent 47ea0f2 commit 6399320e0043275baae2b20c87798c27efd6bf33 committed Feb 27, 2012
Showing with 152 additions and 3 deletions.
1. +4 −1 CHANGES.md
2. +2 −0 THANKS.md
3. +103 −1 lib/Numeric/GSL/Integration.hs
4. +42 −0 lib/Numeric/GSL/gsl-aux.c
5. +1 −1 lib/Numeric/LinearAlgebra/Algorithms.hs
5 CHANGES.md
 @@ -1,7 +1,10 @@ 0.13.2.0 -------- -msadams and msbdf methods for ode +- integration over infinite intervals + +- msadams and msbdf methods for ode + 0.13.0.0 --------
2 THANKS.md
 @@ -107,3 +107,5 @@ module reorganization, monadic mapVectorM, and many other improvements. - Daniel Fischer reported some Haddock markup errors. +- Danny Chan added support for integration over infinite intervals. +
104 lib/Numeric/GSL/Integration.hs
 @@ -17,7 +17,10 @@ Numerical integration routines. module Numeric.GSL.Integration ( integrateQNG, - integrateQAGS + integrateQAGS, + integrateQAGI, + integrateQAGIU, + integrateQAGIL ) where import Foreign.C.Types @@ -94,3 +97,102 @@ integrateQNG prec f a b = unsafePerformIO \$ do foreign import ccall "gsl-aux.h integrate_qng" c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS). +For example: + +@\> let quad = integrateQAGI 1E-9 1000 +\> let f a x = exp(-a * x^2) +\> quad (f 0.5) +(2.5066282746310002,6.229215880648858e-11)@ + +-} + +integrateQAGI :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf) + -> (Double, Double) -- ^ result of the integration and error +integrateQAGI prec n f = unsafePerformIO \$ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagi fp prec (fromIntegral n) r e // check "integrate_qagi" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall "gsl-aux.h integrate_qagi" + c_integrate_qagi :: FunPtr (Double-> Ptr() -> Double) -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf). +For example: + +@\> let quad = integrateQAGIU 1E-9 1000 +\> let f a x = exp(-a * x^2) +\> quad (f 0.5) 0 +(1.2533141373155001,3.114607940324429e-11)@ + +-} + +integrateQAGIU :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) + -> Double -- ^ a + -> (Double, Double) -- ^ result of the integration and error +integrateQAGIU prec n f a = unsafePerformIO \$ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagiu fp a prec (fromIntegral n) r e // check "integrate_qagiu" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall "gsl-aux.h integrate_qagiu" + c_integrate_qagiu :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b). +For example: + +@\> let quad = integrateQAGIL 1E-9 1000 +\> let f a x = exp(-a * x^2) +\> quad (f 0.5) 0 +(1.2533141373155001,3.114607940324429e-11)@ + +-} + +integrateQAGIL :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) + -> Double -- ^ b + -> (Double, Double) -- ^ result of the integration and error +integrateQAGIL prec n f b = unsafePerformIO \$ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagil fp b prec (fromIntegral n) r e // check "integrate_qagil" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall "gsl-aux.h integrate_qagil" + c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt +
42 lib/Numeric/GSL/gsl-aux.c
 @@ -761,6 +761,48 @@ int integrate_qags(double f(double,void*), double a, double b, double prec, int OK } +int integrate_qagi(double f(double,void*), double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagi"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagi (&F, 0, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + +int integrate_qagiu(double f(double,void*), double a, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagiu"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagiu (&F, a, 0, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + +int integrate_qagil(double f(double,void*), double b, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagil"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagil (&F, b, 0, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + int polySolve(KRVEC(a), CVEC(z)) { DEBUGMSG("polySolve"); REQUIRES(an>1,BAD_SIZE);
2 lib/Numeric/LinearAlgebra/Algorithms.hs
 @@ -209,7 +209,7 @@ rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v) | otherwise = let (_,s,v) = svd m in (s,v) --- | Singular values and all right singular vectors. +-- | Singular values and all left singular vectors. leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) leftSV m | vertical m = let (u,s,_) = svd m in (u,s) | otherwise = let (u,s,_) = thinSVD m in (u,s)