Skip to content
Browse files

Add primops for conversion between symmetric/hermitian matrices

in-place. Required for implementation of multiplication of
symmetric matrices
  • Loading branch information...
1 parent 181101b commit 88e6ea1ec6edcd6deb06182e65d8eee187a9458d @Shimuuar committed Aug 30, 2012
Showing with 35 additions and 0 deletions.
  1. +35 −0 Data/Matrix/Symmetric/Mutable.hs
View
35 Data/Matrix/Symmetric/Mutable.hs
@@ -23,6 +23,8 @@ module Data.Matrix.Symmetric.Mutable (
, IsHermitian
-- * Function
, new
+ , symmetricAsDense
+ , hermitianAsDense
, symmIndex
-- * Complex number
, NumberType
@@ -48,6 +50,7 @@ import Foreign.Storable
import Data.Internal
import Data.Vector.Storable.Strided.Mutable
import Data.Matrix.Generic.Mutable
+import Data.Matrix.Dense.Mutable (MMatrix(..))
import Unsafe.Coerce
@@ -151,6 +154,38 @@ new n = do
return $ MSymmetricRaw n n fp
+-- | Convert to dense matrix. Dense matrix will use same buffer as
+-- symmetric
+symmetricAsDense
+ :: (PrimMonad m, Storable a)
+ => MSymmetricRaw IsSymmetric (PrimState m) a
+ -> m (MMatrix (PrimState m) a)
+symmetricAsDense (MSymmetricRaw n lda fp) = do
+ let m = MMatrix n n lda fp
+ forM_ [0 .. n-2] $ \i -> do
+ let row = MVector (n-i-1) lda $ updPtr (`advancePtr` ( i + (i+1) * lda)) fp
+ col = MVector (n-i-1) 1 $ updPtr (`advancePtr` (1+i + i * lda)) fp
+ M.move col row
+ return m
+
+
+-- | Convert to dense matrix. Dense matrix will use same buffer as
+-- symmetric
+hermitianAsDense
+ :: (PrimMonad m, Storable a, Conjugate a)
+ => MSymmetricRaw IsSymmetric (PrimState m) a
+ -> m (MMatrix (PrimState m) a)
+hermitianAsDense (MSymmetricRaw n lda fp) = do
+ let m = MMatrix n n lda fp
+ forM_ [0 .. n-2] $ \i -> do
+ let len = n - i - 1
+ row = MVector (n-i-1) lda $ updPtr (`advancePtr` ( i + (i+1) * lda)) fp
+ col = MVector (n-i-1) 1 $ updPtr (`advancePtr` (1+i + i * lda)) fp
+ forM_ [0 .. len - 1] $ \j -> do
+ M.unsafeWrite col j . conjugateNum =<< M.unsafeRead row j
+ return m
+
+
----------------------------------------------------------------
-- Real/Complex distinction

0 comments on commit 88e6ea1

Please sign in to comment.
Something went wrong with that request. Please try again.