4 files changed, 80 insertions, 10 deletions

diff --git a/lib/Numeric/HMatrix.hs b/lib/Numeric/HMatrix.hs
index f49ea53..09ad518 100644
--- a/lib/Numeric/HMatrix.hs
+++ b/lib/Numeric/HMatrix.hs
@@ -86,7 +86,7 @@ module Numeric.HMatrix (
     geigSH',
 
     -- * QR
-    qr, rq,
+    qr, rq, qrRaw, qrgr,
 
     -- * Cholesky
     chol, cholSH, mbCholSH,
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index de9cb00..8c4b610 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -49,7 +49,7 @@ module Numeric.LinearAlgebra.Algorithms (
     eigenvalues, eigenvaluesSH, eigenvaluesSH',
     geigSH',
 -- ** QR
-    qr, rq,
+    qr, rq, qrRaw, qrgr,
 -- ** Cholesky
     chol, cholSH, mbCholSH,
 -- ** Hessenberg
@@ -115,7 +115,8 @@ class (Product t,
     eigOnlySH    :: Matrix t -> Vector Double
     cholSH'      :: Matrix t -> Matrix t
     mbCholSH'    :: Matrix t -> Maybe (Matrix t)
-    qr'          :: Matrix t -> (Matrix t, Matrix t)
+    qr'          :: Matrix t -> (Matrix t, Vector t)
+    qrgr'        :: Int -> (Matrix t, Vector t) -> Matrix t
     hess'        :: Matrix t -> (Matrix t, Matrix t)
     schur'       :: Matrix t -> (Matrix t, Matrix t)
 
@@ -136,7 +137,8 @@ instance Field Double where
     eigOnlySH = eigOnlyS
     cholSH' = cholS
     mbCholSH' = mbCholS
-    qr' = unpackQR . qrR
+    qr' = qrR
+    qrgr' = qrgrR
     hess' = unpackHess hessR
     schur' = schurR
 
@@ -161,7 +163,8 @@ instance Field (Complex Double) where
     eigOnlySH = eigOnlyH
     cholSH' = cholH
     mbCholSH' = mbCholH
-    qr' = unpackQR . qrC
+    qr' = qrC
+    qrgr' = qrgrC
     hess' = unpackHess hessC
     schur' = schurC
 
@@ -285,7 +288,15 @@ eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m
 --
 -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
 qr :: Field t => Matrix t -> (Matrix t, Matrix t)
-qr = {-# SCC "qr" #-} qr'
+qr = {-# SCC "qr" #-} unpackQR . qr'
+
+qrRaw m = qr' m
+
+{- | generate a matrix with k orthogonal columns from the output of qrRaw
+-}
+qrgr n (a,t)
+    | dim t > min (cols a) (rows a) || n < 0 || n > dim t = error "qrgr expects k <= min(rows,cols)"
+    | otherwise = qrgr' n (a,t)
 
 -- | RQ factorization.
 --
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index a12f6d8..11394a6 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -35,7 +35,7 @@ module Numeric.LinearAlgebra.LAPACK (
     -- * Cholesky
     cholS, cholH, mbCholS, mbCholH,
     -- * QR
-    qrR, qrC,
+    qrR, qrC, qrgrR, qrgrC,
     -- * Hessenberg
     hessR, hessC,
     -- * Schur
@@ -444,9 +444,27 @@ qrAux f st a = unsafePerformIO $ do
     tau <- createVector mn
     app3 f mat a vec tau mat r st
     return (r,tau)
-  where m = rows a
-        n = cols a
-        mn = min m n
+  where
+    m = rows a
+    n = cols a
+    mn = min m n
+
+foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM
+foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM
+
+-- | build rotation from reflectors
+qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double
+qrgrR = qrgrAux dorgqr "qrgrR"
+-- | build rotation from reflectors
+qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double)
+qrgrC = qrgrAux zungqr "qrgrC"
+
+qrgrAux f st n (a, tau) = unsafePerformIO $ do
+    res <- createMatrix ColumnMajor (rows a) n
+    app3 f mat (fmat a) vec (subVector 0 n tau') mat res st
+    return res
+  where
+    tau' = vjoin [tau, constantD 0 n]
 
 -----------------------------------------------------------------------------------
 foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 51972eb..e5e45ef 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -930,6 +930,47 @@ int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
     OK
 }
 
+/* Subroutine */ int dorgqr_(integer *m, integer *n, integer *k, doublereal *
+        a, integer *lda, doublereal *tau, doublereal *work, integer *lwork,
+        integer *info);
+
+int c_dorgqr(KDMAT(a), KDVEC(tau), DMAT(r)) {
+    integer m = ar;
+    integer n = MIN(ac,ar);
+    integer k = taun;
+    DEBUGMSG("c_dorgqr");
+    integer lwork = 8*n; // FIXME
+    double *WORK = (double*)malloc(lwork*sizeof(double));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*k*sizeof(double));
+    integer res;
+    dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
+
+/* Subroutine */ int zungqr_(integer *m, integer *n, integer *k,
+        doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
+        work, integer *lwork, integer *info);
+
+int c_zungqr(KCMAT(a), KCVEC(tau), CMAT(r)) {
+    integer m = ar;
+    integer n = MIN(ac,ar);
+    integer k = taun;
+    DEBUGMSG("z_ungqr");
+    integer lwork = 8*n; // FIXME
+    doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*k*sizeof(doublecomplex));
+    integer res;
+    zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
+
+
 //////////////////// Hessenberg factorization /////////////////////////
 
 /* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi,