matrix 0.3.2.0 → 0.3.3.0
raw patch · 4 files changed
+123/−37 lines, 4 filesPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
+ Data.Matrix: elementwiseUnsafe :: (a -> b -> c) -> (Matrix a -> Matrix b -> Matrix c)
Files
- Data/Matrix.hs +116/−34
- bench/mult.hs +1/−1
- matrix.cabal +1/−1
- test/Main.hs +5/−1
Data/Matrix.hs view
@@ -36,7 +36,7 @@ , (<|>) , (<->) , joinBlocks -- * Matrix operations- , elementwise+ , elementwise, elementwiseUnsafe -- * Matrix multiplication -- ** About matrix multiplication -- $mult@@ -136,7 +136,7 @@ -- Useful when using 'submatrix' from a big matrix. -- forceMatrix :: Matrix a -> Matrix a-forceMatrix m = matrix (nrows m) (ncols m) $ \(i,j) -> m ! (i,j)+forceMatrix m = matrix (nrows m) (ncols m) $ \(i,j) -> unsafeGet i j m ------------------------------------------------------- -------------------------------------------------------@@ -223,7 +223,16 @@ -> ((Int,Int) -> a) -- ^ Generator function -> Matrix a {-# INLINE matrix #-}-matrix n m f = M n m 0 0 m $ V.generate (n*m) $ f . decode m+-- matrix n m f = M n m 0 0 m $ V.generate (n*m) $ f . decode m+matrix n m f = M n m 0 0 m $ V.create $ do+ v <- MV.new $ n * m+ let en = encode m+ sequence_+ [ MV.unsafeWrite v (en (i,j)) (f (i,j))+ | i <- [1 .. n]+ , j <- [1 .. m]+ ]+ return v -- | /O(rows*cols)/. Identity matrix of the given order. --@@ -347,6 +356,11 @@ {-# INLINE (!) #-} m ! (i,j) = getElem i j m +-- | Internal alias for 'unsafeGet'.+(!.) :: Matrix a -> (Int,Int) -> a+{-# INLINE (!.) #-}+m !. (i,j) = unsafeGet i j m+ -- | Variant of 'getElem' that returns Maybe instead of an error. safeGet :: Int -> Int -> Matrix a -> Maybe a safeGet i j a@(M n m _ _ _ _)@@ -572,12 +586,36 @@ ------------------------------------------------------- ---- MATRIX OPERATIONS --- | Perform an operation element-wise. The input matrices are assumed--- to have the same dimensions, but this is not checked.+-- | Perform an operation element-wise.+-- The second matrix must have at least as many rows+-- and columns as the first matrix. If it's bigger,+-- the leftover items will be ignored.+-- If it's smaller, it will cause a run-time error.+-- You may want to use 'elementwiseUnsafe' if you+-- are definitely sure that a run-time error won't+-- arise. elementwise :: (a -> b -> c) -> (Matrix a -> Matrix b -> Matrix c) elementwise f m m' = matrix (nrows m) (ncols m) $- \(i,j) -> f (m ! (i,j)) (m' ! (i,j))+ \k -> f (m ! k) (m' ! k) +-- | Unsafe version of 'elementwise', but faster.+elementwiseUnsafe :: (a -> b -> c) -> (Matrix a -> Matrix b -> Matrix c)+{-# INLINE elementwiseUnsafe #-}+elementwiseUnsafe f m m' = matrix (nrows m) (ncols m) $+ \(i,j) -> f (unsafeGet i j m) (unsafeGet i j m')++infixl 6 +., -.++-- | Internal unsafe addition.+(+.) :: Num a => Matrix a -> Matrix a -> Matrix a+{-# INLINE (+.) #-}+(+.) = elementwiseUnsafe (+)++-- | Internal unsafe substraction.+(-.) :: Num a => Matrix a -> Matrix a -> Matrix a+{-# INLINE (-.) #-}+(-.) = elementwiseUnsafe (-)+ ------------------------------------------------------- ------------------------------------------------------- ---- MATRIX MULTIPLICATION@@ -594,8 +632,8 @@ * 'multStd2': Matrix multiplication following directly the definition. However, using a different definition from 'multStd'.- Some benchmarks show 'multStd2' beats in performance- to 'multStd' when the size of the matrix is around 100x100.+ According to our benchmarks with this version, 'multStd2' is+ around 3 times faster than 'multStd'. * 'multStrassen': Matrix multiplication following the Strassen's algorithm.@@ -643,11 +681,11 @@ multStd_ a@(M 2 2 _ _ _ _) b@(M 2 2 _ _ _ _) = M 2 2 0 0 2 $ let -- A- a11 = a ! (1,1) ; a12 = a ! (1,2)- a21 = a ! (2,1) ; a22 = a ! (2,2)+ a11 = a !. (1,1) ; a12 = a !. (1,2)+ a21 = a !. (2,1) ; a22 = a !. (2,2) -- B- b11 = b ! (1,1) ; b12 = b ! (1,2)- b21 = b ! (2,1) ; b22 = b ! (2,2)+ b11 = b !. (1,1) ; b12 = b !. (1,2)+ b21 = b !. (2,1) ; b22 = b !. (2,2) in V.fromList [ a11*b11 + a12*b21 , a11*b12 + a12*b22 , a21*b11 + a22*b21 , a21*b12 + a22*b22@@ -655,19 +693,19 @@ multStd_ a@(M 3 3 _ _ _ _) b@(M 3 3 _ _ _ _) = M 3 3 0 0 3 $ let -- A- a11 = a ! (1,1) ; a12 = a ! (1,2) ; a13 = a ! (1,3)- a21 = a ! (2,1) ; a22 = a ! (2,2) ; a23 = a ! (2,3)- a31 = a ! (3,1) ; a32 = a ! (3,2) ; a33 = a ! (3,3)+ a11 = a !. (1,1) ; a12 = a !. (1,2) ; a13 = a !. (1,3)+ a21 = a !. (2,1) ; a22 = a !. (2,2) ; a23 = a !. (2,3)+ a31 = a !. (3,1) ; a32 = a !. (3,2) ; a33 = a !. (3,3) -- B- b11 = b ! (1,1) ; b12 = b ! (1,2) ; b13 = b ! (1,3)- b21 = b ! (2,1) ; b22 = b ! (2,2) ; b23 = b ! (2,3)- b31 = b ! (3,1) ; b32 = b ! (3,2) ; b33 = b ! (3,3)+ b11 = b !. (1,1) ; b12 = b !. (1,2) ; b13 = b !. (1,3)+ b21 = b !. (2,1) ; b22 = b !. (2,2) ; b23 = b !. (2,3)+ b31 = b !. (3,1) ; b32 = b !. (3,2) ; b33 = b !. (3,3) in V.fromList [ a11*b11 + a12*b21 + a13*b31 , a11*b12 + a12*b22 + a13*b32 , a11*b13 + a12*b23 + a13*b33 , a21*b11 + a22*b21 + a23*b31 , a21*b12 + a22*b22 + a23*b32 , a21*b13 + a22*b23 + a23*b33 , a31*b11 + a32*b21 + a33*b31 , a31*b12 + a32*b22 + a33*b32 , a31*b13 + a32*b23 + a33*b33 ]-multStd_ a@(M n m _ _ _ _) b@(M _ m' _ _ _ _) = matrix n m' $ \(i,j) -> sum [ a ! (i,k) * b ! (k,j) | k <- [1 .. m] ]+multStd_ a@(M n m _ _ _ _) b@(M _ m' _ _ _ _) = matrix n m' $ \(i,j) -> sum [ a !. (i,k) * b !. (k,j) | k <- [1 .. m] ] multStd__ :: Num a => Matrix a -> Matrix a -> Matrix a {-# INLINE multStd__ #-}@@ -731,7 +769,7 @@ in submatrix 1 n 1 m' $ strassen b1 b2 strmixFactor :: Int-strmixFactor = 2 ^ (9 :: Int)+strmixFactor = 500 -- | Strassen's mixed algorithm. strassenMixed :: Num a => Matrix a -> Matrix a -> Matrix a@@ -744,7 +782,52 @@ a' = setSize 0 r' r' a b' = setSize 0 r' r' b in submatrix 1 r 1 r $ strassenMixed a' b'- | otherwise = joinBlocks (c11,c12,c21,c22)+ | otherwise =+ M r r 0 0 r $ V.create $ do+ v <- MV.unsafeNew (r*r)+ let en = encode r+ n' = n + 1+ -- c11 = p1 + p4 - p5 + p7+ sequence_ [ MV.write v k $+ unsafeGet i j p1+ + unsafeGet i j p4+ - unsafeGet i j p5+ + unsafeGet i j p7+ | i <- [1..n]+ , j <- [1..n]+ , let k = en (i,j)+ ]+ -- c12 = p3 + p5+ sequence_ [ MV.write v k $+ unsafeGet i j' p3+ + unsafeGet i j' p5+ | i <- [1..n]+ , j <- [n'..r]+ , let k = en (i,j)+ , let j' = j - n+ ]+ -- c21 = p2 + p4+ sequence_ [ MV.write v k $+ unsafeGet i' j p2+ + unsafeGet i' j p4+ | i <- [n'..r]+ , j <- [1..n]+ , let k = en (i,j)+ , let i' = i - n+ ]+ -- c22 = p1 - p2 + p3 + p6+ sequence_ [ MV.write v k $+ unsafeGet i' j' p1+ - unsafeGet i' j' p2+ + unsafeGet i' j' p3+ + unsafeGet i' j' p6+ | i <- [n'..r]+ , j <- [n'..r]+ , let k = en (i,j)+ , let i' = i - n+ , let j' = j - n+ ]+ return v where r = nrows a -- Size of the subproblem is halved.@@ -753,18 +836,13 @@ (a11,a12,a21,a22) = splitBlocks n n a (b11,b12,b21,b22) = splitBlocks n n b -- The seven Strassen's products.- p1 = strassenMixed (a11 + a22) (b11 + b22)- p2 = strassenMixed (a21 + a22) b11- p3 = strassenMixed a11 (b12 - b22)- p4 = strassenMixed a22 (b21 - b11)- p5 = strassenMixed (a11 + a12) b22- p6 = strassenMixed (a21 - a11) (b11 + b12)- p7 = strassenMixed (a12 - a22) (b21 + b22)- -- Merging blocks- c11 = p1 + p4 - p5 + p7- c12 = p3 + p5- c21 = p2 + p4- c22 = p1 - p2 + p3 + p6+ p1 = strassenMixed (a11 +. a22) (b11 +. b22)+ p2 = strassenMixed (a21 +. a22) b11+ p3 = strassenMixed a11 (b12 -. b22)+ p4 = strassenMixed a22 (b21 -. b11)+ p5 = strassenMixed (a11 +. a12) b22+ p6 = strassenMixed (a21 -. a11) (b11 +. b12)+ p7 = strassenMixed (a12 -. a22) (b21 +. b22) -- | Mixed Strassen's matrix multiplication. multStrassenMixed :: Num a => Matrix a -> Matrix a -> Matrix a@@ -772,7 +850,7 @@ multStrassenMixed a1@(M n m _ _ _ _) a2@(M n' m' _ _ _ _) | m /= n' = error $ "Multiplication of " ++ sizeStr n m ++ " and " ++ sizeStr n' m' ++ " matrices."- | n < strmixFactor = multStd_ a1 a2+ | n < strmixFactor = multStd__ a1 a2 | otherwise = let mx = maximum [n,m,n',m'] n2 = if even mx then mx else mx+1@@ -897,6 +975,8 @@ -- > ( 1 2 0 ) ( 2 0 2 ) ( 1 0 0 ) ( 0 0 1 ) -- > ( 0 2 1 ) ( 0 2 -1 ) ( 1/2 1 0 ) ( 1 0 0 ) -- > luDecomp ( 2 0 2 ) = ( ( 0 0 2 ) , ( 0 1 1 ) , ( 0 1 0 ) , 1 )+--+-- 'Nothing' is returned if no LU decomposition exists. luDecomp :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a,Matrix a,Matrix a,a) luDecomp a = recLUDecomp a i i 1 1 n where@@ -975,6 +1055,8 @@ -- > ( 1 0 ) ( 2 1 ) ( 1 0 0 ) ( 0 0 1 ) -- > ( 0 2 ) ( 0 2 ) ( 0 1 0 ) ( 0 1 0 ) ( 1 0 ) -- > luDecomp' ( 2 1 ) = ( ( 0 0 ) , ( 1/2 -1/4 1 ) , ( 1 0 0 ) , ( 0 1 ) , -1 , 1 )+--+-- 'Nothing' is returned if no LU decomposition exists. luDecomp' :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a,Matrix a,Matrix a,Matrix a,a,a) luDecomp' a = recLUDecomp' a i i (identity $ ncols a) 1 1 1 n where
bench/mult.hs view
@@ -4,7 +4,7 @@ import Data.Matrix mat :: Int -> Matrix Int -mat n = matrix n n $ \(i,j) -> i - j +mat n = fromList n n [1..] testdef :: Int -> Matrix Int testdef n = multStd (mat n) (mat n)
matrix.cabal view
@@ -1,5 +1,5 @@ Name: matrix -Version: 0.3.2.0 +Version: 0.3.3.0 Author: Daniel Díaz Category: Math Build-type: Simple
test/Main.hs view
@@ -62,8 +62,12 @@ $ \m1 -> forAll (genMatrix b c) $ \m2 -> forAll (genMatrix c d) $ \m3 -> m1 * (m2 * m3) == (m1 * m2) * m3+ , QC.testProperty "multStd a b = multStd2 a b"+ $ \(I a) (I b) (I c) -> forAll (genMatrix a b)+ $ \m1 -> forAll (genMatrix b c)+ $ \m2 -> multStd m1 m2 == multStd2 m1 m2 , QC.testProperty "getMatrixAsVector m = mconcat [ getRow i m | i <- [1 .. nrows m]]"- $ \m -> getMatrixAsVector (m :: Matrix R) == mconcat [ getRow i m | i <- [1 .. nrows m] ]+ $ \m -> getMatrixAsVector (m :: Matrix R) == mconcat [ getRow i m | i <- [1 .. nrows m] ] , QC.testProperty "permMatrix n i j * permMatrix n i j = identity n" $ \(I n) -> forAll (choose (1,n)) $ \i -> forAll (choose (1,n))