diff options
-rw-r--r-- | lib/Numeric/LinearAlgebra.hs | 4 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 205 |
2 files changed, 127 insertions, 82 deletions
diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs index cc2c1d3..b8fd01e 100644 --- a/lib/Numeric/LinearAlgebra.hs +++ b/lib/Numeric/LinearAlgebra.hs | |||
@@ -1,7 +1,7 @@ | |||
1 | ----------------------------------------------------------------------------- | 1 | ----------------------------------------------------------------------------- |
2 | {- | | 2 | {- | |
3 | Module : Numeric.LinearAlgebra | 3 | Module : Numeric.LinearAlgebra |
4 | Copyright : (c) Alberto Ruiz 2006-7 | 4 | Copyright : (c) Alberto Ruiz 2006-9 |
5 | License : GPL-style | 5 | License : GPL-style |
6 | 6 | ||
7 | Maintainer : Alberto Ruiz (aruiz at um dot es) | 7 | Maintainer : Alberto Ruiz (aruiz at um dot es) |
@@ -10,6 +10,8 @@ Portability : uses ffi | |||
10 | 10 | ||
11 | Basic matrix computations implemented by BLAS, LAPACK and GSL. | 11 | Basic matrix computations implemented by BLAS, LAPACK and GSL. |
12 | 12 | ||
13 | This is module reexports the most comon functions (including "Numeric.LinearAlgebra.Instances"). | ||
14 | |||
13 | -} | 15 | -} |
14 | ----------------------------------------------------------------------------- | 16 | ----------------------------------------------------------------------------- |
15 | module Numeric.LinearAlgebra ( | 17 | module Numeric.LinearAlgebra ( |
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index d978fa4..e22246f 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -3,7 +3,7 @@ | |||
3 | ----------------------------------------------------------------------------- | 3 | ----------------------------------------------------------------------------- |
4 | {- | | 4 | {- | |
5 | Module : Numeric.LinearAlgebra.Algorithms | 5 | Module : Numeric.LinearAlgebra.Algorithms |
6 | Copyright : (c) Alberto Ruiz 2006-7 | 6 | Copyright : (c) Alberto Ruiz 2006-9 |
7 | License : GPL-style | 7 | License : GPL-style |
8 | 8 | ||
9 | Maintainer : Alberto Ruiz (aruiz at um dot es) | 9 | Maintainer : Alberto Ruiz (aruiz at um dot es) |
@@ -20,27 +20,33 @@ imported from "Numeric.LinearAlgebra.LAPACK". | |||
20 | ----------------------------------------------------------------------------- | 20 | ----------------------------------------------------------------------------- |
21 | 21 | ||
22 | module Numeric.LinearAlgebra.Algorithms ( | 22 | module Numeric.LinearAlgebra.Algorithms ( |
23 | -- * Linear Systems | 23 | -- * Supported types |
24 | Field(), | ||
25 | -- * Products | ||
24 | multiply, dot, | 26 | multiply, dot, |
27 | outer, kronecker, | ||
28 | -- * Linear Systems | ||
25 | linearSolve, | 29 | linearSolve, |
30 | luSolve, | ||
31 | linearSolveSVD, | ||
26 | inv, pinv, | 32 | inv, pinv, |
27 | pinvTol, det, rank, rcond, | 33 | det, rank, rcond, |
28 | -- * Matrix factorizations | 34 | -- * Matrix factorizations |
29 | -- ** Singular value decomposition | 35 | -- ** Singular value decomposition |
30 | svd, | 36 | svd, |
31 | full, economy, --thin, | 37 | full, economy, --thin, |
32 | -- ** Eigensystems | 38 | -- ** Eigensystems |
33 | eig, eigSH, | 39 | eig, eigSH, eigSH', |
34 | -- ** QR | 40 | -- ** QR |
35 | qr, | 41 | qr, |
36 | -- ** Cholesky | 42 | -- ** Cholesky |
37 | chol, | 43 | chol, cholSH, |
38 | -- ** Hessenberg | 44 | -- ** Hessenberg |
39 | hess, | 45 | hess, |
40 | -- ** Schur | 46 | -- ** Schur |
41 | schur, | 47 | schur, |
42 | -- ** LU | 48 | -- ** LU |
43 | lu, luPacked, luSolve, | 49 | lu, luPacked, |
44 | -- * Matrix functions | 50 | -- * Matrix functions |
45 | expm, | 51 | expm, |
46 | sqrtm, | 52 | sqrtm, |
@@ -53,11 +59,10 @@ module Numeric.LinearAlgebra.Algorithms ( | |||
53 | -- * Misc | 59 | -- * Misc |
54 | ctrans, | 60 | ctrans, |
55 | eps, i, | 61 | eps, i, |
56 | outer, kronecker, | ||
57 | -- * Util | 62 | -- * Util |
58 | haussholder, | 63 | haussholder, |
59 | unpackQR, unpackHess, | 64 | unpackQR, unpackHess, |
60 | Field(linearSolveSVD,eigSH',cholSH) | 65 | pinvTol |
61 | ) where | 66 | ) where |
62 | 67 | ||
63 | 68 | ||
@@ -71,82 +76,120 @@ import Data.List(foldl1') | |||
71 | import Data.Array | 76 | import Data.Array |
72 | 77 | ||
73 | 78 | ||
74 | |||
75 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. | 79 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. |
76 | class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where | 80 | class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where |
77 | -- | Singular value decomposition using lapack's dgesvd or zgesvd. | 81 | svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) |
78 | svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) | 82 | luPacked' :: Matrix t -> (Matrix t, [Int]) |
79 | -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. | 83 | luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t |
80 | luPacked :: Matrix t -> (Matrix t, [Int]) | 84 | linearSolve' :: Matrix t -> Matrix t -> Matrix t |
81 | -- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization | 85 | linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t |
82 | -- obtained by 'luPacked'. | 86 | eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) |
83 | luSolve :: (Matrix t, [Int]) -> Matrix t -> Matrix t | 87 | eigSH'' :: Matrix t -> (Vector Double, Matrix t) |
84 | -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv. | 88 | cholSH' :: Matrix t -> Matrix t |
85 | -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system. | 89 | qr' :: Matrix t -> (Matrix t, Matrix t) |
86 | -- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK". | 90 | hess' :: Matrix t -> (Matrix t, Matrix t) |
87 | linearSolve :: Matrix t -> Matrix t -> Matrix t | 91 | schur' :: Matrix t -> (Matrix t, Matrix t) |
88 | linearSolveSVD :: Matrix t -> Matrix t -> Matrix t | 92 | ctrans' :: Matrix t -> Matrix t |
89 | -- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev. | 93 | multiply' :: Matrix t -> Matrix t -> Matrix t |
90 | -- | 94 | |
91 | -- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ | 95 | |
92 | eig :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) | 96 | -- | Singular value decomposition using lapack's dgesvd or zgesvd. |
93 | -- | Similar to eigSH without checking that the input matrix is hermitian or symmetric. | 97 | svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) |
94 | eigSH' :: Matrix t -> (Vector Double, Matrix t) | 98 | svd = svd' |
95 | -- | Similar to chol without checking that the input matrix is hermitian or symmetric. | 99 | |
96 | cholSH :: Matrix t -> Matrix t | 100 | -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. |
97 | -- | QR factorization using lapack's dgeqr2 or zgeqr2. | 101 | luPacked :: Field t => Matrix t -> (Matrix t, [Int]) |
98 | -- | 102 | luPacked = luPacked' |
99 | -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. | 103 | |
100 | qr :: Matrix t -> (Matrix t, Matrix t) | 104 | -- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization |
101 | -- | Hessenberg factorization using lapack's dgehrd or zgehrd. | 105 | -- obtained by 'luPacked'. |
102 | -- | 106 | luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t |
103 | -- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary | 107 | luSolve = luSolve' |
104 | -- and h is in upper Hessenberg form. | 108 | |
105 | hess :: Matrix t -> (Matrix t, Matrix t) | 109 | -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv. |
106 | -- | Schur factorization using lapack's dgees or zgees. | 110 | -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system. |
107 | -- | 111 | -- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK". |
108 | -- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary | 112 | linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t |
109 | -- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is | 113 | linearSolve = linearSolve' |
110 | -- upper triangular in 2x2 blocks. | 114 | linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t |
111 | -- | 115 | linearSolveSVD = linearSolveSVD' |
112 | -- \"Anything that the Jordan decomposition can do, the Schur decomposition | 116 | |
113 | -- can do better!\" (Van Loan) | 117 | -- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev. |
114 | schur :: Matrix t -> (Matrix t, Matrix t) | 118 | -- |
115 | -- | Conjugate transpose. | 119 | -- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ |
116 | ctrans :: Matrix t -> Matrix t | 120 | eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) |
117 | -- | Matrix product. | 121 | eig = eig' |
118 | multiply :: Matrix t -> Matrix t -> Matrix t | 122 | |
123 | -- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. | ||
124 | eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t) | ||
125 | eigSH' = eigSH'' | ||
126 | |||
127 | -- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric. | ||
128 | cholSH :: Field t => Matrix t -> Matrix t | ||
129 | cholSH = cholSH' | ||
130 | |||
131 | -- | QR factorization using lapack's dgeqr2 or zgeqr2. | ||
132 | -- | ||
133 | -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. | ||
134 | qr :: Field t => Matrix t -> (Matrix t, Matrix t) | ||
135 | qr = qr' | ||
136 | |||
137 | -- | Hessenberg factorization using lapack's dgehrd or zgehrd. | ||
138 | -- | ||
139 | -- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary | ||
140 | -- and h is in upper Hessenberg form. | ||
141 | hess :: Field t => Matrix t -> (Matrix t, Matrix t) | ||
142 | hess = hess' | ||
143 | |||
144 | -- | Schur factorization using lapack's dgees or zgees. | ||
145 | -- | ||
146 | -- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary | ||
147 | -- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is | ||
148 | -- upper triangular in 2x2 blocks. | ||
149 | -- | ||
150 | -- \"Anything that the Jordan decomposition can do, the Schur decomposition | ||
151 | -- can do better!\" (Van Loan) | ||
152 | schur :: Field t => Matrix t -> (Matrix t, Matrix t) | ||
153 | schur = schur' | ||
154 | |||
155 | -- | Generic conjugate transpose. | ||
156 | ctrans :: Field t => Matrix t -> Matrix t | ||
157 | ctrans = ctrans' | ||
158 | |||
159 | -- | Matrix product. | ||
160 | multiply :: Field t => Matrix t -> Matrix t -> Matrix t | ||
161 | multiply = multiply' | ||
119 | 162 | ||
120 | 163 | ||
121 | instance Field Double where | 164 | instance Field Double where |
122 | svd = svdR | 165 | svd' = svdR |
123 | luPacked = luR | 166 | luPacked' = luR |
124 | luSolve (l_u,perm) = lusR l_u perm | 167 | luSolve' (l_u,perm) = lusR l_u perm |
125 | linearSolve = linearSolveR -- (luSolve . luPacked) ?? | 168 | linearSolve' = linearSolveR -- (luSolve . luPacked) ?? |
126 | linearSolveSVD = linearSolveSVDR Nothing | 169 | linearSolveSVD' = linearSolveSVDR Nothing |
127 | ctrans = trans | 170 | ctrans' = trans |
128 | eig = eigR | 171 | eig' = eigR |
129 | eigSH' = eigS | 172 | eigSH'' = eigS |
130 | cholSH = cholS | 173 | cholSH' = cholS |
131 | qr = unpackQR . qrR | 174 | qr' = unpackQR . qrR |
132 | hess = unpackHess hessR | 175 | hess' = unpackHess hessR |
133 | schur = schurR | 176 | schur' = schurR |
134 | multiply = multiplyR | 177 | multiply' = multiplyR |
135 | 178 | ||
136 | instance Field (Complex Double) where | 179 | instance Field (Complex Double) where |
137 | svd = svdC | 180 | svd' = svdC |
138 | luPacked = luC | 181 | luPacked' = luC |
139 | luSolve (l_u,perm) = lusC l_u perm | 182 | luSolve' (l_u,perm) = lusC l_u perm |
140 | linearSolve = linearSolveC | 183 | linearSolve' = linearSolveC |
141 | linearSolveSVD = linearSolveSVDC Nothing | 184 | linearSolveSVD' = linearSolveSVDC Nothing |
142 | ctrans = conj . trans | 185 | ctrans' = conj . trans |
143 | eig = eigC | 186 | eig' = eigC |
144 | eigSH' = eigH | 187 | eigSH'' = eigH |
145 | cholSH = cholH | 188 | cholSH' = cholH |
146 | qr = unpackQR . qrC | 189 | qr' = unpackQR . qrC |
147 | hess = unpackHess hessC | 190 | hess' = unpackHess hessC |
148 | schur = schurC | 191 | schur' = schurC |
149 | multiply = multiplyC | 192 | multiply' = multiplyC |
150 | 193 | ||
151 | 194 | ||
152 | -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev. | 195 | -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev. |
@@ -193,8 +236,8 @@ pinv m = linearSolveSVD m (ident (rows m)) | |||
193 | -- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@. | 236 | -- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@. |
194 | full :: Element t | 237 | full :: Element t |
195 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) | 238 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) |
196 | full svd' m = (u, d ,v) where | 239 | full svdFun m = (u, d ,v) where |
197 | (u,s,v) = svd' m | 240 | (u,s,v) = svdFun m |
198 | d = diagRect s r c | 241 | d = diagRect s r c |
199 | r = rows m | 242 | r = rows m |
200 | c = cols m | 243 | c = cols m |
@@ -204,8 +247,8 @@ full svd' m = (u, d ,v) where | |||
204 | -- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. | 247 | -- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. |
205 | economy :: Element t | 248 | economy :: Element t |
206 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) | 249 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) |
207 | economy svd' m = (u', subVector 0 d s, v') where | 250 | economy svdFun m = (u', subVector 0 d s, v') where |
208 | (u,s,v) = svd' m | 251 | (u,s,v) = svdFun m |
209 | sl@(g:_) = toList s | 252 | sl@(g:_) = toList s |
210 | s' = fromList . filter (>tol) $ sl | 253 | s' = fromList . filter (>tol) $ sl |
211 | t = 1 | 254 | t = 1 |