diagRect, svdC

3 files changed, 23 insertions, 7 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index bd333d4..4383e79 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -315,3 +315,10 @@ diagG v = reshape c $ fromList $ [ l!!(i-1) * delta k i | k <- [1..c], i <- [1..
           l = toList v
           delta i j | i==j      = 1
                     | otherwise = 0
+
+diagRect s r c
+    | dim s < min r c = error "diagRect"
+    | r == c    = diag s
+    | r < c     = joinHoriz [diag s , zeros (r,c-r)]
+    | otherwise = joinVert  [diag s , zeros (r-c,c)]
+    where zeros (r,c) = reshape c $ constant (r*c) 0
diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs
index 67e49af..088b6d7 100644
--- a/lib/LAPACK.hs
+++ b/lib/LAPACK.hs
@@ -14,7 +14,7 @@
 
 module LAPACK (
     --module LAPACK.Internal
-    svdR, svdR',
+    svdR, svdR', svdC, svdC',
     eigC,
     linearSolveLSR
 ) where
diff --git a/lib/LAPACK/Internal.hs b/lib/LAPACK/Internal.hs
index 2569215..e39dd10 100644
--- a/lib/LAPACK/Internal.hs
+++ b/lib/LAPACK/Internal.hs
@@ -31,12 +31,7 @@ foreign import ccall "lapack-aux.h svd_l_R"
 --
 -- @(u,s,v)=svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
 svdR :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double)
-svdR x@M {rows = r, cols = c} = (u, s, v)
-    where (u,s',v) = svdR' x
-          s | r == c    = diag s'
-            | r < c     = joinHoriz [diag s' , zeros (r,c-r)]
-            | otherwise = joinVert  [diag s' , zeros (r-c,c)]
-          zeros (r,c) = reshape c $ constant (r*c) 0
+svdR x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdR' x
 
 svdR' x@M {rows = r, cols = c} = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r r
@@ -55,6 +50,20 @@ foreign import ccall "lapack-aux.h svd_l_Rdd"
 foreign import ccall "lapack-aux.h svd_l_C"
     zgesvd :: (Complex Double) ::> (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
 
+-- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix.
+--
+-- @(u,s,v)=svdC m@ so that @m=u \<\> s \<\> 'trans' v@.
+svdC :: Matrix (Complex Double)
+     -> (Matrix (Complex Double), Matrix Double, Matrix (Complex Double))
+svdC x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdC' x
+
+svdC' x@M {rows = r, cols = c} = unsafePerformIO $ do
+    u <- createMatrix ColumnMajor r r
+    s <- createVector (min r c)
+    v <- createMatrix ColumnMajor c c
+    zgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdC" [fdat x]
+    return (u,s,trans v)
+
 -----------------------------------------------------------------------------
 -- zgeev
 foreign import ccall "lapack-aux.h eig_l_C"