hmatrix 0.7.2.1 → 0.8.0.0
raw patch · 10 files changed
+698/−274 lines, 10 files
Files
- CHANGES +12/−0
- hmatrix.cabal +3/−3
- lib/Numeric/LinearAlgebra/Algorithms.hs +194/−106
- lib/Numeric/LinearAlgebra/Instances.hs +2/−2
- lib/Numeric/LinearAlgebra/LAPACK.hs +207/−95
- lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +152/−48
- lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h +2/−2
- lib/Numeric/LinearAlgebra/Tests.hs +30/−2
- lib/Numeric/LinearAlgebra/Tests/Instances.hs +4/−4
- lib/Numeric/LinearAlgebra/Tests/Properties.hs +92/−12
CHANGES view
@@ -1,3 +1,15 @@+0.8.0.0+=======++- singularValues, 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.++- eigenvalues, eigenvaluesSH++- linearSolveLS, rq+ 0.7.2.0 =======
hmatrix.cabal view
@@ -1,5 +1,5 @@ Name: hmatrix-Version: 0.7.2.1+Version: 0.8.0.0 License: GPL License-file: LICENSE Author: Alberto Ruiz@@ -7,11 +7,11 @@ 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-tested-with: GHC ==6.12.1+tested-with: GHC ==6.10.4, GHC ==6.12.1 cabal-version: >=1.2 build-type: Custom
lib/Numeric/LinearAlgebra/Algorithms.hs view
@@ -10,10 +10,9 @@ 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". -}@@ -28,17 +27,23 @@ -- * Linear Systems linearSolve, luSolve,+ linearSolveLS, linearSolveSVD, inv, pinv, det, rank, rcond, -- * Matrix factorizations -- ** Singular value decomposition svd,- full, economy, --thin,+ fullSVD,+ thinSVD,+ compactSVD,+ singularValues,+ leftSV, rightSV, -- ** Eigensystems eig, eigSH, eigSH',+ eigenvalues, eigenvaluesSH, eigenvaluesSH', -- ** QR- qr,+ qr, rq, -- ** Cholesky chol, cholSH, -- ** Hessenberg@@ -54,6 +59,7 @@ -- * Nullspace nullspacePrec, nullVector,+ nullspaceSVD, -- * Norms Normed(..), NormType(..), -- * Misc@@ -63,8 +69,8 @@ haussholder, unpackQR, unpackHess, pinvTol,- rankSVD, ranksv,- nullspaceSVD+ ranksv,+ full, economy ) where @@ -80,12 +86,17 @@ -- | 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 linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t+ linearSolveLS' :: Matrix t -> Matrix t -> Matrix t eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) eigSH'' :: Matrix t -> (Vector Double, Matrix t)+ eigOnly :: Matrix t -> Vector (Complex Double)+ eigOnlySH :: Matrix t -> Vector Double cholSH' :: Matrix t -> Matrix t qr' :: Matrix t -> (Matrix t, Matrix t) hess' :: Matrix t -> (Matrix t, Matrix t)@@ -94,55 +105,191 @@ 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)+instance Field Double where+ svd' = svdRd+ thinSVD' = thinSVDRd+ sv' = svR+ luPacked' = luR+ luSolve' (l_u,perm) = lusR l_u perm+ linearSolve' = linearSolveR -- (luSolve . luPacked) ??+ linearSolveLS' = linearSolveLSR+ linearSolveSVD' = linearSolveSVDR Nothing+ ctrans' = trans+ eig' = eigR+ eigSH'' = eigS+ eigOnly = eigOnlyR+ eigOnlySH = eigOnlyS+ cholSH' = cholS+ qr' = unpackQR . qrR+ hess' = unpackHess hessR+ schur' = schurR+ multiply' = multiplyR++instance Field (Complex Double) where+ svd' = svdCd+ thinSVD' = thinSVDCd+ sv' = svC+ luPacked' = luC+ luSolve' (l_u,perm) = lusC l_u perm+ linearSolve' = linearSolveC+ linearSolveLS' = linearSolveLSC+ linearSolveSVD' = linearSolveSVDC Nothing+ ctrans' = conj . trans+ eig' = eigC+ eigOnly = eigOnlyC+ eigSH'' = eigH+ eigOnlySH = eigOnlyH+ cholSH' = cholH+ qr' = unpackQR . qrC+ hess' = unpackHess hessC+ schur' = schurC+ multiply' = multiplyC++--------------------------------------------------------------++-- | Full singular value decomposition.+svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) svd = svd' +-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.+--+-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@.+thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)+thinSVD = thinSVD'++-- | Singular values only.+singularValues :: Field t => Matrix t -> Vector Double+singularValues = sv'++-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.+--+-- 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++-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.+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.+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.+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)+++{-# 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++{-# DEPRECATED economy "use compactSVD instead" #-}+economy svdFun m = (u', subVector 0 d s, v') where+ (u,s,v) = svdFun m+ d = rankSVD (1*eps) m s `max` 1+ u' = takeColumns d u+ v' = takeColumns d v+++--------------------------------------------------------------+ -- | 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' --- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization--- obtained by 'luPacked'.+-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'. luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t luSolve = luSolve' --- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv.+-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.--- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK". linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t linearSolve = linearSolve'++-- | Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular valures which are less than 'eps' times the largest singular value. linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t linearSolveSVD = linearSolveSVD' --- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev.++-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.+linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t+linearSolveLS = linearSolveLS'++--------------------------------------------------------------++-- | Eigenvalues and eigenvectors of a general square matrix. -- -- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) eig = eig' +-- | Eigenvalues of a general square matrix.+eigenvalues :: Field t => Matrix t -> Vector (Complex Double)+eigenvalues = eigOnly+ -- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t) eigSH' = eigSH'' --- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.-cholSH :: Field t => Matrix t -> Matrix t-cholSH = cholSH'+-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric.+eigenvaluesSH' :: Field t => Matrix t -> Vector Double+eigenvaluesSH' = eigOnlySH --- | QR factorization using lapack's dgeqr2 or zgeqr2.+-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix. --+-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@+eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)+eigSH m | m `equal` ctrans m = eigSH' m+ | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"++-- | Eigenvalues of a complex hermitian or real symmetric matrix.+eigenvaluesSH :: Field t => Matrix t -> Vector Double+eigenvaluesSH m | m `equal` ctrans m = eigenvaluesSH' m+ | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"++--------------------------------------------------------------++-- | 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 = qr' --- | Hessenberg factorization using lapack's dgehrd or zgehrd.+-- | RQ factorization. --+-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.+rq :: Field t => Matrix t -> (Matrix t, Matrix t)+rq m = (r,q) where+ (q',r') = qr $ trans $ rev1 m+ r = rev2 (trans r')+ q = rev2 (trans q')+ rev1 = flipud . fliprl+ rev2 = fliprl . flipud++-- | Hessenberg factorization.+-- -- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary--- and h is in upper Hessenberg form.+-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal). hess :: Field t => Matrix t -> (Matrix t, Matrix t) hess = hess' --- | Schur factorization using lapack's dgees or zgees.+-- | Schur factorization. -- -- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary -- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is@@ -162,106 +309,52 @@ multiply = multiply' -instance Field Double where- svd' = svdR- luPacked' = luR- luSolve' (l_u,perm) = lusR l_u perm- linearSolve' = linearSolveR -- (luSolve . luPacked) ??- linearSolveSVD' = linearSolveSVDR Nothing- ctrans' = trans- eig' = eigR- eigSH'' = eigS- cholSH' = cholS- qr' = unpackQR . qrR- hess' = unpackHess hessR- schur' = schurR- multiply' = multiplyR--instance Field (Complex Double) where- svd' = svdC- luPacked' = luC- luSolve' (l_u,perm) = lusC l_u perm- linearSolve' = linearSolveC- linearSolveSVD' = linearSolveSVDC Nothing- ctrans' = conj . trans- eig' = eigC- eigSH'' = eigH- cholSH' = cholH- qr' = unpackQR . qrC- hess' = unpackHess hessC- schur' = schurC- multiply' = multiplyC----- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.------ If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@-eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)-eigSH m | m `equal` ctrans m = eigSH' m- | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"+-- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.+cholSH :: Field t => Matrix t -> Matrix t+cholSH = cholSH' --- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf.+-- | Cholesky factorization of a positive definite hermitian or symmetric matrix. -- -- If @c = chol m@ then @m == ctrans c \<> c@. chol :: Field 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" +++ square m = rows m == cols m --- | determinant of a square matrix, computed from the LU decomposition.+-- | Determinant of a square matrix. det :: Field t => Matrix t -> t det m | square m = s * (product $ toList $ takeDiag $ lup) | otherwise = error "det of nonsquare matrix" where (lup,perm) = luPacked m s = signlp (rows m) perm --- | Explicit LU factorization of a general matrix using lapack's dgetrf or zgetrf.+-- | Explicit LU factorization of a general matrix. -- -- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular, -- u is upper triangular, p is a permutation matrix and s is the signature of the permutation. lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) lu = luFact . luPacked --- | Inverse of a square matrix using lapacks' dgesv and zgesv.+-- | Inverse of a square matrix. inv :: Field t => Matrix t -> Matrix t inv m | square m = m `linearSolve` ident (rows m) | otherwise = error "inv of nonsquare matrix" --- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss.+-- | Pseudoinverse of a general matrix. 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)-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.------ 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)-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' = takeColumns d u- v' = takeColumns d v-- -- | 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 +409,12 @@ pnormCV Infinity = vectorMax . mapVector magnitude --pnormCV _ = error "pnormCV not yet defined" -pnormRM PNorm2 m = head (toList s) where (_,s,_) = svdR m+pnormRM PNorm2 m = singularValues 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 = singularValues 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 +447,9 @@ => 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,22 +457,18 @@ _ -> 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'. nullspacePrec :: Field t- => Double -- ^ relative tolerance in 'eps' units+ => Double -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps') -> 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'.+-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision. nullVector :: Field t => Matrix t -> Vector t nullVector = last . nullspacePrec 1 @@ -405,7 +494,7 @@ -} --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 +549,14 @@ -------------------------------------------------------------------------- --- | 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 (singularValues 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 (singularValues m) {- expm' m = case diagonalize (complex m) of
lib/Numeric/LinearAlgebra/Instances.hs view
@@ -190,8 +190,8 @@ --------------------------------------------------------------- -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
lib/Numeric/LinearAlgebra/LAPACK.hs view
@@ -1,4 +1,3 @@-{-# OPTIONS_GHC #-} ----------------------------------------------------------------------------- -- | -- Module : Numeric.LinearAlgebra.LAPACK@@ -9,21 +8,35 @@ -- Stability : provisional -- Portability : portable (uses FFI) ----- Wrappers for a few LAPACK functions (<http://www.netlib.org/lapack>).+-- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>). -- ----------------------------------------------------------------------------- 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',+ eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH,+ -- * LU+ luR, luC,+ -- * Cholesky cholS, cholH,+ -- * QR qrR, qrC,+ -- * Hessenberg hessR, hessC,+ -- * Schur schurR, schurC ) where @@ -61,28 +74,27 @@ 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,59 +104,152 @@ 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+ ------------------------------------------------------------------------------eigAux f st m- | r == 1 = (fromList [flatten m `at` 0], singleton 1)- | otherwise = unsafePerformIO $ do++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 :: CInt -> TMVM+foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: CInt -> TCMVCM++eigAux f st m = unsafePerformIO $ do l <- createVector r v <- createMatrix ColumnMajor r r- dummy <- createMatrix ColumnMajor 1 1- app4 f mat m mat dummy vec l mat v st+ app3 g mat m vec l mat v st return (l,v) where r = rows m+ g ra ca pa = f ra ca pa 0 0 nullPtr -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 ------------------------------------------------------------------------------+eigOnlyAux f st m = unsafePerformIO $ do+ l <- createVector r+ app2 g mat m vec l st+ return l+ where r = rows m+ g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr --- | 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.+-- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'. -- The eigenvalues are not sorted.+eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)+eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat++-- | 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)- s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))+ s' = fixeig1 s v' = toRows $ trans v v'' = fromColumns $ fixeig (toList s') v' r = rows m eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)-eigRaux m- | r == 1 = (fromList [(flatten m `at` 0):+0], singleton 1)- | otherwise = unsafePerformIO $ do+eigRaux m = unsafePerformIO $ do l <- createVector r v <- createMatrix ColumnMajor r r- dummy <- createMatrix ColumnMajor 1 1- app4 dgeev mat m mat dummy vec l mat v "eigR"+ app3 g mat m vec l mat v "eigR" return (l,v) where r = rows m+ g ra ca pa = dgeev ra ca pa 0 0 nullPtr +fixeig1 s = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))+ where r = dim s+ fixeig [] _ = [] fixeig [_] [v] = [comp v] fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)@@ -153,50 +258,59 @@ where scale = vectorMapValR Scale fixeig _ _ = error "fixeig with impossible inputs" ++-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.+-- The eigenvalues are not sorted.+eigOnlyR :: Matrix Double -> Vector (Complex Double)+eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat++ ----------------------------------------------------------------------------- --- | 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@.---+eigSHAux f st m = unsafePerformIO $ do+ l <- createVector r+ v <- createMatrix ColumnMajor r r+ app3 f mat m vec l mat v st+ return (l,v)+ where r = rows m++-- | 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' m- | r == 1 = (fromList [flatten m `at` 0], singleton 1)- | otherwise = unsafePerformIO $ do- l <- createVector r- v <- createMatrix ColumnMajor r r- app3 dsyev mat m vec l mat v "eigS"- return (l,v)- where r = rows m-------------------------------------------------------------------------------+-- | 'eigS' in ascending order+eigS' :: Matrix Double -> (Vector Double, Matrix Double)+eigS' = eigSHAux (dsyev 1) "eigS'" . fmat --- | 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' m- | r == 1 = (fromList [realPart (flatten m `at` 0)], singleton 1)- | otherwise = unsafePerformIO $ do- l <- createVector r- v <- createMatrix ColumnMajor r r- app3 zheev mat m vec l mat v "eigH"- return (l,v)- where r = rows m+-- | 'eigH' in ascending order+eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))+eigH' = eigSHAux (zheev 1) "eigH'" . fmat ++-- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'.+-- The eigenvalues are sorted in descending order.+eigOnlyS :: Matrix Double -> Vector Double+eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat++-- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'.+-- The eigenvalues are sorted in descending order.+eigOnlyH :: Matrix (Complex Double) -> Vector Double+eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat++vrev = flatten . flipud . reshape 1+ ----------------------------------------------------------------------------- foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" zgesv :: TCMCMCM@@ -212,11 +326,11 @@ 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 +348,17 @@ 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 +367,7 @@ 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 +380,11 @@ 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 +398,11 @@ 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 +419,11 @@ 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 +440,11 @@ 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 +459,11 @@ 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 +482,11 @@ 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)
lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c view
@@ -48,7 +48,27 @@ 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 @@ 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 @@ 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 @@ // 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 @@ 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,23 +213,81 @@ 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 //////////// -int eig_l_C(KCMAT(a),CMAT(u), CVEC(s),CMAT(v)) {+int eig_l_C(KCMAT(a), CMAT(u), CVEC(s),CMAT(v)) { integer n = ar;- REQUIRES(n>=2 && ac==n && (ur==1 || (ur==n && uc==n)) && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);+ REQUIRES(ac==n && sn==n, BAD_SIZE);+ REQUIRES(up==NULL || ur==n && uc==n, BAD_SIZE);+ char jobvl = up==NULL?'N':'V';+ REQUIRES(vp==NULL || vr==n && vc==n, BAD_SIZE);+ char jobvr = vp==NULL?'N':'V'; DEBUGMSG("eig_l_C");- double *B = (double*)malloc(2*n*n*sizeof(double));+ doublecomplex *B = (doublecomplex*)malloc(n*n*sizeof(doublecomplex)); CHECK(!B,MEM);- memcpy(B,ap,n*n*2*sizeof(double));-+ memcpy(B,ap,n*n*sizeof(doublecomplex)); double *rwork = (double*) malloc(2*n*sizeof(double)); CHECK(!rwork,MEM); integer lwork = -1;- char jobvl = ur==1?'N':'V';- char jobvr = vr==1?'N':'V'; integer res; // ask for optimal lwork doublecomplex ans;@@ -202,7 +302,7 @@ &res); lwork = ceil(ans.r); //printf("ans = %d\n",lwork);- doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double));+ doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); CHECK(!work,MEM); //printf("zgeev\n"); zgeev_ (&jobvl,&jobvr,@@ -226,14 +326,16 @@ int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) { integer n = ar;- REQUIRES(n>=2 && ac == n && (ur==1 || (ur==n && uc==n)) && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);+ REQUIRES(ac==n && sn==n, BAD_SIZE);+ REQUIRES(up==NULL || ur==n && uc==n, BAD_SIZE);+ char jobvl = up==NULL?'N':'V';+ REQUIRES(vp==NULL || vr==n && vc==n, BAD_SIZE);+ char jobvr = vp==NULL?'N':'V'; DEBUGMSG("eig_l_R"); double *B = (double*)malloc(n*n*sizeof(double)); CHECK(!B,MEM); memcpy(B,ap,n*n*sizeof(double)); integer lwork = -1;- char jobvl = ur==1?'N':'V';- char jobvr = vr==1?'N':'V'; integer res; // ask for optimal lwork double ans;@@ -267,13 +369,14 @@ //////////////////// symmetric real eigensystem //////////// -int eig_l_S(KDMAT(a),DVEC(s),DMAT(v)) {+int eig_l_S(int wantV,KDMAT(a),DVEC(s),DMAT(v)) { integer n = ar;- REQUIRES(n>=2 && ac == n && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);+ REQUIRES(ac==n && sn==n, BAD_SIZE);+ REQUIRES(vr==n && vc==n, BAD_SIZE);+ char jobz = wantV?'V':'N'; DEBUGMSG("eig_l_S"); memcpy(vp,ap,n*n*sizeof(double)); integer lwork = -1;- char jobz = vr==1?'N':'V'; char uplo = 'U'; integer res; // ask for optimal lwork@@ -300,15 +403,16 @@ //////////////////// hermitian complex eigensystem //////////// -int eig_l_H(KCMAT(a),DVEC(s),CMAT(v)) {+int eig_l_H(int wantV,KCMAT(a),DVEC(s),CMAT(v)) { integer n = ar;- REQUIRES(n>=2 && ac==n && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);+ REQUIRES(ac==n && sn==n, BAD_SIZE);+ REQUIRES(vr==n && vc==n, BAD_SIZE);+ char jobz = wantV?'V':'N'; DEBUGMSG("eig_l_H"); memcpy(vp,ap,2*n*n*sizeof(double)); double *rwork = (double*) malloc((3*n-2)*sizeof(double)); CHECK(!rwork,MEM); integer lwork = -1;- char jobz = vr==1?'N':'V'; char uplo = 'U'; integer res; // ask for optimal lwork@@ -322,7 +426,7 @@ &res); lwork = ceil(ans.r); //printf("ans = %d\n",lwork);- doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double));+ doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); CHECK(!work,MEM); zheev_ (&jobz,&uplo, &n,(doublecomplex*)vp,&n,
lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h view
@@ -68,8 +68,8 @@ int eig_l_C(KCMAT(a),CMAT(u),CVEC(s),CMAT(v)); int eig_l_R(KDMAT(a),DMAT(u),CVEC(s),DMAT(v)); -int eig_l_S(KDMAT(a),DVEC(s),DMAT(v));-int eig_l_H(KCMAT(a),DVEC(s),CMAT(v));+int eig_l_S(int,KDMAT(a),DVEC(s),DMAT(v));+int eig_l_H(int,KCMAT(a),DVEC(s),CMAT(v)); int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)); int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x));
lib/Numeric/LinearAlgebra/Tests.hs view
@@ -22,6 +22,7 @@ ) where import Numeric.LinearAlgebra+import Numeric.LinearAlgebra.LAPACK import Numeric.LinearAlgebra.Tests.Instances import Numeric.LinearAlgebra.Tests.Properties import Test.HUnit hiding ((~:),test,Testable)@@ -177,6 +178,11 @@ putStrLn "------ luSolve" test (linearSolveProp (luSolve.luPacked) . rSqWC) test (linearSolveProp (luSolve.luPacked) . cSqWC)+ putStrLn "------ luSolveLS"+ test (linearSolveProp linearSolveLS . rSqWC)+ test (linearSolveProp linearSolveLS . cSqWC)+ test (linearSolveProp2 linearSolveLS . rConsist)+ test (linearSolveProp2 linearSolveLS . cConsist) putStrLn "------ pinv (linearSolveSVD)" test (pinvProp . rM) test (pinvProp . cM)@@ -186,19 +192,41 @@ 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) test (eigProp . rSq) test (eigProp . cSq)+ test (eigSHProp2 . rHer)+ test (eigSHProp2 . cHer)+ test (eigProp2 . rSq)+ test (eigProp2 . cSq) putStrLn "------ nullSpace" test (nullspaceProp . rM) test (nullspaceProp . cM) putStrLn "------ qr" test (qrProp . rM) test (qrProp . cM)+ test (rqProp . rM)+ test (rqProp . cM) putStrLn "------ hess" test (hessProp . rSq) test (hessProp . cSq)
lib/Numeric/LinearAlgebra/Tests/Instances.hs view
@@ -143,8 +143,8 @@ 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 @@ 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)
lib/Numeric/LinearAlgebra/Tests/Properties.hs view
@@ -29,23 +29,26 @@ pinvProp, detProp, nullspaceProp,- svdProp1, svdProp2,- eigProp, eigSHProp,- qrProp,+ svdProp1, svdProp1a, svdProp2, svdProp3, svdProp4,+ svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,+ eigProp, eigSHProp, eigProp2, eigSHProp2,+ qrProp, rqProp, hessProp, schurProp1, schurProp2, cholProp, expmDiagProp, multProp1, multProp2,- linearSolveProp+ linearSolveProp, linearSolveProp2 ) 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 dist a b = r@@ -71,8 +74,11 @@ 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 wellCond m = rcond m > 1/100@@ -119,13 +125,66 @@ 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' && s |~| s'''+ where (u,s,v) = svd m+ (s',v') = rightSV m+ (u',s'') = leftSV m+ s''' = singularValues m++------------------------------------------------------------------+ eigProp m = complex m <> v |~| v <> diag s where (s, v) = eig m @@ -134,9 +193,25 @@ && m |~| v <> real (diag s) <> ctrans v where (s, v) = eigSH m +eigProp2 m = fst (eig m) |~| eigenvalues m++eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m++------------------------------------------------------------------+ qrProp m = q <> r |~| m && unitary q && upperTriang r where (q,r) = qr m +rqProp m = r <> q |~| m && unitary q && utr+ where (r,q) = rq m+ upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1+ where t = f-c+ utr = upptr (rows r) (cols r) * r |~| r++upperTriang' m = rows m == 1 || down |~| z+ where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))+ z = constant 0 (dim down)+ hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h where (p,h) = hess m @@ -162,3 +237,8 @@ multProp2 (a,b) = ctrans (a <> b) |~| ctrans b <> ctrans a linearSolveProp f m = f m m |~| ident (rows m)++linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)+ where q = min (rows a) (cols a)+ b = a <> x+ wc = rank a == q