From 475ab94e1e2210d7f2e49cf42a1dfc7959bbebd3 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 17 Jul 2015 14:44:20 +0200 Subject: QR type --- packages/base/src/Internal/Algorithms.hs | 18 ++++++++++++++---- packages/base/src/Internal/Modular.hs | 2 +- packages/base/src/Numeric/LinearAlgebra.hs | 1 + packages/tests/src/Numeric/LinearAlgebra/Tests.hs | 6 +++--- 4 files changed, 19 insertions(+), 8 deletions(-) diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs index d2f17f4..ee3ddff 100644 --- a/packages/base/src/Internal/Algorithms.hs +++ b/packages/base/src/Internal/Algorithms.hs @@ -470,17 +470,25 @@ eigenvaluesSH (Her m) = eigenvaluesSH' m -------------------------------------------------------------- +-- | QR decomposition of a matrix in compact form. (The orthogonal matrix is not explicitly formed.) +data QR t = QR (Matrix t) (Vector t) + -- | QR factorization. -- -- 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" #-} unpackQR . qr' -qrRaw m = qr' m +-- | Compute the QR decomposition of a matrix in compact form. +qrRaw :: Field t => Matrix t -> QR t +qrRaw m = QR x v + where + (x,v) = qr' m -{- | generate a matrix with k orthogonal columns from the output of qrRaw --} -qrgr n (a,t) +-- | generate a matrix with k orthogonal columns from the compact QR decomposition obtained by 'qrRaw'. +-- +qrgr :: Field t => Int -> QR t -> Matrix t +qrgr n (QR 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) @@ -870,6 +878,8 @@ fixPerm' s = res $ mutable f s0 triang r c h v = (r>=h then v else 1 - v +-- | Compute the explicit LU decomposition from the compact one obtained by 'luPacked'. +luFact :: Numeric t => LU t -> (Matrix t, Matrix t, Matrix t, t) luFact (LU l_u perm) | r <= c = (l ,u ,p, s) | otherwise = (l',u',p, s) diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs index a3421a8..239c742 100644 --- a/packages/base/src/Internal/Modular.hs +++ b/packages/base/src/Internal/Modular.hs @@ -371,7 +371,7 @@ test = (ok, info) checkLU okf t = norm_Inf $ flatten (l <> u <> p - t) where - (l,u,p,_ :: Int) = luFact (LU x' p') + (l,u,p,_) = luFact (LU x' p') where (x',p') = mutable (luST okf) t diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index 7be2600..a1c0158 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -162,6 +162,7 @@ module Numeric.LinearAlgebra ( Transposable, LU(..), LDL(..), + QR(..), CGState(..), Testable(..) ) where diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs index 30480d7..ecda6c1 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs @@ -500,10 +500,10 @@ sliceTest = utest "slice test" $ and , testSlice (test_qrgr 4 tau2) qrr2 ] where - (qrr1,tau1) = qrRaw (rec :: Matrix R) - (qrr2,tau2) = qrRaw (rec :: Matrix C) + QR qrr1 tau1 = qrRaw (rec :: Matrix R) + QR qrr2 tau2 = qrRaw (rec :: Matrix C) - test_qrgr n t x = qrgr n (x,t) + test_qrgr n t x = qrgr n (QR x t) ok_qrgr x = simeq 1E-15 q q' where -- cgit v1.2.3