flag -fvector

3 files changed, 67 insertions, 9 deletions

diff --git a/hmatrix.cabal b/hmatrix.cabal
index e684acc..cbf4a55 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -65,6 +65,10 @@ flag unsafe
     description:    Compile the library with bound checking disabled.
     default:        False
 
+flag vector
+    description:    Use Data.Vector.Storable type from "vector" package.
+    default:        False
+
 library
 
     Build-Depends:      base >= 4 && < 5,
@@ -110,6 +114,10 @@ library
     C-sources:          lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,
                         lib/Numeric/GSL/gsl-aux.c
 
+    if flag(vector)
+       Build-Depends:   vector
+       cpp-options:     -DVECTOR
+
     if flag(tests)
        Build-Depends:   QuickCheck, HUnit
        exposed-modules: Numeric.LinearAlgebra.Tests
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 9ca8e58..a6868d9 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -47,6 +47,22 @@ import GHC.Base
 import GHC.IOBase
 #endif
 
+#ifdef VECTOR
+import qualified Data.Vector.Storable as Vector
+import Data.Vector.Storable(Vector,
+                            unsafeToForeignPtr,
+                            unsafeFromForeignPtr,
+                            unsafeWith)
+#endif
+
+#ifdef VECTOR
+
+-- | Number of elements
+dim :: (Storable t) => Vector t -> Int
+dim = Vector.length
+
+#else
+
 -- | One-dimensional array of objects stored in a contiguous memory block.
 data Vector t =
     V { ioff :: {-# UNPACK #-} !Int              -- ^ offset of first element
@@ -65,11 +81,13 @@ unsafeFromForeignPtr fp i n | n > 0 = V {ioff = i, idim = n, fptr = fp}
 unsafeWith (V i _ fp) m = withForeignPtr fp $ \p -> m (p `advancePtr` i)
 {-# INLINE unsafeWith #-}
 
-
 -- | Number of elements
-dim :: Vector t -> Int
+dim :: (Storable t) => Vector t -> Int
 dim = idim
 
+#endif
+
+
 -- C-Haskell vector adapter
 -- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
 vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
@@ -85,7 +103,7 @@ 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 0 n fp
+    return $ unsafeFromForeignPtr fp 0 n
   where
     --
     -- Use the much cheaper Haskell heap allocated storage
@@ -177,6 +195,12 @@ subVector :: Storable t => Int       -- ^ index of the starting element
                         -> Vector t  -- ^ source
                         -> Vector t  -- ^ result
 
+#ifdef VECTOR
+
+subVector = Vector.slice
+
+#else
+
 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}
@@ -189,6 +213,8 @@ subVectorCopy k l (v@V {idim=n})
         app2 f vec v vec r "subVector"
         return r
 
+#endif
+
 {- | Reads a vector position:
 
 @> fromList [0..9] \@\> 7
@@ -239,16 +265,21 @@ takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (
 
 -- | transforms a complex vector into a real vector with alternating real and imaginary parts 
 asReal :: Vector (Complex Double) -> Vector Double
-asReal v = V { ioff = 2*ioff v, idim = 2*dim v, fptr =  castForeignPtr (fptr v) }
+--asReal v = V { ioff = 2*ioff v, idim = 2*dim v, fptr =  castForeignPtr (fptr v) }
+asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)
+    where (fp,i,n) = unsafeToForeignPtr v
 
 -- | transforms a real vector into a complex vector with alternating real and imaginary parts
 asComplex :: Vector Double -> Vector (Complex Double)
-asComplex v = V { ioff = ioff v `div` 2, 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) }
+asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)
+    where (fp,i,n) = unsafeToForeignPtr v
 
 ----------------------------------------------------------------
 
 cloneVector :: Storable t => Vector t -> IO (Vector t)
-cloneVector (v@V {idim=n}) = do
+cloneVector v = do
+        let n = dim v
         r <- createVector n
         let f _ s _ d =  copyArray d s n >> return 0
         app2 f vec v vec r "cloneVector"
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
index b4a769e..1992db0 100644
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -24,10 +24,11 @@ import Numeric.GSL.Vector
 import Data.Packed.Matrix
 import Data.Complex
 import Data.List(transpose,intersperse)
-import Foreign(Storable)
--- import Data.Monoid
 import Data.Packed.Internal.Vector
--- import Control.Parallel.Strategies
+
+#ifndef VECTOR
+import Foreign(Storable)
+#endif
 
 ------------------------------------------------------------------
 
@@ -43,9 +44,13 @@ dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unw
         pad n str = replicate (n - length str) ' ' ++ str
         unwords' = concat . intersperse ", "
 
+#ifndef VECTOR
+
 instance (Show a, Storable a) => (Show (Vector a)) where
     show v = (show (dim v))++" |> " ++ show (toList v)
 
+#endif
+
 ------------------------------------------------------------------
 
 instance (Element a, Read a) => Read (Matrix a) where
@@ -55,12 +60,23 @@ instance (Element a, Read a) => Read (Matrix a) where
               cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
               rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
 
+#ifdef VECTOR
+
+instance (Element a, Read a) => Read (Vector a) where
+    readsPrec _ s = [(fromList . read $ listnums, rest)]
+        where (thing,trest) = breakAt ']' s
+              (dims,listnums) = breakAt ' ' (dropWhile (==' ') thing)
+              rest = drop 31 trest
+#else
+
 instance (Element a, Read a) => Read (Vector a) where
     readsPrec _ s = [((d |>) . read $ listnums, rest)]
         where (thing,rest) = breakAt ']' s
               (dims,listnums) = breakAt '>' thing
               d = read . init . fst . breakAt '|' $ dims
 
+#endif
+
 breakAt c l = (a++[c],tail b) where
     (a,b) = break (==c) l
 
@@ -71,10 +87,13 @@ adaptScalar f1 f2 f3 x y
     | dim y == 1 = f3 x (y@>0)
     | otherwise = f2 x y
 
+#ifndef VECTOR
 
 instance Linear Vector a => Eq (Vector a) where
     (==) = equal
 
+#endif
+
 instance Num (Vector Double) where
     (+) = adaptScalar addConstant add (flip addConstant)
     negate = scale (-1)