compatible with vector

14 files changed, 110 insertions, 168 deletions

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