generalized diagRect

8 files changed, 44 insertions, 47 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index c0824a3..94b56cf 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -221,13 +221,13 @@ where r is the desired number of rows.)
  , 9.0, 10.0, 11.0, 12.0 ]@
 
 -}
-reshape :: Element t => Int -> Vector t -> Matrix t
+reshape :: Storable t => Int -> Vector t -> Matrix t
 reshape c v = matrixFromVector RowMajor c v
 
 singleton x = reshape 1 (fromList [x])
 
 -- | application of a vector function on the flattened matrix elements
-liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
+liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
 liftMatrix f MC { icols = c, cdat = d } = matrixFromVector RowMajor    c (f d)
 liftMatrix f MF { icols = c, fdat = d } = matrixFromVector ColumnMajor c (f d)
 
@@ -246,7 +246,6 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
 ------------------------------------------------------------------
 
 -- | Supported element types for basic matrix operations.
---class (Storable a, Floating a) => Element a where
 class (Storable a) => Element a where
     subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position 
                -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index b8c309c..ea16748 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -22,7 +22,7 @@ module Data.Packed.Matrix (
     Element,
     Matrix,rows,cols,
     (><),
-    trans, ctrans,
+    trans,
     reshape, flatten,
     fromLists, toLists, buildMatrix,
     (@@>),
@@ -33,7 +33,7 @@ module Data.Packed.Matrix (
     flipud, fliprl,
     subMatrix, takeRows, dropRows, takeColumns, dropColumns,
     extractRows,
-    ident, diag, diagRect, takeDiag,
+    diagRect, takeDiag,
     liftMatrix, liftMatrix2, liftMatrix2Auto,
     dispf, disps, dispcf, vecdisp, latexFormat, format,
     loadMatrix, saveMatrix, fromFile, fileDimensions,
@@ -169,28 +169,19 @@ fliprl m = fromColumns . reverse . toColumns $ m
 
 ------------------------------------------------------------
 
--- | Creates a square matrix with a given diagonal.
-diag :: (Num a, Element a) => Vector a -> Matrix a
-diag v = ST.runSTMatrix $ do
-    let d = dim v
-    m <- ST.newMatrix 0 d d
-    mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]
-    return m
+{- | creates a rectangular diagonal matrix:
 
-{- | creates a rectangular diagonal matrix
-
-@> diagRect (constant 5 3) 3 4 :: Matrix Double
-(3><4)
- [ 5.0, 0.0, 0.0, 0.0
- , 0.0, 5.0, 0.0, 0.0
- , 0.0, 0.0, 5.0, 0.0 ]@
+@> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double
+(4><5)
+ [ 10.0,  7.0,  7.0, 7.0, 7.0
+ ,  7.0, 20.0,  7.0, 7.0, 7.0
+ ,  7.0,  7.0, 30.0, 7.0, 7.0
+ ,  7.0,  7.0,  7.0, 7.0, 7.0 ]@
 -}
-diagRect :: (Element t, Num t) => Vector t -> Int -> Int -> Matrix t
-diagRect v r c
-    | dim v < min r c = error "diagRect called with dim v < min r c"
-    | otherwise = ST.runSTMatrix $ do
-        m <- ST.newMatrix 0 r c
-        let d = min r c
+diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t
+diagRect z v r c = ST.runSTMatrix $ do
+        m <- ST.newMatrix z r c
+        let d = min r c `min` (dim v)
         mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]
         return m
 
@@ -198,10 +189,6 @@ diagRect v r c
 takeDiag :: (Element t) => Matrix t -> Vector t
 takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
 
--- | creates the identity matrix of given dimension
-ident :: (Num a, Element a) => Int -> Matrix a
-ident n = diag (constantD 1 n)
-
 ------------------------------------------------------------
 
 {- | An easy way to create a matrix:
@@ -225,7 +212,7 @@ Example:
  , 4.0, 5.0, 6.0 ]@
 
 -}
-(><) :: (Element a) => Int -> Int -> [a] -> Matrix a
+(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a
 r >< c = f where
     f l | dim v == r*c = matrixFromVector RowMajor c v
         | otherwise    = error $ "inconsistent list size = "
@@ -261,16 +248,13 @@ fromLists :: Element t => [[t]] -> Matrix t
 fromLists = fromRows . map fromList
 
 -- | creates a 1-row matrix from a vector
-asRow :: Element a => Vector a -> Matrix a
+asRow :: Storable a => Vector a -> Matrix a
 asRow v = reshape (dim v) v
 
 -- | creates a 1-column matrix from a vector
-asColumn :: Element a => Vector a -> Matrix a
+asColumn :: Storable a => Vector a -> Matrix a
 asColumn v = reshape 1 v
 
--- | conjugate transpose
-ctrans :: Element e => Matrix e -> Matrix e
-ctrans = liftMatrix conjugateD . trans
 
 
 {- | creates a Matrix of the specified size using the supplied function to
@@ -289,7 +273,7 @@ buildMatrix rc cc f =
 
 -----------------------------------------------------
 
-fromArray2D :: (Element e) => Array (Int, Int) e -> Matrix e
+fromArray2D :: (Storable e) => Array (Int, Int) e -> Matrix e
 fromArray2D m = (r><c) (elems m)
     where ((r0,c0),(r1,c1)) = bounds m
           r = r1-r0+1
diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs
index 48e35b4..652f43e 100644
--- a/lib/Data/Packed/ST.hs
+++ b/lib/Data/Packed/ST.hs
@@ -90,11 +90,11 @@ writeVector :: Storable t => STVector s t -> Int -> t -> ST s ()
 writeVector = safeIndexV unsafeWriteVector
 
 {-# NOINLINE newUndefinedVector #-}
-newUndefinedVector :: Element t => Int -> ST s (STVector s t)
+newUndefinedVector :: Storable t => Int -> ST s (STVector s t)
 newUndefinedVector = unsafeIOToST . fmap STVector . createVector
 
 {-# INLINE newVector #-}
-newVector :: Element t => t -> Int -> ST s (STVector s t)
+newVector :: Storable t => t -> Int -> ST s (STVector s t)
 newVector x n = do
     v <- newUndefinedVector n
     let go (-1) = return v
@@ -164,9 +164,9 @@ writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
 writeMatrix = safeIndexM unsafeWriteMatrix
 
 {-# NOINLINE newUndefinedMatrix #-}
-newUndefinedMatrix :: Element t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t)
+newUndefinedMatrix :: Storable t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t)
 newUndefinedMatrix order r c = unsafeIOToST $ fmap STMatrix $ createMatrix order r c
 
 {-# NOINLINE newMatrix #-}
-newMatrix :: Element t => t -> Int -> Int -> ST s (STMatrix s t)
+newMatrix :: Storable t => t -> Int -> Int -> ST s (STMatrix s t)
 newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c)
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index 83bf44e..1afc5a1 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -22,6 +22,7 @@
 module Numeric.Container (
     -- * Generic operations
     Container(..),
+    ctrans, diag, ident,
     -- * Matrix product and related functions
     Product(..),
     mXm,mXv,vXm,
@@ -221,6 +222,20 @@ instance (Container Vector a) => Container Matrix a where
 
 ----------------------------------------------------
 
+-- | conjugate transpose
+ctrans :: Element e => Matrix e -> Matrix e
+ctrans = liftMatrix conjugateD . trans
+
+-- | Creates a square matrix with a given diagonal.
+diag :: (Num a, Element a) => Vector a -> Matrix a
+diag v = diagRect 0 v n n where n = dim v
+
+-- | creates the identity matrix of given dimension
+ident :: (Num a, Element a) => Int -> Matrix a
+ident n = diag (constantD 1 n)
+
+----------------------------------------------------
+
 
 -- | Matrix product and related functions
 class Element e => Product e where
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 4f6f54d..394a1d7 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -183,7 +183,7 @@ singularValues = {-# SCC "singularValues" #-} sv'
 fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
 fullSVD m = (u,d,v) where
     (u,s,v) = svd m
-    d = diagRect s r c
+    d = diagRect 0 s r c
     r = rows m
     c = cols m
 
@@ -210,7 +210,7 @@ leftSV m  | vertical m = let (u,s,_) = svd m     in (u,s)
 {-# DEPRECATED full "use fullSVD instead" #-}
 full svdFun m = (u, d ,v) where
     (u,s,v) = svdFun m
-    d = diagRect s r c
+    d = diagRect 0 s r c
     r = rows m
     c = cols m
 
@@ -624,10 +624,10 @@ luFact (l_u,perm) | r <= c    = (l ,u ,p, s)
     c = cols l_u
     tu = triang r c 0 1
     tl = triang r c 0 0
-    l = takeColumns r (l_u |*| tl) |+| diagRect (konst 1 r) r r
+    l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst 1 r) r r
     u = l_u |*| tu
     (p,s) = fixPerm r perm
-    l' = (l_u |*| tl) |+| diagRect (konst 1 c) r c
+    l' = (l_u |*| tl) |+| diagRect 0 (konst 1 c) r c
     u' = takeRows c (l_u |*| tu)
     (|+|) = add
     (|*|) = mul
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
index 804c481..771739a 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -150,7 +150,7 @@ instance (ArbitraryField a) => Arbitrary (WC a) where
             c = cols m
             n = min r c
         sv' <- replicateM n (choose (1,100))
-        let s = diagRect (fromList sv') r c
+        let s = diagRect 0 (fromList sv') r c
         return $ WC (u <> real s <> trans v)
 
 #if MIN_VERSION_QuickCheck(2,0,0)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index 623b78c..a35f591 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -138,7 +138,7 @@ svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v
 
 svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where
     (u,s,v) = svdfun m
-    d = diagRect s (rows m) (cols m)
+    d = diagRect 0 s (rows m) (cols m)
 
 svdProp1b svdfun m = unitary u && unitary v where
     (u,_,v) = svdfun m
diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs
index d5c6f44..9260bd5 100644
--- a/lib/Numeric/Matrix.hs
+++ b/lib/Numeric/Matrix.hs
@@ -28,7 +28,6 @@ module Numeric.Matrix (
     -- * Basic functions
     module Data.Packed.Matrix,
     module Numeric.Vector,
-    --module Numeric.Container,
     optimiseMult,
     -- * Operators
     (<>), (<\>)