From b2715e91d7aef5cee1b64b641b8f173167a7145a Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 28 Dec 2009 15:47:26 +0000 Subject: additional svd functions --- CHANGES | 8 + hmatrix.cabal | 4 +- lib/Numeric/LinearAlgebra/Algorithms.hs | 116 ++++++++----- lib/Numeric/LinearAlgebra/Instances.hs | 4 +- lib/Numeric/LinearAlgebra/LAPACK.hs | 225 ++++++++++++++++++-------- lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | 159 ++++++++++++++---- lib/Numeric/LinearAlgebra/Tests.hs | 21 ++- lib/Numeric/LinearAlgebra/Tests/Instances.hs | 8 +- lib/Numeric/LinearAlgebra/Tests/Properties.hs | 78 +++++++-- 9 files changed, 466 insertions(+), 157 deletions(-) diff --git a/CHANGES b/CHANGES index 117acc4..deca1d5 100644 --- a/CHANGES +++ b/CHANGES @@ -1,3 +1,11 @@ +0.7.3.0 +======= + +- added sv, fullSVD, thinSVD, compactSVD, leftSV, rightSV + and complete interface to [d|z]gesdd. + Algorithms based on the SVD of large matrices can now be + significantly faster. + 0.7.2.0 ======= diff --git a/hmatrix.cabal b/hmatrix.cabal index a7662b9..a7dca00 100644 --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -1,5 +1,5 @@ Name: hmatrix -Version: 0.7.2.1 +Version: 0.7.3.0 License: GPL License-file: LICENSE Author: Alberto Ruiz @@ -7,7 +7,7 @@ Maintainer: Alberto Ruiz Stability: provisional Homepage: http://www.hmatrix.googlepages.com Synopsis: Linear algebra and numerical computations -Description: This library provides a purely functional interface to basic linear algebra +Description: Purely functional interface to basic linear algebra and other numerical computations, internally implemented using GSL, BLAS and LAPACK. Category: Math diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index e0cc0a0..1544d7c 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs @@ -10,10 +10,9 @@ Maintainer : Alberto Ruiz (aruiz at um dot es) Stability : provisional Portability : uses ffi -A generic interface for some common functions. Using it we can write higher level algorithms -and testing properties both for real and complex matrices. +Generic interface for the most common functions. Using it we can write higher level algorithms and testing properties for both real and complex matrices. -In any case, the specific functions for particular base types can also be explicitly +Specific functions for particular base types can also be explicitly imported from "Numeric.LinearAlgebra.LAPACK". -} @@ -34,7 +33,11 @@ module Numeric.LinearAlgebra.Algorithms ( -- * Matrix factorizations -- ** Singular value decomposition svd, - full, economy, --thin, + fullSVD, + thinSVD, + compactSVD, + sv, + leftSV, rightSV, -- ** Eigensystems eig, eigSH, eigSH', -- ** QR @@ -54,6 +57,7 @@ module Numeric.LinearAlgebra.Algorithms ( -- * Nullspace nullspacePrec, nullVector, + nullspaceSVD, -- * Norms Normed(..), NormType(..), -- * Misc @@ -63,8 +67,8 @@ module Numeric.LinearAlgebra.Algorithms ( haussholder, unpackQR, unpackHess, pinvTol, - rankSVD, ranksv, - nullspaceSVD + ranksv, + full, economy ) where @@ -80,6 +84,8 @@ import Data.Array -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) + thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) + sv' :: Matrix t -> Vector Double luPacked' :: Matrix t -> (Matrix t, [Int]) luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t linearSolve' :: Matrix t -> Matrix t -> Matrix t @@ -94,10 +100,18 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where multiply' :: Matrix t -> Matrix t -> Matrix t --- | Singular value decomposition using lapack's dgesvd or zgesvd. -svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) +-- | Full singular value decomposition using lapack's dgesdd or zgesdd. +svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) svd = svd' +-- | Thin singular value decomposition using lapack's dgesvd or zgesvd. +thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) +thinSVD = thinSVD' + +-- | Singular values using lapack's dgesvd or zgesvd. +sv :: Field t => Matrix t -> Vector Double +sv = sv' + -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. luPacked :: Field t => Matrix t -> (Matrix t, [Int]) luPacked = luPacked' @@ -163,7 +177,9 @@ multiply = multiply' instance Field Double where - svd' = svdR + svd' = svdRd + thinSVD' = thinSVDRd + sv' = svR luPacked' = luR luSolve' (l_u,perm) = lusR l_u perm linearSolve' = linearSolveR -- (luSolve . luPacked) ?? @@ -178,7 +194,9 @@ instance Field Double where multiply' = multiplyR instance Field (Complex Double) where - svd' = svdC + svd' = svdCd + thinSVD' = thinSVDCd + sv' = svC luPacked' = luC luSolve' (l_u,perm) = lusC l_u perm linearSolve' = linearSolveC @@ -232,36 +250,61 @@ inv m | square m = m `linearSolve` ident (rows m) pinv :: Field t => Matrix t -> Matrix t pinv m = linearSolveSVD m (ident (rows m)) --- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. --- --- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@. -full :: Element t - => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) + +{-# DEPRECATED full "use fullSVD instead" #-} full svdFun m = (u, d ,v) where (u,s,v) = svdFun m d = diagRect s r c r = rows m c = cols m --- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations. +-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. -- --- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. -economy :: Element t - => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) +-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@. +fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) +fullSVD m = (u,d,v) where + (u,s,v) = svd m + d = diagRect s r c + r = rows m + c = cols m + + +{-# DEPRECATED economy "use compactSVD instead" #-} economy svdFun m = (u', subVector 0 d s, v') where - r@(u,s,v) = svdFun m - d = rankSVD (1*eps) m r `max` 1 + (u,s,v) = svdFun m + d = rankSVD (1*eps) m s `max` 1 + u' = takeColumns d u + v' = takeColumns d v + +-- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations. +-- +-- If @(u,s,v) = compactSVD m@ then @m == u \<> diag s \<> trans v@. +compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) +compactSVD m = (u', subVector 0 d s, v') where + (u,s,v) = thinSVD m + d = rankSVD (1*eps) m s `max` 1 u' = takeColumns d u v' = takeColumns d v +vertical m = rows m >= cols m + +-- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd. +rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) +rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v) + | otherwise = let (_,s,v) = svd m in (s,v) + +-- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd. +leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) +leftSV m | vertical m = let (u,s,_) = svd m in (u,s) + | otherwise = let (u,s,_) = thinSVD m in (u,s) -- | Numeric rank of a matrix from the SVD decomposition. rankSVD :: Element t => Double -- ^ numeric zero (e.g. 1*'eps') -> Matrix t -- ^ input matrix m - -> (Matrix t, Vector Double, Matrix t) -- ^ 'svd' of m + -> Vector Double -- ^ 'sv' of m -> Int -- ^ rank of m -rankSVD teps m (_,s,_) = ranksv teps (max (rows m) (cols m)) (toList s) +rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s) -- | Numeric rank of a matrix from its singular values. ranksv :: Double -- ^ numeric zero (e.g. 1*'eps') @@ -316,12 +359,12 @@ pnormCV PNorm1 = norm1 . mapVector magnitude pnormCV Infinity = vectorMax . mapVector magnitude --pnormCV _ = error "pnormCV not yet defined" -pnormRM PNorm2 m = head (toList s) where (_,s,_) = svdR m +pnormRM PNorm2 m = sv m @> 0 pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m) --pnormRM _ _ = error "p norm not yet defined" -pnormCM PNorm2 m = head (toList s) where (_,s,_) = svdC m +pnormCM PNorm2 m = sv m @> 0 pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m) --pnormCM _ _ = error "p norm not yet defined" @@ -354,9 +397,9 @@ nullspaceSVD :: Field t => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'), -- or Right \"theoretical\" matrix rank. -> Matrix t -- ^ input matrix m - -> (Matrix t, Vector Double, Matrix t) -- ^ 'svd' of m + -> (Vector Double, Matrix t) -- ^ 'rightSV' of m -> [Vector t] -- ^ list of unitary vectors spanning the nullspace -nullspaceSVD hint a z@(_,_,v) = vs where +nullspaceSVD hint a (s,v) = vs where r = rows a c = cols a tol = case hint of @@ -364,9 +407,8 @@ nullspaceSVD hint a z@(_,_,v) = vs where _ -> eps k = case hint of Right t -> t - _ -> rankSVD tol a z - nsol = c - k - vs = drop k (toColumns v) + _ -> rankSVD tol a s + vs = drop k $ toRows $ ctrans v -- | The nullspace of a matrix. See also 'nullspaceSVD'. @@ -374,10 +416,7 @@ nullspacePrec :: Field t => Double -- ^ relative tolerance in 'eps' units -> Matrix t -- ^ input matrix -> [Vector t] -- ^ list of unitary vectors spanning the nullspace -nullspacePrec t m = ns where - r@(_,_,v) = svd m - k = rankSVD (t*eps) m r - ns = drop k $ toRows $ ctrans v +nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m) -- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance 1*'eps'. nullVector :: Field t => Matrix t -> Vector t @@ -405,7 +444,7 @@ multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv'). -} --pinvTol :: Double -> Matrix Double -> Matrix Double pinvTol t m = v' `mXm` diag s' `mXm` trans u' where - (u,s,v) = svdR m + (u,s,v) = thinSVDRd m sl@(g:_) = toList s s' = fromList . map rec $ sl rec x = if x < g*tol then 1 else 1/x @@ -460,15 +499,14 @@ uH (pq, tau) = (p,h) -------------------------------------------------------------------------- --- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD. +-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values. rcond :: Field t => Matrix t -> Double rcond m = last s / head s - where (_,s',_) = svd m - s = toList s' + where s = toList (sv m) -- | Number of linearly independent rows or columns. rank :: Field t => Matrix t -> Int -rank m = rankSVD eps m (svd m) +rank m = rankSVD eps m (sv m) {- expm' m = case diagonalize (complex m) of diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs index 04bbd68..a96cf15 100644 --- a/lib/Numeric/LinearAlgebra/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Instances.hs @@ -190,8 +190,8 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati --------------------------------------------------------------- -instance (Storable a) => Monoid (Vector a) where - mempty = V { dim = 0, fptr = undefined } +instance (Storable a, Num (Vector a)) => Monoid (Vector a) where + mempty = 0 { dim = 0 } mappend a b = mconcat [a,b] mconcat = j . filter ((>0).dim) where j [] = mempty diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs index 30a3d3b..2eae9b7 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK.hs +++ b/lib/Numeric/LinearAlgebra/LAPACK.hs @@ -1,4 +1,3 @@ -{-# OPTIONS_GHC #-} ----------------------------------------------------------------------------- -- | -- Module : Numeric.LinearAlgebra.LAPACK @@ -9,21 +8,34 @@ -- Stability : provisional -- Portability : portable (uses FFI) -- --- Wrappers for a few LAPACK functions (). +-- Functional interface to selected LAPACK functions (). -- ----------------------------------------------------------------------------- module Numeric.LinearAlgebra.LAPACK ( + -- * Matrix product multiplyR, multiplyC, - svdR, svdRdd, svdC, - eigC, eigR, eigS, eigH, eigS', eigH', + -- * Linear systems linearSolveR, linearSolveC, + lusR, lusC, linearSolveLSR, linearSolveLSC, linearSolveSVDR, linearSolveSVDC, - luR, luC, lusR, lusC, + -- * SVD + svR, svRd, svC, svCd, + svdR, svdRd, svdC, svdCd, + thinSVDR, thinSVDRd, thinSVDC, thinSVDCd, + rightSVR, rightSVC, leftSVR, leftSVC, + -- * Eigensystems + eigR, eigC, eigS, eigS', eigH, eigH', + -- * LU + luR, luC, + -- * Cholesky cholS, cholH, + -- * QR qrR, qrC, + -- * Hessenberg hessR, hessC, + -- * Schur schurR, schurC ) where @@ -61,28 +73,27 @@ multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Compl multiplyC a b = multiplyAux zgemmc "zgemmc" a b ----------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h svd_l_R" dgesvd :: TMMVM -foreign import ccall "LAPACK/lapack-aux.h svd_l_C" zgesvd :: TCMCMVCM -foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd" dgesdd :: TMMVM +foreign import ccall "svd_l_R" dgesvd :: TMMVM +foreign import ccall "svd_l_C" zgesvd :: TCMCMVCM +foreign import ccall "svd_l_Rdd" dgesdd :: TMMVM +foreign import ccall "svd_l_Cdd" zgesdd :: TCMCMVCM --- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix. --- --- @(u,s,v)=full svdR m@ so that @m=u \<\> s \<\> 'trans' v@. +-- | Full SVD of a real matrix using LAPACK's /dgesvd/. svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) svdR = svdAux dgesvd "svdR" . fmat --- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix. --- --- @(u,s,v)=full svdRdd m@ so that @m=u \<\> s \<\> 'trans' v@. -svdRdd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) -svdRdd = svdAux dgesdd "svdRdd" . fmat +-- | Full SVD of a real matrix using LAPACK's /dgesdd/. +svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +svdRd = svdAux dgesdd "svdRdd" . fmat --- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix. --- --- @(u,s,v)=full svdC m@ so that @m=u \<\> comp s \<\> 'trans' v@. +-- | Full SVD of a complex matrix using LAPACK's /zgesvd/. svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) svdC = svdAux zgesvd "svdC" . fmat +-- | Full SVD of a complex matrix using LAPACK's /zgesdd/. +svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +svdCd = svdAux zgesdd "svdCdd" . fmat + svdAux f st x = unsafePerformIO $ do u <- createMatrix ColumnMajor r r s <- createVector (min r c) @@ -92,7 +103,104 @@ svdAux f st x = unsafePerformIO $ do where r = rows x c = cols x + +-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'. +thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat + +-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'. +thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat + +-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'. +thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat + +-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'. +thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat + +thinSVDAux f st x = unsafePerformIO $ do + u <- createMatrix ColumnMajor r q + s <- createVector q + v <- createMatrix ColumnMajor q c + app4 f mat x mat u vec s mat v st + return (u,s,trans v) + where r = rows x + c = cols x + q = min r c + + +-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'. +svR :: Matrix Double -> Vector Double +svR = svAux dgesvd "svR" . fmat + +-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'. +svC :: Matrix (Complex Double) -> Vector Double +svC = svAux zgesvd "svC" . fmat + +-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'. +svRd :: Matrix Double -> Vector Double +svRd = svAux dgesdd "svRd" . fmat + +-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'. +svCd :: Matrix (Complex Double) -> Vector Double +svCd = svAux zgesdd "svCd" . fmat + +svAux f st x = unsafePerformIO $ do + s <- createVector q + app2 g mat x vec s st + return s + where r = rows x + c = cols x + q = min r c + g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr + + +-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'. +rightSVR :: Matrix Double -> (Vector Double, Matrix Double) +rightSVR = rightSVAux dgesvd "rightSVR" . fmat + +-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'. +rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) +rightSVC = rightSVAux zgesvd "rightSVC" . fmat + +rightSVAux f st x = unsafePerformIO $ do + s <- createVector q + v <- createMatrix ColumnMajor c c + app3 g mat x vec s mat v st + return (s,trans v) + where r = rows x + c = cols x + q = min r c + g ra ca pa = f ra ca pa 0 0 nullPtr + + +-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'. +leftSVR :: Matrix Double -> (Matrix Double, Vector Double) +leftSVR = leftSVAux dgesvd "leftSVR" . fmat + +-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'. +leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double) +leftSVC = leftSVAux zgesvd "leftSVC" . fmat + +leftSVAux f st x = unsafePerformIO $ do + u <- createMatrix ColumnMajor r r + s <- createVector q + app3 g mat x mat u vec s st + return (u,s) + where r = rows x + c = cols x + q = min r c + g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr + ----------------------------------------------------------------------------- + +foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM +foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM +foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM +foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM + eigAux f st m | r == 1 = (fromList [flatten m `at` 0], singleton 1) | otherwise = unsafePerformIO $ do @@ -104,28 +212,13 @@ eigAux f st m where r = rows m -foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM -foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM -foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM -foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM - --- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix: --- --- if @(l,v)=eigC m@ then @m \<\> v = v \<\> diag l@. --- --- The eigenvectors are the columns of v. --- The eigenvalues are not sorted. +-- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/. +-- The eigenvectors are the columns of v. The eigenvalues are not sorted. eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) eigC = eigAux zgeev "eigC" . fmat ------------------------------------------------------------------------------ - --- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix: --- --- if @(l,v)=eigR m@ then @m \<\> v = v \<\> diag l@. --- --- The eigenvectors are the columns of v. --- The eigenvalues are not sorted. +-- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/. +-- The eigenvectors are the columns of v. The eigenvalues are not sorted. eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) eigR m = (s', v'') where (s,v) = eigRaux (fmat m) @@ -155,17 +248,16 @@ fixeig _ _ = error "fixeig with impossible inputs" ----------------------------------------------------------------------------- --- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix: --- --- if @(l,v)=eigSl m@ then @m \<\> v = v \<\> diag l@. --- +-- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/. -- The eigenvectors are the columns of v. --- The eigenvalues are sorted in descending order (use eigS' for ascending order). +-- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order). eigS :: Matrix Double -> (Vector Double, Matrix Double) eigS m = (s', fliprl v) where (s,v) = eigS' (fmat m) s' = fromList . reverse . toList $ s +-- | 'eigS' in ascending order +eigS' :: Matrix Double -> (Vector Double, Matrix Double) eigS' m | r == 1 = (fromList [flatten m `at` 0], singleton 1) | otherwise = unsafePerformIO $ do @@ -177,17 +269,16 @@ eigS' m ----------------------------------------------------------------------------- --- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix: --- --- if @(l,v)=eigH m@ then @m \<\> s v = v \<\> diag l@. --- +-- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/. -- The eigenvectors are the columns of v. --- The eigenvalues are sorted in descending order (use eigH' for ascending order). +-- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order). eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) eigH m = (s', fliprl v) where (s,v) = eigH' (fmat m) s' = fromList . reverse . toList $ s +-- | 'eigH' in ascending order +eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) eigH' m | r == 1 = (fromList [realPart (flatten m `at` 0)], singleton 1) | otherwise = unsafePerformIO $ do @@ -212,11 +303,11 @@ linearSolveSQAux f st a b r = rows b c = cols b --- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition. +-- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'. linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b) --- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition. +-- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'. linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b) @@ -234,17 +325,17 @@ linearSolveAux f st a b = unsafePerformIO $ do n = cols a nrhs = cols b --- | Wrapper for LAPACK's /dgels/, which obtains the least squared error solution of an overconstrained real linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDR'. +-- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'. linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b) --- | Wrapper for LAPACK's /zgels/, which obtains the least squared error solution of an overconstrained complex linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDC'. +-- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'. linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b) --- | Wrapper for LAPACK's /dgelss/, which obtains the minimum norm solution to a real linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. +-- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. linearSolveSVDR :: Maybe Double -- ^ rcond -> Matrix Double -- ^ coefficient matrix -> Matrix Double -- ^ right hand sides (as columns) @@ -253,7 +344,7 @@ linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b) linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b) --- | Wrapper for LAPACK's /zgelss/, which obtains the minimum norm solution to a complex linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. +-- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. linearSolveSVDC :: Maybe Double -- ^ rcond -> Matrix (Complex Double) -- ^ coefficient matrix -> Matrix (Complex Double) -- ^ right hand sides (as columns) @@ -266,13 +357,11 @@ linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b) foreign import ccall "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM --- | Wrapper for LAPACK's /zpotrf/, which computes the Cholesky factorization of a --- complex Hermitian positive definite matrix. +-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/. cholH :: Matrix (Complex Double) -> Matrix (Complex Double) cholH = cholAux zpotrf "cholH" . fmat --- | Wrapper for LAPACK's /dpotrf/, which computes the Cholesky factorization of a --- real symmetric positive definite matrix. +-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/. cholS :: Matrix Double -> Matrix Double cholS = cholAux dpotrf "cholS" . fmat @@ -286,11 +375,11 @@ cholAux f st a = unsafePerformIO $ do foreign import ccall "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM foreign import ccall "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM --- | Wrapper for LAPACK's /dgeqr2/, which computes a QR factorization of a real matrix. +-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/. qrR :: Matrix Double -> (Matrix Double, Vector Double) qrR = qrAux dgeqr2 "qrR" . fmat --- | Wrapper for LAPACK's /zgeqr2/, which computes a QR factorization of a complex matrix. +-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/. qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) qrC = qrAux zgeqr2 "qrC" . fmat @@ -307,11 +396,11 @@ qrAux f st a = unsafePerformIO $ do foreign import ccall "LAPACK/lapack-aux.h hess_l_R" dgehrd :: TMVM foreign import ccall "LAPACK/lapack-aux.h hess_l_C" zgehrd :: TCMCVCM --- | Wrapper for LAPACK's /dgehrd/, which computes a Hessenberg factorization of a square real matrix. +-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/. hessR :: Matrix Double -> (Matrix Double, Vector Double) hessR = hessAux dgehrd "hessR" . fmat --- | Wrapper for LAPACK's /zgehrd/, which computes a Hessenberg factorization of a square complex matrix. +-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/. hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) hessC = hessAux zgehrd "hessC" . fmat @@ -328,11 +417,11 @@ hessAux f st a = unsafePerformIO $ do foreign import ccall "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM --- | Wrapper for LAPACK's /dgees/, which computes a Schur factorization of a square real matrix. +-- | Schur factorization of a square real matrix, using LAPACK's /dgees/. schurR :: Matrix Double -> (Matrix Double, Matrix Double) schurR = schurAux dgees "schurR" . fmat --- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix. +-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/. schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double)) schurC = schurAux zgees "schurC" . fmat @@ -347,11 +436,11 @@ schurAux f st a = unsafePerformIO $ do foreign import ccall "LAPACK/lapack-aux.h lu_l_R" dgetrf :: TMVM foreign import ccall "LAPACK/lapack-aux.h lu_l_C" zgetrf :: TCMVCM --- | Wrapper for LAPACK's /dgetrf/, which computes a LU factorization of a general real matrix. +-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/. luR :: Matrix Double -> (Matrix Double, [Int]) luR = luAux dgetrf "luR" . fmat --- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix. +-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/. luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int]) luC = luAux zgetrf "luC" . fmat @@ -370,11 +459,11 @@ type TQ a = CInt -> CInt -> PC -> a foreign import ccall "LAPACK/lapack-aux.h luS_l_R" dgetrs :: TMVMM foreign import ccall "LAPACK/lapack-aux.h luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt)))) --- | Wrapper for LAPACK's /dgetrs/, which solves a general real linear system (for several right-hand sides) from a precomputed LU decomposition. +-- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/. lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b) --- | Wrapper for LAPACK's /zgetrs/, which solves a general real linear system (for several right-hand sides) from a precomputed LU decomposition. +-- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/. lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b) diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c index f18bbed..db94cc6 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c @@ -48,7 +48,27 @@ int svd_l_R(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) { integer m = ar; integer n = ac; integer q = MIN(m,n); - REQUIRES(ur==m && uc==m && sn==q && vr==n && vc==n,BAD_SIZE); + REQUIRES(sn==q,BAD_SIZE); + REQUIRES(up==NULL || ur==m && (uc==m || uc==q),BAD_SIZE); + char* jobu = "A"; + if (up==NULL) { + jobu = "N"; + } else { + if (uc==q) { + jobu = "S"; + } + } + REQUIRES(vp==NULL || vc==n && (vr==n || vr==q),BAD_SIZE); + char* jobvt = "A"; + integer ldvt = n; + if (vp==NULL) { + jobvt = "N"; + } else { + if (vr==q) { + jobvt = "S"; + ldvt = q; + } + } DEBUGMSG("svd_l_R"); double *B = (double*)malloc(m*n*sizeof(double)); CHECK(!B,MEM); @@ -57,25 +77,21 @@ int svd_l_R(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) { integer res; // ask for optimal lwork double ans; - //printf("ask zgesvd\n"); - char* job = "A"; - dgesvd_ (job,job, + dgesvd_ (jobu,jobvt, &m,&n,B,&m, sp, up,&m, - vp,&n, + vp,&ldvt, &ans, &lwork, &res); lwork = ceil(ans); - //printf("ans = %d\n",lwork); double * work = (double*)malloc(lwork*sizeof(double)); CHECK(!work,MEM); - //printf("dgesdd\n"); - dgesvd_ (job,job, + dgesvd_ (jobu,jobvt, &m,&n,B,&m, sp, up,&m, - vp,&n, + vp,&ldvt, work, &lwork, &res); CHECK(res,res); @@ -90,25 +106,36 @@ int svd_l_Rdd(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) { integer m = ar; integer n = ac; integer q = MIN(m,n); - REQUIRES(ur==m && uc==m && sn==q && vr==n && vc==n,BAD_SIZE); + REQUIRES(sn==q,BAD_SIZE); + REQUIRES(up == NULL && vp == NULL + || ur==m && vc==n + && (uc == q && vr == q + || uc == m && vc==n),BAD_SIZE); + char* jobz = "A"; + integer ldvt = n; + if (up==NULL) { + jobz = "N"; + } else { + if (uc==q && vr == q) { + jobz = "S"; + ldvt = q; + } + } DEBUGMSG("svd_l_Rdd"); double *B = (double*)malloc(m*n*sizeof(double)); CHECK(!B,MEM); memcpy(B,ap,m*n*sizeof(double)); - integer* iwk = (integer*) malloc(8*q*sizeof(int)); + integer* iwk = (integer*) malloc(8*q*sizeof(integer)); CHECK(!iwk,MEM); integer lwk = -1; integer res; // ask for optimal lwk double ans; - //printf("ask dgesdd\n"); - dgesdd_ ("A",&m,&n,B,&m,sp,up,&m,vp,&n,&ans,&lwk,iwk,&res); - lwk = 2*ceil(ans); // ????? otherwise 50x100 rejects lwk - //printf("lwk = %d\n",lwk); + dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res); + lwk = ans; double * workv = (double*)malloc(lwk*sizeof(double)); CHECK(!workv,MEM); - //printf("dgesdd\n"); - dgesdd_ ("A",&m,&n,B,&m,sp,up,&m,vp,&n,workv,&lwk,iwk,&res); + dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res); CHECK(res,res); free(iwk); free(workv); @@ -120,17 +147,36 @@ int svd_l_Rdd(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) { // not in clapack.h -int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n, - doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, - integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, +int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n, + doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, + integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, integer *lwork, doublereal *rwork, integer *info); int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) { integer m = ar; integer n = ac; integer q = MIN(m,n); - REQUIRES(ur==m && uc==m && sn==q && vr==n && vc==n,BAD_SIZE); - DEBUGMSG("svd_l_C"); + REQUIRES(sn==q,BAD_SIZE); + REQUIRES(up==NULL || ur==m && (uc==m || uc==q),BAD_SIZE); + char* jobu = "A"; + if (up==NULL) { + jobu = "N"; + } else { + if (uc==q) { + jobu = "S"; + } + } + REQUIRES(vp==NULL || vc==n && (vr==n || vr==q),BAD_SIZE); + char* jobvt = "A"; + integer ldvt = n; + if (vp==NULL) { + jobvt = "N"; + } else { + if (vr==q) { + jobvt = "S"; + ldvt = q; + } + }DEBUGMSG("svd_l_C"); double *B = (double*)malloc(2*m*n*sizeof(double)); CHECK(!B,MEM); memcpy(B,ap,m*n*2*sizeof(double)); @@ -141,26 +187,22 @@ int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) { integer res; // ask for optimal lwork doublecomplex ans; - //printf("ask zgesvd\n"); - char* job = "A"; - zgesvd_ (job,job, + zgesvd_ (jobu,jobvt, &m,&n,(doublecomplex*)B,&m, sp, (doublecomplex*)up,&m, - (doublecomplex*)vp,&n, + (doublecomplex*)vp,&ldvt, &ans, &lwork, rwork, &res); lwork = ceil(ans.r); - //printf("ans = %d\n",lwork); doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double)); CHECK(!work,MEM); - //printf("zgesvd\n"); - zgesvd_ (job,job, + zgesvd_ (jobu,jobvt, &m,&n,(doublecomplex*)B,&m, sp, (doublecomplex*)up,&m, - (doublecomplex*)vp,&n, + (doublecomplex*)vp,&ldvt, work, &lwork, rwork, &res); @@ -171,7 +213,64 @@ int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) { OK } +int zgesdd_ (char *jobz, integer *m, integer *n, + doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, + integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, + integer *lwork, doublereal *rwork, integer* iwork, integer *info); +int svd_l_Cdd(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) { + //printf("entro\n"); + integer m = ar; + integer n = ac; + integer q = MIN(m,n); + REQUIRES(sn==q,BAD_SIZE); + REQUIRES(up == NULL && vp == NULL + || ur==m && vc==n + && (uc == q && vr == q + || uc == m && vc==n),BAD_SIZE); + char* jobz = "A"; + integer ldvt = n; + if (up==NULL) { + jobz = "N"; + } else { + if (uc==q && vr == q) { + jobz = "S"; + ldvt = q; + } + } + DEBUGMSG("svd_l_Cdd"); + doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); + CHECK(!B,MEM); + memcpy(B,ap,m*n*sizeof(doublecomplex)); + integer* iwk = (integer*) malloc(8*q*sizeof(integer)); + CHECK(!iwk,MEM); + int lrwk; + if (0 && *jobz == 'N') { + lrwk = 5*q; // does not work, crash at free below + } else { + lrwk = 5*q*q + 7*q; + } + double *rwk = (double*)malloc(lrwk*sizeof(double));; + CHECK(!rwk,MEM); + //printf("%s %ld %d\n",jobz,q,lrwk); + integer lwk = -1; + integer res; + // ask for optimal lwk + doublecomplex ans; + zgesdd_ (jobz,&m,&n,B,&m,sp,(doublecomplex*)up,&m,(doublecomplex*)vp,&ldvt,&ans,&lwk,rwk,iwk,&res); + lwk = ans.r; + //printf("lwk = %ld\n",lwk); + doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex)); + CHECK(!workv,MEM); + zgesdd_ (jobz,&m,&n,B,&m,sp,(doublecomplex*)up,&m,(doublecomplex*)vp,&ldvt,workv,&lwk,rwk,iwk,&res); + //printf("res = %ld\n",res); + CHECK(res,res); + free(workv); // printf("freed workv\n"); + free(rwk); // printf("freed rwk\n"); + free(iwk); // printf("freed iwk\n"); + free(B); // printf("freed B, salgo\n"); + OK +} //////////////////// general complex eigensystem //////////// diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs index e591285..efc8c40 100644 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ b/lib/Numeric/LinearAlgebra/Tests.hs @@ -22,6 +22,7 @@ module Numeric.LinearAlgebra.Tests( ) where import Numeric.LinearAlgebra +import Numeric.LinearAlgebra.LAPACK import Numeric.LinearAlgebra.Tests.Instances import Numeric.LinearAlgebra.Tests.Properties import Test.HUnit hiding ((~:),test,Testable) @@ -186,8 +187,24 @@ runTests n = do putStrLn "------ svd" test (svdProp1 . rM) test (svdProp1 . cM) - test (svdProp2 . rM) - test (svdProp2 . cM) + test (svdProp1a svdR) + test (svdProp1a svdC) + test (svdProp1a svdRd) + test (svdProp1a svdCd) + test (svdProp2 thinSVDR) + test (svdProp2 thinSVDC) + test (svdProp2 thinSVDRd) + test (svdProp2 thinSVDCd) + test (svdProp3 . rM) + test (svdProp3 . cM) + test (svdProp4 . rM) + test (svdProp4 . cM) + test (svdProp5a) + test (svdProp5b) + test (svdProp6a) + test (svdProp6b) + test (svdProp7 . rM) + test (svdProp7 . cM) putStrLn "------ eig" test (eigSHProp . rHer) test (eigSHProp . cHer) diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs index 9b18513..4995e39 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs @@ -143,8 +143,8 @@ instance (Field a, Arbitrary a) => Arbitrary (WC a) where r = rows m c = cols m n = min r c - sv <- replicateM n (choose (1,100)) - let s = diagRect (fromList sv) r c + sv' <- replicateM n (choose (1,100)) + let s = diagRect (fromList sv') r c return $ WC (u <> real s <> trans v) #if MIN_VERSION_QuickCheck(2,0,0) @@ -160,8 +160,8 @@ instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where Sq m <- arbitrary let (u,_,v) = svd m n = rows m - sv <- replicateM n (choose (1,100)) - let s = diag (fromList sv) + sv' <- replicateM n (choose (1,100)) + let s = diag (fromList sv') return $ SqWC (u <> real s <> trans v) #if MIN_VERSION_QuickCheck(2,0,0) diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs index d4c2770..d4dff34 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs @@ -29,7 +29,8 @@ module Numeric.LinearAlgebra.Tests.Properties ( pinvProp, detProp, nullspaceProp, - svdProp1, svdProp2, + svdProp1, svdProp1a, svdProp2, svdProp3, svdProp4, + svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, eigProp, eigSHProp, qrProp, hessProp, @@ -41,10 +42,12 @@ module Numeric.LinearAlgebra.Tests.Properties ( ) where import Numeric.LinearAlgebra +import Numeric.LinearAlgebra.LAPACK +import Debug.Trace #include "quickCheckCompat.h" --- import Debug.Trace --- debug x = trace (show x) x + +debug x = trace (show x) x -- relative error dist :: (Normed t, Num t) => t -> t -> Double @@ -71,7 +74,10 @@ a :~n~: b = dist a b < 10^^(-n) square m = rows m == cols m -unitary m = square m && m <> ctrans m |~| ident (rows m) +-- orthonormal columns +orthonormal m = ctrans m <> m |~| ident (cols m) + +unitary m = square m && orthonormal m hermitian m = square m && m |~| ctrans m @@ -119,12 +125,64 @@ nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)) r = rows m c = cols m - rank m -svdProp1 m = u <> real d <> trans v |~| m - && unitary u && unitary v - where (u,d,v) = full svd m - -svdProp2 m = (m |~| 0) `trivial` ((m |~| 0) || u <> real (diag s) <> trans v |~| m) - where (u,s,v) = economy svd m +------------------------------------------------------------------ + +-- fullSVD +svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v + where (u,d,v) = fullSVD m + +svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where + (u,s,v) = svdfun m + d = diagRect s (rows m) (cols m) + +-- thinSVD +svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) + where (u,s,v) = thinSVDfun m + +-- compactSVD +svdProp3 m = (m |~| u <> real (diag s) <> trans v + && orthonormal u && orthonormal v) + where (u,s,v) = compactSVD m + +svdProp4 m' = m |~| u <> real (diag s) <> trans v + && orthonormal u && orthonormal v + && (dim s == r || r == 0 && dim s == 1) + where (u,s,v) = compactSVD m + m = m' <-> m' + r = rank m' + +svdProp5a m = and (map (s1|~|) [s2,s3,s4,s5,s6]) where + s1 = svR m + s2 = svRd m + (_,s3,_) = svdR m + (_,s4,_) = svdRd m + (_,s5,_) = thinSVDR m + (_,s6,_) = thinSVDRd m + +svdProp5b m = and (map (s1|~|) [s2,s3,s4,s5,s6]) where + s1 = svC m + s2 = svCd m + (_,s3,_) = svdC m + (_,s4,_) = svdCd m + (_,s5,_) = thinSVDC m + (_,s6,_) = thinSVDCd m + +svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' + where (u,s,v) = svdR m + (s',v') = rightSVR m + (u',s'') = leftSVR m + +svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' + where (u,s,v) = svdC m + (s',v') = rightSVC m + (u',s'') = leftSVC m + +svdProp7 m = s |~| s' && u |~| u' && v |~| v' + where (u,s,v) = svd m + (s',v') = rightSV m + (u',s'') = leftSV m + +------------------------------------------------------------------ eigProp m = complex m <> v |~| v <> diag s where (s, v) = eig m -- cgit v1.2.3