lapack QR (unpacked)

3 files changed, 69 insertions, 1 deletions

diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs
index e84647b..43a313b 100644
--- a/lib/LAPACK.hs
+++ b/lib/LAPACK.hs
@@ -19,7 +19,8 @@ module LAPACK (
     linearSolveR, linearSolveC,
     linearSolveLSR, linearSolveLSC,
     linearSolveSVDR, linearSolveSVDC,
-    cholS, cholH
+    cholS, cholH,
+    qrR, qrC
 ) where
 
 import Data.Packed.Internal
@@ -304,3 +305,31 @@ cholS a = unsafePerformIO $ do
     return r
   where n = rows a
 
+-----------------------------------------------------------------------------------
+foreign import ccall "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM
+
+-- | Wrapper for LAPACK's /dgeqr2/,which computes a QR factorization of a real matrix.
+qrR :: Matrix Double -> (Matrix Double, Vector Double)
+qrR a = unsafePerformIO $ do
+    r <- createMatrix ColumnMajor m n
+    tau <- createVector mn
+    dgeqr2 // mat fdat a // vec tau // mat dat r // check "qrR" [fdat a]
+    return (r,tau)
+  where m = rows a
+        n = cols a
+        mn = min m n
+
+-----------------------------------------------------------------------------------
+foreign import ccall "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM
+
+-- | Wrapper for LAPACK's /zgeqr2/,which computes a QR factorization of a complex matrix.
+qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
+qrC a = unsafePerformIO $ do
+    r <- createMatrix ColumnMajor m n
+    tau <- createVector mn
+    zgeqr2 // mat fdat a // vec tau // mat dat r // check "qrC" [fdat a]
+    return (r,tau)
+  where m = rows a
+        n = cols a
+        mn = min m n
+
diff --git a/lib/LAPACK/lapack-aux.c b/lib/LAPACK/lapack-aux.c
index 5781cd1..4ef9a6e 100644
--- a/lib/LAPACK/lapack-aux.c
+++ b/lib/LAPACK/lapack-aux.c
@@ -619,3 +619,38 @@ int chol_l_S(KDMAT(a),DMAT(l)) {
     }
     OK
 }
+
+//////////////////// QR factorization /////////////////////////
+// TO DO: unpack
+
+int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE);
+    DEBUGMSG("qr_l_R");
+    double *WORK = (double*)malloc(m*sizeof(double));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*n*sizeof(double));
+    integer res;
+    dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
+
+int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE);
+    DEBUGMSG("qr_l_C");
+    doublecomplex *WORK = (doublecomplex*)malloc(m*sizeof(doublecomplex));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*n*sizeof(doublecomplex));
+    integer res;
+    zgeqr2_ (&m,&n,(doublecomplex*)rp,&m,(doublecomplex*)taup,WORK,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
diff --git a/lib/LAPACK/lapack-aux.h b/lib/LAPACK/lapack-aux.h
index c5bf3cc..ea71847 100644
--- a/lib/LAPACK/lapack-aux.h
+++ b/lib/LAPACK/lapack-aux.h
@@ -40,3 +40,7 @@ int linearSolveSVDC_l(double,KCMAT(a),KCMAT(b),CMAT(x));
 int chol_l_H(KCMAT(a),CMAT(r));
 
 int chol_l_S(KDMAT(a),DMAT(r));
+
+int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r));
+
+int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r));