diff options
Diffstat (limited to 'lib/Numeric/GSL/Matrix.hs')
-rw-r--r-- | lib/Numeric/GSL/Matrix.hs | 311 |
1 files changed, 0 insertions, 311 deletions
diff --git a/lib/Numeric/GSL/Matrix.hs b/lib/Numeric/GSL/Matrix.hs deleted file mode 100644 index f62bb82..0000000 --- a/lib/Numeric/GSL/Matrix.hs +++ /dev/null | |||
@@ -1,311 +0,0 @@ | |||
1 | ----------------------------------------------------------------------------- | ||
2 | -- | | ||
3 | -- Module : Numeric.GSL.Matrix | ||
4 | -- Copyright : (c) Alberto Ruiz 2007 | ||
5 | -- License : GPL-style | ||
6 | -- | ||
7 | -- Maintainer : Alberto Ruiz <aruiz@um.es> | ||
8 | -- Stability : provisional | ||
9 | -- Portability : portable (uses FFI) | ||
10 | -- | ||
11 | -- A few linear algebra computations based on GSL. | ||
12 | -- | ||
13 | ----------------------------------------------------------------------------- | ||
14 | -- #hide | ||
15 | |||
16 | module Numeric.GSL.Matrix( | ||
17 | eigSg, eigHg, | ||
18 | svdg, | ||
19 | qr, qrPacked, unpackQR, | ||
20 | cholR, cholC, | ||
21 | luSolveR, luSolveC, | ||
22 | luR, luC | ||
23 | ) where | ||
24 | |||
25 | import Data.Packed.Internal | ||
26 | import Data.Packed.Matrix(ident) | ||
27 | import Numeric.GSL.Vector | ||
28 | import Foreign | ||
29 | import Complex | ||
30 | |||
31 | {- | eigendecomposition of a real symmetric matrix using /gsl_eigen_symmv/. | ||
32 | |||
33 | > > let (l,v) = eigS $ 'fromLists' [[1,2],[2,1]] | ||
34 | > > l | ||
35 | > 3.000 -1.000 | ||
36 | > | ||
37 | > > v | ||
38 | > 0.707 -0.707 | ||
39 | > 0.707 0.707 | ||
40 | > | ||
41 | > > v <> diag l <> trans v | ||
42 | > 1.000 2.000 | ||
43 | > 2.000 1.000 | ||
44 | |||
45 | -} | ||
46 | eigSg :: Matrix Double -> (Vector Double, Matrix Double) | ||
47 | eigSg = eigSg' . cmat | ||
48 | |||
49 | eigSg' m | ||
50 | | r == 1 = (fromList [cdat m `at` 0], singleton 1) | ||
51 | | otherwise = unsafePerformIO $ do | ||
52 | l <- createVector r | ||
53 | v <- createMatrix RowMajor r r | ||
54 | app3 c_eigS mat m vec l mat v "eigSg" | ||
55 | return (l,v) | ||
56 | where r = rows m | ||
57 | foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM | ||
58 | |||
59 | ------------------------------------------------------------------ | ||
60 | |||
61 | |||
62 | |||
63 | {- | eigendecomposition of a complex hermitian matrix using /gsl_eigen_hermv/ | ||
64 | |||
65 | > > let (l,v) = eigH $ 'fromLists' [[1,2+i],[2-i,3]] | ||
66 | > | ||
67 | > > l | ||
68 | > 4.449 -0.449 | ||
69 | > | ||
70 | > > v | ||
71 | > -0.544 0.839 | ||
72 | > (-0.751,0.375) (-0.487,0.243) | ||
73 | > | ||
74 | > > v <> diag l <> (conjTrans) v | ||
75 | > 1.000 (2.000,1.000) | ||
76 | > (2.000,-1.000) 3.000 | ||
77 | |||
78 | -} | ||
79 | eigHg :: Matrix (Complex Double)-> (Vector Double, Matrix (Complex Double)) | ||
80 | eigHg = eigHg' . cmat | ||
81 | |||
82 | eigHg' m | ||
83 | | r == 1 = (fromList [realPart $ cdat m `at` 0], singleton 1) | ||
84 | | otherwise = unsafePerformIO $ do | ||
85 | l <- createVector r | ||
86 | v <- createMatrix RowMajor r r | ||
87 | app3 c_eigH mat m vec l mat v "eigHg" | ||
88 | return (l,v) | ||
89 | where r = rows m | ||
90 | foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM | ||
91 | |||
92 | |||
93 | {- | Singular value decomposition of a real matrix, using /gsl_linalg_SV_decomp_mod/: | ||
94 | |||
95 | |||
96 | -} | ||
97 | svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) | ||
98 | svdg x = if rows x >= cols x | ||
99 | then svd' (cmat x) | ||
100 | else (v, s, u) where (u,s,v) = svd' (cmat (trans x)) | ||
101 | |||
102 | svd' x = unsafePerformIO $ do | ||
103 | u <- createMatrix RowMajor r c | ||
104 | s <- createVector c | ||
105 | v <- createMatrix RowMajor c c | ||
106 | app4 c_svd mat x mat u vec s mat v "svdg" | ||
107 | return (u,s,v) | ||
108 | where r = rows x | ||
109 | c = cols x | ||
110 | foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM | ||
111 | |||
112 | {- | QR decomposition of a real matrix using /gsl_linalg_QR_decomp/ and /gsl_linalg_QR_unpack/. | ||
113 | |||
114 | -} | ||
115 | qr :: Matrix Double -> (Matrix Double, Matrix Double) | ||
116 | qr = qr' . cmat | ||
117 | |||
118 | qr' x = unsafePerformIO $ do | ||
119 | q <- createMatrix RowMajor r r | ||
120 | rot <- createMatrix RowMajor r c | ||
121 | app3 c_qr mat x mat q mat rot "qr" | ||
122 | return (q,rot) | ||
123 | where r = rows x | ||
124 | c = cols x | ||
125 | foreign import ccall "gsl-aux.h QR" c_qr :: TMMM | ||
126 | |||
127 | qrPacked :: Matrix Double -> (Matrix Double, Vector Double) | ||
128 | qrPacked = qrPacked' . cmat | ||
129 | |||
130 | qrPacked' x = unsafePerformIO $ do | ||
131 | qrp <- createMatrix RowMajor r c | ||
132 | tau <- createVector (min r c) | ||
133 | app3 c_qrPacked mat x mat qrp vec tau "qrUnpacked" | ||
134 | return (qrp,tau) | ||
135 | where r = rows x | ||
136 | c = cols x | ||
137 | foreign import ccall "gsl-aux.h QRpacked" c_qrPacked :: TMMV | ||
138 | |||
139 | unpackQR :: (Matrix Double, Vector Double) -> (Matrix Double, Matrix Double) | ||
140 | unpackQR (qrp,tau) = unpackQR' (cmat qrp, tau) | ||
141 | |||
142 | unpackQR' (qrp,tau) = unsafePerformIO $ do | ||
143 | q <- createMatrix RowMajor r r | ||
144 | res <- createMatrix RowMajor r c | ||
145 | app4 c_qrUnpack mat qrp vec tau mat q mat res "qrUnpack" | ||
146 | return (q,res) | ||
147 | where r = rows qrp | ||
148 | c = cols qrp | ||
149 | foreign import ccall "gsl-aux.h QRunpack" c_qrUnpack :: TMVMM | ||
150 | |||
151 | {- | Cholesky decomposition of a symmetric positive definite real matrix using /gsl_linalg_cholesky_decomp/. | ||
152 | |||
153 | @\> chol $ (2><2) [1,2, | ||
154 | 2,9::Double] | ||
155 | (2><2) | ||
156 | [ 1.0, 0.0 | ||
157 | , 2.0, 2.23606797749979 ]@ | ||
158 | |||
159 | -} | ||
160 | cholR :: Matrix Double -> Matrix Double | ||
161 | cholR = cholR' . cmat | ||
162 | |||
163 | cholR' x = unsafePerformIO $ do | ||
164 | r <- createMatrix RowMajor n n | ||
165 | app2 c_cholR mat x mat r "cholR" | ||
166 | return r | ||
167 | where n = rows x | ||
168 | foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM | ||
169 | |||
170 | cholC :: Matrix (Complex Double) -> Matrix (Complex Double) | ||
171 | cholC = cholC' . cmat | ||
172 | |||
173 | cholC' x = unsafePerformIO $ do | ||
174 | r <- createMatrix RowMajor n n | ||
175 | app2 c_cholC mat x mat r "cholC" | ||
176 | return r | ||
177 | where n = rows x | ||
178 | foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM | ||
179 | |||
180 | |||
181 | -------------------------------------------------------- | ||
182 | |||
183 | {- -| efficient multiplication by the inverse of a matrix (for real matrices) | ||
184 | -} | ||
185 | luSolveR :: Matrix Double -> Matrix Double -> Matrix Double | ||
186 | luSolveR a b = luSolveR' (cmat a) (cmat b) | ||
187 | |||
188 | luSolveR' a b | ||
189 | | n1==n2 && n1==r = unsafePerformIO $ do | ||
190 | s <- createMatrix RowMajor r c | ||
191 | app3 c_luSolveR mat a mat b mat s "luSolveR" | ||
192 | return s | ||
193 | | otherwise = error "luSolveR of nonsquare matrix" | ||
194 | where n1 = rows a | ||
195 | n2 = cols a | ||
196 | r = rows b | ||
197 | c = cols b | ||
198 | foreign import ccall "gsl-aux.h luSolveR" c_luSolveR :: TMMM | ||
199 | |||
200 | {- -| efficient multiplication by the inverse of a matrix (for complex matrices). | ||
201 | -} | ||
202 | luSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
203 | luSolveC a b = luSolveC' (cmat a) (cmat b) | ||
204 | |||
205 | luSolveC' a b | ||
206 | | n1==n2 && n1==r = unsafePerformIO $ do | ||
207 | s <- createMatrix RowMajor r c | ||
208 | app3 c_luSolveC mat a mat b mat s "luSolveC" | ||
209 | return s | ||
210 | | otherwise = error "luSolveC of nonsquare matrix" | ||
211 | where n1 = rows a | ||
212 | n2 = cols a | ||
213 | r = rows b | ||
214 | c = cols b | ||
215 | foreign import ccall "gsl-aux.h luSolveC" c_luSolveC :: TCMCMCM | ||
216 | |||
217 | {- | lu decomposition of real matrix (packed as a vector including l, u, the permutation and sign) | ||
218 | -} | ||
219 | luRaux :: Matrix Double -> Vector Double | ||
220 | luRaux = luRaux' . cmat | ||
221 | |||
222 | luRaux' x = unsafePerformIO $ do | ||
223 | res <- createVector (r*r+r+1) | ||
224 | app2 c_luRaux mat x vec res "luRaux" | ||
225 | return res | ||
226 | where r = rows x | ||
227 | foreign import ccall "gsl-aux.h luRaux" c_luRaux :: TMV | ||
228 | |||
229 | {- | lu decomposition of complex matrix (packed as a vector including l, u, the permutation and sign) | ||
230 | -} | ||
231 | luCaux :: Matrix (Complex Double) -> Vector (Complex Double) | ||
232 | luCaux = luCaux' . cmat | ||
233 | |||
234 | luCaux' x = unsafePerformIO $ do | ||
235 | res <- createVector (r*r+r+1) | ||
236 | app2 c_luCaux mat x vec res "luCaux" | ||
237 | return res | ||
238 | where r = rows x | ||
239 | foreign import ccall "gsl-aux.h luCaux" c_luCaux :: TCMCV | ||
240 | |||
241 | {- | The LU decomposition of a square matrix. Is based on /gsl_linalg_LU_decomp/ and /gsl_linalg_complex_LU_decomp/ as described in <http://www.gnu.org/software/Numeric.GSL/manual/Numeric.GSL-ref_13.html#SEC223>. | ||
242 | |||
243 | @\> let m = 'fromLists' [[1,2,-3],[2+3*i,-7,0],[1,-i,2*i]] | ||
244 | \> let (l,u,p,s) = luR m@ | ||
245 | |||
246 | L is the lower triangular: | ||
247 | |||
248 | @\> l | ||
249 | 1. 0. 0. | ||
250 | 0.154-0.231i 1. 0. | ||
251 | 0.154-0.231i 0.624-0.522i 1.@ | ||
252 | |||
253 | U is the upper triangular: | ||
254 | |||
255 | @\> u | ||
256 | 2.+3.i -7. 0. | ||
257 | 0. 3.077-1.615i -3. | ||
258 | 0. 0. 1.873+0.433i@ | ||
259 | |||
260 | p is a permutation: | ||
261 | |||
262 | @\> p | ||
263 | [1,0,2]@ | ||
264 | |||
265 | L \* U obtains a permuted version of the original matrix: | ||
266 | |||
267 | @\> extractRows p m | ||
268 | 2.+3.i -7. 0. | ||
269 | 1. 2. -3. | ||
270 | 1. -1.i 2.i | ||
271 | \ -- CPP | ||
272 | \> l \<\> u | ||
273 | 2.+3.i -7. 0. | ||
274 | 1. 2. -3. | ||
275 | 1. -1.i 2.i@ | ||
276 | |||
277 | s is the sign of the permutation, required to obtain sign of the determinant: | ||
278 | |||
279 | @\> s * product ('toList' $ 'takeDiag' u) | ||
280 | (-18.0) :+ (-16.000000000000004) | ||
281 | \> 'LinearAlgebra.Algorithms.det' m | ||
282 | (-18.0) :+ (-16.000000000000004)@ | ||
283 | |||
284 | -} | ||
285 | luR :: Matrix Double -> (Matrix Double, Matrix Double, [Int], Double) | ||
286 | luR m = (l,u,p, fromIntegral s') where | ||
287 | r = rows m | ||
288 | v = luRaux m | ||
289 | lu = reshape r $ subVector 0 (r*r) v | ||
290 | s':p = map round . toList . subVector (r*r) (r+1) $ v | ||
291 | u = triang r r 0 1`mul` lu | ||
292 | l = (triang r r 0 0 `mul` lu) `add` ident r | ||
293 | add = liftMatrix2 $ vectorZipR Add | ||
294 | mul = liftMatrix2 $ vectorZipR Mul | ||
295 | |||
296 | -- | Complex version of 'luR'. | ||
297 | luC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double), [Int], Complex Double) | ||
298 | luC m = (l,u,p, fromIntegral s') where | ||
299 | r = rows m | ||
300 | v = luCaux m | ||
301 | lu = reshape r $ subVector 0 (r*r) v | ||
302 | s':p = map (round.realPart) . toList . subVector (r*r) (r+1) $ v | ||
303 | u = triang r r 0 1 `mul` lu | ||
304 | l = (triang r r 0 0 `mul` lu) `add` liftMatrix comp (ident r) | ||
305 | add = liftMatrix2 $ vectorZipC Add | ||
306 | mul = liftMatrix2 $ vectorZipC Mul | ||
307 | |||
308 | {- auxiliary function to get triangular matrices | ||
309 | -} | ||
310 | triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]] | ||
311 | where el i j = if j-i>=h then v else 1 - v | ||