hmatrix 0.13.1.0 → 0.14.0.0

raw patch · 27 files changed

+489/−285 lines, 27 filesdep +randomPVP ok

version bump matches the API change (PVP)

Dependencies added: random

API changes (from Hackage documentation)

Files

• CHANGES.md +13/−0
• INSTALL.md +2/−2
• THANKS.md +6/−0
• examples/Real.hs +0/−104
• examples/devel/ej1/wrappers.hs +2/−2
• examples/devel/ej2/wrappers.hs +1/−1
• examples/ode.hs +16/−0
• hmatrix.cabal +5/−5
• lib/Data/Packed/Internal/Common.hs +2/−2
• lib/Data/Packed/Internal/Matrix.hs +11/−11
• lib/Data/Packed/Internal/Vector.hs +12/−12
• lib/Numeric/Container.hs +11/−5
• lib/Numeric/GSL.hs +1/−1
• lib/Numeric/GSL/Differentiation.hs +2/−2
• lib/Numeric/GSL/Fitting.hs +1/−1
• lib/Numeric/GSL/Fourier.hs +1/−1
• lib/Numeric/GSL/Integration.hs +106/−4
• lib/Numeric/GSL/Internal.hs +5/−5
• lib/Numeric/GSL/Minimization.hs +2/−2
• lib/Numeric/GSL/ODE.hs +38/−13
• lib/Numeric/GSL/Polynomials.hs +1/−1
• lib/Numeric/GSL/Root.hs +2/−2
• lib/Numeric/GSL/Vector.hs +29/−29
• lib/Numeric/GSL/gsl-aux.c +74/−30
• lib/Numeric/LinearAlgebra/Algorithms.hs +10/−18
• lib/Numeric/LinearAlgebra/LAPACK.hs +32/−32
• lib/Numeric/LinearAlgebra/Util.hs +104/−0

@@ -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 -------- 

@@ -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 ############################################### 

@@ -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+

@@ -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-

@@ -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

@@ -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

@@ -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++

@@ -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 <aruiz@um.es> 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 

@@ -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)

@@ -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  ---------------------------------------------------------------------- 

@@ -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 

@@ -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  --------------------------------------------------------

@@ -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 ()

@@ -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)) 

@@ -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  -------------------------------------------------------

@@ -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.

@@ -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+

@@ -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

@@ -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

@@ -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  -------------------------------------------------------

@@ -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

@@ -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  -------------------------------------------------------

@@ -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

@@ -32,7 +32,7 @@ #include <gsl/gsl_complex_math.h> #include <gsl/gsl_rng.h> #include <gsl/gsl_randist.h>-#include <gsl/gsl_odeiv.h>+#include <gsl/gsl_odeiv2.h> #include <gsl/gsl_multifit_nlin.h> #include <string.h> #include <stdio.h>@@ -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<xin; j++) {                solp[i*xin + j] = y[j];            }          }      free(y);-    gsl_odeiv_evolve_free (e);-    gsl_odeiv_control_free (c);-    gsl_odeiv_step_free (s);-    return 0;+    gsl_odeiv2_driver_free (d);+    +    return status; }+

@@ -64,6 +64,7 @@     nullspacePrec,     nullVector,     nullspaceSVD,+    orth, -- * Norms     Normed(..), NormType(..),     relativeError,@@ -73,8 +74,7 @@     haussholder,     unpackQR, unpackHess,     pinvTol,-    ranksv,-    full, economy+    ranksv ) where  @@ -209,27 +209,12 @@ 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)  -{-# 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)  ------------------------------------------------------------------------ 

@@ -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

@@ -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++