diff options
-rw-r--r-- | examples/tests.hs | 1 | ||||
-rw-r--r-- | lib/Numeric/GSL/Minimization.hs | 2 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 39 |
3 files changed, 29 insertions, 13 deletions
diff --git a/examples/tests.hs b/examples/tests.hs index ab1b485..974ec24 100644 --- a/examples/tests.hs +++ b/examples/tests.hs | |||
@@ -204,6 +204,7 @@ eigTest m = complex m <> v |~| v <> diag s | |||
204 | 204 | ||
205 | eigTestSH m = m <> v |~| v <> real (diag s) | 205 | eigTestSH m = m <> v |~| v <> real (diag s) |
206 | && orthonormal v | 206 | && orthonormal v |
207 | && m |~| v <> real (diag s) <> ctrans v | ||
207 | where (s, v) = eigSH m | 208 | where (s, v) = eigSH m |
208 | 209 | ||
209 | rank m | m |~| 0 = 0 | 210 | rank m | m |~| 0 = 0 |
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs index 23ca51b..e93f8cb 100644 --- a/lib/Numeric/GSL/Minimization.hs +++ b/lib/Numeric/GSL/Minimization.hs | |||
@@ -150,11 +150,9 @@ minimizeConjugateGradient istep minimpar tol maxit f df xi = unsafePerformIO $ d | |||
150 | df' = (fromList . df . toList) | 150 | df' = (fromList . df . toList) |
151 | fp <- mkVecfun (iv f') | 151 | fp <- mkVecfun (iv f') |
152 | dfp <- mkVecVecfun (aux_vTov df') | 152 | dfp <- mkVecVecfun (aux_vTov df') |
153 | print "entro" | ||
154 | rawpath <- createMIO maxit (n+2) | 153 | rawpath <- createMIO maxit (n+2) |
155 | (c_minimizeConjugateGradient fp dfp istep minimpar tol maxit // vec xiv) | 154 | (c_minimizeConjugateGradient fp dfp istep minimpar tol maxit // vec xiv) |
156 | "minimizeDerivV" [xiv] | 155 | "minimizeDerivV" [xiv] |
157 | print "salgo" | ||
158 | let it = round (rawpath @@> (maxit-1,0)) | 156 | let it = round (rawpath @@> (maxit-1,0)) |
159 | path = takeRows it rawpath | 157 | path = takeRows it rawpath |
160 | sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path | 158 | sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path |
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index 5492e07..3757f33 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -9,11 +9,11 @@ Maintainer : Alberto Ruiz (aruiz at um dot es) | |||
9 | Stability : provisional | 9 | Stability : provisional |
10 | Portability : uses ffi | 10 | Portability : uses ffi |
11 | 11 | ||
12 | A generic interface for a number of essential functions. Using it some higher level algorithms | 12 | A generic interface for some common functions. Using it we can write higher level algorithms |
13 | and testing properties can be written for both real and complex matrices. | 13 | and testing properties both for real and complex matrices. |
14 | 14 | ||
15 | In any case, the specific functions for particular base types can also be explicitly | 15 | In any case, the specific functions for particular base types can also be explicitly |
16 | imported from the LAPACK and GSL.Matrix modules. | 16 | imported from "Numeric.LinearAlgebra.LAPACK". |
17 | 17 | ||
18 | -} | 18 | -} |
19 | ----------------------------------------------------------------------------- | 19 | ----------------------------------------------------------------------------- |
@@ -29,16 +29,20 @@ module Numeric.LinearAlgebra.Algorithms ( | |||
29 | full, economy, | 29 | full, economy, |
30 | -- ** Eigensystems | 30 | -- ** Eigensystems |
31 | eig, eigSH, | 31 | eig, eigSH, |
32 | -- ** Other | 32 | -- ** QR |
33 | Numeric.LinearAlgebra.Algorithms.qr, chol, | 33 | qr, |
34 | -- ** Cholesky | ||
35 | chol, | ||
34 | -- * Nullspace | 36 | -- * Nullspace |
35 | nullspacePrec, | 37 | nullspacePrec, |
36 | nullVector, | 38 | nullVector, |
39 | -- * Norms | ||
40 | Normed(..), NormType(..), | ||
37 | -- * Misc | 41 | -- * Misc |
38 | eps, i, | ||
39 | ctrans, | 42 | ctrans, |
40 | Normed(..), NormType(..), | 43 | eps, i, |
41 | GenMat(linearSolveSVD,lu,eigSH'), unpackQR | 44 | -- * Util |
45 | GenMat(linearSolveSVD,lu,eigSH',cholSH), unpackQR, haussholder | ||
42 | ) where | 46 | ) where |
43 | 47 | ||
44 | 48 | ||
@@ -53,15 +57,22 @@ import Data.List(foldl1') | |||
53 | 57 | ||
54 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. | 58 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. |
55 | class (Linear Matrix t) => GenMat t where | 59 | class (Linear Matrix t) => GenMat t where |
60 | -- | Singular value decomposition using lapack's dgesvd or zgesvd. | ||
56 | svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) | 61 | svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) |
57 | lu :: Matrix t -> (Matrix t, Matrix t, [Int], t) | 62 | lu :: Matrix t -> (Matrix t, Matrix t, [Int], t) |
63 | -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv and zgesv. | ||
64 | -- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK". | ||
58 | linearSolve :: Matrix t -> Matrix t -> Matrix t | 65 | linearSolve :: Matrix t -> Matrix t -> Matrix t |
59 | 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@ | ||
60 | 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. | ||
61 | 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. | ||
62 | cholSH :: Matrix t -> Matrix t | 72 | cholSH :: Matrix t -> Matrix t |
73 | -- | QR factorization using lapack's dgeqr2 or zgeqr2. | ||
63 | qr :: Matrix t -> (Matrix t, Matrix t) | 74 | qr :: Matrix t -> (Matrix t, Matrix t) |
64 | -- | conjugate transpose | 75 | -- | conjugate transpose. |
65 | ctrans :: Matrix t -> Matrix t | 76 | ctrans :: Matrix t -> Matrix t |
66 | 77 | ||
67 | instance GenMat Double where | 78 | instance GenMat Double where |
@@ -86,12 +97,12 @@ instance GenMat (Complex Double) where | |||
86 | cholSH = cholH | 97 | cholSH = cholH |
87 | qr = unpackQR . qrC | 98 | qr = unpackQR . qrC |
88 | 99 | ||
89 | -- | eigensystem of complex hermitian or real symmetric matrix | 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@ |
90 | eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) | 101 | eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) |
91 | eigSH m | m `equal` ctrans m = eigSH' m | 102 | eigSH m | m `equal` ctrans m = eigSH' m |
92 | | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" | 103 | | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" |
93 | 104 | ||
94 | -- | Cholesky factorization of a positive definite hermitian or symmetric matrix | 105 | -- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf. |
95 | chol :: GenMat t => Matrix t -> Matrix t | 106 | chol :: GenMat t => Matrix t -> Matrix t |
96 | chol m | m `equal` ctrans m = cholSH m | 107 | chol m | m `equal` ctrans m = cholSH m |
97 | | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" | 108 | | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" |
@@ -103,13 +114,16 @@ det m | square m = s * (product $ toList $ takeDiag $ u) | |||
103 | | otherwise = error "det of nonsquare matrix" | 114 | | otherwise = error "det of nonsquare matrix" |
104 | where (_,u,_,s) = lu m | 115 | where (_,u,_,s) = lu m |
105 | 116 | ||
117 | -- | Inverse of a square matrix using lapacks' dgesv and zgesv. | ||
106 | inv :: GenMat t => Matrix t -> Matrix t | 118 | inv :: GenMat t => Matrix t -> Matrix t |
107 | inv m | square m = m `linearSolve` ident (rows m) | 119 | inv m | square m = m `linearSolve` ident (rows m) |
108 | | otherwise = error "inv of nonsquare matrix" | 120 | | otherwise = error "inv of nonsquare matrix" |
109 | 121 | ||
122 | -- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss. | ||
110 | pinv :: GenMat t => Matrix t -> Matrix t | 123 | pinv :: GenMat t => Matrix t -> Matrix t |
111 | pinv m = linearSolveSVD m (ident (rows m)) | 124 | pinv m = linearSolveSVD m (ident (rows m)) |
112 | 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@. | ||
113 | full :: Field t | 127 | full :: Field t |
114 | => (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) |
115 | full svd m = (u, d ,v) where | 129 | full svd m = (u, d ,v) where |
@@ -118,6 +132,7 @@ full svd m = (u, d ,v) where | |||
118 | r = rows m | 132 | r = rows m |
119 | c = cols m | 133 | c = cols m |
120 | 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@. | ||
121 | economy :: Field t | 136 | economy :: Field t |
122 | => (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) |
123 | economy svd m = (u', subVector 0 d s, v') where | 138 | economy svd m = (u', subVector 0 d s, v') where |
@@ -263,9 +278,11 @@ pinvTol t m = v' `mXm` diag s' `mXm` trans u' where | |||
263 | 278 | ||
264 | -- many thanks, quickcheck! | 279 | -- many thanks, quickcheck! |
265 | 280 | ||
281 | haussholder :: (GenMat a) => a -> Vector a -> Matrix a | ||
266 | haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) | 282 | haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) |
267 | where w = asColumn v | 283 | where w = asColumn v |
268 | 284 | ||
285 | unpackQR :: (GenMat t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) | ||
269 | unpackQR (pq, tau) = (q,r) | 286 | unpackQR (pq, tau) = (q,r) |
270 | where cs = toColumns pq | 287 | where cs = toColumns pq |
271 | m = rows pq | 288 | m = rows pq |