added Schur factorization

4 files changed, 85 insertions, 2 deletions

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));