diff options
Diffstat (limited to 'lib/LAPACK.hs')
-rw-r--r-- | lib/LAPACK.hs | 261 |
1 files changed, 259 insertions, 2 deletions
diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs index 0019fbe..e6249a1 100644 --- a/lib/LAPACK.hs +++ b/lib/LAPACK.hs | |||
@@ -1,3 +1,4 @@ | |||
1 | {-# OPTIONS_GHC -fglasgow-exts #-} | ||
1 | ----------------------------------------------------------------------------- | 2 | ----------------------------------------------------------------------------- |
2 | -- | | 3 | -- | |
3 | -- Module : LAPACK | 4 | -- Module : LAPACK |
@@ -13,11 +14,267 @@ | |||
13 | ----------------------------------------------------------------------------- | 14 | ----------------------------------------------------------------------------- |
14 | 15 | ||
15 | module LAPACK ( | 16 | module LAPACK ( |
16 | svdR, svdR', svdC, svdC', | 17 | svdR, svdR', svdRdd, svdRdd', svdC, svdC', |
17 | eigC, eigR, eigS, eigH, | 18 | eigC, eigR, eigS, eigH, |
18 | linearSolveR, linearSolveC, | 19 | linearSolveR, linearSolveC, |
19 | linearSolveLSR, linearSolveLSC, | 20 | linearSolveLSR, linearSolveLSC, |
20 | linearSolveSVDR, linearSolveSVDC, | 21 | linearSolveSVDR, linearSolveSVDC, |
21 | ) where | 22 | ) where |
22 | 23 | ||
23 | import LAPACK.Internal | 24 | import Data.Packed.Internal.Vector |
25 | import Data.Packed.Internal.Matrix | ||
26 | import Complex | ||
27 | import Foreign | ||
28 | |||
29 | ----------------------------------------------------------------------------- | ||
30 | -- dgesvd | ||
31 | foreign import ccall "LAPACK/lapack-aux.h svd_l_R" | ||
32 | dgesvd :: Double ::> Double ::> (Double :> Double ::> IO Int) | ||
33 | |||
34 | -- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix. | ||
35 | -- | ||
36 | -- @(u,s,v)=svdR m@ so that @m=u \<\> s \<\> 'trans' v@. | ||
37 | svdR :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double) | ||
38 | svdR x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdR' x | ||
39 | |||
40 | svdR' x@M {rows = r, cols = c} = unsafePerformIO $ do | ||
41 | u <- createMatrix ColumnMajor r r | ||
42 | s <- createVector (min r c) | ||
43 | v <- createMatrix ColumnMajor c c | ||
44 | dgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdR" [fdat x] | ||
45 | return (u,s,trans v) | ||
46 | |||
47 | ----------------------------------------------------------------------------- | ||
48 | -- dgesdd | ||
49 | foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd" | ||
50 | dgesdd :: Double ::> Double ::> (Double :> Double ::> IO Int) | ||
51 | |||
52 | -- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix. | ||
53 | -- | ||
54 | -- @(u,s,v)=svdRdd m@ so that @m=u \<\> s \<\> 'trans' v@. | ||
55 | svdRdd :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double) | ||
56 | svdRdd x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdRdd' x | ||
57 | |||
58 | svdRdd' x@M {rows = r, cols = c} = unsafePerformIO $ do | ||
59 | u <- createMatrix ColumnMajor r r | ||
60 | s <- createVector (min r c) | ||
61 | v <- createMatrix ColumnMajor c c | ||
62 | dgesdd // mat fdat x // mat dat u // vec s // mat dat v // check "svdRdd" [fdat x] | ||
63 | return (u,s,trans v) | ||
64 | |||
65 | ----------------------------------------------------------------------------- | ||
66 | -- zgesvd | ||
67 | foreign import ccall "LAPACK/lapack-aux.h svd_l_C" | ||
68 | zgesvd :: (Complex Double) ::> (Complex Double) ::> (Double :> (Complex Double) ::> IO Int) | ||
69 | |||
70 | -- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix. | ||
71 | -- | ||
72 | -- @(u,s,v)=svdC m@ so that @m=u \<\> s \<\> 'trans' v@. | ||
73 | svdC :: Matrix (Complex Double) | ||
74 | -> (Matrix (Complex Double), Matrix Double, Matrix (Complex Double)) | ||
75 | svdC x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdC' x | ||
76 | |||
77 | svdC' x@M {rows = r, cols = c} = unsafePerformIO $ do | ||
78 | u <- createMatrix ColumnMajor r r | ||
79 | s <- createVector (min r c) | ||
80 | v <- createMatrix ColumnMajor c c | ||
81 | zgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdC" [fdat x] | ||
82 | return (u,s,trans v) | ||
83 | |||
84 | ----------------------------------------------------------------------------- | ||
85 | -- zgeev | ||
86 | foreign import ccall "LAPACK/lapack-aux.h eig_l_C" | ||
87 | zgeev :: (Complex Double) ::> (Complex Double) ::> ((Complex Double) :> (Complex Double) ::> IO Int) | ||
88 | |||
89 | -- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix: | ||
90 | -- | ||
91 | -- if @(l,v)=eigC m@ then @m \<\> v = v \<\> diag l@. | ||
92 | -- | ||
93 | -- The eigenvectors are the columns of v. | ||
94 | -- The eigenvalues are not sorted. | ||
95 | eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) | ||
96 | eigC (m@M {rows = r}) | ||
97 | | r == 1 = (fromList [cdat m `at` 0], singleton 1) | ||
98 | | otherwise = unsafePerformIO $ do | ||
99 | l <- createVector r | ||
100 | v <- createMatrix ColumnMajor r r | ||
101 | dummy <- createMatrix ColumnMajor 1 1 | ||
102 | zgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigC" [fdat m] | ||
103 | return (l,v) | ||
104 | |||
105 | ----------------------------------------------------------------------------- | ||
106 | -- dgeev | ||
107 | foreign import ccall "LAPACK/lapack-aux.h eig_l_R" | ||
108 | dgeev :: Double ::> Double ::> ((Complex Double) :> Double ::> IO Int) | ||
109 | |||
110 | -- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix: | ||
111 | -- | ||
112 | -- if @(l,v)=eigR m@ then @m \<\> v = v \<\> diag l@. | ||
113 | -- | ||
114 | -- The eigenvectors are the columns of v. | ||
115 | -- The eigenvalues are not sorted. | ||
116 | eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) | ||
117 | eigR (m@M {rows = r}) = (s', v'') | ||
118 | where (s,v) = eigRaux m | ||
119 | s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s)) | ||
120 | v' = toRows $ trans v | ||
121 | v'' = fromColumns $ fixeig (toList s') v' | ||
122 | |||
123 | eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double) | ||
124 | eigRaux (m@M {rows = r}) | ||
125 | | r == 1 = (fromList [(cdat m `at` 0):+0], singleton 1) | ||
126 | | otherwise = unsafePerformIO $ do | ||
127 | l <- createVector r | ||
128 | v <- createMatrix ColumnMajor r r | ||
129 | dummy <- createMatrix ColumnMajor 1 1 | ||
130 | dgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigR" [fdat m] | ||
131 | return (l,v) | ||
132 | |||
133 | fixeig [] _ = [] | ||
134 | fixeig [r] [v] = [comp v] | ||
135 | fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) | ||
136 | | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs | ||
137 | | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs) | ||
138 | |||
139 | scale r v = fromList [r] `outer` v | ||
140 | |||
141 | ----------------------------------------------------------------------------- | ||
142 | -- dsyev | ||
143 | foreign import ccall "LAPACK/lapack-aux.h eig_l_S" | ||
144 | dsyev :: Double ::> (Double :> Double ::> IO Int) | ||
145 | |||
146 | -- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix: | ||
147 | -- | ||
148 | -- if @(l,v)=eigSl m@ then @m \<\> v = v \<\> diag l@. | ||
149 | -- | ||
150 | -- The eigenvectors are the columns of v. | ||
151 | -- The eigenvalues are sorted in descending order (use eigS' for ascending order). | ||
152 | eigS :: Matrix Double -> (Vector Double, Matrix Double) | ||
153 | eigS m = (s', fliprl v) | ||
154 | where (s,v) = eigS' m | ||
155 | s' = fromList . reverse . toList $ s | ||
156 | |||
157 | eigS' (m@M {rows = r}) | ||
158 | | r == 1 = (fromList [cdat m `at` 0], singleton 1) | ||
159 | | otherwise = unsafePerformIO $ do | ||
160 | l <- createVector r | ||
161 | v <- createMatrix ColumnMajor r r | ||
162 | dsyev // mat fdat m // vec l // mat dat v // check "eigS" [fdat m] | ||
163 | return (l,v) | ||
164 | |||
165 | ----------------------------------------------------------------------------- | ||
166 | -- zheev | ||
167 | foreign import ccall "LAPACK/lapack-aux.h eig_l_H" | ||
168 | zheev :: (Complex Double) ::> (Double :> (Complex Double) ::> IO Int) | ||
169 | |||
170 | -- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix: | ||
171 | -- | ||
172 | -- if @(l,v)=eigH m@ then @m \<\> s v = v \<\> diag l@. | ||
173 | -- | ||
174 | -- The eigenvectors are the columns of v. | ||
175 | -- The eigenvalues are sorted in descending order. | ||
176 | eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | ||
177 | eigH m = (s', fliprl v) | ||
178 | where (s,v) = eigH' m | ||
179 | s' = fromList . reverse . toList $ s | ||
180 | |||
181 | eigH' (m@M {rows = r}) | ||
182 | | r == 1 = (fromList [realPart (cdat m `at` 0)], singleton 1) | ||
183 | | otherwise = unsafePerformIO $ do | ||
184 | l <- createVector r | ||
185 | v <- createMatrix ColumnMajor r r | ||
186 | zheev // mat fdat m // vec l // mat dat v // check "eigH" [fdat m] | ||
187 | return (l,v) | ||
188 | |||
189 | ----------------------------------------------------------------------------- | ||
190 | -- dgesv | ||
191 | foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" | ||
192 | dgesv :: Double ::> Double ::> Double ::> IO Int | ||
193 | |||
194 | -- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition. | ||
195 | linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double | ||
196 | linearSolveR a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c}) | ||
197 | | n1==n2 && n1==r = unsafePerformIO $ do | ||
198 | s <- createMatrix ColumnMajor r c | ||
199 | dgesv // mat fdat a // mat fdat b // mat dat s // check "linearSolveR" [fdat a, fdat b] | ||
200 | return s | ||
201 | | otherwise = error "linearSolveR of nonsquare matrix" | ||
202 | |||
203 | ----------------------------------------------------------------------------- | ||
204 | -- zgesv | ||
205 | foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" | ||
206 | zgesv :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int | ||
207 | |||
208 | -- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition. | ||
209 | linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
210 | linearSolveC a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c}) | ||
211 | | n1==n2 && n1==r = unsafePerformIO $ do | ||
212 | s <- createMatrix ColumnMajor r c | ||
213 | zgesv // mat fdat a // mat fdat b // mat dat s // check "linearSolveC" [fdat a, fdat b] | ||
214 | return s | ||
215 | | otherwise = error "linearSolveC of nonsquare matrix" | ||
216 | |||
217 | ----------------------------------------------------------------------------------- | ||
218 | -- dgels | ||
219 | foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l" | ||
220 | dgels :: Double ::> Double ::> Double ::> IO Int | ||
221 | |||
222 | -- | Wrapper for LAPACK's /dgels/, which obtains the least squared error solution of an overconstrained real linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDR'. | ||
223 | linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double | ||
224 | linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSR_l a b | ||
225 | |||
226 | linearSolveLSR_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do | ||
227 | r <- createMatrix ColumnMajor (max m n) nrhs | ||
228 | dgels // mat fdat a // mat fdat b // mat dat r // check "linearSolveLSR" [fdat a, fdat b] | ||
229 | return r | ||
230 | |||
231 | ----------------------------------------------------------------------------------- | ||
232 | -- zgels | ||
233 | foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l" | ||
234 | zgels :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int | ||
235 | |||
236 | -- | Wrapper for LAPACK's /zgels/, which obtains the least squared error solution of an overconstrained complex linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDC'. | ||
237 | linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
238 | linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSC_l a b | ||
239 | |||
240 | linearSolveLSC_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do | ||
241 | r <- createMatrix ColumnMajor (max m n) nrhs | ||
242 | zgels // mat fdat a // mat fdat b // mat dat r // check "linearSolveLSC" [fdat a, fdat b] | ||
243 | return r | ||
244 | |||
245 | ----------------------------------------------------------------------------------- | ||
246 | -- dgelss | ||
247 | foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l" | ||
248 | dgelss :: Double -> Double ::> Double ::> Double ::> IO Int | ||
249 | |||
250 | -- | Wrapper for LAPACK's /dgelss/, which obtains the minimum norm solution to a real linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. | ||
251 | linearSolveSVDR :: Maybe Double -- ^ rcond | ||
252 | -> Matrix Double -- ^ coefficient matrix | ||
253 | -> Matrix Double -- ^ right hand sides (as columns) | ||
254 | -> Matrix Double -- ^ solution vectors (as columns) | ||
255 | linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDR_l rcond a b | ||
256 | linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) a b | ||
257 | |||
258 | linearSolveSVDR_l rcond a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do | ||
259 | r <- createMatrix ColumnMajor (max m n) nrhs | ||
260 | dgelss rcond // mat fdat a // mat fdat b // mat dat r // check "linearSolveSVDR" [fdat a, fdat b] | ||
261 | return r | ||
262 | |||
263 | ----------------------------------------------------------------------------------- | ||
264 | -- zgelss | ||
265 | foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l" | ||
266 | zgelss :: Double -> (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int | ||
267 | |||
268 | -- | Wrapper for LAPACK's /zgelss/, which obtains the minimum norm solution to a complex linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. | ||
269 | linearSolveSVDC :: Maybe Double -- ^ rcond | ||
270 | -> Matrix (Complex Double) -- ^ coefficient matrix | ||
271 | -> Matrix (Complex Double) -- ^ right hand sides (as columns) | ||
272 | -> Matrix (Complex Double) -- ^ solution vectors (as columns) | ||
273 | linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDC_l rcond a b | ||
274 | linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) a b | ||
275 | |||
276 | linearSolveSVDC_l rcond a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do | ||
277 | r <- createMatrix ColumnMajor (max m n) nrhs | ||
278 | zgelss rcond // mat fdat a // mat fdat b // mat dat r // check "linearSolveSVDC" [fdat a, fdat b] | ||
279 | return r | ||
280 | |||