Permalink
Browse files

Add multiplication of symmetric and hermitian matrices

  • Loading branch information...
1 parent 88e6ea1 commit 8bf9badb43b1d81ddac2788e6b5a9dbaf3094532 @Shimuuar committed Aug 30, 2012
Showing with 114 additions and 3 deletions.
  1. +57 −1 Numeric/BLAS/Mutable.hs
  2. +57 −2 Numeric/BLAS/Mutable/Unsafe.hs
View
@@ -40,11 +40,14 @@ module Numeric.BLAS.Mutable (
, crossHVV
-- * Level 3 BLAS (Matrix-matrix operations)
, multMM
+ , multSymMM
+ , multHerMM
-- * BLAS type classes and data types
, BLAS1
, BLAS2
, BLAS3
, Trans(..)
+ , Side(..)
-- * Type classes and helpers
, MVectorBLAS(..)
, colsT
@@ -55,6 +58,8 @@ import Control.Monad.Primitive
import qualified Data.Matrix.Generic.Mutable as M
import qualified Data.Matrix.Dense.Mutable as MD
+import Data.Matrix.Dense.Mutable (MMatrix)
+import Data.Matrix.Symmetric.Mutable (MSymmetric,MHermitian,Conjugate)
import Numeric.BLAS.Mutable.Unsafe
@@ -193,7 +198,7 @@ crossHVV a v u m
-- Level 3 BLAS
----------------------------------------------------------------
--- | Unsafe multiplication of dense matrices:
+-- | Multiplication of dense matrices:
--
-- > C ← α·op(A)·op(B) + β·C
multMM :: (PrimMonad m, BLAS3 a)
@@ -215,3 +220,54 @@ multMM a ta ma tb mb b mc
rowA = rowsT ta ma ; colA = colsT ta ma
rowB = rowsT tb mb ; colB = colsT tb mb
rowC = M.rows mc ; colC = M.cols mc
+
+
+-- | Multiplication of dense matrix /B/ by symmetric matrix
+-- /A/. It could be one of operations.
+--
+-- > C ← α·A·B + β·C
+-- > C ← α·B·A + β·C
+multSymMM :: (PrimMonad m, BLAS3 a)
+ => Side -- ^ On which side matrix /A/ appears
+ -> a -- ^ /α/
+ -> MSymmetric (PrimState m) a -- ^ /A/
+ -> MMatrix (PrimState m) a -- ^ /B/
+ -> a -- ^ /β/
+ -> MMatrix (PrimState m) a -- ^ /C/
+ -> m ()
+{-# INLINE multSymMM #-}
+multSymMM side α ma mb β mc
+ | M.rows ma /= M.rows mc = error "!"
+ | M.cols ma /= M.cols mc = error "!"
+ | not okB = error "!"
+ | otherwise = unsafeMultSymMM side α ma mb β mc
+ where
+ ordB = M.cols mb
+ okB = case side of
+ RightSide -> ordB == M.rows ma
+ LeftSide -> ordB == M.cols ma
+
+-- | Multiplication of dense matrix /B/ by hermitian matrix
+-- /A/. It could be one of operations.
+--
+-- > C ← α·A·B + β·C
+-- > C ← α·B·A + β·C
+multHerMM :: (PrimMonad m, BLAS3 a, Conjugate a)
+ => Side -- ^ On which side matrix /A/ appears
+ -> a -- ^ /α/
+ -> MHermitian (PrimState m) a -- ^ /A/
+ -> MMatrix (PrimState m) a -- ^ /B/
+ -> a -- ^ /β/
+ -> MMatrix (PrimState m) a -- ^ /C/
+ -> m ()
+{-# INLINE multHerMM #-}
+multHerMM side α ma mb β mc
+ | M.rows ma /= M.rows mc = error "!"
+ | M.cols ma /= M.cols mc = error "!"
+ | not okB = error "!"
+ | otherwise = unsafeMultHerMM side α ma mb β mc
+ where
+ ordB = M.cols mb
+ okB = case side of
+ RightSide -> ordB == M.rows ma
+ LeftSide -> ordB == M.cols ma
@@ -38,11 +38,14 @@ module Numeric.BLAS.Mutable.Unsafe (
, unsafeCrossHVV
-- * Level 3 BLAS (Matrix-matrix operations)
, unsafeMultMM
+ , unsafeMultSymMM
+ , unsafeMultHerMM
-- * BLAS type classes and data types
, BLAS1
, BLAS2
, BLAS3
, Trans(..)
+ , Side(..)
-- * Type classes and helpers
, MVectorBLAS(..)
, colsT
@@ -56,7 +59,7 @@ import Foreign.ForeignPtr
import qualified Numeric.BLAS.Bindings as BLAS
import Numeric.BLAS.Bindings (BLAS1,BLAS2,BLAS3,RealType,
- RowOrder(..), Trans(..), Uplo(..)
+ RowOrder(..), Trans(..), Uplo(..), Side
)
import qualified Data.Vector.Storable as S
@@ -65,7 +68,7 @@ import qualified Data.Matrix.Generic.Mutable as M
-- Concrete matrices
import Data.Matrix.Dense.Mutable (MMatrix(..))
import Data.Matrix.Symmetric.Mutable
- (MSymmetricRaw(..),IsSymmetric,IsHermitian,Conjugate(..), NumberType, IsReal)
+ (MSymmetric,MHermitian,MSymmetricRaw(..),IsSymmetric,IsHermitian,Conjugate(..),NumberType,IsReal)
@@ -293,6 +296,58 @@ unsafeMultMM a ta ma@(MMatrix _ _ lda fpa) tb (MMatrix _ _ ldb fpb) b (MMatrix r
a pa lda pb ldb
b pc ldc
+-- | Unsafe multiplication of dense matrix /B/ by symmetric matrix
+-- /A/. It could be one of operations.
+--
+-- > C ← α·A·B + β·C
+-- > C ← α·B·A + β·C
+unsafeMultSymMM :: (PrimMonad m, BLAS3 a)
+ => Side -- ^ On which side matrix /A/ appears
+ -> a -- ^ /α/
+ -> MSymmetric (PrimState m) a -- ^ /A/
+ -> MMatrix (PrimState m) a -- ^ /B/
+ -> a -- ^ /β/
+ -> MMatrix (PrimState m) a -- ^ /C/
+ -> m ()
+{-# INLINE unsafeMultSymMM #-}
+unsafeMultSymMM side α (MSymmetricRaw _ lda fpa) (MMatrix _ _ ldb fpb) β m@(MMatrix _ _ ldc fpc)
+ = unsafePrimToPrim
+ $ withForeignPtr fpa $ \pa ->
+ withForeignPtr fpb $ \pb ->
+ withForeignPtr fpc $ \pc ->
+ BLAS.symm ColMajor side Upper
+ (M.rows m) (M.cols m)
+ α pa lda
+ pb ldb
+ β pc ldc
+
+-- | Unsafe multiplication of dense matrix /B/ by hermitian matrix
+-- /A/. It could be one of operations.
+--
+-- > C ← α·A·B + β·C
+-- > C ← α·B·A + β·C
+unsafeMultHerMM :: (PrimMonad m, BLAS3 a)
+ => Side -- ^ On which side matrix /A/ appears
+ -> a -- ^ /α/
+ -> MHermitian (PrimState m) a -- ^ /A/
+ -> MMatrix (PrimState m) a -- ^ /B/
+ -> a -- ^ /β/
+ -> MMatrix (PrimState m) a -- ^ /C/
+ -> m ()
+{-# INLINE unsafeMultHerMM #-}
+unsafeMultHerMM side α (MSymmetricRaw _ lda fpa) (MMatrix _ _ ldb fpb) β m@(MMatrix _ _ ldc fpc)
+ = unsafePrimToPrim
+ $ withForeignPtr fpa $ \pa ->
+ withForeignPtr fpb $ \pb ->
+ withForeignPtr fpc $ \pc ->
+ BLAS.hemm ColMajor side Upper
+ (M.rows m) (M.cols m)
+ α pa lda
+ pb ldb
+ β pc ldc
+
+
+
----------------------------------------------------------------

0 comments on commit 8bf9bad

Please sign in to comment.