diff --git a/CHANGES.md b/CHANGES.md --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,16 @@ +0.14.0.0 +-------- + +- integration over infinite intervals + +- msadams and msbdf methods for ode + +- Numeric.LinearAlgebra.Util + +- (<\>) extended to multiple right-hand sides + +- orth + 0.13.0.0 -------- diff --git a/INSTALL.md b/INSTALL.md --- a/INSTALL.md +++ b/INSTALL.md @@ -61,8 +61,8 @@ [install]: http://code.haskell.org/hmatrix/INSTALL [install2]: http://patch-tag.com/r/aruiz/hmatrix/snapshot/current/content/pretty/INSTALL -[winpack]: http://perception.inf.um.es/hmatrix/gsl-lapack-windows.zip - +[winpack2]: http://perception.inf.um.es/hmatrix/gsl-lapack-windows.zip +[winpack]: https://github.com/downloads/AlbertoRuiz/hmatrix/gsl-lapack-windows.zip ## Tests ############################################### diff --git a/THANKS.md b/THANKS.md --- a/THANKS.md +++ b/THANKS.md @@ -107,3 +107,9 @@ - Daniel Fischer reported some Haddock markup errors. +- Danny Chan added support for integration over infinite intervals. + +- Clark Gaebel removed superfluous thread safety. + +- Jeffrey Burdges reported a glpk link problem on OS/X + diff --git a/examples/Real.hs b/examples/Real.hs deleted file mode 100644 --- a/examples/Real.hs +++ /dev/null @@ -1,104 +0,0 @@ - --- Alternative interface and utilities for creation of real arrays, useful to work in interactive mode. - -module Real( - module Numeric.LinearAlgebra, - (<>), (*>), (<*), (<\>), (\>), - vector, - eye, - zeros, ones, - diagl, - row, - col, - (#),(&), (//), blocks, - rand -) where - -import Numeric.LinearAlgebra hiding ((<>), (<|>), (<->), (<\>), (.*), (*/)) -import System.Random(randomIO) - -infixl 7 <> --- | Matrix product ('multiply') -(<>) :: Field t => Matrix t -> Matrix t -> Matrix t -(<>) = multiply - -infixl 7 *> --- | matrix x vector -(*>) :: Field t => Matrix t -> Vector t -> Vector t -m *> v = flatten $ m <> (asColumn v) - -infixl 7 <* --- | vector x matrix -(<*) :: Field t => Vector t -> Matrix t -> Vector t -v <* m = flatten $ (asRow v) <> m - - --- | Least squares solution of a linear system for several right-hand sides, similar to the \\ operator of Matlab\/Octave. (\<\\\>) = 'linearSolveSVD'. -(<\>) :: (Field a) => Matrix a -> Matrix a -> Matrix a -infixl 7 <\> -(<\>) = linearSolveSVD - --- | Least squares solution of a linear system for a single right-hand side. See '(\<\\\>)'. -(\>) :: (Field a) => Matrix a -> Vector a -> Vector a -infixl 7 \> -m \> v = flatten (m <\> reshape 1 v) - --- | Pseudorandom matrix with uniform elements between 0 and 1. -rand :: Int -- ^ rows - -> Int -- ^ columns - -> IO (Matrix Double) -rand r c = do - seed <- randomIO - return (reshape c $ randomVector seed Uniform (r*c)) - --- | Real identity matrix. -eye :: Int -> Matrix Double -eye = ident - --- | Create a real vector from a list. -vector :: [Double] -> Vector Double -vector = fromList - --- | Create a real diagonal matrix from a list. -diagl :: [Double] -> Matrix Double -diagl = diag . vector - --- | Create a matrix or zeros. -zeros :: Int -- ^ rows - -> Int -- ^ columns - -> Matrix Double -zeros r c = konst 0 (r,c) - --- | Create a matrix or ones. -ones :: Int -- ^ rows - -> Int -- ^ columns - -> Matrix Double -ones r c = konst 1 (r,c) - --- | Concatenation of real vectors. -infixl 9 # -(#) :: Vector Double -> Vector Double -> Vector Double -a # b = join [a,b] - --- | Horizontal concatenation of real matrices. -infixl 8 & -(&) :: Matrix Double -> Matrix Double -> Matrix Double -a & b = fromBlocks [[a,b]] - --- | Vertical concatenation of real matrices. -infixl 7 // -(//) :: Matrix Double -> Matrix Double -> Matrix Double -a // b = fromBlocks [[a],[b]] - --- | Real block matrix from a rectangular list of lists. -blocks :: [[Matrix Double]] -> Matrix Double -blocks = fromBlocks - --- | A real matrix with a single row, create from a list of elements. -row :: [Double] -> Matrix Double -row = asRow . vector - --- | A real matrix with a single column, created from a list of elements. -col :: [Double] -> Matrix Double -col = asColumn . vector - diff --git a/examples/devel/ej1/wrappers.hs b/examples/devel/ej1/wrappers.hs --- a/examples/devel/ej1/wrappers.hs +++ b/examples/devel/ej1/wrappers.hs @@ -15,7 +15,7 @@ ----------------------------------------------------- -foreign import ccall "c_scale_vector" +foreign import ccall unsafe "c_scale_vector" cScaleVector :: Double -- scale -> CInt -> Ptr Double -- argument -> CInt -> Ptr Double -- result @@ -29,7 +29,7 @@ ----------------------------------------------------- -- forcing row order -foreign import ccall "c_diag" +foreign import ccall unsafe "c_diag" cDiag :: CInt -> CInt -> Ptr Double -- argument -> CInt -> Ptr Double -- result1 -> CInt -> CInt -> Ptr Double -- result2 diff --git a/examples/devel/ej2/wrappers.hs b/examples/devel/ej2/wrappers.hs --- a/examples/devel/ej2/wrappers.hs +++ b/examples/devel/ej2/wrappers.hs @@ -15,7 +15,7 @@ ----------------------------------------------------- -- arbitrary data order -foreign import ccall "c_diag" +foreign import ccall unsafe "c_diag" cDiag :: CInt -- matrix order -> CInt -> CInt -> Ptr Double -- argument -> CInt -> Ptr Double -- result1 diff --git a/examples/ode.hs b/examples/ode.hs --- a/examples/ode.hs +++ b/examples/ode.hs @@ -1,6 +1,9 @@ +{-# LANGUAGE ViewPatterns #-} import Numeric.GSL.ODE import Numeric.LinearAlgebra import Graphics.Plot +import Debug.Trace(trace) +debug x = trace (show x) x vanderpol mu = do let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)] @@ -32,3 +35,16 @@ harmonic 1 0.1 kepler 0.3 60 kepler 0.4 70 + vanderpol' 2 + +-- example of odeSolveV with jacobian +vanderpol' mu = do + let xdot mu t (toList->[x,v]) = fromList [v, -x + mu * v * (1-x^2)] + jac t (toList->[x,v]) = (2><2) [ 0 , 1 + , -1-2*x*v*mu, mu*(1-x**2) ] + ts = linspace 1000 (0,50) + hi = (ts@>1 - ts@>0)/100 + sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts + mplot sol + + diff --git a/hmatrix.cabal b/hmatrix.cabal --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -1,11 +1,11 @@ Name: hmatrix -Version: 0.13.1.0 +Version: 0.14.0.0 License: GPL License-file: LICENSE Author: Alberto Ruiz Maintainer: Alberto Ruiz Stability: provisional -Homepage: http://perception.inf.um.es/hmatrix +Homepage: https://github.com/albertoruiz/hmatrix Synopsis: Linear algebra and numerical computation Description: Purely functional interface to basic linear algebra and other numerical computations, internally implemented using @@ -49,7 +49,6 @@ examples/devel/ej1/functions.c examples/devel/ej2/wrappers.hs examples/devel/ej2/functions.c - examples/Real.hs examples/vector.hs examples/monadic.hs examples/bool.hs @@ -87,7 +86,7 @@ Build-Depends: base >= 4 && < 5, array, storable-complex, - process, + process, random, vector >= 0.8, binary @@ -111,6 +110,7 @@ Numeric.LinearAlgebra, Numeric.LinearAlgebra.LAPACK, Numeric.LinearAlgebra.Algorithms, + Numeric.LinearAlgebra.Util, Graphics.Plot, Data.Packed.ST, Data.Packed.Development @@ -194,7 +194,7 @@ source-repository head type: git - location: https://github.com/AlbertoRuiz/hmatrix + location: https://github.com/albertoruiz/hmatrix -- The tests are in package hmatrix-tests diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs --- a/lib/Data/Packed/Internal/Common.hs +++ b/lib/Data/Packed/Internal/Common.hs @@ -140,7 +140,7 @@ -- | clear the fpu -foreign import ccall "asm_finit" finit :: IO () +foreign import ccall unsafe "asm_finit" finit :: IO () -- | check the error code check :: String -> IO CInt -> IO () @@ -158,7 +158,7 @@ return () -- | description of GSL error codes -foreign import ccall "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar) +foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar) -- | Error capture and conversion to Maybe mbCatch :: IO x -> IO (Maybe x) diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs @@ -331,11 +331,11 @@ sz = sizeOf (d @> 0) noneed = r1 == 1 || c1 == 1 -foreign import ccall "transF" ctransF :: TFMFM -foreign import ccall "transR" ctransR :: TMM -foreign import ccall "transQ" ctransQ :: TQMQM -foreign import ccall "transC" ctransC :: TCMCM -foreign import ccall "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt +foreign import ccall unsafe "transF" ctransF :: TFMFM +foreign import ccall unsafe "transR" ctransR :: TMM +foreign import ccall unsafe "transQ" ctransQ :: TQMQM +foreign import ccall unsafe "transC" ctransC :: TCMCM +foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt ---------------------------------------------------------------------- @@ -358,19 +358,19 @@ constantF :: Float -> Int -> Vector Float constantF = constantAux cconstantF -foreign import ccall "constantF" cconstantF :: Ptr Float -> TF +foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF constantR :: Double -> Int -> Vector Double constantR = constantAux cconstantR -foreign import ccall "constantR" cconstantR :: Ptr Double -> TV +foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV constantQ :: Complex Float -> Int -> Vector (Complex Float) constantQ = constantAux cconstantQ -foreign import ccall "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV +foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV constantC :: Complex Double -> Int -> Vector (Complex Double) constantC = constantAux cconstantC -foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV +foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV constantP :: Storable a => a -> Int -> Vector a constantP a n = unsafePerformIO $ do @@ -381,7 +381,7 @@ poke k a cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP" return v -foreign import ccall "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt +foreign import ccall unsafe "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt ---------------------------------------------------------------------- @@ -427,7 +427,7 @@ free charname free charfmt -foreign import ccall "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM +foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM ---------------------------------------------------------------------- diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs --- a/lib/Data/Packed/Internal/Vector.hs +++ b/lib/Data/Packed/Internal/Vector.hs @@ -308,8 +308,8 @@ return r -foreign import ccall "float2double" c_float2double:: TFV -foreign import ccall "double2float" c_double2float:: TVF +foreign import ccall unsafe "float2double" c_float2double:: TFV +foreign import ccall unsafe "double2float" c_double2float:: TVF --------------------------------------------------------------- @@ -325,8 +325,8 @@ app2 c_stepD vec v vec r "stepD" return r -foreign import ccall "stepF" c_stepF :: TFF -foreign import ccall "stepD" c_stepD :: TVV +foreign import ccall unsafe "stepF" c_stepF :: TFF +foreign import ccall unsafe "stepD" c_stepD :: TVV --------------------------------------------------------------- @@ -342,8 +342,8 @@ app6 c_condD vec x vec y vec l vec e vec g vec r "condD" return r -foreign import ccall "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF -foreign import ccall "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV +foreign import ccall unsafe "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF +foreign import ccall unsafe "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV -------------------------------------------------------------------------------- @@ -354,11 +354,11 @@ conjugateQ :: Vector (Complex Float) -> Vector (Complex Float) conjugateQ = conjugateAux c_conjugateQ -foreign import ccall "conjugateQ" c_conjugateQ :: TQVQV +foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV conjugateC :: Vector (Complex Double) -> Vector (Complex Double) conjugateC = conjugateAux c_conjugateC -foreign import ccall "conjugateC" c_conjugateC :: TCVCV +foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV -------------------------------------------------------------------------------- @@ -547,7 +547,7 @@ free charname return res -foreign import ccall "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV +foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV -- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file. fprintfVector :: FilePath -> String -> Vector Double -> IO () @@ -558,7 +558,7 @@ free charname free charfmt -foreign import ccall "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV +foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV -- | Loads a vector from a binary file (the number of elements must be known in advance). freadVector :: FilePath -> Int -> IO (Vector Double) @@ -569,7 +569,7 @@ free charname return res -foreign import ccall "vector_fread" gsl_vector_fread:: Ptr CChar -> TV +foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV -- | Saves the elements of a vector to a binary file. fwriteVector :: FilePath -> Vector Double -> IO () @@ -578,5 +578,5 @@ app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite" free charname -foreign import ccall "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV +foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs --- a/lib/Numeric/Container.hs +++ b/lib/Numeric/Container.hs @@ -36,7 +36,7 @@ -- * Matrix product Product(..), optimiseMult, - mXm,mXv,vXm,(<.>),Mul(..),(<\>), + mXm,mXv,vXm,(<.>),Mul(..),LSDiv(..), outer, kronecker, -- * Random numbers RandDist(..), @@ -120,9 +120,15 @@ -------------------------------------------------------- --- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD). -(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a -infixl 7 <\> -m <\> v = flatten (linearSolveSVD m (reshape 1 v)) +class LSDiv b c | b -> c, c->b where + infixl 7 <\> + -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) + (<\>) :: Field t => Matrix t -> b t -> c t + +instance LSDiv Vector Vector where + m <\> v = flatten (linearSolveSVD m (reshape 1 v)) + +instance LSDiv Matrix Matrix where + (<\>) = linearSolveSVD -------------------------------------------------------- diff --git a/lib/Numeric/GSL.hs b/lib/Numeric/GSL.hs --- a/lib/Numeric/GSL.hs +++ b/lib/Numeric/GSL.hs @@ -40,4 +40,4 @@ -- | This action removes the GSL default error handler (which aborts the program), so that -- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort. -foreign import ccall "GSL/gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO () +foreign import ccall unsafe "GSL/gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO () diff --git a/lib/Numeric/GSL/Differentiation.hs b/lib/Numeric/GSL/Differentiation.hs --- a/lib/Numeric/GSL/Differentiation.hs +++ b/lib/Numeric/GSL/Differentiation.hs @@ -48,7 +48,7 @@ freeHaskellFunPtr fp return result -foreign import ccall "gsl-aux.h deriv" +foreign import ccall safe "gsl-aux.h deriv" c_deriv :: CInt -> FunPtr (Double -> Ptr () -> Double) -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt @@ -84,4 +84,4 @@ {- | conversion of Haskell functions into function pointers that can be used in the C side -} -foreign import ccall "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) +foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) diff --git a/lib/Numeric/GSL/Fitting.hs b/lib/Numeric/GSL/Fitting.hs --- a/lib/Numeric/GSL/Fitting.hs +++ b/lib/Numeric/GSL/Fitting.hs @@ -90,7 +90,7 @@ freeHaskellFunPtr jp return (subVector 2 p sol, path) -foreign import ccall "nlfit" +foreign import ccall safe "nlfit" c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM ------------------------------------------------------- diff --git a/lib/Numeric/GSL/Fourier.hs b/lib/Numeric/GSL/Fourier.hs --- a/lib/Numeric/GSL/Fourier.hs +++ b/lib/Numeric/GSL/Fourier.hs @@ -30,7 +30,7 @@ app2 (c_fft code) vec v vec r "fft" return r -foreign import ccall "gsl-aux.h fft" c_fft :: CInt -> TCVCV +foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCVCV {- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave. diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs --- a/lib/Numeric/GSL/Integration.hs +++ b/lib/Numeric/GSL/Integration.hs @@ -17,7 +17,10 @@ module Numeric.GSL.Integration ( integrateQNG, - integrateQAGS + integrateQAGS, + integrateQAGI, + integrateQAGIU, + integrateQAGIL ) where import Foreign.C.Types @@ -29,7 +32,7 @@ {- | conversion of Haskell functions into function pointers that can be used in the C side -} -foreign import ccall "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) +foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) -------------------------------------------------------------------- {- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example: @@ -60,7 +63,7 @@ freeHaskellFunPtr fp return result -foreign import ccall "gsl-aux.h integrate_qags" +foreign import ccall safe "gsl-aux.h integrate_qags" c_integrate_qags :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt @@ -91,6 +94,105 @@ freeHaskellFunPtr fp return result -foreign import ccall "gsl-aux.h integrate_qng" +foreign import ccall safe "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 safe "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 safe "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 safe "gsl-aux.h integrate_qagil" + c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt + diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs --- a/lib/Numeric/GSL/Internal.hs +++ b/lib/Numeric/GSL/Internal.hs @@ -26,14 +26,14 @@ return 0 -- | conversion of Haskell functions into function pointers that can be used in the C side -foreign import ccall "wrapper" +foreign import ccall safe "wrapper" mkVecfun :: (CInt -> Ptr Double -> Double) -> IO( FunPtr (CInt -> Ptr Double -> Double)) -foreign import ccall "wrapper" +foreign import ccall safe "wrapper" mkVecVecfun :: TVV -> IO (FunPtr TVV) -foreign import ccall "wrapper" +foreign import ccall safe "wrapper" mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV)) aux_vTov :: (Vector Double -> Vector Double) -> TVV @@ -46,10 +46,10 @@ g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr) return 0 -foreign import ccall "wrapper" +foreign import ccall safe "wrapper" mkVecMatfun :: TVM -> IO (FunPtr TVM) -foreign import ccall "wrapper" +foreign import ccall safe "wrapper" mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM)) aux_vTom :: (Vector Double -> Matrix Double) -> TVM diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs --- a/lib/Numeric/GSL/Minimization.hs +++ b/lib/Numeric/GSL/Minimization.hs @@ -122,7 +122,7 @@ return (sol, path) -foreign import ccall "gsl-aux.h minimize" +foreign import ccall safe "gsl-aux.h minimize" c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM ---------------------------------------------------------------------------------- @@ -179,7 +179,7 @@ freeHaskellFunPtr dfp return (sol,path) -foreign import ccall "gsl-aux.h minimizeD" +foreign import ccall safe "gsl-aux.h minimizeD" c_minimizeD :: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> FunPtr TVV diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs --- a/lib/Numeric/GSL/ODE.hs +++ b/lib/Numeric/GSL/ODE.hs @@ -19,7 +19,7 @@ xdot t [x,v] = [v, -0.95*x - 0.1*v] -ts = linspace 100 (0,20) +ts = linspace 100 (0,20 :: Double) sol = odeSolve xdot [10,0] ts @@ -29,7 +29,7 @@ ----------------------------------------------------------------------------- module Numeric.GSL.ODE ( - odeSolve, odeSolveV, ODEMethod(..) + odeSolve, odeSolveV, ODEMethod(..), Jacobian ) where import Data.Packed.Internal @@ -41,26 +41,29 @@ ------------------------------------------------------------------------- +type Jacobian = Double -> Vector Double -> Matrix Double + -- | Stepping functions data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method. - | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use 'RKf45'. + | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use the embedded methods. | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator. | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method. | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method. - | RK2imp -- ^ Implicit 2nd order Runge-Kutta at Gaussian points. - | RK4imp -- ^ Implicit 4th order Runge-Kutta at Gaussian points. - | BSimp -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian. - | Gear1 -- ^ M=1 implicit Gear method. - | Gear2 -- ^ M=2 implicit Gear method. - deriving (Enum,Eq,Show,Bounded) + | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points. + | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points. + | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems. + | RK1imp Jacobian -- ^ Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method. + | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12. + | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems. + -- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. odeSolve :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x) -> [Double] -- ^ initial conditions -> Vector Double -- ^ desired solution times -> Matrix Double -- ^ solution -odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) Nothing (fromList xi) ts +odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts where hi = (ts@>1 - ts@>0)/100 epsAbs = 1.49012e-08 epsRel = 1.49012e-08 @@ -73,11 +76,33 @@ -> Double -- ^ absolute tolerance for the state vector -> Double -- ^ relative tolerance for the state vector -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) + -> Vector Double -- ^ initial conditions + -> Vector Double -- ^ desired solution times + -> Matrix Double -- ^ solution +odeSolveV RK2 = odeSolveV' 0 Nothing +odeSolveV RK4 = odeSolveV' 1 Nothing +odeSolveV RKf45 = odeSolveV' 2 Nothing +odeSolveV RKck = odeSolveV' 3 Nothing +odeSolveV RK8pd = odeSolveV' 4 Nothing +odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) +odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) +odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) +odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac) +odeSolveV MSAdams = odeSolveV' 9 Nothing +odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac) + + +odeSolveV' + :: CInt -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian + -> Double -- ^ initial step size + -> Double -- ^ absolute tolerance for the state vector + -> Double -- ^ relative tolerance for the state vector + -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) -> Vector Double -- ^ initial conditions -> Vector Double -- ^ desired solution times -> Matrix Double -- ^ solution -odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do +odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do let n = dim xiv fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) jp <- case mbjac of @@ -86,12 +111,12 @@ sol <- vec xiv $ \xiv' -> vec (checkTimes ts) $ \ts' -> createMIO (dim ts) n - (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' ) + (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) "ode" freeHaskellFunPtr fp return sol -foreign import ccall "ode" +foreign import ccall safe "ode" ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM ------------------------------------------------------- diff --git a/lib/Numeric/GSL/Polynomials.hs b/lib/Numeric/GSL/Polynomials.hs --- a/lib/Numeric/GSL/Polynomials.hs +++ b/lib/Numeric/GSL/Polynomials.hs @@ -55,4 +55,4 @@ return r | otherwise = error "polySolve on a polynomial of degree zero" -foreign import ccall "gsl-aux.h polySolve" c_polySolve:: TVCV +foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TVCV diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs --- a/lib/Numeric/GSL/Root.hs +++ b/lib/Numeric/GSL/Root.hs @@ -91,7 +91,7 @@ return (take n $ drop 1 sol, path) -foreign import ccall "root" +foreign import ccall safe "root" c_root:: CInt -> FunPtr TVV -> Double -> CInt -> TVM ------------------------------------------------------------------------- @@ -130,7 +130,7 @@ return (take n $ drop 1 sol, path) -foreign import ccall "rootj" +foreign import ccall safe "rootj" c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM ------------------------------------------------------- diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs --- a/lib/Numeric/GSL/Vector.hs +++ b/lib/Numeric/GSL/Vector.hs @@ -109,10 +109,10 @@ app2 c_sumC vec x vec r "sumC" return $ r @> 0 -foreign import ccall safe "gsl-aux.h sumF" c_sumF :: TFF -foreign import ccall safe "gsl-aux.h sumR" c_sumR :: TVV -foreign import ccall safe "gsl-aux.h sumQ" c_sumQ :: TQVQV -foreign import ccall safe "gsl-aux.h sumC" c_sumC :: TCVCV +foreign import ccall unsafe "gsl-aux.h sumF" c_sumF :: TFF +foreign import ccall unsafe "gsl-aux.h sumR" c_sumR :: TVV +foreign import ccall unsafe "gsl-aux.h sumQ" c_sumQ :: TQVQV +foreign import ccall unsafe "gsl-aux.h sumC" c_sumC :: TCVCV -- | product of elements prodF :: Vector Float -> Float @@ -142,10 +142,10 @@ app2 c_prodC vec x vec r "prodC" return $ r @> 0 -foreign import ccall safe "gsl-aux.h prodF" c_prodF :: TFF -foreign import ccall safe "gsl-aux.h prodR" c_prodR :: TVV -foreign import ccall safe "gsl-aux.h prodQ" c_prodQ :: TQVQV -foreign import ccall safe "gsl-aux.h prodC" c_prodC :: TCVCV +foreign import ccall unsafe "gsl-aux.h prodF" c_prodF :: TFF +foreign import ccall unsafe "gsl-aux.h prodR" c_prodR :: TVV +foreign import ccall unsafe "gsl-aux.h prodQ" c_prodQ :: TQVQV +foreign import ccall unsafe "gsl-aux.h prodC" c_prodC :: TCVCV -- | dot product dotF :: Vector Float -> Vector Float -> Float @@ -175,10 +175,10 @@ app3 c_dotC vec x vec y vec r "dotC" return $ r @> 0 -foreign import ccall safe "gsl-aux.h dotF" c_dotF :: TFFF -foreign import ccall safe "gsl-aux.h dotR" c_dotR :: TVVV -foreign import ccall safe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV -foreign import ccall safe "gsl-aux.h dotC" c_dotC :: TCVCVCV +foreign import ccall unsafe "gsl-aux.h dotF" c_dotF :: TFFF +foreign import ccall unsafe "gsl-aux.h dotR" c_dotR :: TVVV +foreign import ccall unsafe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV +foreign import ccall unsafe "gsl-aux.h dotC" c_dotC :: TCVCVCV ------------------------------------------------------------------ @@ -210,25 +210,25 @@ toScalarR :: FunCodeS -> Vector Double -> Double toScalarR oper = toScalarAux c_toScalarR (fromei oper) -foreign import ccall safe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV +foreign import ccall unsafe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV -- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc. toScalarF :: FunCodeS -> Vector Float -> Float toScalarF oper = toScalarAux c_toScalarF (fromei oper) -foreign import ccall safe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF +foreign import ccall unsafe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF -- | obtains different functions of a vector: only norm1, norm2 toScalarC :: FunCodeS -> Vector (Complex Double) -> Double toScalarC oper = toScalarAux c_toScalarC (fromei oper) -foreign import ccall safe "gsl-aux.h toScalarC" c_toScalarC :: CInt -> TCVV +foreign import ccall unsafe "gsl-aux.h toScalarC" c_toScalarC :: CInt -> TCVV -- | obtains different functions of a vector: only norm1, norm2 toScalarQ :: FunCodeS -> Vector (Complex Float) -> Float toScalarQ oper = toScalarAux c_toScalarQ (fromei oper) -foreign import ccall safe "gsl-aux.h toScalarQ" c_toScalarQ :: CInt -> TQVF +foreign import ccall unsafe "gsl-aux.h toScalarQ" c_toScalarQ :: CInt -> TQVF ------------------------------------------------------------------ @@ -236,25 +236,25 @@ vectorMapR :: FunCodeV -> Vector Double -> Vector Double vectorMapR = vectorMapAux c_vectorMapR -foreign import ccall safe "gsl-aux.h mapR" c_vectorMapR :: CInt -> TVV +foreign import ccall unsafe "gsl-aux.h mapR" c_vectorMapR :: CInt -> TVV -- | map of complex vectors with given function vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double) vectorMapC oper = vectorMapAux c_vectorMapC (fromei oper) -foreign import ccall safe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV +foreign import ccall unsafe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV -- | map of real vectors with given function vectorMapF :: FunCodeV -> Vector Float -> Vector Float vectorMapF = vectorMapAux c_vectorMapF -foreign import ccall safe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF +foreign import ccall unsafe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF -- | map of real vectors with given function vectorMapQ :: FunCodeV -> Vector (Complex Float) -> Vector (Complex Float) vectorMapQ = vectorMapAux c_vectorMapQ -foreign import ccall safe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV +foreign import ccall unsafe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV ------------------------------------------------------------------- @@ -262,25 +262,25 @@ vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromei oper) -foreign import ccall safe "gsl-aux.h mapValR" c_vectorMapValR :: CInt -> Ptr Double -> TVV +foreign import ccall unsafe "gsl-aux.h mapValR" c_vectorMapValR :: CInt -> Ptr Double -> TVV -- | map of complex vectors with given function vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double) vectorMapValC = vectorMapValAux c_vectorMapValC -foreign import ccall safe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV +foreign import ccall unsafe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV -- | map of real vectors with given function vectorMapValF :: FunCodeSV -> Float -> Vector Float -> Vector Float vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper) -foreign import ccall safe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF +foreign import ccall unsafe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF -- | map of complex vectors with given function vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float) vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper) -foreign import ccall safe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV +foreign import ccall unsafe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV ------------------------------------------------------------------- @@ -288,25 +288,25 @@ vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double vectorZipR = vectorZipAux c_vectorZipR -foreign import ccall safe "gsl-aux.h zipR" c_vectorZipR :: CInt -> TVVV +foreign import ccall unsafe "gsl-aux.h zipR" c_vectorZipR :: CInt -> TVVV -- | elementwise operation on complex vectors vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) vectorZipC = vectorZipAux c_vectorZipC -foreign import ccall safe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV +foreign import ccall unsafe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV -- | elementwise operation on real vectors vectorZipF :: FunCodeVV -> Vector Float -> Vector Float -> Vector Float vectorZipF = vectorZipAux c_vectorZipF -foreign import ccall safe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF +foreign import ccall unsafe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF -- | elementwise operation on complex vectors vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float) vectorZipQ = vectorZipAux c_vectorZipQ -foreign import ccall safe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV +foreign import ccall unsafe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV ----------------------------------------------------------------------- @@ -324,4 +324,4 @@ app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" return r -foreign import ccall safe "random_vector" c_random_vector :: CInt -> CInt -> TV +foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -32,7 +32,7 @@ #include #include #include -#include +#include #include #include #include @@ -761,6 +761,48 @@ 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); @@ -1314,38 +1356,38 @@ int jac(double, int, const double*, int, int, double*), KRVEC(xi), KRVEC(ts), RMAT(sol)) { - const gsl_odeiv_step_type * T; + const gsl_odeiv2_step_type * T; switch(method) { - case 0 : {T = gsl_odeiv_step_rk2; break; } - case 1 : {T = gsl_odeiv_step_rk4; break; } - case 2 : {T = gsl_odeiv_step_rkf45; break; } - case 3 : {T = gsl_odeiv_step_rkck; break; } - case 4 : {T = gsl_odeiv_step_rk8pd; break; } - case 5 : {T = gsl_odeiv_step_rk2imp; break; } - case 6 : {T = gsl_odeiv_step_rk4imp; break; } - case 7 : {T = gsl_odeiv_step_bsimp; break; } - case 8 : {T = gsl_odeiv_step_gear1; break; } - case 9 : {T = gsl_odeiv_step_gear2; break; } + case 0 : {T = gsl_odeiv2_step_rk2; break; } + case 1 : {T = gsl_odeiv2_step_rk4; break; } + case 2 : {T = gsl_odeiv2_step_rkf45; break; } + case 3 : {T = gsl_odeiv2_step_rkck; break; } + case 4 : {T = gsl_odeiv2_step_rk8pd; break; } + case 5 : {T = gsl_odeiv2_step_rk2imp; break; } + case 6 : {T = gsl_odeiv2_step_rk4imp; break; } + case 7 : {T = gsl_odeiv2_step_bsimp; break; } + case 8 : {T = gsl_odeiv2_step_rk1imp; break; } + case 9 : {T = gsl_odeiv2_step_msadams; break; } + case 10: {T = gsl_odeiv2_step_msbdf; break; } default: ERROR(BAD_CODE); } - - gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); - gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel); - gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); - Tode P; P.f = f; P.j = jac; P.n = xin; - gsl_odeiv_system sys = {odefunc, odejac, xin, &P}; + gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; + gsl_odeiv2_driver * d = + gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); + double t = tsp[0]; double* y = (double*)calloc(xin,sizeof(double)); int i,j; + int status; for(i=0; i< xin; i++) { y[i] = xip[i]; solp[i] = xip[i]; @@ -1354,22 +1396,24 @@ for (i = 1; i < tsn ; i++) { double ti = tsp[i]; - while (t < ti) - { - gsl_odeiv_evolve_apply (e, c, s, - &sys, - &t, ti, &h, - y); - // if (h < hmin) h = hmin; - } + + status = gsl_odeiv2_driver_apply (d, &t, ti, y); + + if (status != GSL_SUCCESS) { + printf ("error in ode, return value=%d\n", status); + break; + } + +// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); + for(j=0; j 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) -{-# DEPRECATED full "use fullSVD instead" #-} -full svdFun m = (u, d ,v) where - (u,s,v) = svdFun m - d = diagRect 0 s r c - r = rows m - c = cols m - -{-# DEPRECATED economy "use compactSVD instead" #-} -economy svdFun m = (u', subVector 0 d s, v') where - (u,s,v) = svdFun m - d = rankSVD (1*eps) m s `max` 1 - u' = takeColumns d u - v' = takeColumns d v - - -------------------------------------------------------------- -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. @@ -444,6 +429,13 @@ -- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision. nullVector :: Field t => Matrix t -> Vector t nullVector = last . nullspacePrec 1 + +orth :: Field t => Matrix t -> [Vector t] +-- ^ Return an orthonormal basis of the range space of a matrix +orth m = take r $ toColumns u + where + (u,s,_) = compactSVD m + r = ranksv eps (max (rows m) (cols m)) (toList s) ------------------------------------------------------------------------ diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs --- a/lib/Numeric/LinearAlgebra/LAPACK.hs +++ b/lib/Numeric/LinearAlgebra/LAPACK.hs @@ -53,10 +53,10 @@ ----------------------------------------------------------------------------------- -foreign import ccall "multiplyR" dgemmc :: CInt -> CInt -> TMMM -foreign import ccall "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM -foreign import ccall "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM -foreign import ccall "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM +foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM +foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM +foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM +foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM isT Matrix{order = ColumnMajor} = 0 isT Matrix{order = RowMajor} = 1 @@ -88,10 +88,10 @@ multiplyQ a b = multiplyAux cgemmc "cgemmc" a b ----------------------------------------------------------------------------- -foreign import ccall "svd_l_R" dgesvd :: TMMVM -foreign import ccall "svd_l_C" zgesvd :: TCMCMVCM -foreign import ccall "svd_l_Rdd" dgesdd :: TMMVM -foreign import ccall "svd_l_Cdd" zgesdd :: TCMCMVCM +foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM +foreign import ccall unsafe "svd_l_C" zgesvd :: TCMCMVCM +foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TMMVM +foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TCMCMVCM -- | Full SVD of a real matrix using LAPACK's /dgesvd/. svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) @@ -211,10 +211,10 @@ ----------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM -foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM -foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: CInt -> TMVM -foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: CInt -> TCMVCM +foreign import ccall unsafe "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM +foreign import ccall unsafe "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM +foreign import ccall unsafe "LAPACK/lapack-aux.h eig_l_S" dsyev :: CInt -> TMVM +foreign import ccall unsafe "LAPACK/lapack-aux.h eig_l_H" zheev :: CInt -> TCMVCM eigAux f st m = unsafePerformIO $ do l <- createVector r @@ -325,10 +325,10 @@ vrev = flatten . flipud . reshape 1 ----------------------------------------------------------------------------- -foreign import ccall "linearSolveR_l" dgesv :: TMMM -foreign import ccall "linearSolveC_l" zgesv :: TCMCMCM -foreign import ccall "cholSolveR_l" dpotrs :: TMMM -foreign import ccall "cholSolveC_l" zpotrs :: TCMCMCM +foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM +foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM +foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM +foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM linearSolveSQAux f st a b | n1==n2 && n1==r = unsafePerformIO $ do @@ -359,10 +359,10 @@ cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b) ----------------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l" dgels :: TMMM -foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l" zgels :: TCMCMCM -foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l" dgelss :: Double -> TMMM -foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l" zgelss :: Double -> TCMCMCM +foreign import ccall unsafe "LAPACK/lapack-aux.h linearSolveLSR_l" dgels :: TMMM +foreign import ccall unsafe "LAPACK/lapack-aux.h linearSolveLSC_l" zgels :: TCMCMCM +foreign import ccall unsafe "LAPACK/lapack-aux.h linearSolveSVDR_l" dgelss :: Double -> TMMM +foreign import ccall unsafe "LAPACK/lapack-aux.h linearSolveSVDC_l" zgelss :: Double -> TCMCMCM linearSolveAux f st a b = unsafePerformIO $ do r <- createMatrix ColumnMajor (max m n) nrhs @@ -401,8 +401,8 @@ linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b) ----------------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM -foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM +foreign import ccall unsafe "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM +foreign import ccall unsafe "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM cholAux f st a = do r <- createMatrix ColumnMajor n n @@ -427,8 +427,8 @@ mbCholS = unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat ----------------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM -foreign import ccall "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM +foreign import ccall unsafe "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM +foreign import ccall unsafe "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM -- | QR factorization of a real matrix, using LAPACK's /dgeqr2/. qrR :: Matrix Double -> (Matrix Double, Vector Double) @@ -448,8 +448,8 @@ mn = min m n ----------------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h hess_l_R" dgehrd :: TMVM -foreign import ccall "LAPACK/lapack-aux.h hess_l_C" zgehrd :: TCMCVCM +foreign import ccall unsafe "LAPACK/lapack-aux.h hess_l_R" dgehrd :: TMVM +foreign import ccall unsafe "LAPACK/lapack-aux.h hess_l_C" zgehrd :: TCMCVCM -- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/. hessR :: Matrix Double -> (Matrix Double, Vector Double) @@ -469,8 +469,8 @@ mn = min m n ----------------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM -foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM +foreign import ccall unsafe "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM +foreign import ccall unsafe "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM -- | Schur factorization of a square real matrix, using LAPACK's /dgees/. schurR :: Matrix Double -> (Matrix Double, Matrix Double) @@ -488,8 +488,8 @@ where n = rows a ----------------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h lu_l_R" dgetrf :: TMVM -foreign import ccall "LAPACK/lapack-aux.h lu_l_C" zgetrf :: TCMVCM +foreign import ccall unsafe "LAPACK/lapack-aux.h lu_l_R" dgetrf :: TMVM +foreign import ccall unsafe "LAPACK/lapack-aux.h lu_l_C" zgetrf :: TCMVCM -- | LU factorization of a general real matrix, using LAPACK's /dgetrf/. luR :: Matrix Double -> (Matrix Double, [Int]) @@ -511,8 +511,8 @@ type TW a = CInt -> PD -> a type TQ a = CInt -> CInt -> PC -> a -foreign import ccall "LAPACK/lapack-aux.h luS_l_R" dgetrs :: TMVMM -foreign import ccall "LAPACK/lapack-aux.h luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt)))) +foreign import ccall unsafe "LAPACK/lapack-aux.h luS_l_R" dgetrs :: TMVMM +foreign import ccall unsafe "LAPACK/lapack-aux.h luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt)))) -- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/. lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double diff --git a/lib/Numeric/LinearAlgebra/Util.hs b/lib/Numeric/LinearAlgebra/Util.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/LinearAlgebra/Util.hs @@ -0,0 +1,104 @@ +----------------------------------------------------------------------------- +{- | +Module : Numeric.LinearAlgebra.Util +Copyright : (c) Alberto Ruiz 2012 +License : GPL + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional + +-} +----------------------------------------------------------------------------- + +module Numeric.LinearAlgebra.Util( + disp, + zeros, ones, + diagl, + row, + col, + (&),(!), (#), + rand, randn, + cross, + norm +) where + +import Numeric.LinearAlgebra +import System.Random(randomIO) + + +disp :: Int -> Matrix Double -> IO () +-- ^ show a matrix with given number of digits after the decimal point +disp n = putStrLn . dispf n + +-- | pseudorandom matrix with uniform elements between 0 and 1 +randm :: RandDist + -> Int -- ^ rows + -> Int -- ^ columns + -> IO (Matrix Double) +randm d r c = do + seed <- randomIO + return (reshape c $ randomVector seed d (r*c)) + +-- | pseudorandom matrix with uniform elements between 0 and 1 +rand :: Int -> Int -> IO (Matrix Double) +rand = randm Uniform + +-- | pseudorandom matrix with normal elements +randn :: Int -> Int -> IO (Matrix Double) +randn = randm Gaussian + +-- | create a real diagonal matrix from a list +diagl :: [Double] -> Matrix Double +diagl = diag . fromList + +-- | a real matrix of zeros +zeros :: Int -- ^ rows + -> Int -- ^ columns + -> Matrix Double +zeros r c = konst 0 (r,c) + +-- | a real matrix of ones +ones :: Int -- ^ rows + -> Int -- ^ columns + -> Matrix Double +ones r c = konst 1 (r,c) + +-- | concatenation of real vectors +infixl 3 & +(&) :: Vector Double -> Vector Double -> Vector Double +a & b = join [a,b] + +-- | horizontal concatenation of real matrices +infixl 3 ! +(!) :: Matrix Double -> Matrix Double -> Matrix Double +a ! b = fromBlocks [[a,b]] + +-- | vertical concatenation of real matrices +(#) :: Matrix Double -> Matrix Double -> Matrix Double +infixl 2 # +a # b = fromBlocks [[a],[b]] + +-- | create a single row real matrix from a list +row :: [Double] -> Matrix Double +row = asRow . fromList + +-- | create a single column real matrix from a list +col :: [Double] -> Matrix Double +col = asColumn . fromList + +cross :: Vector Double -> Vector Double -> Vector Double +-- ^ cross product of dimension 3 real vectors +cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3] + | otherwise = error $ "cross ("++show x++") ("++show y++")" + where + [x1,x2,x3] = toList x + [y1,y2,y3] = toList y + z1 = x2*y3-x3*y2 + z2 = x3*y1-x1*y3 + z3 = x1*y2-x2*y1 + +norm :: Vector Double -> Double +-- ^ 2-norm of real vectors +norm = pnorm PNorm2 + +