lapack real and complex unpacked QR

6 files changed, 88 insertions, 12 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));
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 3513b18..5492e07 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -38,17 +38,18 @@ module Numeric.LinearAlgebra.Algorithms (
     eps, i,
     ctrans,
     Normed(..), NormType(..),
-    GenMat(linearSolveSVD,lu,eigSH')
+    GenMat(linearSolveSVD,lu,eigSH'), unpackQR
 ) where
 
 
 import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj)
 import Data.Packed
-import Numeric.GSL.Matrix(luR,luC,qr)
+import qualified Numeric.GSL.Matrix as GSL
 import Numeric.GSL.Vector
 import Numeric.LinearAlgebra.LAPACK as LAPACK
 import Complex
 import Numeric.LinearAlgebra.Linear
+import Data.List(foldl1')
 
 -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
 class (Linear Matrix t) => GenMat t where
@@ -59,28 +60,31 @@ class (Linear Matrix t) => GenMat t where
     eig         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
     eigSH'      :: Matrix t -> (Vector Double, Matrix t)
     cholSH      :: Matrix t -> Matrix t
+    qr          :: Matrix t -> (Matrix t, Matrix t)
     -- | conjugate transpose
     ctrans :: Matrix t -> Matrix t
 
 instance GenMat Double where
     svd = svdR
-    lu  = luR
+    lu  = GSL.luR
     linearSolve = linearSolveR
     linearSolveSVD = linearSolveSVDR Nothing
     ctrans = trans
     eig = eigR
     eigSH' = eigS
     cholSH = cholS
+    qr = GSL.unpackQR . qrR
 
 instance GenMat (Complex Double) where
     svd = svdC
-    lu  = luC
+    lu  = GSL.luC
     linearSolve = linearSolveC
     linearSolveSVD = linearSolveSVDC Nothing
     ctrans = conjTrans
     eig = eigC
     eigSH' = eigH
     cholSH = cholH
+    qr = unpackQR . qrC
 
 -- | eigensystem of complex hermitian or real symmetric matrix
 eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t)
@@ -92,9 +96,6 @@ chol :: GenMat t => Matrix t ->  Matrix t
 chol m | m `equal` ctrans m = cholSH m
        | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
 
-qr :: Matrix Double -> (Matrix Double, Matrix Double)
-qr = Numeric.GSL.Matrix.qr
-
 square m = rows m == cols m
 
 det :: GenMat t => Matrix t -> t
@@ -257,3 +258,26 @@ pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
     d = dim s
     u' = takeColumns d u
     v' = takeColumns d v
+
+---------------------------------------------------------------------
+
+-- many thanks, quickcheck!
+
+haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w))
+    where w = asColumn v
+
+unpackQR (pq, tau) = (q,r)
+    where cs = toColumns pq
+          m = rows pq
+          n = cols pq
+          mn = min m n
+          r = fromColumns $ zipWith zt ([m-1, m-2 .. 1] ++ repeat 0) cs
+          vs = zipWith zh [1..mn] cs
+          hs = zipWith haussholder (toList tau) vs
+          q = foldl1' mXm hs
+
+          zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs)
+              where xs = toList v
+
+          zt 0 v = v
+          zt k v = join [subVector 0 (dim v - k) v, constant 0 k]
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 4ef9a6e..bf0cfba 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -621,7 +621,6 @@ int chol_l_S(KDMAT(a),DMAT(l)) {
 }
 
 //////////////////// QR factorization /////////////////////////
-// TO DO: unpack
 
 int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
     integer m = ar;
@@ -629,7 +628,7 @@ int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
     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));
+    double *WORK = (double*)malloc(n*sizeof(double));
     CHECK(!WORK,MEM);
     memcpy(rp,ap,m*n*sizeof(double));
     integer res;
@@ -645,7 +644,7 @@ int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
     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));
+    doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex));
     CHECK(!WORK,MEM);
     memcpy(rp,ap,m*n*sizeof(doublecomplex));
     integer res;