From 42bec1ac9911131b552f66779203eb599a86563d Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Tue, 2 Oct 2007 18:59:50 +0000 Subject: lapack real and complex unpacked QR --- lib/Numeric/GSL/Matrix.hs | 26 +++++++++++++++++- lib/Numeric/GSL/Vector.hs | 1 - lib/Numeric/GSL/gsl-aux.c | 27 +++++++++++++++++++ lib/Numeric/GSL/gsl-aux.h | 3 +++ lib/Numeric/LinearAlgebra/Algorithms.hs | 38 ++++++++++++++++++++++----- lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | 5 ++-- 6 files changed, 88 insertions(+), 12 deletions(-) (limited to 'lib/Numeric') 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; -- cgit v1.2.3