summaryrefslogtreecommitdiff
path: root/lib/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs31
1 files changed, 13 insertions, 18 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 3757f33..425532d 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -64,15 +64,15 @@ class (Linear Matrix t) => GenMat t where
64 -- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK". 64 -- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".
65 linearSolve :: Matrix t -> Matrix t -> Matrix t 65 linearSolve :: Matrix t -> Matrix t -> Matrix t
66 linearSolveSVD :: Matrix t -> Matrix t -> Matrix t 66 linearSolveSVD :: Matrix t -> Matrix t -> Matrix t
67 -- | eigensystem of general square matrix using lapack's dgeev or zgeev. If @(s,v) = eig m@ then @m <> v =~= v <> diag s@ 67 -- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev. If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
68 eig :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) 68 eig :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
69 -- | similar to eigSH without checking that the input matrix is hermitian or symmetric. 69 -- | Similar to eigSH without checking that the input matrix is hermitian or symmetric.
70 eigSH' :: Matrix t -> (Vector Double, Matrix t) 70 eigSH' :: Matrix t -> (Vector Double, Matrix t)
71 -- | similar to chol without checking that the input matrix is hermitian or symmetric. 71 -- | Similar to chol without checking that the input matrix is hermitian or symmetric.
72 cholSH :: Matrix t -> Matrix t 72 cholSH :: Matrix t -> Matrix t
73 -- | QR factorization using lapack's dgeqr2 or zgeqr2. 73 -- | QR factorization using lapack's dgeqr2 or zgeqr2.
74 qr :: Matrix t -> (Matrix t, Matrix t) 74 qr :: Matrix t -> (Matrix t, Matrix t)
75 -- | conjugate transpose. 75 -- | Conjugate transpose.
76 ctrans :: Matrix t -> Matrix t 76 ctrans :: Matrix t -> Matrix t
77 77
78instance GenMat Double where 78instance GenMat Double where
@@ -97,7 +97,7 @@ instance GenMat (Complex Double) where
97 cholSH = cholH 97 cholSH = cholH
98 qr = unpackQR . qrC 98 qr = unpackQR . qrC
99 99
100-- | eigensystem of complex hermitian or real symmetric matrix using lapack's dsyev or zheev. If @(s,v) = eigSH m@ then @m =~= v <> diag s <> ctrans v@ 100-- | 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@
101eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) 101eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t)
102eigSH m | m `equal` ctrans m = eigSH' m 102eigSH m | m `equal` ctrans m = eigSH' m
103 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" 103 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
@@ -123,7 +123,7 @@ inv m | square m = m `linearSolve` ident (rows m)
123pinv :: GenMat t => Matrix t -> Matrix t 123pinv :: GenMat t => Matrix t -> Matrix t
124pinv m = linearSolveSVD m (ident (rows m)) 124pinv m = linearSolveSVD m (ident (rows m))
125 125
126-- | 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@. 126-- | 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@.
127full :: Field t 127full :: Field t
128 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) 128 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)
129full svd m = (u, d ,v) where 129full svd m = (u, d ,v) where
@@ -132,7 +132,7 @@ full svd m = (u, d ,v) where
132 r = rows m 132 r = rows m
133 c = cols m 133 c = cols m
134 134
135-- | 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@. 135-- | 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@.
136economy :: Field t 136economy :: Field t
137 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) 137 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t)
138economy svd m = (u', subVector 0 d s, v') where 138economy svd m = (u', subVector 0 d s, v') where
@@ -149,24 +149,24 @@ economy svd m = (u', subVector 0 d s, v') where
149 v' = takeColumns d v 149 v' = takeColumns d v
150 150
151 151
152-- | The machine precision of a Double: @eps == 2.22044604925031e-16@ (the value used by GNU-Octave). 152-- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave).
153eps :: Double 153eps :: Double
154eps = 2.22044604925031e-16 154eps = 2.22044604925031e-16
155 155
156-- | The imaginary unit: @i == 0.0 :+ 1.0@ 156-- | The imaginary unit: @i = 0.0 :+ 1.0@
157i :: Complex Double 157i :: Complex Double
158i = 0:+1 158i = 0:+1
159 159
160 160
161-- | matrix product 161-- matrix product
162mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t 162mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t
163mXm = multiply 163mXm = multiply
164 164
165-- | matrix - vector product 165-- matrix - vector product
166mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t 166mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t
167mXv m v = flatten $ m `mXm` (asColumn v) 167mXv m v = flatten $ m `mXm` (asColumn v)
168 168
169-- | vector - matrix product 169-- vector - matrix product
170vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t 170vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t
171vXm v m = flatten $ (asRow v) `mXm` m 171vXm v m = flatten $ (asRow v) `mXm` m
172 172
@@ -201,10 +201,6 @@ pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (liftVector
201pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` constant 1 (cols m) 201pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` constant 1 (cols m)
202--pnormCM _ _ = error "p norm not yet defined" 202--pnormCM _ _ = error "p norm not yet defined"
203 203
204-- -- | computes the p-norm of a matrix or vector (with the same definitions as GNU-octave). pnorm 0 denotes \\inf-norm. See also 'norm'.
205--pnorm :: (Container t, GenMat a) => Int -> t a -> Double
206--pnorm = pnormG
207
208class Normed t where 204class Normed t where
209 pnorm :: NormType -> t -> Double 205 pnorm :: NormType -> t -> Double
210 norm :: t -> Double 206 norm :: t -> Double
@@ -226,7 +222,7 @@ instance Normed (Matrix (Complex Double)) where
226 222
227-- | The nullspace of a matrix from its SVD decomposition. 223-- | The nullspace of a matrix from its SVD decomposition.
228nullspacePrec :: GenMat t 224nullspacePrec :: GenMat t
229 => Double -- ^ relative tolerance in 'eps' units 225 => Double -- ^ relative tolerance in 'eps' units
230 -> Matrix t -- ^ input matrix 226 -> Matrix t -- ^ input matrix
231 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace 227 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
232nullspacePrec t m = ns where 228nullspacePrec t m = ns where
@@ -234,7 +230,6 @@ nullspacePrec t m = ns where
234 sl@(g:_) = toList s 230 sl@(g:_) = toList s
235 tol = (fromIntegral (max (rows m) (cols m)) * g * t * eps) 231 tol = (fromIntegral (max (rows m) (cols m)) * g * t * eps)
236 rank = length (filter (> g*tol) sl) 232 rank = length (filter (> g*tol) sl)
237-- ns = drop rank (toColumns v)
238 ns = drop rank $ toRows $ ctrans v 233 ns = drop rank $ toRows $ ctrans v
239 234
240-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@). 235-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@).