added Schur factorization

5 files changed, 106 insertions, 4 deletions

diff --git a/examples/tests.hs b/examples/tests.hs
index e91b9f1..8224255 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -193,6 +193,14 @@ unitary m = square m && m <> ctrans m |~| ident (rows m)
 
 hermitian m = m |~| ctrans m
 
+upperTriang m = rows m == 1 || down == z
+    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
+upperHessenberg m = rows m < 3 || down == z
+    where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
 svdTest svd m = u <> real d <> trans v |~| m
           && unitary u && unitary v
     where (u,d,v) = full svd m
@@ -274,16 +282,24 @@ cholCTest = chol ((2><2) [1,2,2,9::Complex Double]) == (2><2) [1,2,0,2.236067977
 
 ---------------------------------------------------------------------
 
-qrTest qr m = q <> r |~| m && unitary q
+qrTest qr m = q <> r |~| m && unitary q && upperTriang r
     where (q,r) = qr m
 
 ---------------------------------------------------------------------
 
-hessTest m = m |~| p <> h <> ctrans p && unitary p
+hessTest m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
     where (p,h) = hess m
 
 ---------------------------------------------------------------------
 
+schurTest1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s
+    where (u,s) = schur m
+
+schurTest2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme
+    where (u,s) = schur m
+
+---------------------------------------------------------------------
+
 asFortran m = (rows m >|< cols m) $ toList (fdat m)
 asC m = (rows m >< cols m) $ toList (cdat m)
 
@@ -346,6 +362,9 @@ tests = do
     putStrLn "--------- hess --------"
     quickCheck (hessTest . sqm ::SqM Double->Bool)
     quickCheck (hessTest . sqm ::SqM (Complex Double) -> Bool)
+    putStrLn "--------- schur --------"
+    quickCheck (schurTest2 . sqm ::SqM Double->Bool)
+    quickCheck (schurTest1 . sqm ::SqM (Complex Double) -> Bool)
     putStrLn "--------- nullspace ------"
     quickCheck (nullspaceTest :: RM -> Bool)
     quickCheck (nullspaceTest :: CM -> Bool)
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 1345975..7f8e84c 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -35,6 +35,8 @@ module Numeric.LinearAlgebra.Algorithms (
     chol,
 -- ** Hessenberg
     hess,
+-- ** Schur
+    schur,
 -- * Nullspace
     nullspacePrec,
     nullVector,
@@ -79,6 +81,8 @@ class (Linear Matrix t) => GenMat t where
     qr          :: Matrix t -> (Matrix t, Matrix t)
     -- | Hessenberg factorization using lapack's dgehrd or zgehrd.
     hess        :: Matrix t -> (Matrix t, Matrix t)
+    -- | Schur factorization using lapack's dgees or zgees.
+    schur       :: Matrix t -> (Matrix t, Matrix t)
     -- | Conjugate transpose.
     ctrans :: Matrix t -> Matrix t
 
@@ -93,6 +97,7 @@ instance GenMat Double where
     cholSH = cholS
     qr = GSL.unpackQR . qrR
     hess = unpackHess hessR
+    schur = schurR
 
 instance GenMat (Complex Double) where
     svd = svdC
@@ -105,6 +110,7 @@ instance GenMat (Complex Double) where
     cholSH = cholH
     qr = unpackQR . qrC
     hess = unpackHess hessC
+    schur = schurC
 
 -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev. If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
 eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t)
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index a84a17e..628d4f8 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -21,7 +21,8 @@ module Numeric.LinearAlgebra.LAPACK (
     linearSolveSVDR, linearSolveSVDC,
     cholS, cholH,
     qrR, qrC,
-    hessR, hessC
+    hessR, hessC,
+    schurR, schurC
 ) where
 
 import Data.Packed.Internal
@@ -351,7 +352,7 @@ hessR a = unsafePerformIO $ do
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h hess_l_C" zgehrd :: TCMCVCM
 
--- | Wrapper for LAPACK's /zgeqr2/, which computes a Hessenberg factorization of a square complex matrix.
+-- | Wrapper for LAPACK's /zgehrd/, which computes a Hessenberg factorization of a square complex matrix.
 hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
 hessC a = unsafePerformIO $ do
     r <- createMatrix ColumnMajor m n
@@ -362,3 +363,26 @@ hessC a = unsafePerformIO $ do
         n = cols a
         mn = min m n
 
+-----------------------------------------------------------------------------------
+foreign import ccall "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM
+
+-- | Wrapper for LAPACK's /dgees/, which computes a Schur factorization of a square real matrix.
+schurR :: Matrix Double -> (Matrix Double, Matrix Double)
+schurR a = unsafePerformIO $ do
+    u <- createMatrix ColumnMajor n n
+    s <- createMatrix ColumnMajor n n
+    dgees // mat fdat a // mat dat u // mat dat s // check "schurR" [fdat a]
+    return (u,s)
+  where n = rows a
+
+-----------------------------------------------------------------------------------
+foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM
+
+-- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix.
+schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
+schurC a = unsafePerformIO $ do
+    u <- createMatrix ColumnMajor n n
+    s <- createMatrix ColumnMajor n n
+    zgees // mat fdat a // mat dat u // mat dat s // check "schurC" [fdat a]
+    return (u,s)
+  where n = rows a
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 04ef416..cab3f5b 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -701,3 +701,52 @@ int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
     free(WORK);
     OK
 }
+
+//////////////////// Schur factorization /////////////////////////
+
+int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) {
+    integer m = ar;
+    integer n = ac;
+    REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
+    DEBUGMSG("schur_l_R");
+    memcpy(sp,ap,n*n*sizeof(double));
+    integer lwork = 6*n; // fixme
+    double *WORK = (double*)malloc(lwork*sizeof(double));
+    double *WR = (double*)malloc(n*sizeof(double));
+    double *WI = (double*)malloc(n*sizeof(double));
+    // WR and WI not really required in this call
+    logical *BWORK = (logical*)malloc(n*sizeof(logical));
+    integer res;
+    integer sdim;
+    dgees_ ("V","N",NULL,&n,sp,&n,&sdim,WR,WI,up,&n,WORK,&lwork,BWORK,&res);
+    CHECK(res,res);
+    free(WR);
+    free(WI);
+    free(BWORK);
+    free(WORK);
+    OK
+}
+
+int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) {
+    integer m = ar;
+    integer n = ac;
+    REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
+    DEBUGMSG("schur_l_C");
+    memcpy(sp,ap,n*n*sizeof(doublecomplex));
+    integer lwork = 6*n; // fixme
+    doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex));
+    // W not really required in this call
+    logical *BWORK = (logical*)malloc(n*sizeof(logical));
+    double *RWORK = (double*)malloc(n*sizeof(double));
+    integer res;
+    integer sdim;
+    zgees_ ("V","N",NULL,&n,(doublecomplex*)sp,&n,&sdim,W,
+                            (doublecomplex*)up,&n,
+                            WORK,&lwork,RWORK,BWORK,&res);
+    CHECK(res,res);
+    free(W);
+    free(BWORK);
+    free(WORK);
+    OK
+}
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index 52ac41e..e5d74d7 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -48,3 +48,7 @@ int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r));
 int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r));
 
 int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r));
+
+int schur_l_R(KDMAT(a), DMAT(u), DMAT(s));
+
+int schur_l_C(KCMAT(a), CMAT(u), CMAT(s));