diff --git a/Data/Matrix.hs b/Data/Matrix.hs --- a/Data/Matrix.hs +++ b/Data/Matrix.hs @@ -10,13 +10,15 @@ , forceMatrix -- * Builders , matrix - , fromList , fromLists , rowVector , colVector -- ** Special matrices , zero , identity , permMatrix + -- * List conversions + , fromList , fromLists + , toList , toLists -- * Accessing , getElem , (!) , unsafeGet , safeGet , getRow , getCol @@ -66,6 +68,7 @@ -- Classes import Control.DeepSeq import Control.Monad (forM_) +import Control.Loop (numLoop,numLoopFold) import Data.Foldable (Foldable (..)) import Data.Monoid import Data.Traversable @@ -223,15 +226,12 @@ -> ((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.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] - ] + numLoop 1 n $ + \i -> numLoop 1 m $ + \j -> MV.unsafeWrite v (en (i,j)) (f (i,j)) return v -- | /O(rows*cols)/. Identity matrix of the given order. @@ -262,6 +262,25 @@ {-# INLINE fromList #-} fromList n m = M n m 0 0 m . V.fromListN (n*m) +-- | Get the elements of a matrix stored in a list. +-- +-- > ( 1 2 3 ) +-- > ( 4 5 6 ) +-- > toList ( 7 8 9 ) = [1,2,3,4,5,6,7,8,9] +-- +toList :: Matrix a -> [a] +toList m = [ unsafeGet i j m | i <- [1 .. nrows m] , j <- [1 .. ncols m] ] + +-- | Get the elements of a matrix stored in a list of lists, +-- where each list contains the elements of a single row. +-- +-- > ( 1 2 3 ) [ [1,2,3] +-- > ( 4 5 6 ) , [4,5,6] +-- > toLists ( 7 8 9 ) = , [7,8,9] ] +-- +toLists :: Matrix a -> [[a]] +toLists m = [ [ unsafeGet i j m | j <- [1 .. ncols m] ] | i <- [1 .. nrows m] ] + -- | Create a matrix from a non-empty list of non-empty lists. -- /Each list must have at least as many elements as the first list/. -- Examples: @@ -545,9 +564,25 @@ -- > (bl <|> br) joinBlocks :: (Matrix a,Matrix a,Matrix a,Matrix a) -> Matrix a {-# INLINE[1] joinBlocks #-} -joinBlocks (tl,tr,bl,br) = (tl <|> tr) - <-> - (bl <|> br) +joinBlocks (tl,tr,bl,br) = + let n = nrows tl + nb = nrows bl + n' = n + nb + m = ncols tl + mr = ncols tr + m' = m + mr + en = encode m' + in M n' m' 0 0 m' $ V.create $ do + v <- MV.new (n'*m') + let wr = MV.write v + numLoop 1 n $ \i -> do + numLoop 1 m $ \j -> wr (en (i ,j )) $ tl ! (i,j) + numLoop 1 mr $ \j -> wr (en (i ,j+m)) $ tr ! (i,j) + numLoop 1 nb $ \i -> do + let i' = i+n + numLoop 1 m $ \j -> wr (en (i',j )) $ bl ! (i,j) + numLoop 1 mr $ \j -> wr (en (i',j+m)) $ br ! (i,j) + return v {-# RULES "matrix/splitAndJoin" @@ -718,16 +753,21 @@ dotProduct :: Num a => V.Vector a -> V.Vector a -> a {-# INLINE dotProduct #-} +dotProduct v1 v2 = numLoopFold 0 (V.length v1 - 1) 0 $ + \r i -> V.unsafeIndex v1 i * V.unsafeIndex v2 i + r + +{- dotProduct v1 v2 = go (V.length v1 - 1) 0 where go (-1) a = a go i a = go (i-1) $ (V.unsafeIndex v1 i) * (V.unsafeIndex v2 i) + a +-} first :: (a -> Bool) -> [a] -> a first f = go where go (x:xs) = if f x then x else go xs - go [] = error "first: no element match the condition." + go _ = error "first: no element match the condition." -- | Strassen's algorithm over square matrices of order @2^n@. strassen :: Num a => Matrix a -> Matrix a -> Matrix a @@ -769,7 +809,7 @@ in submatrix 1 n 1 m' $ strassen b1 b2 strmixFactor :: Int -strmixFactor = 500 +strmixFactor = 300 -- | Strassen's mixed algorithm. strassenMixed :: Num a => Matrix a -> Matrix a -> Matrix a @@ -926,7 +966,7 @@ -> Matrix a -- ^ Original matrix. -> Matrix a -- ^ Matrix with rows 1 and 2 switched. switchRows r1 r2 (M n m ro co w vs) = M n m ro co w $ V.modify (\mv -> do - forM_ [1..m] $ \j -> + numLoop 1 m $ \j -> MV.swap mv (encode w (r1+ro,j+co)) (encode w (r2+ro,j+co))) vs -- | Switch two coumns of a matrix. @@ -940,7 +980,7 @@ -> Matrix a -- ^ Original matrix. -> Matrix a -- ^ Matrix with cols 1 and 2 switched. switchCols c1 c2 (M n m ro co w vs) = M n m ro co w $ V.modify (\mv -> do - forM_ [1..n] $ \j -> + numLoop 1 n $ \j -> MV.swap mv (encode m (j+ro,c1+co)) (encode m (j+ro,c2+co))) vs ------------------------------------------------------- @@ -1000,16 +1040,16 @@ i = maximumBy (\x y -> compare (abs $ u ! (x,k)) (abs $ u ! (y,k))) [ k .. n ] -- Switching to place pivot in current row. u' = switchRows k i u - -- l' = switchRows k i l l' = let lw = vcols l + en = encode lw lro = rowOffset l lco = colOffset l in if i == k then l else M (nrows l) (ncols l) lro lco lw $ V.modify (\mv -> forM_ [1 .. k-1] $ - \j -> MV.swap mv (encode lw (i+lro,j+lco)) - (encode lw (k+lro,j+lco)) + \j -> MV.swap mv (en (i+lro,j+lco)) + (en (k+lro,j+lco)) ) $ mvect l p' = switchRows k i p -- Permutation determinant diff --git a/bench/mult.hs b/bench/mult.hs --- a/bench/mult.hs +++ b/bench/mult.hs @@ -20,11 +20,13 @@ bmat :: Int -> Benchmark bmat n = bgroup ("mult" ++ show n) - [ bench "Definition" $ nf testdef n - , bench "Definition 2" $ nf testdef2 n + [ -- testdef seems to be the slowest in every situation. + -- We comment it to speed up benchmarking. + -- bench "Definition" $ nf testdef n + bench "Definition 2" $ nf testdef2 n -- , bench "Strassen" $ nf teststr n , bench "Strassen mixed" $ nf teststrm n ] main :: IO () -main = defaultMain $ fmap bmat [10,25,100,150,250,500] +main = defaultMain $ fmap bmat [10,25,100,150,250,400,500] diff --git a/matrix.cabal b/matrix.cabal --- a/matrix.cabal +++ b/matrix.cabal @@ -1,5 +1,5 @@ Name: matrix -Version: 0.3.3.0 +Version: 0.3.4.0 Author: Daniel Díaz Category: Math Build-type: Simple @@ -33,8 +33,9 @@ , vector >= 0.10 && < 0.11 , deepseq >= 1.3.0.0 && < 1.4 , primitive >= 0.5 && < 0.6 + , loop >= 0.2 Exposed-modules: Data.Matrix - GHC-Options: -Wall -O2 + GHC-Options: -Wall -O2 -fstatic-argument-transformation Benchmark matrix-mult type: exitcode-stdio-1.0 @@ -43,7 +44,7 @@ build-depends: base == 4.* , matrix , criterion - ghc-options: -O2 + ghc-options: -O2 -fstatic-argument-transformation Test-Suite matrix-test type: exitcode-stdio-1.0 @@ -54,3 +55,10 @@ , tasty , QuickCheck , tasty-quickcheck + +Test-Suite matrix-examples + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: Examples.hs + build-depends: base == 4.* + , matrix diff --git a/test/Examples.hs b/test/Examples.hs new file mode 100644 --- /dev/null +++ b/test/Examples.hs @@ -0,0 +1,111 @@ + +import Data.Matrix +import System.Exit (exitFailure) +import System.IO (hFlush,stdout) + +-- We flush stdout explictly to get output printed +-- in real-time in Windows systems. +putStr' :: String -> IO () +putStr' str = putStr str >> hFlush stdout + +testExample :: String -> Bool -> IO () +testExample n b = do + putStr' $ "Test (" ++ n ++ "): " + if b then putStrLn "OK." + else putStrLn "Failed." >> exitFailure + +testEquality :: Eq a => String -> (a,a) -> IO () +testEquality n (x,y) = testExample n $ x == y + +main :: IO () +main = sequence_ + [ testEquality "matrix" + ( matrix 4 4 $ \(i,j) -> 2*i - j + , fromList 4 4 + [ 1 , 0 , -1 , -2 + , 3 , 2 , 1 , 0 + , 5 , 4 , 3 , 2 + , 7 , 6 , 5 , 4 ] + ) + , testEquality "fromList" + ( fromList 3 3 [1..] + , fromList 3 3 + [ 1 , 2 , 3 + , 4 , 5 , 6 + , 7 , 8 , 9 ] + ) + , testEquality "fromLists (1)" + ( fromLists [ [1,2,3] , [4,5,6] , [7,8,9] ] + , fromList 3 3 + [ 1 , 2 , 3 + , 4 , 5 , 6 + , 7 , 8 , 9 ] + ) + , testEquality "fromLists (2)" + ( fromLists [ [1,2,3] , [4,5,6,7] , [8,9,0] ] + , fromList 3 3 + [ 1 , 2 , 3 + , 4 , 5 , 6 + , 8 , 9 , 0 ] + ) + , testEquality "identity" + ( identity 3 , fromList 3 3 [1,0,0 , 0,1,0 , 0,0,1] + ) + , testEquality "transpose" + ( transpose $ fromList 3 3 [1..9] + , fromList 3 3 [1,4,7 , 2,5,8 , 3,6,9] + ) + , testEquality "extendTo" + ( extendTo 0 4 5 $ fromList 3 3 [1..9] + , fromList 4 5 + [ 1 , 2 , 3 , 0 , 0 + , 4 , 5 , 6 , 0 , 0 + , 7 , 8 , 9 , 0 , 0 + , 0 , 0 , 0 , 0 , 0 + ] + ) + , testEquality "mapRow" + ( mapRow (\_ x -> x + 1) 2 $ fromList 3 3 [1..9] + , fromList 3 3 [1,2,3 , 5,6,7 , 7,8,9] + ) + , testEquality "mapCol" + ( mapCol (\_ x -> x + 1) 2 $ fromList 3 3 [1..9] + , fromList 3 3 [1,3,3 , 4,6,6 , 7,9,9] + ) + , testEquality "submatrix" + ( submatrix 1 2 2 3 $ fromList 3 3 [1..9] + , fromList 2 2 [2,3 , 5,6] + ) + , testEquality "minorMatrix" + ( minorMatrix 2 2 $ fromList 3 3 [1..9] + , fromList 2 2 [1,3 , 7,9] + ) + , testEquality "scaleMatrix" + ( scaleMatrix 2 $ fromList 3 3 [1..9] + , fromList 3 3 [2,4,6 , 8,10,12 , 14,16,18] + ) + , testEquality "scaleRow" + ( scaleRow 2 2 $ fromList 3 3 [1..9] + , fromList 3 3 [1,2,3 , 8,10,12 , 7,8,9] + ) + , testEquality "combineRows" + ( combineRows 2 2 1 $ fromList 3 3 [1..9] + , fromList 3 3 [1,2,3 , 6,9,12 , 7,8,9] + ) + , testEquality "switchRows" + ( switchRows 1 2 $ fromList 3 3 [1..9] + , fromList 3 3 [4,5,6 , 1,2,3 , 7,8,9] + ) + , testEquality "switchCols" + ( switchCols 1 2 $ fromList 3 3 [1..9] + , fromList 3 3 [2,1,3 , 5,4,6 , 8,7,9] + ) + , testEquality "toList" + ( toList $ fromList 3 3 [1..9] + , [1..9] + ) + , testEquality "toLists" + ( toLists $ fromList 3 3 [1..9] + , [ [1,2,3] , [4,5,6] , [7,8,9] ] + ) + ] diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -55,7 +55,9 @@ main :: IO () main = defaultMain $ testGroup "matrix tests" [ - QC.testProperty "identity * m = m * identity = m" + QC.testProperty "zero n m = matrix n m (const 0)" + $ \(I n) (I m) -> zero n m == matrix n m (const 0) + , QC.testProperty "identity * m = m * identity = m" $ \(Sq m) -> let n = nrows m in identity n * m == m && m * identity n == m , QC.testProperty "a * (b * c) = (a * b) * c" $ \(I a) (I b) (I c) (I d) -> forAll (genMatrix a b) @@ -68,6 +70,8 @@ $ \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] ] + , QC.testProperty "fmap id = id" + $ \m -> fmap id m == (m :: Matrix R) , QC.testProperty "permMatrix n i j * permMatrix n i j = identity n" $ \(I n) -> forAll (choose (1,n)) $ \i -> forAll (choose (1,n)) @@ -91,6 +95,8 @@ $ \m -> forAll (choose (1,nrows m)) $ \i -> forAll (choose (1,ncols m)) $ \j -> joinBlocks (splitBlocks i j m) == (m :: Matrix R) + , QC.testProperty "scaleMatrix k m = fmap (*k) m" + $ \k m -> scaleMatrix k m == fmap (*k) (m :: Matrix R) , QC.testProperty "(+) = elementwise (+)" $ \m1 -> forAll (genMatrix (nrows m1) (ncols m1)) $ \m2 -> m1 + m2 == elementwise (+) m1 m2 @@ -115,4 +121,8 @@ $ \(Sq m) -> let n = nrows m in forAll (choose (1,n)) $ \i -> forAll (choose (1,n)) $ \j -> detLU (switchRows i j m) == detLU (permMatrix n i j) * detLU m + , QC.testProperty "fromList n m . toList = id" + $ \m -> fromList (nrows m) (ncols m) (toList m) == (m :: Matrix R) + , QC.testProperty "fromLists . toLists = id" + $ \m -> fromLists (toLists m) == (m :: Matrix R) ]