Add Fast modules to DLA (#59)

* Added a Fast module to DLA, modified benchmarks accordingly

Author: Magalame

dense-linear-algebra/bench/ChronosBench.hs

@@ -1,6 +1,8 @@
 module Main where
 
 import qualified Statistics.Matrix as M
+import qualified Statistics.Matrix.Fast as F
+import qualified Statistics.Matrix.Fast.Algorithms as FA
 import           Statistics.Matrix (Matrix (..))
 import qualified Statistics.Matrix.Algorithms as A
 
@@ -31,9 +33,11 @@ runtimelight v a = do
       v2 = U.take n v
 
     C.defaultMainWith (C.defaultConfig {C.timeout = Just 3}) [
-                   C.bench "norm" M.norm v,
+                   C.bench "norm" M.norm v2,
+                   C.bench "Fast.norm" F.norm v2,
 
                    C.bench "multiplyV" (M.multiplyV a) (v2),
+                   C.bench "Fast.multiplyV" (F.multiplyV a) (v2),
 
                    C.bench "transpose" M.transpose a ,
                    C.bench "ident" M.ident n,
@@ -45,7 +49,9 @@ runtimeheavy a b = do
 
     C.defaultMainWith (C.defaultConfig {C.timeout = Just 1}) [
                    C.bench "multiply" (M.multiply a) b,
-                   C.bench "qr" A.qr a
+                   C.bench "Fast.multiply" (F.multiply a) b,
+                   C.bench "qr" A.qr a,
+                   C.bench "Fast.qr" FA.qr a
                    ]
 
 
@@ -60,4 +66,4 @@ main = do
     -- we split heavy and light, we lose some precision in the bar plots from chronos
     runtimelight v a
     putStrLn "---Benchmarking heavy operations---"
-    runtimeheavy a b
+    runtimeheavy a b

dense-linear-algebra/bench/WeighBench.hs

@@ -1,6 +1,8 @@
 module Main where
 
 import qualified Statistics.Matrix as M
+import qualified Statistics.Matrix.Fast as F
+import qualified Statistics.Matrix.Fast.Algorithms as FA
 import           Statistics.Matrix (Matrix (..))
 import qualified Statistics.Matrix.Algorithms as A
 
@@ -32,14 +34,21 @@ weight v a b = do
       v2 = U.take n v
     W.mainWith (do 
         W.func "norm" M.norm v2
+        W.func "Fast.norm" F.norm v2
 
-        W.func "multiply" (M.multiply a) b
         W.func "multiplyV" (M.multiplyV a) (v2)
-        W.func "qr" A.qr a
-
+        W.func "Fast.multiplyV" (F.multiplyV a) (v2)
         W.func "transpose" M.transpose a
         W.func "ident" M.ident n
-        W.func "diag" M.diag v2)
+        W.func "diag" M.diag v2
+
+        W.func "multiply" (M.multiply a) b
+        W.func "Fast.multiply" (F.multiply a) b
+        
+        W.func "qr" A.qr a
+        W.func "Fast.qr" FA.qr a
+
+        )
 
 
 

dense-linear-algebra/dense-linear-algebra.cabal

@@ -26,6 +26,8 @@ library
                        Statistics.Matrix.Function
                        Statistics.Matrix.Mutable
                        Statistics.Matrix.Types
+                       Statistics.Matrix.Fast
+                       Statistics.Matrix.Fast.Algorithms
   build-depends:       base                    >= 4.5 && < 5
                      , deepseq                 >= 1.1.0.2
                      , math-functions          >= 0.1.7

dense-linear-algebra/src/Statistics/Matrix/Fast.hs

@@ -0,0 +1,50 @@
+{-# LANGUAGE BangPatterns #-}
+
+module Statistics.Matrix.Fast (
+    multiply,
+    norm,
+    multiplyV
+    ) where
+
+import Prelude hiding (exponent, map)
+import Control.Monad.ST
+import qualified Data.Vector.Unboxed as U
+
+
+import Statistics.Matrix (row)
+import Statistics.Matrix.Function
+import Statistics.Matrix.Types
+import Statistics.Matrix.Mutable  (unsafeNew,unsafeWrite,unsafeFreeze)
+
+-- | Matrix-matrix multiplication in a more imperative fashion. Matrices must be of compatible
+-- sizes (/note: not checked/). Faster but less accurate than Statistics.Matrix.multiply
+multiply :: Matrix -> Matrix -> Matrix
+multiply m1@(Matrix r1 _ _) m2@(Matrix _ c2 _) = runST $ do
+  m3 <- unsafeNew r1 c2
+  for 0 c2 $ \j -> do
+    for 0 r1 $ \i -> do
+      let 
+        z = accum i m1 j m2 
+      unsafeWrite m3 i j z
+  unsafeFreeze m3
+
+accum :: Int -> Matrix -> Int -> Matrix -> Double
+accum ithrow (Matrix r1 c1 v1) jthcol (Matrix _ c2 v2) = sub 0 0
+  where sub !acc !ij | ij == r1 = acc
+                     | otherwise = sub ( valRow*valCol + acc ) (ij+1)
+                                   where 
+                                    valRow = U.unsafeIndex v1 (ithrow*c1 + ij)
+                                    valCol = U.unsafeIndex v2 (ij*c2+jthcol)
+
+-- | Matrix-vector multiplication, with better performances but not as accurate as
+-- Statistics.Matrix.multiplyV
+multiplyV :: Matrix -> Vector -> Vector
+multiplyV m v
+  | cols m == c = U.generate (rows m) (U.sum . U.zipWith (*) v . row m)
+  | otherwise   = error $ "matrix/vector unconformable " ++ show (cols m,c)
+  where c = U.length v
+
+-- | Norm of a vector. Faster but less accurate than Statistics.Matrix.norm
+norm :: Vector -> Double
+norm = sqrt . U.sum . U.map square
+

dense-linear-algebra/src/Statistics/Matrix/Fast/Algorithms.hs

@@ -0,0 +1,59 @@
+{-# LANGUAGE BangPatterns #-}
+
+-- |
+-- Module    : Statistics.Matrix.Fast.Algorithms
+-- Copyright : 2019 Magalame
+-- License   : BSD3
+--
+-- Useful matrix functions.
+
+module Statistics.Matrix.Fast.Algorithms
+    (
+      qr
+    ) where
+
+import Control.Applicative ((<$>), (<*>))
+import Control.Monad.ST (ST, runST)
+import Prelude hiding (replicate)
+import Statistics.Matrix (Matrix (..),dimension, for)
+import qualified Statistics.Matrix.Mutable as M
+import qualified Data.Vector.Unboxed as U
+
+-- | /O(r*c)/ Compute the QR decomposition of a matrix.
+-- The result returned is the matrices (/q/,/r/).
+qr :: Matrix -> (Matrix, Matrix)
+qr mat = runST $ do
+  let (m,n) = dimension mat
+
+  r <- M.replicate n n 0
+  a <- M.thaw mat
+  for 0 n $ \j -> do
+    cn <- M.immutably a $ \aa -> sqrt $ normCol j aa
+    M.unsafeWrite r j j cn
+    for 0 m $ \i -> M.unsafeModify a i j (/ cn)
+    for (j+1) n $ \jj -> do
+      p <- innerProduct a j jj
+      M.unsafeWrite r j jj p
+      for 0 m $ \i -> do
+        aij <- M.unsafeRead a i j
+        M.unsafeModify a i jj $ subtract (p * aij)
+  (,) <$> M.unsafeFreeze a <*> M.unsafeFreeze r
+
+normCol :: Int -> Matrix -> Double
+normCol jthcol (Matrix r c v) = sub 0 0
+  where sub !acc !ij | ij == r = acc
+                     | otherwise = sub ( valCol*valCol + acc ) (ij+1)
+                                   where
+                                    valCol = U.unsafeIndex v (ij*c+jthcol)
+
+innerProduct :: M.MMatrix s -> Int -> Int -> ST s Double
+innerProduct mmat j k = M.immutably mmat $ \mat ->
+  dotCol j mat k mat
+
+dotCol :: Int -> Matrix -> Int -> Matrix -> Double
+dotCol jthcol (Matrix r1 c1 v1) kthcol (Matrix _ c2 v2) = sub 0 0
+  where sub !acc !ij | ij == r1 = acc
+                     | otherwise = sub ( valColj*valColk + acc ) (ij+1)
+                                   where
+                                    valColk = U.unsafeIndex v2 (ij*c2+kthcol)
+                                    valColj = U.unsafeIndex v1 (ij*c1+jthcol)

dense-linear-algebra/stack.yaml

@@ -9,6 +9,7 @@ extra-deps:
 - torsor-0.1
 - chronos-1.0.5
 - chronos-bench-0.2.0.2
+- primitive-0.6.4.0
 
 # Override default flag values for local packages and extra-deps
 # flags: {}