diff options
Diffstat (limited to 'lib/Numeric')
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 31 |
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 | ||
78 | instance GenMat Double where | 78 | instance 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@ |
101 | eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) | 101 | eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) |
102 | eigSH m | m `equal` ctrans m = eigSH' m | 102 | eigSH 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) | |||
123 | pinv :: GenMat t => Matrix t -> Matrix t | 123 | pinv :: GenMat t => Matrix t -> Matrix t |
124 | pinv m = linearSolveSVD m (ident (rows m)) | 124 | pinv 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@. |
127 | full :: Field t | 127 | full :: 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) |
129 | full svd m = (u, d ,v) where | 129 | full 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@. |
136 | economy :: Field t | 136 | economy :: 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) |
138 | economy svd m = (u', subVector 0 d s, v') where | 138 | economy 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). |
153 | eps :: Double | 153 | eps :: Double |
154 | eps = 2.22044604925031e-16 | 154 | eps = 2.22044604925031e-16 |
155 | 155 | ||
156 | -- | The imaginary unit: @i == 0.0 :+ 1.0@ | 156 | -- | The imaginary unit: @i = 0.0 :+ 1.0@ |
157 | i :: Complex Double | 157 | i :: Complex Double |
158 | i = 0:+1 | 158 | i = 0:+1 |
159 | 159 | ||
160 | 160 | ||
161 | -- | matrix product | 161 | -- matrix product |
162 | mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t | 162 | mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t |
163 | mXm = multiply | 163 | mXm = multiply |
164 | 164 | ||
165 | -- | matrix - vector product | 165 | -- matrix - vector product |
166 | mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t | 166 | mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t |
167 | mXv m v = flatten $ m `mXm` (asColumn v) | 167 | mXv m v = flatten $ m `mXm` (asColumn v) |
168 | 168 | ||
169 | -- | vector - matrix product | 169 | -- vector - matrix product |
170 | vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t | 170 | vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t |
171 | vXm v m = flatten $ (asRow v) `mXm` m | 171 | vXm v m = flatten $ (asRow v) `mXm` m |
172 | 172 | ||
@@ -201,10 +201,6 @@ pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (liftVector | |||
201 | pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` constant 1 (cols m) | 201 | pnormCM 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 | |||
208 | class Normed t where | 204 | class 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. |
228 | nullspacePrec :: GenMat t | 224 | nullspacePrec :: 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 |
232 | nullspacePrec t m = ns where | 228 | nullspacePrec 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@). |