11 files changed, 2235 insertions, 0 deletions

diff --git a/packages/hmatrix/src/Data/Packed/Development.hs b/packages/hmatrix/src/Data/Packed/Development.hs
new file mode 100644
index 0000000..471e560
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Development.hs
@@ -0,0 +1,32 @@
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Development
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- The library can be easily extended with additional foreign functions
+-- using the tools in this module. Illustrative usage examples can be found
+-- in the @examples\/devel@ folder included in the package.
+--
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+module Data.Packed.Development (
+    createVector, createMatrix,
+    vec, mat,
+    app1, app2, app3, app4,
+    app5, app6, app7, app8, app9, app10,
+    MatrixOrder(..), orderOf, cmat, fmat,
+    matrixFromVector,
+    unsafeFromForeignPtr,
+    unsafeToForeignPtr,
+    check, (//),
+    at', atM'
+) where
+
+import Data.Packed.Internal
diff --git a/packages/hmatrix/src/Data/Packed/Foreign.hs b/packages/hmatrix/src/Data/Packed/Foreign.hs
new file mode 100644
index 0000000..1ec3694
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Foreign.hs
@@ -0,0 +1,100 @@
+{-# LANGUAGE MagicHash, UnboxedTuples #-}
+-- | FFI and hmatrix helpers.
+--
+-- Sample usage, to upload a perspective matrix to a shader.
+--
+-- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) 
+-- @
+--
+{-# OPTIONS_HADDOCK hide #-}
+module Data.Packed.Foreign 
+    ( app
+    , appVector, appVectorLen
+    , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen
+    , unsafeMatrixToVector, unsafeMatrixToForeignPtr
+    ) where
+import Data.Packed.Internal
+import qualified Data.Vector.Storable as S
+import Foreign (Ptr, ForeignPtr, Storable)
+import Foreign.C.Types (CInt)
+import GHC.Base (IO(..), realWorld#)
+
+{-# INLINE unsafeInlinePerformIO #-}
+-- | If we use unsafePerformIO, it may not get inlined, so in a function that returns IO (which are all safe uses of app* in this module), there would be
+-- unecessary calls to unsafePerformIO or its internals.
+unsafeInlinePerformIO :: IO a -> a
+unsafeInlinePerformIO (IO f) = case f realWorld# of
+    (# _, x #) -> x
+
+{-# INLINE app #-}
+-- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative.
+-- e.g.
+--
+-- @
+-- someFunction
+--     \`appMatrixLen\` m
+--     \`appVectorLen\` v
+--     \`app\` other
+--     \`app\` arguments
+--     \`app\` go here
+-- @
+--
+-- One could also write:
+--
+-- @
+-- (someFunction 
+--     \`appMatrixLen\` m
+--     \`appVectorLen\` v) 
+--     other 
+--     arguments 
+--     (go here)
+-- @
+--
+app :: (a -> b) -> a -> b
+app f = f
+
+{-# INLINE appVector #-}
+appVector :: Storable a => (Ptr a -> b) -> Vector a -> b
+appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))
+
+{-# INLINE appVectorLen #-}
+appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
+appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x))))
+
+{-# INLINE appMatrix #-}
+appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b
+appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f))
+
+{-# INLINE appMatrixLen #-}
+appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
+appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c))
+  where
+    r = fromIntegral (rows x)
+    c = fromIntegral (cols x)
+
+{-# INLINE appMatrixRaw #-}
+appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b
+appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f))
+
+{-# INLINE appMatrixRawLen #-}
+appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
+appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c))
+  where
+    r = fromIntegral (rows x)
+    c = fromIntegral (cols x)
+
+infixl 1 `app`
+infixl 1 `appVector`
+infixl 1 `appMatrix`
+infixl 1 `appMatrixRaw`
+
+{-# INLINE unsafeMatrixToVector #-}
+-- | This will disregard the order of the matrix, and simply return it as-is. 
+-- If the order of the matrix is RowMajor, this function is identical to 'flatten'.
+unsafeMatrixToVector :: Matrix a -> Vector a
+unsafeMatrixToVector = xdat
+
+{-# INLINE unsafeMatrixToForeignPtr #-}
+unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int)
+unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m)
+
diff --git a/packages/hmatrix/src/Data/Packed/Internal.hs b/packages/hmatrix/src/Data/Packed/Internal.hs
new file mode 100644
index 0000000..537e51e
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Internal.hs
@@ -0,0 +1,26 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Internal
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Reexports all internal modules
+--
+-----------------------------------------------------------------------------
+-- #hide
+
+module Data.Packed.Internal (
+    module Data.Packed.Internal.Common,
+    module Data.Packed.Internal.Signatures,
+    module Data.Packed.Internal.Vector,
+    module Data.Packed.Internal.Matrix,
+) where
+
+import Data.Packed.Internal.Common
+import Data.Packed.Internal.Signatures
+import Data.Packed.Internal.Vector
+import Data.Packed.Internal.Matrix
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Common.hs b/packages/hmatrix/src/Data/Packed/Internal/Common.hs
new file mode 100644
index 0000000..edef3c2
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Internal/Common.hs
@@ -0,0 +1,171 @@
+{-# LANGUAGE CPP #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Internal.Common
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable (uses FFI)
+--
+-- Development utilities.
+--
+-----------------------------------------------------------------------------
+-- #hide
+
+module Data.Packed.Internal.Common(
+  Adapt,
+  app1, app2, app3, app4,
+  app5, app6, app7, app8, app9, app10,
+  (//), check, mbCatch,
+  splitEvery, common, compatdim,
+  fi,
+  table
+) where
+
+import Foreign
+import Control.Monad(when)
+import Foreign.C.String(peekCString)
+import Foreign.C.Types
+import Foreign.Storable.Complex()
+import Data.List(transpose,intersperse)
+import Control.Exception as E
+
+-- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@
+splitEvery :: Int -> [a] -> [[a]]
+splitEvery _ [] = []
+splitEvery k l = take k l : splitEvery k (drop k l)
+
+-- | obtains the common value of a property of a list
+common :: (Eq a) => (b->a) -> [b] -> Maybe a
+common f = commonval . map f where
+    commonval :: (Eq a) => [a] -> Maybe a
+    commonval [] = Nothing
+    commonval [a] = Just a
+    commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
+
+-- | common value with \"adaptable\" 1
+compatdim :: [Int] -> Maybe Int
+compatdim [] = Nothing
+compatdim [a] = Just a
+compatdim (a:b:xs)
+    | a==b = compatdim (b:xs)
+    | a==1 = compatdim (b:xs)
+    | b==1 = compatdim (a:xs)
+    | otherwise = Nothing
+
+-- | Formatting tool
+table :: String -> [[String]] -> String
+table sep as = unlines . map unwords' $ transpose mtp where 
+    mt = transpose as
+    longs = map (maximum . map length) mt
+    mtp = zipWith (\a b -> map (pad a) b) longs mt
+    pad n str = replicate (n - length str) ' ' ++ str
+    unwords' = concat . intersperse sep
+
+-- | postfix function application (@flip ($)@)
+(//) :: x -> (x -> y) -> y
+infixl 0 //
+(//) = flip ($)
+
+-- | specialized fromIntegral
+fi :: Int -> CInt
+fi = fromIntegral
+
+-- hmm..
+ww2 w1 o1 w2 o2 f = w1 o1 $ w2 o2 . f
+ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ ww2 w2 o2 w3 o3 . f
+ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ ww3 w2 o2 w3 o3 w4 o4 . f
+ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 f = w1 o1 $ ww4 w2 o2 w3 o3 w4 o4 w5 o5 . f
+ww6 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 f = w1 o1 $ ww5 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 . f
+ww7 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 f = w1 o1 $ ww6 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 . f
+ww8 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 f = w1 o1 $ ww7 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 . f
+ww9 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 f = w1 o1 $ ww8 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 . f
+ww10 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 f = w1 o1 $ ww9 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 . f
+
+type Adapt f t r = t -> ((f -> r) -> IO()) -> IO()
+
+type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO()
+type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2
+type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3
+type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4
+type Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5 = Adapt f t1 r1 -> t1 -> Adapt4 r1 t2 r2 t3 r3 t4 r4 t5
+type Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 = Adapt f t1 r1 -> t1 -> Adapt5 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6
+type Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 = Adapt f t1 r1 -> t1 -> Adapt6 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7
+type Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 = Adapt f t1 r1 -> t1 -> Adapt7 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8
+type Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 = Adapt f t1 r1 -> t1 -> Adapt8 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9
+type Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 = Adapt f t1 r1 -> t1 -> Adapt9 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10
+
+app1 :: f -> Adapt1 f t1
+app2 :: f -> Adapt2 f t1 r1 t2
+app3 :: f -> Adapt3 f t1 r1 t2 r2 t3
+app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4
+app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5
+app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6
+app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7
+app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8
+app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9
+app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10
+
+app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s
+app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s
+app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $
+     \a1 a2 a3 -> f // a1 // a2 // a3 // check s
+app4 f w1 o1 w2 o2 w3 o3 w4 o4 s = ww4 w1 o1 w2 o2 w3 o3 w4 o4 $
+     \a1 a2 a3 a4 -> f // a1 // a2 // a3 // a4 // check s
+app5 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 s = ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 $
+     \a1 a2 a3 a4 a5 -> f // a1 // a2 // a3 // a4 // a5 // check s
+app6 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 s = ww6 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 $
+     \a1 a2 a3 a4 a5 a6 -> f // a1 // a2 // a3 // a4 // a5 // a6 // check s
+app7 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 s = ww7 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 $
+     \a1 a2 a3 a4 a5 a6 a7 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // check s
+app8 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 s = ww8 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 $
+     \a1 a2 a3 a4 a5 a6 a7 a8 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // check s
+app9 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 s = ww9 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 $
+     \a1 a2 a3 a4 a5 a6 a7 a8 a9 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // check s
+app10 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 s = ww10 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 $
+     \a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // a10 // check s
+
+
+
+-- GSL error codes are <= 1024
+-- | error codes for the auxiliary functions required by the wrappers
+errorCode :: CInt -> String
+errorCode 2000 = "bad size"
+errorCode 2001 = "bad function code"
+errorCode 2002 = "memory problem"
+errorCode 2003 = "bad file"
+errorCode 2004 = "singular"
+errorCode 2005 = "didn't converge"
+errorCode 2006 = "the input matrix is not positive definite"
+errorCode 2007 = "not yet supported in this OS"
+errorCode n    = "code "++show n
+
+
+-- | clear the fpu
+foreign import ccall unsafe "asm_finit" finit :: IO ()
+
+-- | check the error code
+check :: String -> IO CInt -> IO ()
+check msg f = do
+#if FINIT
+    finit
+#endif
+    err <- f
+    when (err/=0) $ if err > 1024
+                      then (error (msg++": "++errorCode err)) -- our errors
+                      else do                                 -- GSL errors
+                        ps <- gsl_strerror err
+                        s <- peekCString ps
+                        error (msg++": "++s)
+    return ()
+
+-- | description of GSL error codes
+foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
+
+-- | Error capture and conversion to Maybe
+mbCatch :: IO x -> IO (Maybe x)
+mbCatch act = E.catch (Just `fmap` act) f
+    where f :: SomeException -> IO (Maybe x)
+          f _ = return Nothing
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs b/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs
new file mode 100644
index 0000000..9719fc0
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs
@@ -0,0 +1,491 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+{-# LANGUAGE FlexibleContexts         #-}
+{-# LANGUAGE FlexibleInstances        #-}
+{-# LANGUAGE BangPatterns             #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Internal.Matrix
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable (uses FFI)
+--
+-- Internal matrix representation
+--
+-----------------------------------------------------------------------------
+-- #hide
+
+module Data.Packed.Internal.Matrix(
+    Matrix(..), rows, cols, cdat, fdat,
+    MatrixOrder(..), orderOf,
+    createMatrix, mat,
+    cmat, fmat,
+    toLists, flatten, reshape,
+    Element(..),
+    trans,
+    fromRows, toRows, fromColumns, toColumns,
+    matrixFromVector,
+    subMatrix,
+    liftMatrix, liftMatrix2,
+    (@@>), atM',
+    saveMatrix,
+    singleton,
+    emptyM,
+    size, shSize, conformVs, conformMs, conformVTo, conformMTo
+) where
+
+import Data.Packed.Internal.Common
+import Data.Packed.Internal.Signatures
+import Data.Packed.Internal.Vector
+
+import Foreign.Marshal.Alloc(alloca, free)
+import Foreign.Marshal.Array(newArray)
+import Foreign.Ptr(Ptr, castPtr)
+import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf)
+import Data.Complex(Complex)
+import Foreign.C.Types
+import Foreign.C.String(newCString)
+import System.IO.Unsafe(unsafePerformIO)
+import Control.DeepSeq
+
+-----------------------------------------------------------------
+
+{- Design considerations for the Matrix Type
+   -----------------------------------------
+
+- we must easily handle both row major and column major order,
+  for bindings to LAPACK and GSL/C
+
+- we'd like to simplify redundant matrix transposes:
+   - Some of them arise from the order requirements of some functions
+   - some functions (matrix product) admit transposed arguments
+
+- maybe we don't really need this kind of simplification:
+   - more complex code
+   - some computational overhead
+   - only appreciable gain in code with a lot of redundant transpositions
+     and cheap matrix computations
+
+- we could carry both the matrix and its (lazily computed) transpose.
+  This may save some transpositions, but it is necessary to keep track of the
+  data which is actually computed to be used by functions like the matrix product
+  which admit both orders.
+
+- but if we need the transposed data and it is not in the structure, we must make
+  sure that we touch the same foreignptr that is used in the computation.
+
+- a reasonable solution is using two constructors for a matrix. Transposition just
+  "flips" the constructor. Actual data transposition is not done if followed by a
+  matrix product or another transpose.
+
+-}
+
+data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
+
+transOrder RowMajor = ColumnMajor
+transOrder ColumnMajor = RowMajor
+{- | Matrix representation suitable for GSL and LAPACK computations.
+
+The elements are stored in a continuous memory array.
+
+-}
+
+data Matrix t = Matrix { irows :: {-# UNPACK #-} !Int
+                       , icols :: {-# UNPACK #-} !Int
+                       , xdat :: {-# UNPACK #-} !(Vector t)
+                       , order :: !MatrixOrder }
+-- RowMajor: preferred by C, fdat may require a transposition
+-- ColumnMajor: preferred by LAPACK, cdat may require a transposition
+
+cdat = xdat
+fdat = xdat
+
+rows :: Matrix t -> Int
+rows = irows
+
+cols :: Matrix t -> Int
+cols = icols
+
+orderOf :: Matrix t -> MatrixOrder
+orderOf = order
+
+
+-- | Matrix transpose.
+trans :: Matrix t -> Matrix t
+trans Matrix {irows = r, icols = c, xdat = d, order = o } = Matrix { irows = c, icols = r, xdat = d, order = transOrder o}
+
+cmat :: (Element t) => Matrix t -> Matrix t
+cmat m@Matrix{order = RowMajor} = m
+cmat Matrix {irows = r, icols = c, xdat = d, order = ColumnMajor } = Matrix { irows = r, icols = c, xdat = transdata r d c, order = RowMajor}
+
+fmat :: (Element t) => Matrix t -> Matrix t
+fmat m@Matrix{order = ColumnMajor} = m
+fmat Matrix {irows = r, icols = c, xdat = d, order = RowMajor } = Matrix { irows = r, icols = c, xdat = transdata c d r, order = ColumnMajor}
+
+-- C-Haskell matrix adapter
+-- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r
+
+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
+
+{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
+
+>>> flatten (ident 3)
+fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
+
+-}
+flatten :: Element t => Matrix t -> Vector t
+flatten = xdat . cmat
+
+{-
+type Mt t s = Int -> Int -> Ptr t -> s
+
+infixr 6 ::>
+type t ::> s = Mt t s
+-}
+
+-- | the inverse of 'Data.Packed.Matrix.fromLists'
+toLists :: (Element t) => Matrix t -> [[t]]
+toLists m = splitEvery (cols m) . toList . flatten $ m
+
+-- | Create a matrix from a list of vectors.
+-- All vectors must have the same dimension,
+-- or dimension 1, which is are automatically expanded.
+fromRows :: Element t => [Vector t] -> Matrix t
+fromRows [] = emptyM 0 0
+fromRows vs = case compatdim (map dim vs) of
+    Nothing -> error $ "fromRows expects vectors with equal sizes (or singletons), given: " ++ show (map dim vs)
+    Just 0  -> emptyM r 0
+    Just c  -> matrixFromVector RowMajor r c . vjoin . map (adapt c) $ vs
+  where
+    r = length vs
+    adapt c v
+        | c == 0 = fromList[]
+        | dim v == c = v
+        | otherwise = constantD (v@>0) c
+
+-- | extracts the rows of a matrix as a list of vectors
+toRows :: Element t => Matrix t -> [Vector t]
+toRows m
+    | c == 0    = replicate r (fromList[])
+    | otherwise = toRows' 0
+  where
+    v = flatten m
+    r = rows m
+    c = cols m
+    toRows' k | k == r*c  = []
+              | otherwise = subVector k c v : toRows' (k+c)
+
+-- | Creates a matrix from a list of vectors, as columns
+fromColumns :: Element t => [Vector t] -> Matrix t
+fromColumns m = trans . fromRows $ m
+
+-- | Creates a list of vectors from the columns of a matrix
+toColumns :: Element t => Matrix t -> [Vector t]
+toColumns m = toRows . trans $ m
+
+-- | Reads a matrix position.
+(@@>) :: Storable t => Matrix t -> (Int,Int) -> t
+infixl 9 @@>
+m@Matrix {irows = r, icols = c} @@> (i,j)
+    | safe      = if i<0 || i>=r || j<0 || j>=c
+                    then error "matrix indexing out of range"
+                    else atM' m i j
+    | otherwise = atM' m i j
+{-# INLINE (@@>) #-}
+
+--  Unsafe matrix access without range checking
+atM' Matrix {icols = c, xdat = v, order = RowMajor} i j = v `at'` (i*c+j)
+atM' Matrix {irows = r, xdat = v, order = ColumnMajor} i j = v `at'` (j*r+i)
+{-# INLINE atM' #-}
+
+------------------------------------------------------------------
+
+matrixFromVector o r c v
+    | r * c == dim v = m
+    | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m
+  where
+    m = Matrix { irows = r, icols = c, xdat = v, order = o }
+
+-- allocates memory for a new matrix
+createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a)
+createMatrix ord r c = do
+    p <- createVector (r*c)
+    return (matrixFromVector ord r c p)
+
+{- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = trans . reshape r@
+where r is the desired number of rows.)
+
+>>> reshape 4 (fromList [1..12])
+(3><4)
+ [ 1.0,  2.0,  3.0,  4.0
+ , 5.0,  6.0,  7.0,  8.0
+ , 9.0, 10.0, 11.0, 12.0 ]
+
+-}
+reshape :: Storable t => Int -> Vector t -> Matrix t
+reshape 0 v = matrixFromVector RowMajor 0 0 v
+reshape c v = matrixFromVector RowMajor (dim v `div` c) c v
+
+singleton x = reshape 1 (fromList [x])
+
+-- | application of a vector function on the flattened matrix elements
+liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
+liftMatrix f Matrix { irows = r, icols = c, xdat = d, order = o } = matrixFromVector o r c (f d)
+
+-- | application of a vector function on the flattened matrices elements
+liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
+liftMatrix2 f m1 m2
+    | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2"
+    | otherwise = case orderOf m1 of
+        RowMajor    -> matrixFromVector RowMajor    (rows m1) (cols m1) (f (xdat m1) (flatten m2))
+        ColumnMajor -> matrixFromVector ColumnMajor (rows m1) (cols m1) (f (xdat m1) ((xdat.fmat) m2))
+
+
+compat :: Matrix a -> Matrix b -> Bool
+compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
+
+------------------------------------------------------------------
+
+{- | Supported matrix elements.
+
+    This class provides optimized internal
+    operations for selected element types.
+    It provides unoptimised defaults for any 'Storable' type,
+    so you can create instances simply as:
+    @instance Element Foo@.
+-}
+class (Storable a) => Element a where
+    subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position 
+               -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
+               -> Matrix a -> Matrix a
+    subMatrixD = subMatrix'
+    transdata :: Int -> Vector a -> Int -> Vector a
+    transdata = transdataP -- transdata'
+    constantD  :: a -> Int -> Vector a
+    constantD = constantP -- constant'
+
+
+instance Element Float where
+    transdata  = transdataAux ctransF
+    constantD  = constantAux cconstantF
+
+instance Element Double where
+    transdata  = transdataAux ctransR
+    constantD  = constantAux cconstantR
+
+instance Element (Complex Float) where
+    transdata  = transdataAux ctransQ
+    constantD  = constantAux cconstantQ
+
+instance Element (Complex Double) where
+    transdata  = transdataAux ctransC
+    constantD  = constantAux cconstantC
+
+-------------------------------------------------------------------
+
+transdata' :: Storable a => Int -> Vector a -> Int -> Vector a
+transdata' c1 v c2 =
+    if noneed
+        then v
+        else unsafePerformIO $ do
+                w <- createVector (r2*c2)
+                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)
+                                          pokeElemOff      q (j*c2+i) x
+                                          go i (j-1)
+                        go (r1-1) (c1-1)
+                return w
+  where r1 = dim v `div` c1
+        r2 = dim v `div` c2
+        noneed = dim v == 0 || r1 == 1 || c1 == 1
+
+-- {-# SPECIALIZE transdata' :: Int -> Vector Double -> Int ->  Vector Double #-}
+-- {-# SPECIALIZE transdata' :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double) #-}
+
+-- I don't know how to specialize...
+-- The above pragmas only seem to work on top level defs
+-- Fortunately everything seems to work using the above class
+
+-- C versions, still a little faster:
+
+transdataAux fun c1 d c2 =
+    if noneed
+        then d
+        else unsafePerformIO $ do
+            v <- createVector (dim d)
+            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
+        r2 = dim d `div` c2
+        noneed = dim d == 0 || r1 == 1 || c1 == 1
+
+transdataP :: Storable a => Int -> Vector a -> Int -> Vector a
+transdataP c1 d c2 =
+    if noneed
+       then d
+       else unsafePerformIO $ do
+          v <- createVector (dim d)
+          unsafeWith d $ \pd ->
+              unsafeWith v $ \pv ->
+                  ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataP"
+          return v
+   where r1 = dim d `div` c1
+         r2 = dim d `div` c2
+         sz = sizeOf (d @> 0)
+         noneed = dim d == 0 || r1 == 1 || c1 == 1
+
+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
+
+----------------------------------------------------------------------
+
+constant' v n = unsafePerformIO $ do
+    w <- createVector n
+    unsafeWith w $ \p -> do
+        let go (-1) = return ()
+            go !k = pokeElemOff p k v >> go (k-1)
+        go (n-1)
+    return w
+
+-- C versions
+
+constantAux fun x n = unsafePerformIO $ do
+    v <- createVector n
+    px <- newArray [x]
+    app1 (fun px) vec v "constantAux"
+    free px
+    return v
+
+constantF :: Float -> Int -> Vector Float
+constantF = constantAux cconstantF
+foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF
+
+constantR :: Double -> Int -> Vector Double
+constantR = constantAux cconstantR
+foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV
+
+constantQ :: Complex Float -> Int -> Vector (Complex Float)
+constantQ = constantAux cconstantQ
+foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
+
+constantC :: Complex Double -> Int -> Vector (Complex Double)
+constantC = constantAux cconstantC
+foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV
+
+constantP :: Storable a => a -> Int -> Vector a
+constantP a n = unsafePerformIO $ do
+    let sz = sizeOf a
+    v <- createVector n
+    unsafeWith v $ \p -> do
+       alloca $ \k -> do
+                      poke k a
+                      cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP"
+    return v
+foreign import ccall unsafe "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt
+
+----------------------------------------------------------------------
+
+-- | Extracts a submatrix from a matrix.
+subMatrix :: Element a
+          => (Int,Int) -- ^ (r0,c0) starting position 
+          -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
+          -> Matrix a -- ^ input matrix
+          -> Matrix a -- ^ result
+subMatrix (r0,c0) (rt,ct) m
+    | 0 <= r0 && 0 <= rt && r0+rt <= (rows m) &&
+      0 <= c0 && 0 <= ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m
+    | otherwise = error $ "wrong subMatrix "++
+                          show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m)
+
+subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do
+    w <- createVector (rt*ct)
+    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)
+                              pokeElemOff      q (i*ct+j) x
+                              go i (j-1)
+            go (rt-1) (ct-1)
+    return w
+
+subMatrix' (r0,c0) (rt,ct) (Matrix { icols = c, xdat = v, order = RowMajor}) = Matrix rt ct (subMatrix'' (r0,c0) (rt,ct) c v) RowMajor
+subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m)
+
+--------------------------------------------------------------------------
+
+-- | Saves a matrix as 2D ASCII table.
+saveMatrix :: FilePath
+           -> String     -- ^ format (%f, %g, %e)
+           -> Matrix Double
+           -> IO ()
+saveMatrix filename fmt m = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    let o = if orderOf m == RowMajor then 1 else 0
+    app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
+
+----------------------------------------------------------------------
+
+maxZ xs = if minimum xs == 0 then 0 else maximum xs
+
+conformMs ms = map (conformMTo (r,c)) ms
+  where
+    r = maxZ (map rows ms)
+    c = maxZ (map cols ms)
+    
+
+conformVs vs = map (conformVTo n) vs
+  where
+    n = maxZ (map dim vs)
+
+conformMTo (r,c) m
+    | size m == (r,c) = m
+    | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
+    | size m == (r,1) = repCols c m
+    | size m == (1,c) = repRows r m
+    | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")"
+
+conformVTo n v
+    | dim v == n = v
+    | dim v == 1 = constantD (v@>0) n
+    | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n
+
+repRows n x = fromRows (replicate n (flatten x))
+repCols n x = fromColumns (replicate n (flatten x))
+
+size m = (rows m, cols m)
+
+shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"
+
+emptyM r c = matrixFromVector RowMajor r c (fromList[])
+
+----------------------------------------------------------------------
+
+instance (Storable t, NFData t) => NFData (Matrix t)
+  where
+    rnf m | d > 0     = rnf (v @> 0)
+          | otherwise = ()
+      where
+        d = dim v
+        v = xdat m
+
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs b/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs
new file mode 100644
index 0000000..2835720
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs
@@ -0,0 +1,72 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Internal.Signatures
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable (uses FFI)
+--
+-- Signatures of the C functions.
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Internal.Signatures where
+
+import Foreign.Ptr(Ptr)
+import Data.Complex(Complex)
+import Foreign.C.Types(CInt)
+
+type PF = Ptr Float                             --
+type PD = Ptr Double                            --
+type PQ = Ptr (Complex Float)                   --
+type PC = Ptr (Complex Double)                  --
+type TF = CInt -> PF -> IO CInt                 --
+type TFF = CInt -> PF -> TF                     --
+type TFV = CInt -> PF -> TV                     --
+type TVF = CInt -> PD -> TF                     --
+type TFFF = CInt -> PF -> TFF                   --
+type TV = CInt -> PD -> IO CInt                 --
+type TVV = CInt -> PD -> TV                     --
+type TVVV = CInt -> PD -> TVV                   --
+type TFM = CInt -> CInt -> PF -> IO CInt        --
+type TFMFM =  CInt -> CInt -> PF -> TFM         --
+type TFMFMFM =  CInt -> CInt -> PF -> TFMFM     --
+type TM = CInt -> CInt -> PD -> IO CInt         --
+type TMM =  CInt -> CInt -> PD -> TM            --
+type TVMM = CInt -> PD -> TMM                   --
+type TMVMM = CInt -> CInt -> PD -> TVMM         --
+type TMMM =  CInt -> CInt -> PD -> TMM          --
+type TVM = CInt -> PD -> TM                     --
+type TVVM = CInt -> PD -> TVM                   --
+type TMV = CInt -> CInt -> PD -> TV             --
+type TMMV = CInt -> CInt -> PD -> TMV           --
+type TMVM = CInt -> CInt -> PD -> TVM           --
+type TMMVM = CInt -> CInt -> PD -> TMVM         --
+type TCM = CInt -> CInt -> PC -> IO CInt        --
+type TCVCM = CInt -> PC -> TCM                  --
+type TCMCVCM = CInt -> CInt -> PC -> TCVCM      --
+type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM   --
+type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM  --
+type TCMCM = CInt -> CInt -> PC -> TCM          --
+type TVCM = CInt -> PD -> TCM                   --
+type TCMVCM = CInt -> CInt -> PC -> TVCM        --
+type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM    --
+type TCMCMCM = CInt -> CInt -> PC -> TCMCM      --
+type TCV = CInt -> PC -> IO CInt                --
+type TCVCV = CInt -> PC -> TCV                  --
+type TCVCVCV = CInt -> PC -> TCVCV              --
+type TCVV = CInt -> PC -> TV                    --
+type TQV = CInt -> PQ -> IO CInt                --
+type TQVQV = CInt -> PQ -> TQV                  --
+type TQVQVQV = CInt -> PQ -> TQVQV              --
+type TQVF = CInt -> PQ -> TF                    --
+type TQM = CInt -> CInt -> PQ -> IO CInt        --
+type TQMQM = CInt -> CInt -> PQ -> TQM          --
+type TQMQMQM = CInt -> CInt -> PQ -> TQMQM      --
+type TCMCV = CInt -> CInt -> PC -> TCV          --
+type TVCV = CInt -> PD -> TCV                   --
+type TCVM = CInt -> PC -> TM                    --
+type TMCVM = CInt -> CInt -> PD -> TCVM         --
+type TMMCVM = CInt -> CInt -> PD -> TMCVM       --
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Vector.hs b/packages/hmatrix/src/Data/Packed/Internal/Vector.hs
new file mode 100644
index 0000000..6d03438
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Internal/Vector.hs
@@ -0,0 +1,521 @@
+{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Internal.Vector
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable (uses FFI)
+--
+-- Vector implementation
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Internal.Vector (
+    Vector, dim,
+    fromList, toList, (|>),
+    vjoin, (@>), safe, at, at', subVector, takesV,
+    mapVector, mapVectorWithIndex, zipVectorWith, unzipVectorWith,
+    mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
+    foldVector, foldVectorG, foldLoop, foldVectorWithIndex,
+    createVector, vec,
+    asComplex, asReal, float2DoubleV, double2FloatV,
+    stepF, stepD, condF, condD,
+    conjugateQ, conjugateC,
+    fwriteVector, freadVector, fprintfVector, fscanfVector,
+    cloneVector,
+    unsafeToForeignPtr,
+    unsafeFromForeignPtr,
+    unsafeWith
+) where
+
+import Data.Packed.Internal.Common
+import Data.Packed.Internal.Signatures
+import Foreign.Marshal.Alloc(free)
+import Foreign.Marshal.Array(peekArray, copyArray, advancePtr)
+import Foreign.ForeignPtr(ForeignPtr, castForeignPtr)
+import Foreign.Ptr(Ptr)
+import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf)
+import Foreign.C.String
+import Foreign.C.Types
+import Data.Complex
+import Control.Monad(when)
+import System.IO.Unsafe(unsafePerformIO)
+
+#if __GLASGOW_HASKELL__ >= 605
+import GHC.ForeignPtr           (mallocPlainForeignPtrBytes)
+#else
+import Foreign.ForeignPtr       (mallocForeignPtrBytes)
+#endif
+
+import GHC.Base
+#if __GLASGOW_HASKELL__ < 612
+import GHC.IOBase hiding (liftIO)
+#endif
+
+import qualified Data.Vector.Storable as Vector
+import Data.Vector.Storable(Vector,
+                            fromList,
+                            unsafeToForeignPtr,
+                            unsafeFromForeignPtr,
+                            unsafeWith)
+
+
+-- | Number of elements
+dim :: (Storable t) => Vector t -> Int
+dim = Vector.length
+
+
+-- 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
+vec x f = unsafeWith x $ \p -> do
+    let v g = do
+        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 negative dim: "++show n)
+    fp <- doMalloc undefined
+    return $ unsafeFromForeignPtr fp 0 n
+  where
+    --
+    -- Use the much cheaper Haskell heap allocated storage
+    -- for foreign pointer space we control
+    --
+    doMalloc :: Storable b => b -> IO (ForeignPtr b)
+    doMalloc dummy = do
+#if __GLASGOW_HASKELL__ >= 605
+        mallocPlainForeignPtrBytes (n * sizeOf dummy)
+#else
+        mallocForeignPtrBytes      (n * sizeOf dummy)
+#endif
+
+{- | creates a Vector from a list:
+
+@> fromList [2,3,5,7]
+4 |> [2.0,3.0,5.0,7.0]@
+
+-}
+
+safeRead v = inlinePerformIO . unsafeWith v
+{-# INLINE safeRead #-}
+
+inlinePerformIO :: IO a -> a
+inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
+{-# INLINE inlinePerformIO #-}
+
+{- | extracts the Vector elements to a list
+
+>>> toList (linspace 5 (1,10))
+[1.0,3.25,5.5,7.75,10.0]
+
+-}
+toList :: Storable a => Vector a -> [a]
+toList v = safeRead v $ peekArray (dim v)
+
+{- | Create a vector from a list of elements and explicit dimension. The input
+     list is explicitly truncated if it is too long, so it may safely
+     be used, for instance, with infinite lists.
+
+>>> 5 |> [1..]
+fromList [1.0,2.0,3.0,4.0,5.0]
+
+-}
+(|>) :: (Storable a) => Int -> [a] -> Vector a
+infixl 9 |>
+n |> l = if length l' == n
+            then fromList l'
+            else error "list too short for |>"
+  where l' = take n l
+
+
+-- | access to Vector elements without range checking
+at' :: Storable a => Vector a -> Int -> a
+at' v n = safeRead v $ flip peekElemOff n
+{-# INLINE at' #-}
+
+--
+-- turn off bounds checking with -funsafe at configure time.
+-- ghc will optimise away the salways true case at compile time.
+--
+#if defined(UNSAFE)
+safe :: Bool
+safe = False
+#else
+safe = True
+#endif
+
+-- | access to Vector elements with range checking.
+at :: Storable a => Vector a -> Int -> a
+at v n
+    | safe      = if n >= 0 && n < dim v
+                    then at' v n
+                    else error "vector index out of range"
+    | otherwise = at' v n
+{-# INLINE at #-}
+
+{- | takes a number of consecutive elements from a Vector
+
+>>> subVector 2 3 (fromList [1..10])
+fromList [3.0,4.0,5.0]
+
+-}
+subVector :: Storable t => Int       -- ^ index of the starting element
+                        -> Int       -- ^ number of elements to extract
+                        -> Vector t  -- ^ source
+                        -> Vector t  -- ^ result
+subVector = Vector.slice
+
+
+{- | Reads a vector position:
+
+>>> fromList [0..9] @> 7
+7.0
+
+-}
+(@>) :: Storable t => Vector t -> Int -> t
+infixl 9 @>
+(@>) = at
+
+
+{- | concatenate a list of vectors
+
+>>> vjoin [fromList [1..5::Double], konst 1 3]
+fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]
+
+-}
+vjoin :: Storable t => [Vector t] -> Vector t
+vjoin [] = fromList []
+vjoin [v] = v
+vjoin as = unsafePerformIO $ do
+    let tot = sum (map dim as)
+    r <- createVector tot
+    unsafeWith r $ \ptr ->
+        joiner as tot ptr
+    return r
+  where joiner [] _ _ = return ()
+        joiner (v:cs) _ p = do
+            let n = dim v
+            unsafeWith v $ \pb -> copyArray p pb n
+            joiner cs 0 (advancePtr p n)
+
+
+{- | Extract consecutive subvectors of the given sizes.
+
+>>> takesV [3,4] (linspace 10 (1,10::Double))
+[fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]]
+
+-}
+takesV :: Storable t => [Int] -> Vector t -> [Vector t]
+takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (show $ dim w)
+            | otherwise = go ms w
+    where go [] _ = []
+          go (n:ns) v = subVector 0 n v
+                      : go ns (subVector n (dim v - n) v)
+
+---------------------------------------------------------------
+
+-- | transforms a complex vector into a real vector with alternating real and imaginary parts 
+asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a
+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 :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a)
+asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)
+    where (fp,i,n) = unsafeToForeignPtr v
+
+---------------------------------------------------------------
+
+float2DoubleV :: Vector Float -> Vector Double
+float2DoubleV v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 c_float2double vec v vec r "float2double"
+    return r
+
+double2FloatV :: Vector Double -> Vector Float
+double2FloatV v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 c_double2float vec v vec r "double2float2"
+    return r
+
+
+foreign import ccall unsafe "float2double" c_float2double:: TFV
+foreign import ccall unsafe "double2float" c_double2float:: TVF
+
+---------------------------------------------------------------
+
+stepF :: Vector Float -> Vector Float
+stepF v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 c_stepF vec v vec r "stepF"
+    return r
+
+stepD :: Vector Double -> Vector Double
+stepD v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 c_stepD vec v vec r "stepD"
+    return r
+
+foreign import ccall unsafe "stepF" c_stepF :: TFF
+foreign import ccall unsafe "stepD" c_stepD :: TVV
+
+---------------------------------------------------------------
+
+condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float
+condF x y l e g = unsafePerformIO $ do
+    r <- createVector (dim x)
+    app6 c_condF vec x vec y vec l vec e vec g vec r "condF"
+    return r
+
+condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double
+condD x y l e g = unsafePerformIO $ do
+    r <- createVector (dim x)
+    app6 c_condD vec x vec y vec l vec e vec g vec r "condD"
+    return r
+
+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
+
+--------------------------------------------------------------------------------
+
+conjugateAux fun x = unsafePerformIO $ do
+    v <- createVector (dim x)
+    app2 fun vec x vec v "conjugateAux"
+    return v
+
+conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
+conjugateQ = conjugateAux c_conjugateQ
+foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV
+
+conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
+conjugateC = conjugateAux c_conjugateC
+foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV
+
+--------------------------------------------------------------------------------
+
+cloneVector :: Storable t => Vector t -> IO (Vector t)
+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"
+        return r
+
+------------------------------------------------------------------
+
+-- | map on Vectors
+mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
+mapVector f v = unsafePerformIO $ do
+    w <- createVector (dim v)
+    unsafeWith v $ \p ->
+        unsafeWith w $ \q -> do
+            let go (-1) = return ()
+                go !k = do x <- peekElemOff p k
+                           pokeElemOff      q k (f x)
+                           go (k-1)
+            go (dim v -1)
+    return w
+{-# INLINE mapVector #-}
+
+-- | zipWith for Vectors
+zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c
+zipVectorWith f u v = unsafePerformIO $ do
+    let n = min (dim u) (dim v)
+    w <- createVector n
+    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
+                               pokeElemOff      pw k (f x y)
+                               go (k-1)
+                go (n -1)
+    return w
+{-# INLINE zipVectorWith #-}
+
+-- | unzipWith for Vectors
+unzipVectorWith :: (Storable (a,b), Storable c, Storable d) 
+                   => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)
+unzipVectorWith f u = unsafePerformIO $ do
+      let n = dim u
+      v <- createVector n
+      w <- createVector n
+      unsafeWith u $ \pu ->
+          unsafeWith v $ \pv ->
+              unsafeWith w $ \pw -> do
+                  let go (-1) = return ()
+                      go !k   = do z <- peekElemOff pu k
+                                   let (x,y) = f z 
+                                   pokeElemOff      pv k x
+                                   pokeElemOff      pw k y
+                                   go (k-1)
+                  go (n-1)
+      return (v,w)
+{-# INLINE unzipVectorWith #-}
+
+foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b
+foldVector f x v = unsafePerformIO $
+    unsafeWith v $ \p -> do
+        let go (-1) s = return s
+            go !k !s = do y <- peekElemOff p k
+                          go (k-1::Int) (f y s)
+        go (dim v -1) x
+{-# INLINE foldVector #-}
+
+-- the zero-indexed index is passed to the folding function
+foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b
+foldVectorWithIndex f x v = unsafePerformIO $
+    unsafeWith v $ \p -> do
+        let go (-1) s = return s
+            go !k !s = do y <- peekElemOff p k
+                          go (k-1::Int) (f k y s)
+        go (dim v -1) x
+{-# INLINE foldVectorWithIndex #-}
+
+foldLoop f s0 d = go (d - 1) s0
+     where
+       go 0 s = f (0::Int) s
+       go !j !s = go (j - 1) (f j s)
+
+foldVectorG f s0 v = foldLoop g s0 (dim v)
+    where g !k !s = f k (at' v) s
+          {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)
+{-# INLINE foldVectorG #-}
+
+-------------------------------------------------------------------
+
+-- | monadic map over Vectors
+--    the monad @m@ must be strict
+mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b)
+mapVectorM f v = do
+    w <- return $! unsafePerformIO $! createVector (dim v)
+    mapVectorM' w 0 (dim v -1)
+    return w
+    where mapVectorM' w' !k !t
+              | k == t               = do
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       y <- f x
+                                       return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
+              | otherwise            = do
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       y <- f x
+                                       _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
+                                       mapVectorM' w' (k+1) t
+{-# INLINE mapVectorM #-}
+
+-- | monadic map over Vectors
+mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m ()
+mapVectorM_ f v = do
+    mapVectorM' 0 (dim v -1)
+    where mapVectorM' !k !t
+              | k == t            = do
+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
+                                    f x
+              | otherwise         = do
+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                    _ <- f x
+                                    mapVectorM' (k+1) t
+{-# INLINE mapVectorM_ #-}
+
+-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
+--    the monad @m@ must be strict
+mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b)
+mapVectorWithIndexM f v = do
+    w <- return $! unsafePerformIO $! createVector (dim v)
+    mapVectorM' w 0 (dim v -1)
+    return w
+    where mapVectorM' w' !k !t
+              | k == t               = do
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       y <- f k x
+                                       return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
+              | otherwise            = do
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       y <- f k x
+                                       _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
+                                       mapVectorM' w' (k+1) t
+{-# INLINE mapVectorWithIndexM #-}
+
+-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
+mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m ()
+mapVectorWithIndexM_ f v = do
+    mapVectorM' 0 (dim v -1)
+    where mapVectorM' !k !t
+              | k == t            = do
+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
+                                    f k x
+              | otherwise         = do
+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                    _ <- f k x
+                                    mapVectorM' (k+1) t
+{-# INLINE mapVectorWithIndexM_ #-}
+
+
+mapVectorWithIndex :: (Storable a, Storable b) => (Int -> a -> b) -> Vector a -> Vector b
+--mapVectorWithIndex g = head . mapVectorWithIndexM (\a b -> [g a b])
+mapVectorWithIndex f v = unsafePerformIO $ do
+    w <- createVector (dim v)
+    unsafeWith v $ \p ->
+        unsafeWith w $ \q -> do
+            let go (-1) = return ()
+                go !k = do x <- peekElemOff p k
+                           pokeElemOff      q k (f k x)
+                           go (k-1)
+            go (dim v -1)
+    return w
+{-# INLINE mapVectorWithIndex #-}
+
+-------------------------------------------------------------------
+
+
+-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
+fscanfVector :: FilePath -> Int -> IO (Vector Double)
+fscanfVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
+    free charname
+    return res
+
+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 ()
+fprintfVector filename fmt v = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
+    free charname
+    free charfmt
+
+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)
+freadVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
+    free charname
+    return res
+
+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 ()
+fwriteVector filename v = do
+    charname <- newCString filename
+    app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
+    free charname
+
+foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
+
diff --git a/packages/hmatrix/src/Data/Packed/Matrix.hs b/packages/hmatrix/src/Data/Packed/Matrix.hs
new file mode 100644
index 0000000..d94d167
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Matrix.hs
@@ -0,0 +1,490 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE CPP #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Matrix
+-- Copyright   :  (c) Alberto Ruiz 2007-10
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+--
+-- A Matrix representation suitable for numerical computations using LAPACK and GSL.
+--
+-- This module provides basic functions for manipulation of structure.
+
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+module Data.Packed.Matrix (
+    Matrix,
+    Element,
+    rows,cols,
+    (><),
+    trans,
+    reshape, flatten,
+    fromLists, toLists, buildMatrix,
+    (@@>),
+    asRow, asColumn,
+    fromRows, toRows, fromColumns, toColumns,
+    fromBlocks, diagBlock, toBlocks, toBlocksEvery,
+    repmat,
+    flipud, fliprl,
+    subMatrix, takeRows, dropRows, takeColumns, dropColumns,
+    extractRows, extractColumns,
+    diagRect, takeDiag,
+    mapMatrix, mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_,
+    liftMatrix, liftMatrix2, liftMatrix2Auto,fromArray2D
+) where
+
+import Data.Packed.Internal
+import qualified Data.Packed.ST as ST
+import Data.Array
+
+import Data.List(transpose,intersperse)
+import Foreign.Storable(Storable)
+import Control.Monad(liftM)
+
+-------------------------------------------------------------------
+
+#ifdef BINARY
+
+import Data.Binary
+import Control.Monad(replicateM)
+
+instance (Binary a, Element a, Storable a) => Binary (Matrix a) where
+    put m = do
+            let r = rows m
+            let c = cols m
+            put r
+            put c
+            mapM_ (\i -> mapM_ (\j -> put $ m @@> (i,j)) [0..(c-1)]) [0..(r-1)]
+    get = do
+          r <- get
+          c <- get
+          xs <- replicateM r $ replicateM c get
+          return $ fromLists xs
+
+#endif
+
+-------------------------------------------------------------------
+
+instance (Show a, Element a) => (Show (Matrix a)) where
+    show m | rows m == 0 || cols m == 0 = sizes m ++" []"
+    show m = (sizes m++) . dsp . map (map show) . toLists $ m
+
+sizes m = "("++show (rows m)++"><"++show (cols m)++")\n"
+
+dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
+    where
+        mt = transpose as
+        longs = map (maximum . map length) mt
+        mtp = zipWith (\a b -> map (pad a) b) longs mt
+        pad n str = replicate (n - length str) ' ' ++ str
+        unwords' = concat . intersperse ", "
+
+------------------------------------------------------------------
+
+instance (Element a, Read a) => Read (Matrix a) where
+    readsPrec _ s = [((rs><cs) . read $ listnums, rest)]
+        where (thing,rest) = breakAt ']' s
+              (dims,listnums) = breakAt ')' thing
+              cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
+              rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
+
+
+breakAt c l = (a++[c],tail b) where
+    (a,b) = break (==c) l
+
+------------------------------------------------------------------
+
+-- | creates a matrix from a vertical list of matrices
+joinVert :: Element t => [Matrix t] -> Matrix t
+joinVert [] = emptyM 0 0
+joinVert ms = case common cols ms of
+    Nothing -> error "(impossible) joinVert on matrices with different number of columns"
+    Just c  -> matrixFromVector RowMajor (sum (map rows ms)) c $ vjoin (map flatten ms)
+
+-- | creates a matrix from a horizontal list of matrices
+joinHoriz :: Element t => [Matrix t] -> Matrix t
+joinHoriz ms = trans. joinVert . map trans $ ms
+
+{- | Create a matrix from blocks given as a list of lists of matrices.
+
+Single row-column components are automatically expanded to match the
+corresponding common row and column:
+
+@
+disp = putStr . dispf 2
+@
+
+>>> disp $ fromBlocks [[ident 5, 7, row[10,20]], [3, diagl[1,2,3], 0]]
+8x10
+1  0  0  0  0  7  7  7  10  20
+0  1  0  0  0  7  7  7  10  20
+0  0  1  0  0  7  7  7  10  20
+0  0  0  1  0  7  7  7  10  20
+0  0  0  0  1  7  7  7  10  20
+3  3  3  3  3  1  0  0   0   0
+3  3  3  3  3  0  2  0   0   0
+3  3  3  3  3  0  0  3   0   0
+
+-}
+fromBlocks :: Element t => [[Matrix t]] -> Matrix t
+fromBlocks = fromBlocksRaw . adaptBlocks
+
+fromBlocksRaw mms = joinVert . map joinHoriz $ mms
+
+adaptBlocks ms = ms' where
+    bc = case common length ms of
+          Just c -> c
+          Nothing -> error "fromBlocks requires rectangular [[Matrix]]"
+    rs = map (compatdim . map rows) ms
+    cs = map (compatdim . map cols) (transpose ms)
+    szs = sequence [rs,cs]
+    ms' = splitEvery bc $ zipWith g szs (concat ms)
+
+    g [Just nr,Just nc] m
+                | nr == r && nc == c = m
+                | r == 1 && c == 1 = matrixFromVector RowMajor nr nc (constantD x (nr*nc))
+                | r == 1 = fromRows (replicate nr (flatten m))
+                | otherwise = fromColumns (replicate nc (flatten m))
+      where
+        r = rows m
+        c = cols m
+        x = m@@>(0,0)
+    g _ _ = error "inconsistent dimensions in fromBlocks"
+
+
+--------------------------------------------------------------------------------
+
+{- | create a block diagonal matrix
+
+>>>  disp 2 $ diagBlock [konst 1 (2,2), konst 2 (3,5), col [5,7]]
+7x8
+1  1  0  0  0  0  0  0
+1  1  0  0  0  0  0  0
+0  0  2  2  2  2  2  0
+0  0  2  2  2  2  2  0
+0  0  2  2  2  2  2  0
+0  0  0  0  0  0  0  5
+0  0  0  0  0  0  0  7
+
+>>> diagBlock [(0><4)[], konst 2 (2,3)]  :: Matrix Double
+(2><7)
+ [ 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0
+ , 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0 ]
+
+-}
+diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t
+diagBlock ms = fromBlocks $ zipWith f ms [0..]
+  where
+    f m k = take n $ replicate k z ++ m : repeat z
+    n = length ms
+    z = (1><1) [0]
+
+--------------------------------------------------------------------------------
+
+
+-- | Reverse rows
+flipud :: Element t => Matrix t -> Matrix t
+flipud m = extractRows [r-1,r-2 .. 0] $ m
+  where
+    r = rows m
+
+-- | Reverse columns
+fliprl :: Element t => Matrix t -> Matrix t
+fliprl m = extractColumns [c-1,c-2 .. 0] $ m
+  where
+    c = cols m
+
+------------------------------------------------------------
+
+{- | creates a rectangular diagonal matrix:
+
+>>> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double
+(4><5)
+ [ 10.0,  7.0,  7.0, 7.0, 7.0
+ ,  7.0, 20.0,  7.0, 7.0, 7.0
+ ,  7.0,  7.0, 30.0, 7.0, 7.0
+ ,  7.0,  7.0,  7.0, 7.0, 7.0 ]
+
+-}
+diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t
+diagRect z v r c = ST.runSTMatrix $ do
+        m <- ST.newMatrix z r c
+        let d = min r c `min` (dim v)
+        mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]
+        return m
+
+-- | extracts the diagonal from a rectangular matrix
+takeDiag :: (Element t) => Matrix t -> Vector t
+takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
+
+------------------------------------------------------------
+
+{- | An easy way to create a matrix:
+
+>>> (2><3)[2,4,7,-3,11,0]
+(2><3)
+ [  2.0,  4.0, 7.0
+ , -3.0, 11.0, 0.0 ]
+
+This is the format produced by the instances of Show (Matrix a), which
+can also be used for input.
+
+The input list is explicitly truncated, so that it can
+safely be used with lists that are too long (like infinite lists).
+
+>>> (2><3)[1..]
+(2><3)
+ [ 1.0, 2.0, 3.0
+ , 4.0, 5.0, 6.0 ]
+
+
+-}
+(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a
+r >< c = f where
+    f l | dim v == r*c = matrixFromVector RowMajor r c v
+        | otherwise    = error $ "inconsistent list size = "
+                                 ++show (dim v) ++" in ("++show r++"><"++show c++")"
+        where v = fromList $ take (r*c) l
+
+----------------------------------------------------------------
+
+-- | Creates a matrix with the first n rows of another matrix
+takeRows :: Element t => Int -> Matrix t -> Matrix t
+takeRows n mt = subMatrix (0,0) (n, cols mt) mt
+-- | Creates a copy of a matrix without the first n rows
+dropRows :: Element t => Int -> Matrix t -> Matrix t
+dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt
+-- |Creates a matrix with the first n columns of another matrix
+takeColumns :: Element t => Int -> Matrix t -> Matrix t
+takeColumns n mt = subMatrix (0,0) (rows mt, n) mt
+-- | Creates a copy of a matrix without the first n columns
+dropColumns :: Element t => Int -> Matrix t -> Matrix t
+dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
+
+----------------------------------------------------------------
+
+{- | Creates a 'Matrix' from a list of lists (considered as rows).
+
+>>> fromLists [[1,2],[3,4],[5,6]]
+(3><2)
+ [ 1.0, 2.0
+ , 3.0, 4.0
+ , 5.0, 6.0 ]
+
+-}
+fromLists :: Element t => [[t]] -> Matrix t
+fromLists = fromRows . map fromList
+
+-- | creates a 1-row matrix from a vector
+--
+-- >>> asRow (fromList [1..5])
+--  (1><5)
+--   [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
+--
+asRow :: Storable a => Vector a -> Matrix a
+asRow v = reshape (dim v) v
+
+-- | creates a 1-column matrix from a vector
+--
+-- >>> asColumn (fromList [1..5])
+-- (5><1)
+--  [ 1.0
+--  , 2.0
+--  , 3.0
+--  , 4.0
+--  , 5.0 ]
+--
+asColumn :: Storable a => Vector a -> Matrix a
+asColumn = trans . asRow
+
+
+
+{- | creates a Matrix of the specified size using the supplied function to
+     to map the row\/column position to the value at that row\/column position.
+
+@> buildMatrix 3 4 (\\(r,c) -> fromIntegral r * fromIntegral c)
+(3><4)
+ [ 0.0, 0.0, 0.0, 0.0, 0.0
+ , 0.0, 1.0, 2.0, 3.0, 4.0
+ , 0.0, 2.0, 4.0, 6.0, 8.0]@
+
+Hilbert matrix of order N:
+
+@hilb n = buildMatrix n n (\\(i,j)->1/(fromIntegral i + fromIntegral j +1))@
+
+-}
+buildMatrix :: Element a => Int -> Int -> ((Int, Int) -> a) -> Matrix a
+buildMatrix rc cc f =
+    fromLists $ map (map f)
+        $ map (\ ri -> map (\ ci -> (ri, ci)) [0 .. (cc - 1)]) [0 .. (rc - 1)]
+
+-----------------------------------------------------
+
+fromArray2D :: (Storable e) => Array (Int, Int) e -> Matrix e
+fromArray2D m = (r><c) (elems m)
+    where ((r0,c0),(r1,c1)) = bounds m
+          r = r1-r0+1
+          c = c1-c0+1
+
+
+-- | rearranges the rows of a matrix according to the order given in a list of integers.
+extractRows :: Element t => [Int] -> Matrix t -> Matrix t
+extractRows [] m = emptyM 0 (cols m)
+extractRows l m = fromRows $ extract (toRows m) l
+  where
+    extract l' is = [l'!!i | i<- map verify is]
+    verify k
+        | k >= 0 && k < rows m = k
+        | otherwise = error $ "can't extract row "
+                           ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
+
+-- | rearranges the rows of a matrix according to the order given in a list of integers.
+extractColumns :: Element t => [Int] -> Matrix t -> Matrix t
+extractColumns l m = trans . extractRows (map verify l) . trans $ m
+  where
+    verify k
+        | k >= 0 && k < cols m = k
+        | otherwise = error $ "can't extract column "
+                           ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
+
+
+
+{- | creates matrix by repetition of a matrix a given number of rows and columns
+
+>>> repmat (ident 2) 2 3
+(4><6)
+ [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
+ , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
+ , 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
+ , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]
+
+-}
+repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t
+repmat m r c
+    | r == 0 || c == 0 = emptyM (r*rows m) (c*cols m)
+    | otherwise = fromBlocks $ replicate r $ replicate c $ m
+
+-- | A version of 'liftMatrix2' which automatically adapt matrices with a single row or column to match the dimensions of the other matrix.
+liftMatrix2Auto :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
+liftMatrix2Auto f m1 m2
+    | compat' m1 m2 = lM f m1  m2
+    | ok            = lM f m1' m2'
+    | otherwise = error $ "nonconformable matrices in liftMatrix2Auto: " ++ shSize m1 ++ ", " ++ shSize m2
+  where
+    (r1,c1) = size m1
+    (r2,c2) = size m2
+    r = max r1 r2
+    c = max c1 c2
+    r0 = min r1 r2
+    c0 = min c1 c2
+    ok = r0 == 1 || r1 == r2 && c0 == 1 || c1 == c2
+    m1' = conformMTo (r,c) m1
+    m2' = conformMTo (r,c) m2
+
+-- FIXME do not flatten if equal order
+lM f m1 m2 = matrixFromVector
+                RowMajor
+                (max (rows m1) (rows m2))
+                (max (cols m1) (cols m2))
+                (f (flatten m1) (flatten m2))
+
+compat' :: Matrix a -> Matrix b -> Bool
+compat' m1 m2 = s1 == (1,1) || s2 == (1,1) || s1 == s2
+  where
+    s1 = size m1
+    s2 = size m2
+
+------------------------------------------------------------
+
+toBlockRows [r] m | r == rows m = [m]
+toBlockRows rs m = map (reshape (cols m)) (takesV szs (flatten m))
+    where szs = map (* cols m) rs
+
+toBlockCols [c] m | c == cols m = [m]
+toBlockCols cs m = map trans . toBlockRows cs . trans $ m
+
+-- | Partition a matrix into blocks with the given numbers of rows and columns.
+-- The remaining rows and columns are discarded.
+toBlocks :: (Element t) => [Int] -> [Int] -> Matrix t -> [[Matrix t]]
+toBlocks rs cs m = map (toBlockCols cs) . toBlockRows rs $ m
+
+-- | Fully partition a matrix into blocks of the same size. If the dimensions are not
+-- a multiple of the given size the last blocks will be smaller.
+toBlocksEvery :: (Element t) => Int -> Int -> Matrix t -> [[Matrix t]]
+toBlocksEvery r c m
+    | r < 1 || c < 1 = error $ "toBlocksEvery expects block sizes > 0, given "++show r++" and "++ show c
+    | otherwise = toBlocks rs cs m
+  where
+    (qr,rr) = rows m `divMod` r
+    (qc,rc) = cols m `divMod` c
+    rs = replicate qr r ++ if rr > 0 then [rr] else []
+    cs = replicate qc c ++ if rc > 0 then [rc] else []
+
+-------------------------------------------------------------------
+
+-- Given a column number and a function taking matrix indexes, returns
+-- a function which takes vector indexes (that can be used on the
+-- flattened matrix).
+mk :: Int -> ((Int, Int) -> t) -> (Int -> t)
+mk c g = \k -> g (divMod k c)
+
+{- |
+
+>>> mapMatrixWithIndexM_ (\(i,j) v -> printf "m[%d,%d] = %.f\n" i j v :: IO()) ((2><3)[1 :: Double ..])
+m[0,0] = 1
+m[0,1] = 2
+m[0,2] = 3
+m[1,0] = 4
+m[1,1] = 5
+m[1,2] = 6
+
+-}
+mapMatrixWithIndexM_
+  :: (Element a, Num a, Monad m) =>
+      ((Int, Int) -> a -> m ()) -> Matrix a -> m ()
+mapMatrixWithIndexM_ g m = mapVectorWithIndexM_ (mk c g) . flatten $ m
+  where
+    c = cols m
+
+{- |
+
+>>> mapMatrixWithIndexM (\(i,j) v -> Just $ 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
+Just (3><3)
+ [ 100.0,   1.0,   2.0
+ ,  10.0, 111.0,  12.0
+ ,  20.0,  21.0, 122.0 ]
+
+-}
+mapMatrixWithIndexM
+  :: (Element a, Storable b, Monad m) =>
+      ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
+mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m 
+    where
+      c = cols m
+
+{- |
+
+>>> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
+(3><3)
+ [ 100.0,   1.0,   2.0
+ ,  10.0, 111.0,  12.0
+ ,  20.0,  21.0, 122.0 ]
+
+ -}
+mapMatrixWithIndex
+  :: (Element a, Storable b) =>
+      ((Int, Int) -> a -> b) -> Matrix a -> Matrix b
+mapMatrixWithIndex g m = reshape c . mapVectorWithIndex (mk c g) . flatten $ m
+    where
+      c = cols m
+
+mapMatrix :: (Storable a, Storable b) => (a -> b) -> Matrix a -> Matrix b
+mapMatrix f = liftMatrix (mapVector f)
diff --git a/packages/hmatrix/src/Data/Packed/Random.hs b/packages/hmatrix/src/Data/Packed/Random.hs
new file mode 100644
index 0000000..e8b0268
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Random.hs
@@ -0,0 +1,57 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Vector
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+--
+-- Random vectors and matrices.
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Random (
+    Seed,
+    RandDist(..),
+    randomVector,
+    gaussianSample,
+    uniformSample
+) where
+
+import Numeric.GSL.Vector
+import Data.Packed
+import Numeric.ContainerBoot
+import Numeric.LinearAlgebra.Algorithms
+
+
+type Seed = Int
+
+-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
+-- Gaussian distribution.
+gaussianSample :: Seed
+               -> Int -- ^ number of rows
+               -> Vector Double -- ^ mean vector
+               -> Matrix Double -- ^ covariance matrix
+               -> Matrix Double -- ^ result
+gaussianSample seed n med cov = m where
+    c = dim med
+    meds = konst' 1 n `outer` med
+    rs = reshape c $ randomVector seed Gaussian (c * n)
+    m = rs `mXm` cholSH cov `add` meds
+
+-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
+-- uniform distribution.
+uniformSample :: Seed
+               -> Int -- ^ number of rows
+               -> [(Double,Double)] -- ^ ranges for each column
+               -> Matrix Double -- ^ result
+uniformSample seed n rgs = m where
+    (as,bs) = unzip rgs
+    a = fromList as
+    cs = zipWith subtract as bs
+    d = dim a
+    dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
+    am = konst' 1 n `outer` a
+    m = fromColumns (zipWith scale cs dat) `add` am
+
diff --git a/packages/hmatrix/src/Data/Packed/ST.hs b/packages/hmatrix/src/Data/Packed/ST.hs
new file mode 100644
index 0000000..1cef296
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/ST.hs
@@ -0,0 +1,179 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE Rank2Types    #-}
+{-# LANGUAGE BangPatterns  #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.ST
+-- Copyright   :  (c) Alberto Ruiz 2008
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- In-place manipulation inside the ST monad.
+-- See examples/inplace.hs in the distribution.
+--
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+module Data.Packed.ST (
+    -- * Mutable Vectors
+    STVector, newVector, thawVector, freezeVector, runSTVector,
+    readVector, writeVector, modifyVector, liftSTVector,
+    -- * Mutable Matrices
+    STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix,
+    readMatrix, writeMatrix, modifyMatrix, liftSTMatrix,
+    -- * Unsafe functions
+    newUndefinedVector,
+    unsafeReadVector, unsafeWriteVector,
+    unsafeThawVector, unsafeFreezeVector,
+    newUndefinedMatrix,
+    unsafeReadMatrix, unsafeWriteMatrix,
+    unsafeThawMatrix, unsafeFreezeMatrix
+) where
+
+import Data.Packed.Internal
+
+import Control.Monad.ST(ST, runST)
+import Foreign.Storable(Storable, peekElemOff, pokeElemOff)
+
+#if MIN_VERSION_base(4,4,0)
+import Control.Monad.ST.Unsafe(unsafeIOToST)
+#else
+import Control.Monad.ST(unsafeIOToST)
+#endif
+
+{-# INLINE ioReadV #-}
+ioReadV :: Storable t => Vector t -> Int -> IO t
+ioReadV v k = unsafeWith v $ \s -> peekElemOff s k
+
+{-# INLINE ioWriteV #-}
+ioWriteV :: Storable t => Vector t -> Int -> t -> IO ()
+ioWriteV v k x = unsafeWith v $ \s -> pokeElemOff s k x
+
+newtype STVector s t = STVector (Vector t)
+
+thawVector :: Storable t => Vector t -> ST s (STVector s t)
+thawVector = unsafeIOToST . fmap STVector . cloneVector
+
+unsafeThawVector :: Storable t => Vector t -> ST s (STVector s t)
+unsafeThawVector = unsafeIOToST . return . STVector
+
+runSTVector :: Storable t => (forall s . ST s (STVector s t)) -> Vector t
+runSTVector st = runST (st >>= unsafeFreezeVector)
+
+{-# INLINE unsafeReadVector #-}
+unsafeReadVector :: Storable t => STVector s t -> Int -> ST s t
+unsafeReadVector   (STVector x) = unsafeIOToST . ioReadV x
+
+{-# INLINE unsafeWriteVector #-}
+unsafeWriteVector :: Storable t => STVector s t -> Int -> t -> ST s ()
+unsafeWriteVector  (STVector x) k = unsafeIOToST . ioWriteV x k
+
+{-# INLINE modifyVector #-}
+modifyVector :: (Storable t) => STVector s t -> Int -> (t -> t) -> ST s ()
+modifyVector x k f = readVector x k >>= return . f >>= unsafeWriteVector x k
+
+liftSTVector :: (Storable t) => (Vector t -> a) -> STVector s1 t -> ST s2 a
+liftSTVector f (STVector x) = unsafeIOToST . fmap f . cloneVector $ x
+
+freezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
+freezeVector v = liftSTVector id v
+
+unsafeFreezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
+unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x
+
+{-# INLINE safeIndexV #-}
+safeIndexV f (STVector v) k
+    | k < 0 || k>= dim v = error $ "out of range error in vector (dim="
+                                   ++show (dim v)++", pos="++show k++")"
+    | otherwise = f (STVector v) k
+
+{-# INLINE readVector #-}
+readVector :: Storable t => STVector s t -> Int -> ST s t
+readVector = safeIndexV unsafeReadVector
+
+{-# INLINE writeVector #-}
+writeVector :: Storable t => STVector s t -> Int -> t -> ST s ()
+writeVector = safeIndexV unsafeWriteVector
+
+newUndefinedVector :: Storable t => Int -> ST s (STVector s t)
+newUndefinedVector = unsafeIOToST . fmap STVector . createVector
+
+{-# INLINE newVector #-}
+newVector :: Storable t => t -> Int -> ST s (STVector s t)
+newVector x n = do
+    v <- newUndefinedVector n
+    let go (-1) = return v
+        go !k = unsafeWriteVector v k x >> go (k-1 :: Int)
+    go (n-1)
+
+-------------------------------------------------------------------------
+
+{-# INLINE ioReadM #-}
+ioReadM :: Storable t => Matrix t -> Int -> Int -> IO t
+ioReadM (Matrix _ nc cv RowMajor) r c = ioReadV cv (r*nc+c)
+ioReadM (Matrix nr _ fv ColumnMajor) r c = ioReadV fv (c*nr+r)
+
+{-# INLINE ioWriteM #-}
+ioWriteM :: Storable t => Matrix t -> Int -> Int -> t -> IO ()
+ioWriteM (Matrix _ nc cv RowMajor) r c val = ioWriteV cv (r*nc+c) val
+ioWriteM (Matrix nr _ fv ColumnMajor) r c val = ioWriteV fv (c*nr+r) val
+
+newtype STMatrix s t = STMatrix (Matrix t)
+
+thawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)
+thawMatrix = unsafeIOToST . fmap STMatrix . cloneMatrix
+
+unsafeThawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)
+unsafeThawMatrix = unsafeIOToST . return . STMatrix
+
+runSTMatrix :: Storable t => (forall s . ST s (STMatrix s t)) -> Matrix t
+runSTMatrix st = runST (st >>= unsafeFreezeMatrix)
+
+{-# INLINE unsafeReadMatrix #-}
+unsafeReadMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t
+unsafeReadMatrix   (STMatrix x) r = unsafeIOToST . ioReadM x r
+
+{-# INLINE unsafeWriteMatrix #-}
+unsafeWriteMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
+unsafeWriteMatrix  (STMatrix x) r c = unsafeIOToST . ioWriteM x r c
+
+{-# INLINE modifyMatrix #-}
+modifyMatrix :: (Storable t) => STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
+modifyMatrix x r c f = readMatrix x r c >>= return . f >>= unsafeWriteMatrix x r c
+
+liftSTMatrix :: (Storable t) => (Matrix t -> a) -> STMatrix s1 t -> ST s2 a
+liftSTMatrix f (STMatrix x) = unsafeIOToST . fmap f . cloneMatrix $ x
+
+unsafeFreezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
+unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x
+
+freezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
+freezeMatrix m = liftSTMatrix id m
+
+cloneMatrix (Matrix r c d o) = cloneVector d >>= return . (\d' -> Matrix r c d' o)
+
+{-# INLINE safeIndexM #-}
+safeIndexM f (STMatrix m) r c
+    | r<0 || r>=rows m ||
+      c<0 || c>=cols m = error $ "out of range error in matrix (size="
+                                 ++show (rows m,cols m)++", pos="++show (r,c)++")"
+    | otherwise = f (STMatrix m) r c
+
+{-# INLINE readMatrix #-}
+readMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t
+readMatrix = safeIndexM unsafeReadMatrix
+
+{-# INLINE writeMatrix #-}
+writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
+writeMatrix = safeIndexM unsafeWriteMatrix
+
+newUndefinedMatrix :: Storable t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t)
+newUndefinedMatrix ord r c = unsafeIOToST $ fmap STMatrix $ createMatrix ord r c
+
+{-# NOINLINE newMatrix #-}
+newMatrix :: Storable t => t -> Int -> Int -> ST s (STMatrix s t)
+newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c)
diff --git a/packages/hmatrix/src/Data/Packed/Vector.hs b/packages/hmatrix/src/Data/Packed/Vector.hs
new file mode 100644
index 0000000..b5a4318
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Vector.hs
@@ -0,0 +1,96 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE CPP #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Vector
+-- Copyright   :  (c) Alberto Ruiz 2007-10
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+--
+-- 1D arrays suitable for numeric computations using external libraries.
+--
+-- This module provides basic functions for manipulation of structure.
+--
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+module Data.Packed.Vector (
+    Vector,
+    fromList, (|>), toList, buildVector,
+    dim, (@>),
+    subVector, takesV, vjoin, join,
+    mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith,
+    mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
+    foldLoop, foldVector, foldVectorG, foldVectorWithIndex
+) where
+
+import Data.Packed.Internal.Vector
+import Foreign.Storable
+
+-------------------------------------------------------------------
+
+#ifdef BINARY
+
+import Data.Binary
+import Control.Monad(replicateM)
+
+-- a 64K cache, with a Double taking 13 bytes in Bytestring,
+-- implies a chunk size of 5041
+chunk :: Int
+chunk = 5000
+
+chunks :: Int -> [Int]
+chunks d = let c = d `div` chunk
+               m = d `mod` chunk
+           in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk)  
+
+putVector v = do
+              let d = dim v
+              mapM_ (\i -> put $ v @> i) [0..(d-1)]
+
+getVector d = do
+              xs <- replicateM d get
+              return $! fromList xs
+
+instance (Binary a, Storable a) => Binary (Vector a) where
+    put v = do
+            let d = dim v
+            put d
+            mapM_ putVector $! takesV (chunks d) v
+    get = do
+          d <- get
+          vs <- mapM getVector $ chunks d
+          return $! vjoin vs
+
+#endif
+
+-------------------------------------------------------------------
+
+{- | creates a Vector of the specified length using the supplied function to
+     to map the index to the value at that index.
+
+@> buildVector 4 fromIntegral
+4 |> [0.0,1.0,2.0,3.0]@
+
+-}
+buildVector :: Storable a => Int -> (Int -> a) -> Vector a
+buildVector len f =
+    fromList $ map f [0 .. (len - 1)]
+
+
+-- | zip for Vectors
+zipVector :: (Storable a, Storable b, Storable (a,b)) => Vector a -> Vector b -> Vector (a,b)
+zipVector = zipVectorWith (,)
+
+-- | unzip for Vectors
+unzipVector :: (Storable a, Storable b, Storable (a,b)) => Vector (a,b) -> (Vector a,Vector b)
+unzipVector = unzipVectorWith id
+
+-------------------------------------------------------------------
+
+{-# DEPRECATED join "use vjoin or Data.Vector.concat" #-}
+join ::  Storable t => [Vector t] -> Vector t
+join = vjoin
+