From 283f3033f86fabde2290bb28a59e7d87fd0754f5 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 1 Mar 2010 11:15:22 +0000 Subject: compatible with vector --- CHANGES | 11 +++ examples/vector.hs | 7 +- hmatrix.cabal | 4 +- lib/Data/Packed/Convert.hs | 101 -------------------------- lib/Data/Packed/Internal/Matrix.hs | 24 +++--- lib/Data/Packed/Internal/Vector.hs | 75 +++++++++++-------- lib/Data/Packed/ST.hs | 4 +- lib/Numeric/GSL/Internal.hs | 8 +- lib/Numeric/GSL/Minimization.hs | 4 +- lib/Numeric/GSL/ODE.hs | 4 +- lib/Numeric/GSL/Root.hs | 4 +- lib/Numeric/LinearAlgebra/Instances.hs | 14 ++-- lib/Numeric/LinearAlgebra/Tests.hs | 14 ++++ lib/Numeric/LinearAlgebra/Tests/Properties.hs | 4 + 14 files changed, 110 insertions(+), 168 deletions(-) delete mode 100644 lib/Data/Packed/Convert.hs diff --git a/CHANGES b/CHANGES index 32d34c0..fc089e0 100644 --- a/CHANGES +++ b/CHANGES @@ -1,3 +1,14 @@ +0.9.1.0 +======= + +- GSL special functions moved to separate package hmatrix-special. + +- Added offset of Vector, allowing fast, noncopy subVector (slice). + Vector is now identical to Roman Leshchinskiy's Data.Vector.Storable.Vector, + so we can convert from/to them in O(1). + +- Removed Data.Packed.Convert, see examples/vector.hs + 0.8.3.0 ======= diff --git a/examples/vector.hs b/examples/vector.hs index 12cbc42..855c6b4 100644 --- a/examples/vector.hs +++ b/examples/vector.hs @@ -12,7 +12,7 @@ import qualified Data.Vector.Storable as V fromVector :: Storable t => V.Vector t -> H.Vector t fromVector v = unsafeFromForeignPtr p i n where - (p,i,n) = V.unsafeToForeignPtr (V.copy v) + (p,i,n) = V.unsafeToForeignPtr v toVector :: H.Vector t -> V.Vector t toVector v = V.unsafeFromForeignPtr p i n where @@ -22,8 +22,11 @@ toVector v = V.unsafeFromForeignPtr p i n where v = V.slice 5 10 (V.fromList [1 .. 10::Double] V.++ V.replicate 10 7) -w = linspace 5 (0,2) +w = subVector 2 3 (linspace 10 (0,2)) main = do + print v print $ fromVector v + print w print $ toVector w + diff --git a/hmatrix.cabal b/hmatrix.cabal index d587a73..392b301 100644 --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -1,5 +1,5 @@ Name: hmatrix -Version: 0.9.0.0 +Version: 0.9.1.0 License: GPL License-file: LICENSE Author: Alberto Ruiz @@ -90,7 +90,7 @@ library Numeric.LinearAlgebra.Interface, Numeric.LinearAlgebra.Algorithms, Graphics.Plot, - Data.Packed.Convert, + -- Data.Packed.Convert, Data.Packed.ST, Data.Packed.Development, Data.Packed.Random diff --git a/lib/Data/Packed/Convert.hs b/lib/Data/Packed/Convert.hs deleted file mode 100644 index 8f0eef6..0000000 --- a/lib/Data/Packed/Convert.hs +++ /dev/null @@ -1,101 +0,0 @@ -{-# OPTIONS -XTypeOperators -XRank2Types -XFlexibleContexts #-} - ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Convert --- Copyright : (c) Alberto Ruiz 2008 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable --- --- Conversion of Vectors and Matrices to and from the standard Haskell arrays. --- (provisional) --- ------------------------------------------------------------------------------ - -module Data.Packed.Convert ( - arrayFromVector, vectorFromArray, - mArrayFromVector, vectorFromMArray, - vectorToStorableArray, storableArrayToVector, - arrayFromMatrix, matrixFromArray, - mArrayFromMatrix, matrixFromMArray, --- matrixToStorableArray, storableArrayToMatrix -) where - -import Data.Packed.Internal -import Data.Array.Storable -import Foreign -import Control.Monad.ST -import Data.Array.ST -import Data.Array.Unboxed - --- | Creates a StorableArray indexed from 0 to dim -1. --- (Memory is efficiently copied, so you can then freely modify the obtained array) -vectorToStorableArray :: Storable t => Vector t -> IO (StorableArray Int t) -vectorToStorableArray v = do - r <- cloneVector v - unsafeForeignPtrToStorableArray (fptr r) (0,dim r -1) - --- | Creates a Vector from a StorableArray. --- (Memory is efficiently copied, so posterior changes in the array will not affect the result) -storableArrayToVector :: Storable t => StorableArray Int t -> IO (Vector t) -storableArrayToVector s = do - (a,b) <- getBounds s - let n = (b-a+1) - r <- createVector n - withStorableArray s $ \p -> do - let f _ d = copyArray d p n >> return 0 - app1 f vec r "storableArrayToVector" - return r - - -unsafeVectorToStorableArray :: Storable t => Vector t -> IO (StorableArray Int t) -unsafeVectorToStorableArray v = unsafeForeignPtrToStorableArray (fptr v) (0,dim v -1) - ---unsafeStorableArrayToVector :: Storable t => StorableArray Int t -> IO (Vector t) ---unsafeStorableArrayToVector s = undefined -- the foreign ptr of Storable Array is not available? - ------------------------------------------------------------------ --- provisional, we need Unboxed arrays for Complex Double - - -unsafeFreeze' :: (MArray a e m, Ix i) => a i e -> m (Array i e) -unsafeFreeze' = unsafeFreeze - --- | creates an immutable Array from an hmatrix Vector (to do: unboxed) -arrayFromVector :: (Storable t) => Vector t -> Array Int t -arrayFromVector x = runSTArray (mArrayFromVector x) - --- | creates a mutable array from an hmatrix Vector (to do: unboxed) -mArrayFromVector :: (MArray b t (ST s), Storable t) => Vector t -> ST s (b Int t) -mArrayFromVector v = unsafeThaw =<< unsafeIOToST ( unsafeFreeze' =<< (vectorToStorableArray $ v)) - - --- (creates an hmatrix Vector from an immutable array (to do: unboxed)) -vectorFromArray :: (Storable t) => Array Int t -> Vector t -vectorFromArray a = unsafePerformIO $ storableArrayToVector =<< unsafeThaw a - --- | creates a mutable Array from an hmatrix Vector for manipulation with runSTUArray (to do: unboxed) -vectorFromMArray :: (Storable t, MArray a t (ST s)) => a Int t -> ST s (Vector t) -vectorFromMArray x = fmap vectorFromArray (unsafeFreeze' x) - --------------------------------------------------------------------- --- provisional - -matrixFromArray :: UArray (Int, Int) Double -> Matrix Double -matrixFromArray m = reshape c . fromList . elems $ m - where ((r1,c1),(r2,c2)) = bounds m - _r = r2-r1+1 - c = c2-c1+1 - -arrayFromMatrix :: Matrix Double -> UArray (Int, Int) Double -arrayFromMatrix m = listArray ((0,0),(rows m -1, cols m -1)) (toList $ flatten m) - - -mArrayFromMatrix :: (MArray b Double m) => Matrix Double -> m (b (Int, Int) Double) -mArrayFromMatrix = unsafeThaw . arrayFromMatrix - -matrixFromMArray :: (MArray a Double (ST s)) => a (Int,Int) Double -> ST s (Matrix Double) -matrixFromMArray x = fmap matrixFromArray (unsafeFreeze x) diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index a6b6215..b8e7ab0 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs @@ -18,7 +18,7 @@ module Data.Packed.Internal.Matrix( Matrix(..), rows, cols, MatrixOrder(..), orderOf, - createMatrix, withMatrix, mat, + createMatrix, mat, cmat, fmat, toLists, flatten, reshape, Element(..), @@ -115,11 +115,11 @@ fmat m@MF{} = m fmat MC {irows = r, icols = c, cdat = d } = MF {irows = r, icols = c, fdat = transdata c d r} -- C-Haskell matrix adapter -mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r -mat = withMatrix +-- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r -withMatrix a f = - withForeignPtr (fptr (xdat a)) $ \p -> do +mat :: (Storable t) => Matrix t -> (((CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b +mat a f = + unsafeWith (xdat a) $ \p -> do let m g = do g (fi (rows a)) (fi (cols a)) p f m @@ -273,8 +273,8 @@ transdata' c1 v c2 = then v else unsafePerformIO $ do w <- createVector (r2*c2) - withForeignPtr (fptr v) $ \p -> - withForeignPtr (fptr w) $ \q -> do + unsafeWith v $ \p -> + unsafeWith w $ \q -> do let go (-1) _ = return () go !i (-1) = go (i-1) (c1-1) go !i !j = do x <- peekElemOff p (i*c1+j) @@ -300,8 +300,8 @@ transdataAux fun c1 d c2 = then d else unsafePerformIO $ do v <- createVector (dim d) - withForeignPtr (fptr d) $ \pd -> - withForeignPtr (fptr v) $ \pv -> + unsafeWith d $ \pd -> + unsafeWith v $ \pv -> fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux" return v where r1 = dim d `div` c1 @@ -314,7 +314,7 @@ foreign import ccall "transC" ctransC :: TCMCM constant' v n = unsafePerformIO $ do w <- createVector n - withForeignPtr (fptr w) $ \p -> do + unsafeWith w $ \p -> do let go (-1) = return () go !k = pokeElemOff p k v >> go (k-1) go (n-1) @@ -352,8 +352,8 @@ subMatrix (r0,c0) (rt,ct) m subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do w <- createVector (rt*ct) - withForeignPtr (fptr v) $ \p -> - withForeignPtr (fptr w) $ \q -> do + unsafeWith v $ \p -> + unsafeWith w $ \q -> do let go (-1) _ = return () go !i (-1) = go (i-1) (ct-1) go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0) diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs index 51c0d92..6bbd207 100644 --- a/lib/Data/Packed/Internal/Vector.hs +++ b/lib/Data/Packed/Internal/Vector.hs @@ -12,20 +12,20 @@ -- Vector implementation -- ----------------------------------------------------------------------------- --- #hide module Data.Packed.Internal.Vector ( - Vector(..), dim, + Vector, dim, fromList, toList, (|>), join, (@>), safe, at, at', subVector, mapVector, zipVector, foldVector, foldVectorG, foldLoop, - createVector, withVector, vec, + createVector, vec, asComplex, asReal, fwriteVector, freadVector, fprintfVector, fscanfVector, cloneVector, unsafeToForeignPtr, - unsafeFromForeignPtr + unsafeFromForeignPtr, + unsafeWith ) where import Data.Packed.Internal.Common @@ -47,39 +47,45 @@ import GHC.Base import GHC.IOBase #endif --- | A one-dimensional array of objects stored in a contiguous memory block. +-- | One-dimensional array of objects stored in a contiguous memory block. data Vector t = - V { idim :: {-# UNPACK #-} !Int -- ^ number of elements - , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block + V { ioff :: {-# UNPACK #-} !Int -- ^ offset of first element + , idim :: {-# UNPACK #-} !Int -- ^ number of elements + , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block } unsafeToForeignPtr :: Vector a -> (ForeignPtr a, Int, Int) -unsafeToForeignPtr v = (fptr v, 0, idim v) +unsafeToForeignPtr v = (fptr v, ioff v, idim v) -- | Same convention as in Roman Leshchinskiy's vector package. unsafeFromForeignPtr :: ForeignPtr a -> Int -> Int -> Vector a -unsafeFromForeignPtr fp i n | i == 0 = V {idim = n, fptr = fp} - | otherwise = error "unsafeFromForeignPtr with nonzero offset" +unsafeFromForeignPtr fp i n | n > 0 = V {ioff = i, idim = n, fptr = fp} + | otherwise = error "unsafeFromForeignPtr with dim < 1" + +unsafeWith (V i _ fp) m = withForeignPtr fp $ \p -> m (p `advancePtr` i) +{-# INLINE unsafeWith #-} + -- | Number of elements dim :: Vector t -> Int dim = idim -- C-Haskell vector adapter -vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r -vec = withVector - -withVector (V n fp) f = withForeignPtr fp $ \p -> do +-- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r +vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b +vec x f = unsafeWith x $ \p -> do let v g = do - g (fi n) p + g (fi $ dim x) p f v +{-# INLINE vec #-} + -- allocates memory for a new vector createVector :: Storable a => Int -> IO (Vector a) createVector n = do when (n <= 0) $ error ("trying to createVector of dim "++show n) fp <- doMalloc undefined - return $ V n fp + return $ V 0 n fp where -- -- Use the much cheaper Haskell heap allocated storage @@ -102,11 +108,10 @@ createVector n = do fromList :: Storable a => [a] -> Vector a fromList l = unsafePerformIO $ do v <- createVector (length l) - let f _ p = pokeArray p l >> return 0 - app1 f vec v "fromList" + unsafeWith v $ \ p -> pokeArray p l return v -safeRead v = inlinePerformIO . withForeignPtr (fptr v) +safeRead v = inlinePerformIO . unsafeWith v {-# INLINE safeRead #-} inlinePerformIO :: IO a -> a @@ -171,7 +176,12 @@ subVector :: Storable t => Int -- ^ index of the starting element -> Int -- ^ number of elements to extract -> Vector t -- ^ source -> Vector t -- ^ result -subVector k l (v@V {idim=n}) + +subVector k l v@V{idim = n, ioff = i} + | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range" + | otherwise = v {idim = l, ioff = i+k} + +subVectorCopy k l (v@V {idim=n}) | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range" | otherwise = unsafePerformIO $ do r <- createVector l @@ -201,23 +211,24 @@ join [] = error "joining zero vectors" join [v] = v join as = unsafePerformIO $ do let tot = sum (map dim as) - r@V {fptr = p} <- createVector tot - withForeignPtr p $ \ptr -> + r <- createVector tot + unsafeWith r $ \ptr -> joiner as tot ptr return r where joiner [] _ _ = return () - joiner (V {idim = n, fptr = b} : cs) _ p = do - withForeignPtr b $ \pb -> copyArray p pb n + joiner (v:cs) _ p = do + let n = dim v + unsafeWith v $ \pb -> copyArray p pb n joiner cs 0 (advancePtr p n) -- | transforms a complex vector into a real vector with alternating real and imaginary parts asReal :: Vector (Complex Double) -> Vector Double -asReal v = V { idim = 2*dim v, fptr = castForeignPtr (fptr v) } +asReal v = V { ioff = 2*ioff v, idim = 2*dim v, fptr = castForeignPtr (fptr v) } -- | transforms a real vector into a complex vector with alternating real and imaginary parts asComplex :: Vector Double -> Vector (Complex Double) -asComplex v = V { idim = dim v `div` 2, fptr = castForeignPtr (fptr v) } +asComplex v = V { ioff = ioff v `div` 2, idim = dim v `div` 2, fptr = castForeignPtr (fptr v) } ---------------------------------------------------------------- @@ -234,8 +245,8 @@ cloneVector (v@V {idim=n}) = do mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b mapVector f v = unsafePerformIO $ do w <- createVector (dim v) - withForeignPtr (fptr v) $ \p -> - withForeignPtr (fptr w) $ \q -> do + unsafeWith v $ \p -> + unsafeWith w $ \q -> do let go (-1) = return () go !k = do x <- peekElemOff p k pokeElemOff q k (f x) @@ -249,9 +260,9 @@ zipVector :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> zipVector f u v = unsafePerformIO $ do let n = min (dim u) (dim v) w <- createVector n - withForeignPtr (fptr u) $ \pu -> - withForeignPtr (fptr v) $ \pv -> - withForeignPtr (fptr w) $ \pw -> do + unsafeWith u $ \pu -> + unsafeWith v $ \pv -> + unsafeWith w $ \pw -> do let go (-1) = return () go !k = do x <- peekElemOff pu k y <- peekElemOff pv k @@ -262,7 +273,7 @@ zipVector f u v = unsafePerformIO $ do {-# INLINE zipVector #-} foldVector f x v = unsafePerformIO $ - withForeignPtr (fptr (v::Vector Double)) $ \p -> do + unsafeWith (v::Vector Double) $ \p -> do let go (-1) s = return s go !k !s = do y <- peekElemOff p k go (k-1::Int) (f y s) diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs index a47da05..48e35b4 100644 --- a/lib/Data/Packed/ST.hs +++ b/lib/Data/Packed/ST.hs @@ -37,11 +37,11 @@ import Foreign {-# INLINE ioReadV #-} ioReadV :: Storable t => Vector t -> Int -> IO t -ioReadV v k = withForeignPtr (fptr v) $ \s -> peekElemOff s k +ioReadV v k = unsafeWith v $ \s -> peekElemOff s k {-# INLINE ioWriteV #-} ioWriteV :: Storable t => Vector t -> Int -> t -> IO () -ioWriteV v k x = withForeignPtr (fptr v) $ \s -> pokeElemOff s k x +ioWriteV v k x = unsafeWith v $ \s -> pokeElemOff s k x newtype STVector s t = STVector (Vector t) diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs index 37bcc1b..2ccc3ef 100644 --- a/lib/Numeric/GSL/Internal.hs +++ b/lib/Numeric/GSL/Internal.hs @@ -35,12 +35,12 @@ foreign import ccall "wrapper" aux_vTov :: (Vector Double -> Vector Double) -> TVV aux_vTov f n p nr r = g where - V {fptr = pr} = f x + v = f x x = createV (fromIntegral n) copy "aux_vTov" copy n' q = do copyArray q p (fromIntegral n') return 0 - g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral nr) + g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr) return 0 foreign import ccall "wrapper" @@ -51,12 +51,12 @@ foreign import ccall "wrapper" aux_vTom :: (Vector Double -> Matrix Double) -> TVM aux_vTom f n p rr cr r = g where - V {fptr = pr} = flatten $ f x + v = flatten $ f x x = createV (fromIntegral n) copy "aux_vTov" copy n' q = do copyArray q p (fromIntegral n') return 0 - g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral $ rr*cr) + g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral $ rr*cr) return 0 createV n fun msg = unsafePerformIO $ do diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs index 3240ddf..3bbb2b3 100644 --- a/lib/Numeric/GSL/Minimization.hs +++ b/lib/Numeric/GSL/Minimization.hs @@ -109,7 +109,7 @@ ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 minimizeV method eps maxit szv f xiv = unsafePerformIO $ do let n = dim xiv fp <- mkVecfun (iv f) - rawpath <- ww2 withVector xiv withVector szv $ \xiv' szv' -> + rawpath <- ww2 vec xiv vec szv $ \xiv' szv' -> createMIO maxit (n+3) (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv') "minimize" @@ -166,7 +166,7 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do df' = (checkdim1 n . df) fp <- mkVecfun (iv f') dfp <- mkVecVecfun (aux_vTov df') - rawpath <- withVector xiv $ \xiv' -> + rawpath <- vec xiv $ \xiv' -> createMIO maxit (n+2) (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv') "minimizeD" diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs index 9d890fc..f6f11f9 100644 --- a/lib/Numeric/GSL/ODE.hs +++ b/lib/Numeric/GSL/ODE.hs @@ -83,8 +83,8 @@ odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do jp <- case mbjac of Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) Nothing -> return nullFunPtr - sol <- withVector xiv $ \xiv' -> - withVector (checkTimes ts) $ \ts' -> + 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" diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs index 840e8ee..ddb090d 100644 --- a/lib/Numeric/GSL/Root.hs +++ b/lib/Numeric/GSL/Root.hs @@ -79,7 +79,7 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do let xiv = fromList xi n = dim xiv fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) - rawpath <- withVector xiv $ \xiv' -> + rawpath <- vec xiv $ \xiv' -> createMIO maxit (2*n+1) (c_root m fp epsabs (fi maxit) // xiv') "root" @@ -117,7 +117,7 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do n = dim xiv fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList)) - rawpath <- withVector xiv $ \xiv' -> + rawpath <- vec xiv $ \xiv' -> createMIO maxit (2*n+1) (c_rootj m fp jp epsabs (fi maxit) // xiv') "root" diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs index 67496f2..b17da5a 100644 --- a/lib/Numeric/LinearAlgebra/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Instances.hs @@ -25,7 +25,7 @@ import Data.Packed.Matrix import Complex import Data.List(transpose,intersperse) import Foreign(Storable) -import Data.Monoid +-- import Data.Monoid import Data.Packed.Internal.Vector -- import Control.Parallel.Strategies @@ -182,12 +182,12 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati --------------------------------------------------------------- -instance (Storable a, Num (Vector a)) => Monoid (Vector a) where - mempty = 0 { idim = 0 } - mappend a b = mconcat [a,b] - mconcat = j . filter ((>0).dim) - where j [] = mempty - j l = join l +-- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where +-- mempty = 0 { idim = 0 } +-- mappend a b = mconcat [a,b] +-- mconcat = j . filter ((>0).dim) +-- where j [] = mempty +-- j l = join l --------------------------------------------------------------- diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs index 016b9a1..3f7c847 100644 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ b/lib/Numeric/LinearAlgebra/Tests.hs @@ -180,6 +180,9 @@ runTests n = do test (multProp1 . cConsist) test (multProp2 . rConsist) test (multProp2 . cConsist) + putStrLn "------ sub-trans" + test (subProp . rM) + test (subProp . cM) putStrLn "------ lu" test (luProp . rM) test (luProp . cM) @@ -305,6 +308,7 @@ makeUnitary v | realPart n > 1 = v / scalar n -- | Performance measurements. runBenchmarks :: IO () runBenchmarks = do + subBench multBench svdBench eigBench @@ -335,6 +339,16 @@ multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double)) -------------------------------- +subBench = do + putStrLn "" + let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v)) + time "0.1M subVector " (g (constant 1 (1+10^5) :: Vector Double) @> 0) + let f = foldl1' (.) (replicate (10^5) (fromRows.toRows)) + time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0)) + time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0)) + +-------------------------------- + multBench = do let a = ident 1000 :: Matrix Double let b = ident 2000 :: Matrix Double diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs index 93175f2..6d176fa 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs @@ -38,6 +38,7 @@ module Numeric.LinearAlgebra.Tests.Properties ( cholProp, expmDiagProp, multProp1, multProp2, + subProp, linearSolveProp, linearSolveProp2 ) where @@ -243,3 +244,6 @@ linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b) where q = min (rows a) (cols a) b = a <> x wc = rank a == q + +subProp m = m == (trans . fromColumns . toRows) m + -- cgit v1.2.3