4 files changed, 55 insertions, 2 deletions

diff --git a/lib/Numeric/GSL/Matrix.hs b/lib/Numeric/GSL/Matrix.hs
index eb1931a..5a5c19e 100644
--- a/lib/Numeric/GSL/Matrix.hs
+++ b/lib/Numeric/GSL/Matrix.hs
@@ -16,7 +16,7 @@
 module Numeric.GSL.Matrix(
     eigSg, eigHg,
     svdg,
-    qr,
+    qr, qrPacked, unpackQR,
     cholR, -- cholC,
     luSolveR, luSolveC,
     luR, luC
@@ -149,6 +149,30 @@ qr x = unsafePerformIO $ do
         c = cols x
 foreign import ccall "gsl-aux.h QR" c_qr :: TMMM
 
+qrPacked :: Matrix Double -> (Matrix Double, Vector Double)
+qrPacked x = unsafePerformIO $ do
+    qr <- createMatrix RowMajor r c
+    tau <- createVector (min r c)
+    c_qrPacked // mat cdat x // mat dat qr // vec tau // check "qrUnpacked" [cdat x]
+    return (qr,tau)
+  where r = rows x
+        c = cols x
+foreign import ccall "gsl-aux.h QRpacked" c_qrPacked :: TMMV
+
+unpackQR :: (Matrix Double, Vector Double) -> (Matrix Double, Matrix Double)
+unpackQR (qr,tau) = unsafePerformIO $ do
+    q <- createMatrix RowMajor r r
+    rot <- createMatrix RowMajor r c
+    c_qrUnpack // mat cdat qr // vec tau // mat dat q // mat dat rot // check "qrUnpack" [cdat qr,tau]
+    return (q,rot)
+  where r = rows qr
+        c = cols qr
+foreign import ccall "gsl-aux.h QRunpack" c_qrUnpack :: TMVMM
+
+
+type TMMV = Int -> Int -> PD -> TMV
+type TMVMM = Int -> Int -> PD -> Int -> PD -> TMM
+
 {- | Cholesky decomposition of a symmetric positive definite real matrix using /gsl_linalg_cholesky_decomp/.
 
 @\> chol $ (2><2) [1,2,
diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index ef3d5e8..41efdc0 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -12,7 +12,6 @@
 -- Vector operations
 --
 -----------------------------------------------------------------------------
--- #hide
 
 module Numeric.GSL.Vector (
     FunCodeS(..), toScalarR,
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index c602d5e..83c3bf8 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -471,6 +471,33 @@ int QR(KRMAT(x),RMAT(q),RMAT(r)) {
     OK
 }
 
+int QRpacked(KRMAT(x),RMAT(qr),RVEC(tau)) {
+    //REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE);
+    DEBUGMSG("QRpacked");
+    KDMVIEW(x);
+    DMVIEW(qr);
+    DVVIEW(tau);
+    //gsl_matrix * a = gsl_matrix_alloc(xr,xc);
+    //gsl_vector * tau = gsl_vector_alloc(MIN(xr,xc));
+    gsl_matrix_memcpy(M(qr),M(x));
+    int res = gsl_linalg_QR_decomp(M(qr),V(tau));
+    CHECK(res,res);
+    OK
+}
+
+
+int QRunpack(KRMAT(xqr),KRVEC(tau),RMAT(q),RMAT(r)) {
+    //REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE);
+    DEBUGMSG("QRunpack");
+    KDMVIEW(xqr);
+    KDVVIEW(tau);
+    DMVIEW(q);
+    DMVIEW(r);
+    gsl_linalg_QR_unpack(M(xqr),V(tau),M(q),M(r));
+    OK
+}
+
+
 int cholR(KRMAT(x),RMAT(l)) {
     REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE);
     DEBUGMSG("cholR");
diff --git a/lib/Numeric/GSL/gsl-aux.h b/lib/Numeric/GSL/gsl-aux.h
index 3ccac25..cf22024 100644
--- a/lib/Numeric/GSL/gsl-aux.h
+++ b/lib/Numeric/GSL/gsl-aux.h
@@ -38,6 +38,9 @@ int eigensystemC(KCMAT(x),RVEC(l),CMAT(v));
 
 int QR(KRMAT(x),RMAT(q),RMAT(r));
 
+int QRpacked(KRMAT(x),RMAT(qr),RVEC(tau));
+int QRunpack(KRMAT(qr),KRVEC(tau),RMAT(q),RMAT(r));
+
 int cholR(KRMAT(x),RMAT(l));
 
 int cholC(KCMAT(x),CMAT(l));