packages feed

matrices 0.3.2 → 0.4.0

raw patch · 16 files changed

+838/−571 lines, 16 filesdep −vector-binary-instancesdep ~basedep ~vector

Dependencies removed: vector-binary-instances

Dependency ranges changed: base, vector

Files

matrices.cabal view
@@ -2,17 +2,17 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                matrices-version:             0.3.2+version:             0.4.0 synopsis:            native matrix based on vector description:         This library provide the APIs for creating, indexing,-                     modifying matrices (2d arrays). The underling data-                     structure is vector. It's not intended to be a linear-                     algebra library.+                     modifying matrices (2d arrays), including dense and+                     sparse representations. The underling data+                     structure is vectors. license:             BSD3 license-file:        LICENSE author:              Kai Zhang maintainer:          kai@kzhang.org-copyright:           (c) 2014 Kai Zhang+copyright:           (c) 2015 Kai Zhang category:            Data build-type:          Simple -- extra-source-files:@@ -20,16 +20,19 @@  library   exposed-modules:-      Data.Matrix-    , Data.Matrix.Mutable-    , Data.Matrix.Generic-    , Data.Matrix.Generic.Mutable-    , Data.Matrix.Generic.Base-    , Data.Matrix.Generic.Types-    , Data.Matrix.Storable-    , Data.Matrix.Storable.Mutable-    , Data.Matrix.Unboxed-    , Data.Matrix.Unboxed.Mutable+    Data.Matrix+    Data.Matrix.Mutable+    Data.Matrix.Storable+    Data.Matrix.Storable.Mutable+    Data.Matrix.Unboxed+    Data.Matrix.Unboxed.Mutable+    Data.Matrix.Generic+    Data.Matrix.Generic.Mutable+    Data.Matrix.Dense.Generic+    Data.Matrix.Dense.Generic.Mutable+    Data.Matrix.Sparse.Generic+    Data.Matrix.Symmetric+    Data.Matrix.Symmetric.Mutable    ghc-options: -Wall -fwarn-tabs -funbox-strict-fields @@ -39,7 +42,6 @@       base >=4.0 && <5     , binary     , vector >=0.9-    , vector-binary-instances     , primitive    hs-source-dirs:      src
src/Data/Matrix.hs view
@@ -1,10 +1,10 @@ module Data.Matrix     ( Matrix-    , module Data.Matrix.Generic.Base-    ) where+    , module Data.Matrix.Dense.Generic+    )where -import qualified Data.Matrix.Generic.Types as MG-import Data.Matrix.Generic.Base+import qualified Data.Matrix.Dense.Generic as MG+import Data.Matrix.Dense.Generic hiding (Matrix) import qualified Data.Vector as V -type Matrix a = MG.Matrix V.Vector a+type Matrix = MG.Matrix V.Vector
+ src/Data/Matrix/Dense/Generic.hs view
@@ -0,0 +1,280 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies #-}+module Data.Matrix.Dense.Generic+    ( +    -- * Immutable Matrix+      Matrix(..)++    -- * Accessors+    -- ** length information+    , MG.dim+    , MG.rows+    , MG.cols++    -- ** Indexing+    , MG.unsafeIndex+    , (MG.!)+    , MG.takeRow+    , MG.takeColumn+    , MG.takeDiag++    -- * Construction+    , MG.unsafeFromVector+    , MG.fromVector+    , MG.matrix+    , MG.fromLists+    , MG.fromRows+    , fromColumns+    , MG.empty++    -- * Conversions+    , MG.flatten+    , MG.toRows+    , MG.toColumns+    , MG.toList+    , MG.toLists++    -- * Different matrix types+    , convert++    , tr+    , subMatrix+    , ident+    , diag+    , diagRect+    , fromBlocks+    , isSymmetric+    , force+    +    , Data.Matrix.Dense.Generic.foldl++    -- * Mapping+    , imap+    , Data.Matrix.Dense.Generic.map+    -- * Monadic mapping+    , mapM+    , mapM_+    , forM+    , forM_++    -- * Monadic sequencing+    , Data.Matrix.Dense.Generic.sequence+    , Data.Matrix.Dense.Generic.sequence_++    , generate++    -- * Mutable matrix+    , MG.thaw+    , MG.unsafeThaw+    , MG.freeze+    , MG.unsafeFreeze+    , MG.create+    ) where++import Prelude hiding (mapM_, mapM)+import qualified Data.Vector.Generic as G+import Control.Arrow ((***), (&&&))+import Control.Monad (liftM, foldM, foldM_)+import qualified Data.Foldable as F+import qualified Data.Vector.Generic.Mutable as GM++import qualified Data.Matrix.Generic as MG+import Data.Matrix.Dense.Generic.Mutable (MMatrix(..))++type instance MG.Mutable Matrix = MMatrix++-- | row-major matrix supporting efficient slice+data Matrix v a = Matrix !Int    -- number of rows+                         !Int    -- number of cols+                         !Int    -- physical row dimension+                         !Int    -- offset+                         !(v a)  -- flat matrix+    deriving (Show)+++instance G.Vector v a => MG.Matrix Matrix v a where+    -- | O(1) Return the size of matrix.+    dim (Matrix r c _ _ _) = (r,c)+    {-# INLINE dim #-}++    -- | O(1) Unsafe indexing without bound check.+    unsafeIndex (Matrix _ _ tda offset vec) (i,j) = vec `G.unsafeIndex` idx+      where+        idx = offset + i * tda + j+    {-# INLINE unsafeIndex #-}++    -- | O(1) Create matrix from vector.+    unsafeFromVector (r,c) = Matrix r c c 0+    {-# INLINE unsafeFromVector #-}++    -- | O(1) Extract a row.+    takeRow (Matrix _ c tda offset vec) i = G.slice i' c vec+      where+        i' = offset + i * tda+    {-# INLINE takeRow #-}++    -- | Create a vector by concatenating rows.+    flatten (Matrix r c tda offset vec)+        | c == tda = G.slice offset (r*c) vec+        | otherwise = G.generate (r*c) $ \i ->+            vec `G.unsafeIndex` (offset + (i `div` c) * tda + (i `mod` c))+    {-# INLINE flatten #-}++    thaw (Matrix r c tda offset v) = MMatrix r c tda offset `liftM` G.thaw v+    {-# INLINE thaw #-}++    unsafeThaw (Matrix r c tda offset v) = MMatrix r c tda offset `liftM` G.unsafeThaw v+    {-# INLINE unsafeThaw #-}++    freeze (MMatrix r c tda offset v) = Matrix r c tda offset `liftM` G.freeze v+    {-# INLINE freeze #-}++    unsafeFreeze (MMatrix r c tda offset v) = Matrix r c tda offset `liftM` G.unsafeFreeze v+    {-# INLINE unsafeFreeze #-}++--reshape :: G.Vector v a => Matrix v a -> (Int, Int) -> Matrix v a++-- | O(m*n) Create matrix from columns+fromColumns :: G.Vector v a => [v a] -> Matrix v a+fromColumns = tr . MG.fromRows+{-# INLINE fromColumns #-}++---- | construct upper triangular matrix from vector+--upperTriangular :: (Num a, G.Vector v a) => Int -> v a -> Matrix v a+--upperTriangular n vec =++-- | O(m*n) Convert different matrix type+convert :: (G.Vector v a, G.Vector w a) => Matrix v a -> Matrix w a+convert (Matrix r c tda offset vec) = Matrix r c tda offset . G.convert $ vec+{-# INLINE convert #-}++-- | O(1) Extract sub matrix+subMatrix :: G.Vector v a+          => (Int, Int)  -- ^ upper left corner of the submatrix+          -> (Int, Int)  -- ^ bottom right corner of the submatrix+          -> Matrix v a -> Matrix v a+subMatrix (i,j) (i',j') (Matrix _ n tda offset vec)+    | m' <= 0 || n' <= 0 = MG.empty+    | otherwise = Matrix m' n' tda offset' vec+  where+    m' = i' - i + 1+    n' = j' - j + 1+    offset' = offset + i * n + j+{-# INLINE subMatrix #-}++-- | O(m*n) Matrix transpose+tr :: G.Vector v a => Matrix v a -> Matrix v a+tr (Matrix r c tda offset vec) = MG.fromVector (c,r) $ G.generate (r*c) f+  where+    f i = vec G.! (offset + i `mod` r * tda + i `div` r)+{-# INLINE tr #-}++-- | O(m*n) Create an identity matrix+ident :: (Num a, G.Vector v a) => Int -> Matrix v a+ident n = diagRect 0 (n,n) $ replicate n 1+{-# INLINE ident #-}++-- | O(m*n) Create a square matrix with given diagonal, other entries default to 0+diag :: (Num a, G.Vector v a, F.Foldable t)+     => t a  -- ^ diagonal+     -> Matrix v a+diag d = diagRect 0 (n,n) d+  where n = length . F.toList $ d+{-# INLINE diag #-}++-- | O(m*n) Create a rectangular matrix with default values and given diagonal+diagRect :: (G.Vector v a, F.Foldable t)+         => a         -- ^ default value+         -> (Int, Int)+         -> t a       -- ^ diagonal+         -> Matrix v a+diagRect z0 (r,c) d = MG.fromVector (r,c) $ G.create $ GM.replicate n z0 >>= go d c+  where+    go xs c' v = F.foldlM f 0 xs >> return v+      where+        f !i x = GM.unsafeWrite v (i*(c'+1)) x >> return (i+1)+    n = r * c+{-# INLINE diagRect #-}++fromBlocks :: G.Vector v a+           => a               -- ^ default value+           -> [[Matrix v a]]+           -> Matrix v a+fromBlocks d ms = MG.fromVector (m,n) $ G.create $ GM.replicate (m*n) d >>= go n ms+  where+    go n' xss v = foldM_ f 0 xss >> return v+      where+        f !cr xs = do (r', _) <- foldM g (0, 0) xs+                      return $ cr + r'+          where+            g (!maxR, !cc) x = do+                let (r,c) = MG.dim x+                    vec = MG.flatten x+                    step i u = do+                        GM.unsafeWrite v ((cr + i `div` c) * n' + i `mod` c + cc) u+                        return (i+1)+                G.foldM'_ step (0::Int) vec+                return (max maxR r, cc + c)+    -- figure out the dimension of the new matrix+    (m, n) = (sum *** maximum) . unzip . Prelude.map ((maximum *** sum) .+                unzip . Prelude.map (MG.rows &&& MG.cols)) $ ms+{-# INLINE fromBlocks #-}++isSymmetric :: (Eq a, G.Vector v a) => Matrix v a -> Bool+isSymmetric m@(Matrix r c _ _ _) | r /= c = False+                                 | otherwise = all f [0 .. r-1]+  where+    f i = all g [i + 1 .. c-1]+      where g j = m MG.! (i,j) == m MG.! (j,i)+{-# INLINE isSymmetric #-}++force :: G.Vector v a => Matrix v a -> Matrix v a+force m@(Matrix r c _ _ _) = MG.fromVector (r,c) . G.force . MG.flatten $ m+{-# INLINE force #-}++imap :: (G.Vector v a, G.Vector v b) => ((Int, Int) -> a -> b) -> Matrix v a -> Matrix v b+imap f m@(Matrix r c _ _ _) = MG.fromVector (r,c) $ G.imap f' . MG.flatten $ m+  where+    f' i = f (i `div` c, i `mod` c)+{-# INLINE imap #-}++map :: (G.Vector v a, G.Vector v b) => (a -> b) -> Matrix v a -> Matrix v b+map f m@(Matrix r c _ _ _) = MG.fromVector (r,c) $ G.map f . MG.flatten $ m+{-# INLINE map #-}++foldl :: G.Vector v b => (a -> b -> a) -> a -> Matrix v b -> a+foldl f acc m = G.foldl f acc . MG.flatten $ m+{-# INLINE foldl #-}++mapM :: (G.Vector v a, G.Vector v b, Monad m) => (a -> m b) -> Matrix v a -> m (Matrix v b)+mapM f m@(Matrix r c _ _ _) = liftM (MG.fromVector (r,c)) . G.mapM f . MG.flatten $ m+{-# INLINE mapM #-}++mapM_ :: (G.Vector v a, Monad m) => (a -> m b) -> Matrix v a -> m ()+mapM_ f = G.mapM_ f . MG.flatten+{-# INLINE mapM_ #-}++forM :: (G.Vector v a, G.Vector v b, Monad m) => Matrix v a -> (a -> m b) -> m (Matrix v b)+forM = flip mapM+{-# INLINE forM #-}++forM_ :: (G.Vector v a, Monad m) => Matrix v a -> (a -> m b) -> m ()+forM_ = flip mapM_+{-# INLINE forM_ #-}++sequence :: (G.Vector v a, G.Vector v (m a), Monad m)+         => Matrix v (m a) -> m (Matrix v a)+sequence (Matrix r c tda offset vec) = liftM (Matrix r c tda offset) . G.sequence $ vec+{-# INLINE sequence #-}++sequence_ :: (G.Vector v (m a), Monad m)+          => Matrix v (m a) -> m ()+sequence_ (Matrix _ _ _ _ vec) = G.sequence_ vec+{-# INLINE sequence_ #-}++generate :: G.Vector v a => (Int, Int) -> ((Int, Int) -> a) -> Matrix v a+generate (r,c) f = MG.fromVector (r,c) . G.generate (r*c) $ \i -> f (i `div` c, i `mod` c)+{-# INLINE generate #-}
+ src/Data/Matrix/Dense/Generic/Mutable.hs view
@@ -0,0 +1,47 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module Data.Matrix.Dense.Generic.Mutable+   ( -- * Mutable Matrix+     MMatrix(..)+   , C.dim+   , takeRow+   , C.write+   , C.unsafeWrite+   , C.read+   , C.unsafeRead+   , C.new+   , C.replicate+   ) where++import Prelude hiding (read, replicate)+import Control.Monad (liftM)+import qualified Data.Vector.Generic.Mutable as GM++import qualified Data.Matrix.Generic.Mutable as C++-- | mutable matrix+data MMatrix v s a = MMatrix !Int !Int !Int !Int !(v s a)++instance GM.MVector v a => C.MMatrix MMatrix v a where+    dim (MMatrix r c _ _ _) = (r,c)+    {-# INLINE dim #-}++    unsafeRead (MMatrix _ _ tda offset v) (i,j) = GM.unsafeRead v idx+      where idx = offset + i * tda + j+    {-# INLINE unsafeRead #-}++    unsafeWrite (MMatrix _ _ tda offset v) (i,j) = GM.unsafeWrite v idx+      where idx = offset + i * tda + j+    {-# INLINE unsafeWrite #-}++    new (r,c) = MMatrix r c c 0 `liftM` GM.new (r*c)+    {-# INLINE new #-}++    replicate (r,c) x = MMatrix r c c 0 `liftM` GM.replicate (r*c) x+    {-# INLINE replicate #-}+    +takeRow :: GM.MVector v a => MMatrix v m a -> Int -> v m a+takeRow (MMatrix _ c tda offset vec) i = GM.slice i' c vec+  where+    i' = offset + i * tda+{-# INLINE takeRow #-}
src/Data/Matrix/Generic.hs view
@@ -1,20 +1,168 @@------------------------------------------------------------------------------------ |--- Module      :  $Header$--- Copyright   :  (c) 2014 Kai Zhang--- License     :  BSD3+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE Rank2Types #-}+module Data.Matrix.Generic+    ( Mutable+    , Matrix(..) --- Maintainer  :  kai@kzhang.org--- Stability   :  experimental--- Portability :  portable+    -- * Derived mothods+    , rows+    , cols+    , (!)+    , fromVector+    , empty+    , toList+    , fromLists+    , matrix+    , fromRows+    , toRows+    , toColumns+    , toLists+    , create+    ) where ------------------------------------------------------------------------------------+import Control.Monad.ST (runST, ST)+import Control.Monad.Primitive (PrimMonad, PrimState)+import qualified Data.Vector.Generic as G -module Data.Matrix.Generic-    ( module Data.Matrix.Generic.Base-    , module Data.Matrix.Generic.Types-    ) where+import qualified Data.Matrix.Generic.Mutable as MM -import Data.Matrix.Generic.Base-import Data.Matrix.Generic.Types+type family Mutable (m :: (* -> *) -> * -> *) :: (* -> * -> *) -> * -> * -> *++class (MM.MMatrix (Mutable m) (G.Mutable v) a, G.Vector v a) => Matrix m v a where+    dim :: m v a -> (Int, Int)++    unsafeIndex :: m v a -> (Int, Int) -> a++    unsafeFromVector :: (Int, Int) -> v a -> m v a++    -- | Default algorithm is O((m*n) * O(unsafeIndex)).+    flatten :: m v a -> v a+    flatten mat = G.generate (r*c) $ \i -> unsafeIndex mat (i `div` c, i `mod` c)+      where+        (r,c) = dim mat+    {-# INLINE flatten #-}++    -- | Extract a row. Default algorithm is O(n * O(unsafeIndex)).+    takeRow :: m v a -> Int -> v a+    takeRow mat i = G.generate c $ \j -> unsafeIndex mat (i,j)+      where+        (_,c) = dim mat+    {-# INLINE takeRow #-}++    -- | Extract a column. Default algorithm is O(m * O(unsafeIndex)).+    takeColumn :: m v a -> Int -> v a+    takeColumn mat j = G.generate r $ \i -> unsafeIndex mat (i,j)+      where+        (r,_) = dim mat+    {-# INLINE takeColumn #-}++    -- | Extract the diagonal. Default algorithm is O(min(m,n) * O(unsafeIndex)).+    takeDiag :: m v a -> v a+    takeDiag mat = G.generate n $ \i -> unsafeIndex mat (i,i)+      where+        n = uncurry min . dim $ mat+    {-# INLINE takeDiag #-}++    thaw :: PrimMonad s => m v a -> s ((Mutable m) (G.Mutable v) (PrimState s) a)++    unsafeThaw :: PrimMonad s+               => m v a -> s ((Mutable m) (G.Mutable v) (PrimState s) a)++    freeze :: PrimMonad s+           => (Mutable m) (G.Mutable v) (PrimState s) a -> s (m v a)++    unsafeFreeze :: PrimMonad s+                 => (Mutable m) (G.Mutable v) (PrimState s) a -> s (m v a)++    {-# MINIMAL dim, unsafeIndex, unsafeFromVector, thaw, unsafeThaw, freeze, unsafeFreeze #-}++-- | Derived methods++-- | Return the number of rows+rows :: Matrix m v a => m v a -> Int+rows = fst . dim+{-# INLINE rows #-}++-- | Return the number of columns+cols :: Matrix m v a => m v a -> Int+cols = snd . dim+{-# INLINE cols #-}++-- | Indexing+(!) :: Matrix m v a => m v a -> (Int, Int) -> a+(!) mat (i,j) | i >= r || j >= c = error "Index out of bounds"+              | otherwise = unsafeIndex mat (i,j)+  where+    (r,c) = dim mat+{-# INLINE (!) #-}++-- | O(m*n) Create a list by concatenating rows+toList :: Matrix m v a => m v a -> [a]+toList = G.toList . flatten+{-# INLINE toList #-}++empty :: Matrix m v a => m v a+empty = fromVector (0,0) G.empty+{-# INLINE empty #-}++fromVector :: Matrix m v a => (Int, Int) -> v a -> m v a+fromVector (r,c) vec | r*c /= G.length vec = error "incorrect length"+                     | otherwise = unsafeFromVector (r,c) vec+{-# INLINE fromVector #-}++-- | O(m*n) Matrix construction+matrix :: Matrix m v a+       => Int  -- ^ number of columns+       -> [a]  -- ^ row list+       -> m v a+matrix ncol xs | n `mod` ncol /= 0 = error "incorrect length"+               | otherwise = unsafeFromVector (nrow,ncol) vec+  where+    vec = G.fromList xs+    nrow = n `div` ncol+    n = G.length vec+{-# INLINE matrix #-}++-- | O(m*n) Create matrix from list of lists, it doesn't check if the list of+-- list is a valid matrix+fromLists :: Matrix m v a => [[a]] -> m v a+fromLists xs | null xs = empty+             | otherwise = fromVector (r,c) . G.fromList . concat $ xs+  where+    r = length xs+    c = length . head $ xs+{-# INLINE fromLists #-}++-- | O(m*n) Create matrix from rows+fromRows :: Matrix m v a => [v a] -> m v a+fromRows xs | null xs = empty+            | otherwise = fromVector (r,c) . G.concat $ xs+  where+    r = length xs+    c = G.length . head $ xs+{-# INLINE fromRows #-}++-- | O(m) Return the rows+toRows :: Matrix m v a => m v a -> [v a]+toRows mat = map (takeRow mat) [0..r-1]+  where+    (r,_) = dim mat+{-# INLINE toRows #-}++-- | O(m*n) Return the columns+toColumns :: Matrix m v a => m v a -> [v a]+toColumns mat = map (takeColumn mat) [0..c-1]+  where+    (_,c) = dim mat+{-# INLINE toColumns #-}++-- | O(m*n) List of lists+toLists :: Matrix m v a => m v a -> [[a]]+toLists = map G.toList . toRows+{-# INLINE toLists #-}++create :: Matrix m v a => (forall s . ST s ((Mutable m) (G.Mutable v) s a)) -> m v a+create m = runST $ unsafeFreeze =<< m+{-# INLINE create #-}
− src/Data/Matrix/Generic/Base.hs 
@@ -1,331 +0,0 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE FlexibleContexts #-}------------------------------------------------------------------------------------ |--- Module      :  $Header$--- Copyright   :  (c) 2014 Kai Zhang--- License     :  BSD3---- Maintainer  :  kai@kzhang.org--- Stability   :  experimental--- Portability :  portable--------------------------------------------------------------------------------------module Data.Matrix.Generic.Base-    ( rows-    , cols-    , dim-    , (!)-    , unsafeIndex-    , empty--    -- * Conversions-    , matrix-    , flatten-    , fromVector-    , toRows-    , toColumns-    , fromRows-    , fromColumns-    , toList-    , toLists-    , fromLists--    -- * Different matrix types-    , convert--    , tr-    , takeRow-    , takeColumn-    , subMatrix-    , ident-    , diag-    , diagRect-    , takeDiag-    , fromBlocks-    , isSymmetric-    , force-    , Data.Matrix.Generic.Base.foldl--    -- * Mapping-    , imap-    , Data.Matrix.Generic.Base.map--    -- * Monadic mapping-    , Data.Matrix.Generic.Base.mapM-    , Data.Matrix.Generic.Base.mapM_-    , Data.Matrix.Generic.Base.forM-    , Data.Matrix.Generic.Base.forM_--    , generate--    ) where--import Control.Arrow ((***), (&&&))-import Control.Monad-import qualified Data.Foldable as F-import qualified Data.Vector.Generic as G-import qualified Data.Vector.Generic.Mutable as GM-import Data.Matrix.Generic.Types---- | O(1) Return the number of rows-rows :: G.Vector v a => Matrix v a -> Int-rows (Matrix m _ _ _ _) = m---- | O(1) Return the number of columns-cols :: G.Vector v a => Matrix v a -> Int-cols (Matrix _ n _ _ _) = n---- | O(1) Return the size of matrix-dim :: G.Vector v a => Matrix v a -> (Int, Int)-dim (Matrix r c _ _ _) = (r,c)-{-# INLINE dim #-}---- | O(1) Indexing-(!) :: G.Vector v a => Matrix v a -> (Int, Int) -> a-(!) (Matrix _ _ tda offset vec) (i, j) = vec G.! idx-  where-    idx = offset + i * tda + j-{-# INLINE (!) #-}---- | O(1) Unsafe indexing without bound check-unsafeIndex :: G.Vector v a => Matrix v a -> (Int, Int) -> a-unsafeIndex (Matrix _ _ tda offset vec) (i,j) = vec `G.unsafeIndex` idx-  where-    idx = offset + i * tda + j-{-# INLINE unsafeIndex #-}----reshape :: G.Vector v a => Matrix v a -> (Int, Int) -> Matrix v a---- | O(1) Return an empty matrix-empty :: G.Vector v a => Matrix v a-empty = Matrix 0 0 0 0 G.empty-{-# INLINE empty #-}---- | O(m*n) Matrix construction-matrix :: G.Vector v a-       => Int  -- ^ number of columns-       -> [a]  -- ^ row list-       -> Matrix v a-matrix ncol xs | n `mod` ncol /= 0 = error "incorrect length"-               | otherwise = fromVector (nrow,ncol) vec-  where-    vec = G.fromList xs-    nrow = n `div` ncol-    n = G.length vec-{-# INLINE matrix #-}---- | Create a vector by concatenating rows-flatten :: G.Vector v a => Matrix v a -> v a-flatten (Matrix m n tda offset vec)-    | n == tda = G.slice offset (m * n) vec-    | otherwise = G.generate (m * n) f-  where-    f i = (G.!) vec $ offset + (i `div` n) * tda + i `mod` n-{-# INLINE flatten #-}---- | O(1) Create matrix from vector-fromVector :: G.Vector v a => (Int, Int) -> v a -> Matrix v a-fromVector (r,c) = Matrix r c c 0-{-# INLINE fromVector #-}---- | O(m*n) Create a list by concatenating rows-toList :: G.Vector v a => Matrix v a -> [a]-toList = G.toList . flatten-{-# INLINE toList #-}---- | O(m) Return the rows-toRows :: G.Vector v a => Matrix v a -> [v a]-toRows (Matrix m n tda offset vec) = loop 0-  where-    loop !i | i < m = G.slice (f i) n vec : loop (i+1)-            | otherwise = []-    f i = offset + i * tda-{-# INLINE toRows #-}---- | O(m*n) Return the columns-toColumns :: G.Vector v a => Matrix v a -> [v a]-toColumns m = Prelude.map (takeColumn m) [0 .. c-1]-  where c = cols m-{-# INLINE toColumns #-}---- | O(m*n) Create matrix from rows-fromRows :: G.Vector v a => [v a] -> Matrix v a-fromRows xs | any (\x -> G.length x /= c) xs = error "inequal length"-            | otherwise = fromVector (r,c) . G.concat $ xs-  where-    r = length xs-    c = G.length . head $ xs-{-# INLINE fromRows #-}---- | O(m*n) Create matrix from columns-fromColumns :: G.Vector v a => [v a] -> Matrix v a-fromColumns = tr . fromRows-{-# INLINE fromColumns #-}---- | O(m*n) List of lists-toLists :: G.Vector v a => Matrix v a -> [[a]]-toLists = Prelude.map G.toList . toRows-{-# INLINE toLists #-}---- | O(m*n) Create matrix from list of lists, it doesn't check if the list of--- list is a valid matrix-fromLists :: G.Vector v a => [[a]] -> Matrix v a-fromLists xs = fromVector (r,c) . G.fromList . concat $ xs-  where-    r = length xs-    c = length .head $ xs-{-# INLINE fromLists #-}------ | construct upper triangular matrix from vector---upperTriangular :: (Num a, G.Vector v a) => Int -> v a -> Matrix v a---upperTriangular n vec =---- | O(m*n) Convert different matrix type-convert :: (G.Vector v a, G.Vector w a) => Matrix v a -> Matrix w a-convert (Matrix r c tda offset vec) = Matrix r c tda offset . G.convert $ vec-{-# INLINE convert #-}---- | O(1) Extract a row-takeRow :: G.Vector v a => Matrix v a -> Int -> v a-takeRow (Matrix _ c tda offset vec) i = G.slice i' c vec-  where-    i' = offset + i * tda-{-# INLINE takeRow #-}---- | O(m) Extract a column-takeColumn :: G.Vector v a => Matrix v a -> Int -> v a-takeColumn (Matrix r _ tda offset vec) j = G.create $ GM.new r >>= go idx vec r 0-  where-    go f vec' r' !i v | i >= r' = return v-                      | otherwise = do GM.unsafeWrite v i $ vec' G.! f i-                                       go f vec' r' (i+1) v-    idx i = offset + i * tda + j-{-# INLINE takeColumn #-}---- | O(1) Extract sub matrix-subMatrix :: G.Vector v a-          => (Int, Int)  -- ^ upper left corner of the submatrix-          -> (Int, Int)  -- ^ bottom right corner of the submatrix-          -> Matrix v a -> Matrix v a-subMatrix (i,j) (i',j') (Matrix _ n tda offset vec)-    | m' <= 0 || n' <= 0 = empty-    | otherwise = Matrix m' n' tda offset' vec-  where-    m' = i' - i + 1-    n' = j' - j + 1-    offset' = offset + i * n + j-{-# INLINE subMatrix #-}---- | O(m*n) Matrix transpose-tr :: G.Vector v a => Matrix v a -> Matrix v a-tr (Matrix r c tda offset vec) = fromVector (c,r) $ G.generate (r*c) f-  where-    f i = vec G.! (offset + i `mod` r * tda + i `div` r)-{-# INLINE tr #-}---- | O(m*n) Create an identity matrix-ident :: (Num a, G.Vector v a) => Int -> Matrix v a-ident n = diagRect 0 (n,n) $ replicate n 1-{-# INLINE ident #-}---- | O(m*n) Create a square matrix with given diagonal, other entries default to 0-diag :: (Num a, G.Vector v a, F.Foldable t)-     => t a  -- ^ diagonal-     -> Matrix v a-diag d = diagRect 0 (n,n) d-  where n = length . F.toList $ d-{-# INLINE diag #-}---- | O(m*n) Create a rectangular matrix with default values and given diagonal-diagRect :: (G.Vector v a, F.Foldable t)-         => a         -- ^ default value-         -> (Int, Int)-         -> t a       -- ^ diagonal-         -> Matrix v a-diagRect z0 (r,c) d = fromVector (r,c) $ G.create $ GM.replicate n z0 >>= go d c-  where-    go xs c' v = F.foldlM f 0 xs >> return v-      where-        f !i x = GM.unsafeWrite v (i*(c'+1)) x >> return (i+1)-    n = r * c-{-# INLINE diagRect #-}---- | extracts the diagonal from a rectangular matrix-takeDiag :: G.Vector v a => Matrix v a -> v a-takeDiag mat@(Matrix r c _ _ _) = G.generate n $ \i -> unsafeIndex mat (i,i)-  where-    n = min r c-{-# INLINE takeDiag #-}--fromBlocks :: G.Vector v a-           => a               -- ^ default value-           -> [[Matrix v a]]-           -> Matrix v a-fromBlocks d ms = fromVector (m,n) $ G.create $ GM.replicate (m*n) d >>= go n ms-  where-    go n' xss v = foldM_ f 0 xss >> return v-      where-        f !cr xs = do (r', _) <- foldM g (0, 0) xs-                      return $ cr + r'-          where-            g (!maxR, !cc) x = do-                let c = cols x-                    r = rows x-                    vec = flatten x-                    step i u = do-                        GM.unsafeWrite v ((cr + i `div` c) * n' + i `mod` c + cc) u-                        return (i+1)-                G.foldM'_ step (0::Int) vec-                return (max maxR r, cc + c)-    -- figure out the dimension of the new matrix-    (m, n) = (sum *** maximum) . unzip . Prelude.map ((maximum *** sum) .-                unzip . Prelude.map (rows &&& cols)) $ ms-{-# INLINE fromBlocks #-}--isSymmetric :: (Eq a, G.Vector v a) => Matrix v a -> Bool-isSymmetric m@(Matrix r c _ _ _) | r /= c = False-                                 | otherwise = all f [0 .. r-1]-  where-    f i = all g [i + 1 .. c-1]-      where g j = m ! (i,j) == m ! (j,i)-{-# INLINE isSymmetric #-}--force :: G.Vector v a => Matrix v a -> Matrix v a-force m@(Matrix r c _ _ _) = fromVector (r,c) . G.force . flatten $ m-{-# INLINE force #-}--imap :: (G.Vector v a, G.Vector v b) => ((Int, Int) -> a -> b) -> Matrix v a -> Matrix v b-imap f m@(Matrix r c _ _ _) = fromVector (r,c) $ G.imap f' . flatten $ m-  where-    f' i = f (i `div` c, i `mod` c)-{-# INLINE imap #-}--map :: (G.Vector v a, G.Vector v b) => (a -> b) -> Matrix v a -> Matrix v b-map f m@(Matrix r c _ _ _) = fromVector (r,c) $ G.map f . flatten $ m-{-# INLINE map #-}--foldl :: G.Vector v b => (a -> b -> a) -> a -> Matrix v b -> a-foldl f acc m = G.foldl f acc . flatten $ m-{-# INLINE foldl #-}--mapM :: (G.Vector v a, G.Vector v b, Monad m) => (a -> m b) -> Matrix v a -> m (Matrix v b)-mapM f m@(Matrix r c _ _ _) = liftM (fromVector (r,c)) . G.mapM f . flatten $ m-{-# INLINE mapM #-}--mapM_ :: (G.Vector v a, Monad m) => (a -> m b) -> Matrix v a -> m ()-mapM_ f = G.mapM_ f . flatten-{-# INLINE mapM_ #-}--forM :: (G.Vector v a, G.Vector v b, Monad m) => Matrix v a -> (a -> m b) -> m (Matrix v b)-forM = flip Data.Matrix.Generic.Base.mapM-{-# INLINE forM #-}--forM_ :: (G.Vector v a, Monad m) => Matrix v a -> (a -> m b) -> m ()-forM_ = flip Data.Matrix.Generic.Base.mapM_-{-# INLINE forM_ #-}--generate :: G.Vector v a => (Int, Int) -> ((Int, Int) -> a) -> Matrix v a-generate (r,c) f = fromVector (r,c) . G.generate (r*c) $ \i -> f (i `div` c, i `mod` c)-{-# INLINE generate #-}
src/Data/Matrix/Generic/Mutable.hs view
@@ -1,112 +1,51 @@-{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies #-}+ module Data.Matrix.Generic.Mutable-   ( fromMVector-   , dim-   , flatten-   , takeRow-   , thaw-   , unsafeThaw-   , freeze-   , unsafeFreeze-   , write-   , unsafeWrite-   , read-   , unsafeRead-   , replicate-   , new-   , create-   ) where+    ( MMatrix(..)+    , write+    , read+    ) where -import Prelude hiding (read, replicate)-import Control.Monad-import Control.Monad.ST-import Data.Matrix.Generic.Types-import qualified Data.Vector.Generic as G+import Prelude hiding (read)+import Control.Monad.Primitive (PrimMonad, PrimState) import qualified Data.Vector.Generic.Mutable as GM-import Control.Monad.Primitive --- to be removed in GHC-7.10-(<$>) :: Monad m => (a -> b) -> m a -> m b-(<$>) = liftM--dim :: MMatrix v m a -> (Int, Int)-dim (MMatrix r c _ _ _) = (r,c)-{-# INLINE dim #-}--fromMVector :: GM.MVector v a => (Int, Int) -> v m a -> MMatrix v m a-fromMVector (r,c) = MMatrix r c c 0-{-# INLINE fromMVector #-}+class GM.MVector v a => MMatrix m v a where+    dim ::  m v s a -> (Int, Int) -flatten :: (GM.MVector v a, PrimMonad m)-        => MMatrix v (PrimState m) a -> m (v (PrimState m) a)-flatten (MMatrix m n tda offset vec)-    | n == tda = return $ GM.slice offset (m * n) vec-    | otherwise = do-        vec' <- GM.new (m*n)-        forM_ [0 .. m*n] $ \i ->-            GM.unsafeRead vec (offset + (i `div` n) * tda + i `mod` n) >>=-                GM.unsafeWrite vec' i-        return vec'-{-# INLINE flatten #-}+    unsafeRead :: PrimMonad s => m v (PrimState s) a -> (Int, Int) -> s a -takeRow :: GM.MVector v a => MMatrix v m a -> Int -> v m a-takeRow (MMatrix _ c tda offset vec) i = GM.slice i' c vec-  where-    i' = offset + i * tda-{-# INLINE takeRow #-}+    unsafeWrite :: PrimMonad s => m v (PrimState s) a -> (Int, Int) -> a -> s () -thaw :: (G.Vector v a, PrimMonad m)-     => Matrix v a -> m (MMatrix (G.Mutable v) (PrimState m) a)-thaw (Matrix r c tda offset v) = MMatrix r c tda offset <$> G.thaw v-{-# INLINE thaw #-}+    -- | Create a mutable matrix without initialization+    new :: PrimMonad s => (Int, Int) -> s (m v (PrimState s) a) -unsafeThaw :: (G.Vector v a, PrimMonad m)-           => Matrix v a -> m (MMatrix (G.Mutable v) (PrimState m) a)-unsafeThaw (Matrix r c tda offset v) = MMatrix r c tda offset <$> G.unsafeThaw v-{-# INLINE unsafeThaw #-}+    replicate :: PrimMonad s => (Int, Int) -> a -> s (m v (PrimState s) a) -freeze :: (PrimMonad m, G.Vector v a) => MMatrix (G.Mutable v) (PrimState m) a -> m (Matrix v a)-freeze (MMatrix r c tda offset v) = Matrix r c tda offset <$> G.freeze v-{-# INLINE freeze #-}+    {-# MINIMAL dim, unsafeRead, unsafeWrite, new, replicate #-} -unsafeFreeze :: (PrimMonad m, G.Vector v a) => MMatrix (G.Mutable v) (PrimState m) a -> m (Matrix v a)-unsafeFreeze (MMatrix r c tda offset v) = Matrix r c tda offset <$> G.unsafeFreeze v-{-# INLINE unsafeFreeze #-}+-- | Derived methods -write :: (PrimMonad m, GM.MVector v a)-      => MMatrix v (PrimState m) a -> (Int, Int) -> a -> m ()-write (MMatrix _ _ tda offset v) (i,j) = GM.write v idx-  where idx = offset + i * tda + j+write :: (PrimMonad s, MMatrix m v a)+      => m v (PrimState s) a -> (Int, Int) -> a -> s ()+write mat (i,j) | i >= r || j >= c = error "Index out of bounds"+                | otherwise = unsafeWrite mat (i,j)+  where+    (r,c) = dim mat {-# INLINE write #-} -unsafeWrite :: (PrimMonad m, GM.MVector v a)-            => MMatrix v (PrimState m) a -> (Int, Int) -> a -> m ()-unsafeWrite (MMatrix _ _ tda offset v) (i,j) = GM.unsafeWrite v idx-  where idx = offset + i * tda + j-{-# INLINE unsafeWrite #-}--read :: (PrimMonad m, GM.MVector v a)-     => MMatrix v (PrimState m) a -> (Int, Int) -> m a-read (MMatrix _ _ tda offset v) (i,j) = GM.read v idx-  where idx = offset + i * tda + j+read :: (PrimMonad s, MMatrix m v a)+     => m v (PrimState s) a -> (Int, Int) -> s a+read mat (i,j) | i >= r || j >= c = error "Index out of bounds"+               | otherwise = unsafeRead mat (i,j)+  where+    (r,c) = dim mat {-# INLINE read #-} -unsafeRead :: (PrimMonad m, GM.MVector v a)-           => MMatrix v (PrimState m) a -> (Int, Int) -> m a-unsafeRead (MMatrix _ _ tda offset v) (i,j) = GM.unsafeRead v idx-  where idx = offset + i * tda + j-{-# INLINE unsafeRead #-}--replicate :: (PrimMonad m, GM.MVector v a)-          => (Int, Int) -> a -> m (MMatrix v (PrimState m) a)-replicate (r,c) x = fromMVector (r,c) <$> GM.replicate (r*c) x-{-# INLINE replicate #-}--new :: (PrimMonad m, GM.MVector v a)-    => (Int, Int) -> m (MMatrix v (PrimState m) a)-new (r,c) = fromMVector (r,c) <$> GM.new (r*c)-{-# INLINE new #-}--create :: G.Vector v a => (forall s . ST s (MMatrix (G.Mutable v) s a)) -> Matrix v a-create m = runST $ unsafeFreeze =<< m-{-# INLINE create #-}+{-+replicate :: (PrimMonad s, MMatrix m v a)+          => (Int, Int) -> a -> s (m v (PrimState s) a)+replicate (r,c) x = do+    mat <- new (r,c)+ -}
− src/Data/Matrix/Generic/Types.hs 
@@ -1,85 +0,0 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE OverlappingInstances #-}------------------------------------------------------------------------------------ |--- Module      :  $Header$--- Copyright   :  (c) 2014 Kai Zhang--- License     :  BSD3---- Maintainer  :  kai@kzhang.org--- Stability   :  experimental--- Portability :  portable---------------------------------------------------------------------------------------module Data.Matrix.Generic.Types-    ( Matrix(..)-    , MMatrix(..)-    ) where--import Data.Binary-import qualified Data.Vector as V-import qualified Data.Vector.Unboxed as U-import qualified Data.Vector.Storable as S-import qualified Data.Vector.Binary ()-import qualified Data.Vector.Generic as G---- | row-major matrix supporting efficient slice-data Matrix v a = Matrix !Int    -- number of rows-                         !Int    -- number of cols-                         !Int    -- physical row dimension-                         !Int    -- offset-                         !(v a)  -- flat matrix--instance Binary a => Binary (Matrix V.Vector a) where-    put = putGeneric-    get = getGeneric--instance (U.Unbox a, Binary a) => Binary (Matrix U.Vector a) where-    put = putGeneric-    get = getGeneric--instance (S.Storable a, Binary a) => Binary (Matrix S.Vector a) where-    put = putGeneric-    get = getGeneric--getGeneric :: (Binary (v a), G.Vector v a) => Get (Matrix v a)-getGeneric = do-    r <- get-    c <- get-    tda <- get-    offset <- get-    vec <- get-    return $ Matrix r c tda offset vec--putGeneric :: (Binary (v a), G.Vector v a) => Matrix v a -> Put-putGeneric (Matrix r c tda offset vec) = do-    put r-    put c-    put tda-    put offset-    put vec---- | mutable matrix-data MMatrix v m a = MMatrix !Int !Int !Int !Int !(v m a)--instance (G.Vector v a, Show a) => Show (Matrix v a) where-    show mat = unlines . map (unwords . map show) . toLists $ mat--instance G.Vector v Bool => Show (Matrix v Bool) where-    show mat = unlines . map (unwords . map showBool) . toLists $ mat-      where-        showBool x = if x then "1" else "0"--toRows :: G.Vector v a => Matrix v a -> [v a]-toRows (Matrix m n tda offset vec) = loop 0-    where-    loop !i | i < m = G.slice (f i) n vec : loop (i+1)-            | otherwise = []-    f i = offset + i * tda--toLists :: G.Vector v a => Matrix v a -> [[a]]-toLists = map G.toList . toRows
src/Data/Matrix/Mutable.hs view
@@ -1,10 +1,10 @@ module Data.Matrix.Mutable     ( MMatrix-    , module Data.Matrix.Generic.Mutable+    , module Data.Matrix.Dense.Generic.Mutable     ) where -import qualified Data.Matrix.Generic.Types as MG-import Data.Matrix.Generic.Mutable+import qualified Data.Matrix.Dense.Generic.Mutable as MG+import Data.Matrix.Dense.Generic.Mutable hiding (MMatrix) import qualified Data.Vector.Mutable as VM  type MMatrix a = MG.MMatrix VM.MVector a
+ src/Data/Matrix/Sparse/Generic.hs view
@@ -0,0 +1,125 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies #-}+module Data.Matrix.Sparse.Generic+    ( Zero(..)+    , CSR(..)+    , fromAscAL+    , fromAscStream+    , (!)+    ) where++import Control.Monad (when, forM_)+import Control.Monad.ST (runST)+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as GM+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Fusion.Stream as S+import Data.Bits (shiftR)++import Data.Matrix.Generic+import Data.Matrix.Dense.Generic.Mutable (MMatrix)++class Eq a => Zero a where+    zero :: a++instance Zero Int where+    zero = 0++instance Zero Double where+    zero = 0.0++instance Eq a => Zero ([] a) where+    zero = []++-- | mutable sparse matrix not implemented+type instance Mutable CSR = MMatrix++-- | Compressed Sparse Row (CSR) matrix+data CSR v a = CSR !Int  -- rows+                   !Int  -- cols+                   !(v a)  -- values+                   !(U.Vector Int)  -- column index of values+                   !(U.Vector Int)  -- row pointer+    deriving (Show)++instance (Zero a, G.Vector v a) => Matrix CSR v a where+    dim (CSR r c _ _ _) = (r,c)+    {-# INLINE dim #-}++    unsafeIndex (CSR _ _ vec ci rp) (i,j) =+        case binarySearchByBounds ci j r0 r1 of+            Nothing -> zero+            Just k -> vec `G.unsafeIndex` k+      where+        r0 = rp `U.unsafeIndex` i+        r1 = rp `U.unsafeIndex` (i+1) - 1+    {-# INLINE unsafeIndex #-}++    unsafeFromVector (r,c) vec =+        CSR r c (G.generate n (G.unsafeIndex vec . U.unsafeIndex nz))+                (U.map (`mod` c) nz)+                (U.fromList . g . U.foldr f ((r-1,n-1), [n]) $ nz)+      where+        nz = U.filter (\i -> vec `G.unsafeIndex` i /= zero) . U.enumFromN 0 $ (r*c)+        f i ((!prev,!acc), xs) | stride == 0 = ((prev, acc-1), xs)+                               | otherwise = ((current, acc-1), replicate stride (acc+1) ++ xs)+          where+            stride = prev - current+            current = i `div` c+        g ((a, _), xs) | a == 0 = 0 : xs+                       | otherwise = replicate (a+1) 0 ++ xs+        n = U.length nz+    {-# INLINE unsafeFromVector #-}++    thaw = undefined+    unsafeThaw = undefined+    freeze = undefined+    unsafeFreeze = undefined++type AssocList a = [((Int, Int), a)]++fromAscAL :: G.Vector v a => (Int, Int) -> Int -> AssocList a -> CSR v a+fromAscAL (r,c) n al = fromAscStream (r,c) n . S.fromList $ al+{-# INLINE fromAscAL #-}++fromAscStream :: (GM.MVector (G.Mutable v) a, G.Vector v a) => (Int, Int) -> Int -> S.Stream ((Int,Int),a) -> CSR v a+fromAscStream (r,c) n al = CSR r c values ci rp+  where+    (values, ci, rp) = runST $ do+        v <- GM.new n+        col <- GM.new n+        row <- GM.new (r+1)++        let f (i',acc) ((i,j),x) = do+                GM.write v acc x+                GM.write col acc j++                let stride = i - i'+                when (stride > 0) $ forM_ [0..stride-1] $ \s -> GM.write row (i-s) acc+                +                return (i,acc+1)++        _ <- S.foldM f (0,0) al+        GM.write row r n++        v' <- G.unsafeFreeze v+        col' <- G.unsafeFreeze col+        row' <- G.unsafeFreeze row+        return (v', col', row')+{-# INLINE fromAscStream #-}+    +binarySearchByBounds :: U.Vector Int -> Int -> Int -> Int -> Maybe Int+binarySearchByBounds vec x = loop+  where+    loop !l !u+        | l > u = Nothing+        | x == x' = Just k+        | x < x' = loop l (k-1)+        | otherwise = loop (k+1) u+      where+        k = (u+l) `shiftR` 1+        x' = vec `U.unsafeIndex` k+{-# INLINE binarySearchByBounds #-}
src/Data/Matrix/Storable.hs view
@@ -1,10 +1,10 @@ module Data.Matrix.Storable     ( Matrix-    , module Data.Matrix.Generic.Base+    , module Data.Matrix.Dense.Generic     ) where -import qualified Data.Matrix.Generic.Types as MG-import Data.Matrix.Generic.Base-import qualified Data.Vector.Storable as S+import qualified Data.Matrix.Dense.Generic as MG+import Data.Matrix.Dense.Generic hiding (Matrix)+import qualified Data.Vector.Storable as V -type Matrix a = MG.Matrix S.Vector a+type Matrix a = MG.Matrix V.Vector a
src/Data/Matrix/Storable/Mutable.hs view
@@ -1,10 +1,10 @@ module Data.Matrix.Storable.Mutable     ( MMatrix-    , module Data.Matrix.Generic.Mutable+    , module Data.Matrix.Dense.Generic.Mutable     ) where -import qualified Data.Matrix.Generic.Types as MG-import Data.Matrix.Generic.Mutable-import qualified Data.Vector.Storable.Mutable as SM+import qualified Data.Matrix.Dense.Generic.Mutable as MG+import Data.Matrix.Dense.Generic.Mutable hiding (MMatrix)+import qualified Data.Vector.Storable.Mutable as VM -type MMatrix a = MG.MMatrix SM.MVector a+type MMatrix a = MG.MMatrix VM.MVector a
+ src/Data/Matrix/Symmetric.hs view
@@ -0,0 +1,96 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE TypeFamilies #-}+module Data.Matrix.Symmetric+    ( SymMatrix(..)+    , dim+    , rows+    , cols+    , unsafeIndex+    , (!)+    , flatten+    , unsafeFromVector+    , fromVector+    , takeRow+    , thaw+    , unsafeThaw+    , freeze+    , unsafeFreeze+    , zip+    , zipWith+    ) where++import Prelude hiding (zip, zipWith)+import Control.Monad (liftM)+import Data.Bits (shiftR)+import Data.Binary+import qualified Data.Vector.Generic as G++import Data.Matrix.Generic+import Data.Matrix.Symmetric.Mutable (SymMMatrix(..))++type instance Mutable SymMatrix = SymMMatrix++-- | Symmetric square matrix+data SymMatrix v a = SymMatrix !Int !(v a)+    deriving (Show)++instance (G.Vector v a, Binary a) => Binary (SymMatrix v a) where+    put (SymMatrix n vec) = do put (52 :: Int)+                               put n+                               G.mapM_ put vec+    get = do+        magic <- get :: Get Int+        if magic == 52+           then do+               n <- get +               vec <- G.replicateM ((n+1)*n `shiftR` 1) get+               return $ SymMatrix n vec+           else error "not a valid file"++instance G.Vector v a => Matrix SymMatrix v a where+    dim (SymMatrix n _) = (n,n)+    {-# INLINE dim #-}++    unsafeIndex (SymMatrix n vec) (i,j) = vec `G.unsafeIndex` idx n i j+    {-# INLINE unsafeIndex #-}++    unsafeFromVector (r,c) vec | r /= c = error "columns /= rows"+                               | otherwise = SymMatrix r . G.concat . map f $ [0..r-1]+      where+        f i = G.slice (i*(c+1)) (c-i) vec+--        n = ((r+1)*r) `shiftR` 1+    {-# INLINE unsafeFromVector #-}++    thaw (SymMatrix n v) = SymMMatrix n `liftM` G.thaw v+    {-# INLINE thaw #-}++    unsafeThaw (SymMatrix n v) = SymMMatrix n `liftM` G.thaw v+    {-# INLINE unsafeThaw #-}++    freeze (SymMMatrix n v) = SymMatrix n `liftM` G.freeze v+    {-# INLINE freeze #-}++    unsafeFreeze (SymMMatrix n v) = SymMatrix n `liftM` G.unsafeFreeze v+    {-# INLINE unsafeFreeze #-}++-- row major upper triangular indexing+idx :: Int -> Int -> Int -> Int+idx n i j | i <= j = (i * (2 * n - i - 1)) `shiftR` 1 + j+          | otherwise = (j * (2 * n - j - 1)) `shiftR` 1 + i+{-# INLINE idx #-}++zip :: (G.Vector v a, G.Vector v b, G.Vector v (a,b))+    => SymMatrix v a -> SymMatrix v b -> SymMatrix v (a,b)+zip (SymMatrix n1 v1) (SymMatrix n2 v2)+    | n1 /= n2 = error "imcompatible size"+    | otherwise = SymMatrix n1 $ G.zip v1 v2+{-# INLINE zip #-}++zipWith :: (G.Vector v a, G.Vector v b, G.Vector v c)+        => (a -> b -> c) -> SymMatrix v a -> SymMatrix v b -> SymMatrix v c+zipWith f (SymMatrix n1 v1) (SymMatrix n2 v2)+    | n1 /= n2 = error "imcompatible size"+    | otherwise = SymMatrix n1 . G.zipWith f v1 $ v2+{-# INLINE zipWith #-}
+ src/Data/Matrix/Symmetric/Mutable.hs view
@@ -0,0 +1,46 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module Data.Matrix.Symmetric.Mutable+   ( -- * Mutable Matrix+     SymMMatrix(..)+   , C.dim+   , C.write+   , C.unsafeWrite+   , C.read+   , C.unsafeRead+   , C.new+   , C.replicate+   ) where++import Prelude hiding (read, replicate)+import Control.Monad (liftM)+import Data.Bits (shiftR)+import qualified Data.Vector.Generic.Mutable as GM++import qualified Data.Matrix.Generic.Mutable as C++-- | mutable matrix+data SymMMatrix v s a = SymMMatrix !Int !(v s a)++instance GM.MVector v a => C.MMatrix SymMMatrix v a where+    dim (SymMMatrix n _) = (n,n)+    {-# INLINE dim #-}++    unsafeRead (SymMMatrix n v) (i,j) = GM.unsafeRead v (idx n i j)+    {-# INLINE unsafeRead #-}++    unsafeWrite (SymMMatrix n v) (i,j) = GM.unsafeWrite v (idx n i j)+    {-# INLINE unsafeWrite #-}++    new (r,c) | r /= c = error "colmumns /= rows"+              | otherwise = SymMMatrix r `liftM` GM.new ((r*(r+1)) `shiftR` 1)++    replicate (r,c) x+        | r /= c = error "colmumns /= rows"+        | otherwise = SymMMatrix r `liftM` GM.replicate ((r*(r+1)) `shiftR` 1) x++-- row major upper triangular indexing+idx :: Int -> Int -> Int -> Int+idx n i j | i <= j = (i * (2 * n - i - 1)) `shiftR` 1 + j+          | otherwise = (j * (2 * n - j - 1)) `shiftR` 1 + i+{-# INLINE idx #-}
src/Data/Matrix/Unboxed.hs view
@@ -1,10 +1,10 @@ module Data.Matrix.Unboxed     ( Matrix-    , module Data.Matrix.Generic.Base+    , module Data.Matrix.Dense.Generic     ) where -import qualified Data.Matrix.Generic.Types as MG-import Data.Matrix.Generic.Base+import qualified Data.Matrix.Dense.Generic as MG+import Data.Matrix.Dense.Generic hiding (Matrix) import qualified Data.Vector.Unboxed as V  type Matrix a = MG.Matrix V.Vector a
src/Data/Matrix/Unboxed/Mutable.hs view
@@ -1,10 +1,10 @@ module Data.Matrix.Unboxed.Mutable     ( MMatrix-    , module Data.Matrix.Generic.Mutable+    , module Data.Matrix.Dense.Generic.Mutable     ) where -import qualified Data.Matrix.Generic.Types as MG-import Data.Matrix.Generic.Mutable-import qualified Data.Vector.Unboxed.Mutable as UM+import qualified Data.Matrix.Dense.Generic.Mutable as MG+import Data.Matrix.Dense.Generic.Mutable hiding (MMatrix)+import qualified Data.Vector.Unboxed.Mutable as VM -type MMatrix a = MG.MMatrix UM.MVector a+type MMatrix a = MG.MMatrix VM.MVector a