diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-08 12:18:56 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-08 12:18:56 +0200 |
commit | 561a6c0e21bb77c21114ccbbd86d3af5ddb5a3f1 (patch) | |
tree | d49b67d75d63938229f2d5cbed5c49d06dc02bcf /packages/base/src | |
parent | 5992d92357cfd911c8f2e9f5faaa4fd8a323fd9a (diff) |
Conversion, LAPACK -> base
Diffstat (limited to 'packages/base/src')
-rw-r--r-- | packages/base/src/Data/Packed/Development.hs | 2 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/Vector.hs | 13 | ||||
-rw-r--r-- | packages/base/src/Numeric/Conversion.hs | 91 | ||||
-rw-r--r-- | packages/base/src/Numeric/LAPACK.hs | 552 |
4 files changed, 656 insertions, 2 deletions
diff --git a/packages/base/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs index 777b6c5..9350acb 100644 --- a/packages/base/src/Data/Packed/Development.hs +++ b/packages/base/src/Data/Packed/Development.hs | |||
@@ -24,7 +24,7 @@ module Data.Packed.Development ( | |||
24 | unsafeFromForeignPtr, | 24 | unsafeFromForeignPtr, |
25 | unsafeToForeignPtr, | 25 | unsafeToForeignPtr, |
26 | check, (//), | 26 | check, (//), |
27 | at', atM' | 27 | at', atM', fi, table |
28 | ) where | 28 | ) where |
29 | 29 | ||
30 | import Data.Packed.Internal | 30 | import Data.Packed.Internal |
diff --git a/packages/base/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs index b5a4318..653a257 100644 --- a/packages/base/src/Data/Packed/Vector.hs +++ b/packages/base/src/Data/Packed/Vector.hs | |||
@@ -18,7 +18,7 @@ | |||
18 | 18 | ||
19 | module Data.Packed.Vector ( | 19 | module Data.Packed.Vector ( |
20 | Vector, | 20 | Vector, |
21 | fromList, (|>), toList, buildVector, | 21 | fromList, (|>), toList, buildVector, constant, |
22 | dim, (@>), | 22 | dim, (@>), |
23 | subVector, takesV, vjoin, join, | 23 | subVector, takesV, vjoin, join, |
24 | mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith, | 24 | mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith, |
@@ -27,6 +27,7 @@ module Data.Packed.Vector ( | |||
27 | ) where | 27 | ) where |
28 | 28 | ||
29 | import Data.Packed.Internal.Vector | 29 | import Data.Packed.Internal.Vector |
30 | import Data.Packed.Internal.Matrix | ||
30 | import Foreign.Storable | 31 | import Foreign.Storable |
31 | 32 | ||
32 | ------------------------------------------------------------------- | 33 | ------------------------------------------------------------------- |
@@ -94,3 +95,13 @@ unzipVector = unzipVectorWith id | |||
94 | join :: Storable t => [Vector t] -> Vector t | 95 | join :: Storable t => [Vector t] -> Vector t |
95 | join = vjoin | 96 | join = vjoin |
96 | 97 | ||
98 | {- | creates a vector with a given number of equal components: | ||
99 | |||
100 | @> constant 2 7 | ||
101 | 7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ | ||
102 | -} | ||
103 | constant :: Element a => a -> Int -> Vector a | ||
104 | -- constant x n = runSTVector (newVector x n) | ||
105 | constant = constantD-- about 2x faster | ||
106 | |||
107 | |||
diff --git a/packages/base/src/Numeric/Conversion.hs b/packages/base/src/Numeric/Conversion.hs new file mode 100644 index 0000000..a1f9385 --- /dev/null +++ b/packages/base/src/Numeric/Conversion.hs | |||
@@ -0,0 +1,91 @@ | |||
1 | {-# LANGUAGE TypeFamilies #-} | ||
2 | {-# LANGUAGE FlexibleContexts #-} | ||
3 | {-# LANGUAGE FlexibleInstances #-} | ||
4 | {-# LANGUAGE MultiParamTypeClasses #-} | ||
5 | {-# LANGUAGE FunctionalDependencies #-} | ||
6 | {-# LANGUAGE UndecidableInstances #-} | ||
7 | |||
8 | ----------------------------------------------------------------------------- | ||
9 | -- | | ||
10 | -- Module : Numeric.Conversion | ||
11 | -- Copyright : (c) Alberto Ruiz 2010 | ||
12 | -- License : BSD3 | ||
13 | -- Maintainer : Alberto Ruiz | ||
14 | -- Stability : provisional | ||
15 | -- | ||
16 | -- Conversion routines | ||
17 | -- | ||
18 | ----------------------------------------------------------------------------- | ||
19 | {-# OPTIONS_HADDOCK hide #-} | ||
20 | |||
21 | |||
22 | module Numeric.Conversion ( | ||
23 | Complexable(..), RealElement, | ||
24 | module Data.Complex | ||
25 | ) where | ||
26 | |||
27 | import Data.Packed.Internal.Vector | ||
28 | import Data.Packed.Internal.Matrix | ||
29 | import Data.Complex | ||
30 | import Control.Arrow((***)) | ||
31 | |||
32 | ------------------------------------------------------------------- | ||
33 | |||
34 | -- | Supported single-double precision type pairs | ||
35 | class (Element s, Element d) => Precision s d | s -> d, d -> s where | ||
36 | double2FloatG :: Vector d -> Vector s | ||
37 | float2DoubleG :: Vector s -> Vector d | ||
38 | |||
39 | instance Precision Float Double where | ||
40 | double2FloatG = double2FloatV | ||
41 | float2DoubleG = float2DoubleV | ||
42 | |||
43 | instance Precision (Complex Float) (Complex Double) where | ||
44 | double2FloatG = asComplex . double2FloatV . asReal | ||
45 | float2DoubleG = asComplex . float2DoubleV . asReal | ||
46 | |||
47 | -- | Supported real types | ||
48 | class (Element t, Element (Complex t), RealFloat t | ||
49 | -- , RealOf t ~ t, RealOf (Complex t) ~ t | ||
50 | ) | ||
51 | => RealElement t | ||
52 | |||
53 | instance RealElement Double | ||
54 | instance RealElement Float | ||
55 | |||
56 | |||
57 | -- | Structures that may contain complex numbers | ||
58 | class Complexable c where | ||
59 | toComplex' :: (RealElement e) => (c e, c e) -> c (Complex e) | ||
60 | fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e) | ||
61 | comp' :: (RealElement e) => c e -> c (Complex e) | ||
62 | single' :: Precision a b => c b -> c a | ||
63 | double' :: Precision a b => c a -> c b | ||
64 | |||
65 | |||
66 | instance Complexable Vector where | ||
67 | toComplex' = toComplexV | ||
68 | fromComplex' = fromComplexV | ||
69 | comp' v = toComplex' (v,constantD 0 (dim v)) | ||
70 | single' = double2FloatG | ||
71 | double' = float2DoubleG | ||
72 | |||
73 | |||
74 | -- | creates a complex vector from vectors with real and imaginary parts | ||
75 | toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a) | ||
76 | toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] | ||
77 | |||
78 | -- | the inverse of 'toComplex' | ||
79 | fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a) | ||
80 | fromComplexV z = (r,i) where | ||
81 | [r,i] = toColumns $ reshape 2 $ asReal z | ||
82 | |||
83 | |||
84 | instance Complexable Matrix where | ||
85 | toComplex' = uncurry $ liftMatrix2 $ curry toComplex' | ||
86 | fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z | ||
87 | where c = cols z | ||
88 | comp' = liftMatrix comp' | ||
89 | single' = liftMatrix single' | ||
90 | double' = liftMatrix double' | ||
91 | |||
diff --git a/packages/base/src/Numeric/LAPACK.hs b/packages/base/src/Numeric/LAPACK.hs new file mode 100644 index 0000000..08cd759 --- /dev/null +++ b/packages/base/src/Numeric/LAPACK.hs | |||
@@ -0,0 +1,552 @@ | |||
1 | ----------------------------------------------------------------------------- | ||
2 | -- | | ||
3 | -- Module : Numeric.LAPACK | ||
4 | -- Copyright : (c) Alberto Ruiz 2006-14 | ||
5 | -- License : BSD3 | ||
6 | -- | ||
7 | -- Maintainer : Alberto Ruiz | ||
8 | -- Stability : provisional | ||
9 | -- | ||
10 | -- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>). | ||
11 | -- | ||
12 | ----------------------------------------------------------------------------- | ||
13 | |||
14 | module Numeric.LAPACK ( | ||
15 | -- * Matrix product | ||
16 | multiplyR, multiplyC, multiplyF, multiplyQ, | ||
17 | -- * Linear systems | ||
18 | linearSolveR, linearSolveC, | ||
19 | lusR, lusC, | ||
20 | cholSolveR, cholSolveC, | ||
21 | linearSolveLSR, linearSolveLSC, | ||
22 | linearSolveSVDR, linearSolveSVDC, | ||
23 | -- * SVD | ||
24 | svR, svRd, svC, svCd, | ||
25 | svdR, svdRd, svdC, svdCd, | ||
26 | thinSVDR, thinSVDRd, thinSVDC, thinSVDCd, | ||
27 | rightSVR, rightSVC, leftSVR, leftSVC, | ||
28 | -- * Eigensystems | ||
29 | eigR, eigC, eigS, eigS', eigH, eigH', | ||
30 | eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH, | ||
31 | -- * LU | ||
32 | luR, luC, | ||
33 | -- * Cholesky | ||
34 | cholS, cholH, mbCholS, mbCholH, | ||
35 | -- * QR | ||
36 | qrR, qrC, qrgrR, qrgrC, | ||
37 | -- * Hessenberg | ||
38 | hessR, hessC, | ||
39 | -- * Schur | ||
40 | schurR, schurC | ||
41 | ) where | ||
42 | |||
43 | import Data.Packed.Development | ||
44 | import Data.Packed | ||
45 | import Data.Packed.Internal | ||
46 | import Numeric.Conversion | ||
47 | |||
48 | import Foreign.Ptr(nullPtr) | ||
49 | import Foreign.C.Types | ||
50 | import Control.Monad(when) | ||
51 | import System.IO.Unsafe(unsafePerformIO) | ||
52 | |||
53 | ----------------------------------------------------------------------------------- | ||
54 | |||
55 | foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM | ||
56 | foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM | ||
57 | foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM | ||
58 | foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM | ||
59 | |||
60 | isT Matrix{order = ColumnMajor} = 0 | ||
61 | isT Matrix{order = RowMajor} = 1 | ||
62 | |||
63 | tt x@Matrix{order = ColumnMajor} = x | ||
64 | tt x@Matrix{order = RowMajor} = trans x | ||
65 | |||
66 | multiplyAux f st a b = unsafePerformIO $ do | ||
67 | when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ | ||
68 | show (rows a,cols a) ++ " x " ++ show (rows b, cols b) | ||
69 | s <- createMatrix ColumnMajor (rows a) (cols b) | ||
70 | app3 (f (isT a) (isT b)) mat (tt a) mat (tt b) mat s st | ||
71 | return s | ||
72 | |||
73 | -- | Matrix product based on BLAS's /dgemm/. | ||
74 | multiplyR :: Matrix Double -> Matrix Double -> Matrix Double | ||
75 | multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b | ||
76 | |||
77 | -- | Matrix product based on BLAS's /zgemm/. | ||
78 | multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
79 | multiplyC a b = multiplyAux zgemmc "zgemmc" a b | ||
80 | |||
81 | -- | Matrix product based on BLAS's /sgemm/. | ||
82 | multiplyF :: Matrix Float -> Matrix Float -> Matrix Float | ||
83 | multiplyF a b = multiplyAux sgemmc "sgemmc" a b | ||
84 | |||
85 | -- | Matrix product based on BLAS's /cgemm/. | ||
86 | multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float) | ||
87 | multiplyQ a b = multiplyAux cgemmc "cgemmc" a b | ||
88 | |||
89 | ----------------------------------------------------------------------------- | ||
90 | foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM | ||
91 | foreign import ccall unsafe "svd_l_C" zgesvd :: TCMCMVCM | ||
92 | foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TMMVM | ||
93 | foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TCMCMVCM | ||
94 | |||
95 | -- | Full SVD of a real matrix using LAPACK's /dgesvd/. | ||
96 | svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) | ||
97 | svdR = svdAux dgesvd "svdR" . fmat | ||
98 | |||
99 | -- | Full SVD of a real matrix using LAPACK's /dgesdd/. | ||
100 | svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) | ||
101 | svdRd = svdAux dgesdd "svdRdd" . fmat | ||
102 | |||
103 | -- | Full SVD of a complex matrix using LAPACK's /zgesvd/. | ||
104 | svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) | ||
105 | svdC = svdAux zgesvd "svdC" . fmat | ||
106 | |||
107 | -- | Full SVD of a complex matrix using LAPACK's /zgesdd/. | ||
108 | svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) | ||
109 | svdCd = svdAux zgesdd "svdCdd" . fmat | ||
110 | |||
111 | svdAux f st x = unsafePerformIO $ do | ||
112 | u <- createMatrix ColumnMajor r r | ||
113 | s <- createVector (min r c) | ||
114 | v <- createMatrix ColumnMajor c c | ||
115 | app4 f mat x mat u vec s mat v st | ||
116 | return (u,s,trans v) | ||
117 | where r = rows x | ||
118 | c = cols x | ||
119 | |||
120 | |||
121 | -- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'. | ||
122 | thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) | ||
123 | thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat | ||
124 | |||
125 | -- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'. | ||
126 | thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) | ||
127 | thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat | ||
128 | |||
129 | -- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'. | ||
130 | thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) | ||
131 | thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat | ||
132 | |||
133 | -- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'. | ||
134 | thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) | ||
135 | thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat | ||
136 | |||
137 | thinSVDAux f st x = unsafePerformIO $ do | ||
138 | u <- createMatrix ColumnMajor r q | ||
139 | s <- createVector q | ||
140 | v <- createMatrix ColumnMajor q c | ||
141 | app4 f mat x mat u vec s mat v st | ||
142 | return (u,s,trans v) | ||
143 | where r = rows x | ||
144 | c = cols x | ||
145 | q = min r c | ||
146 | |||
147 | |||
148 | -- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'. | ||
149 | svR :: Matrix Double -> Vector Double | ||
150 | svR = svAux dgesvd "svR" . fmat | ||
151 | |||
152 | -- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'. | ||
153 | svC :: Matrix (Complex Double) -> Vector Double | ||
154 | svC = svAux zgesvd "svC" . fmat | ||
155 | |||
156 | -- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'. | ||
157 | svRd :: Matrix Double -> Vector Double | ||
158 | svRd = svAux dgesdd "svRd" . fmat | ||
159 | |||
160 | -- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'. | ||
161 | svCd :: Matrix (Complex Double) -> Vector Double | ||
162 | svCd = svAux zgesdd "svCd" . fmat | ||
163 | |||
164 | svAux f st x = unsafePerformIO $ do | ||
165 | s <- createVector q | ||
166 | app2 g mat x vec s st | ||
167 | return s | ||
168 | where r = rows x | ||
169 | c = cols x | ||
170 | q = min r c | ||
171 | g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr | ||
172 | |||
173 | |||
174 | -- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'. | ||
175 | rightSVR :: Matrix Double -> (Vector Double, Matrix Double) | ||
176 | rightSVR = rightSVAux dgesvd "rightSVR" . fmat | ||
177 | |||
178 | -- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'. | ||
179 | rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | ||
180 | rightSVC = rightSVAux zgesvd "rightSVC" . fmat | ||
181 | |||
182 | rightSVAux f st x = unsafePerformIO $ do | ||
183 | s <- createVector q | ||
184 | v <- createMatrix ColumnMajor c c | ||
185 | app3 g mat x vec s mat v st | ||
186 | return (s,trans v) | ||
187 | where r = rows x | ||
188 | c = cols x | ||
189 | q = min r c | ||
190 | g ra ca pa = f ra ca pa 0 0 nullPtr | ||
191 | |||
192 | |||
193 | -- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'. | ||
194 | leftSVR :: Matrix Double -> (Matrix Double, Vector Double) | ||
195 | leftSVR = leftSVAux dgesvd "leftSVR" . fmat | ||
196 | |||
197 | -- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'. | ||
198 | leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double) | ||
199 | leftSVC = leftSVAux zgesvd "leftSVC" . fmat | ||
200 | |||
201 | leftSVAux f st x = unsafePerformIO $ do | ||
202 | u <- createMatrix ColumnMajor r r | ||
203 | s <- createVector q | ||
204 | app3 g mat x mat u vec s st | ||
205 | return (u,s) | ||
206 | where r = rows x | ||
207 | c = cols x | ||
208 | q = min r c | ||
209 | g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr | ||
210 | |||
211 | ----------------------------------------------------------------------------- | ||
212 | |||
213 | foreign import ccall unsafe "eig_l_R" dgeev :: TMMCVM | ||
214 | foreign import ccall unsafe "eig_l_C" zgeev :: TCMCMCVCM | ||
215 | foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> TMVM | ||
216 | foreign import ccall unsafe "eig_l_H" zheev :: CInt -> TCMVCM | ||
217 | |||
218 | eigAux f st m = unsafePerformIO $ do | ||
219 | l <- createVector r | ||
220 | v <- createMatrix ColumnMajor r r | ||
221 | app3 g mat m vec l mat v st | ||
222 | return (l,v) | ||
223 | where r = rows m | ||
224 | g ra ca pa = f ra ca pa 0 0 nullPtr | ||
225 | |||
226 | |||
227 | -- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/. | ||
228 | -- The eigenvectors are the columns of v. The eigenvalues are not sorted. | ||
229 | eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) | ||
230 | eigC = eigAux zgeev "eigC" . fmat | ||
231 | |||
232 | eigOnlyAux f st m = unsafePerformIO $ do | ||
233 | l <- createVector r | ||
234 | app2 g mat m vec l st | ||
235 | return l | ||
236 | where r = rows m | ||
237 | g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr | ||
238 | |||
239 | -- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'. | ||
240 | -- The eigenvalues are not sorted. | ||
241 | eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double) | ||
242 | eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat | ||
243 | |||
244 | -- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/. | ||
245 | -- The eigenvectors are the columns of v. The eigenvalues are not sorted. | ||
246 | eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) | ||
247 | eigR m = (s', v'') | ||
248 | where (s,v) = eigRaux (fmat m) | ||
249 | s' = fixeig1 s | ||
250 | v' = toRows $ trans v | ||
251 | v'' = fromColumns $ fixeig (toList s') v' | ||
252 | |||
253 | eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double) | ||
254 | eigRaux m = unsafePerformIO $ do | ||
255 | l <- createVector r | ||
256 | v <- createMatrix ColumnMajor r r | ||
257 | app3 g mat m vec l mat v "eigR" | ||
258 | return (l,v) | ||
259 | where r = rows m | ||
260 | g ra ca pa = dgeev ra ca pa 0 0 nullPtr | ||
261 | |||
262 | fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s)) | ||
263 | where r = dim s | ||
264 | |||
265 | fixeig [] _ = [] | ||
266 | fixeig [_] [v] = [comp' v] | ||
267 | fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) | ||
268 | | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1, mapVector negate v2) : fixeig r vs | ||
269 | | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) | ||
270 | fixeig _ _ = error "fixeig with impossible inputs" | ||
271 | |||
272 | |||
273 | -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. | ||
274 | -- The eigenvalues are not sorted. | ||
275 | eigOnlyR :: Matrix Double -> Vector (Complex Double) | ||
276 | eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat | ||
277 | |||
278 | |||
279 | ----------------------------------------------------------------------------- | ||
280 | |||
281 | eigSHAux f st m = unsafePerformIO $ do | ||
282 | l <- createVector r | ||
283 | v <- createMatrix ColumnMajor r r | ||
284 | app3 f mat m vec l mat v st | ||
285 | return (l,v) | ||
286 | where r = rows m | ||
287 | |||
288 | -- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/. | ||
289 | -- The eigenvectors are the columns of v. | ||
290 | -- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order). | ||
291 | eigS :: Matrix Double -> (Vector Double, Matrix Double) | ||
292 | eigS m = (s', fliprl v) | ||
293 | where (s,v) = eigS' (fmat m) | ||
294 | s' = fromList . reverse . toList $ s | ||
295 | |||
296 | -- | 'eigS' in ascending order | ||
297 | eigS' :: Matrix Double -> (Vector Double, Matrix Double) | ||
298 | eigS' = eigSHAux (dsyev 1) "eigS'" . fmat | ||
299 | |||
300 | -- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/. | ||
301 | -- The eigenvectors are the columns of v. | ||
302 | -- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order). | ||
303 | eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | ||
304 | eigH m = (s', fliprl v) | ||
305 | where (s,v) = eigH' (fmat m) | ||
306 | s' = fromList . reverse . toList $ s | ||
307 | |||
308 | -- | 'eigH' in ascending order | ||
309 | eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | ||
310 | eigH' = eigSHAux (zheev 1) "eigH'" . fmat | ||
311 | |||
312 | |||
313 | -- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'. | ||
314 | -- The eigenvalues are sorted in descending order. | ||
315 | eigOnlyS :: Matrix Double -> Vector Double | ||
316 | eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat | ||
317 | |||
318 | -- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'. | ||
319 | -- The eigenvalues are sorted in descending order. | ||
320 | eigOnlyH :: Matrix (Complex Double) -> Vector Double | ||
321 | eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat | ||
322 | |||
323 | vrev = flatten . flipud . reshape 1 | ||
324 | |||
325 | ----------------------------------------------------------------------------- | ||
326 | foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM | ||
327 | foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM | ||
328 | foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM | ||
329 | foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM | ||
330 | |||
331 | linearSolveSQAux f st a b | ||
332 | | n1==n2 && n1==r = unsafePerformIO $ do | ||
333 | s <- createMatrix ColumnMajor r c | ||
334 | app3 f mat a mat b mat s st | ||
335 | return s | ||
336 | | otherwise = error $ st ++ " of nonsquare matrix" | ||
337 | where n1 = rows a | ||
338 | n2 = cols a | ||
339 | r = rows b | ||
340 | c = cols b | ||
341 | |||
342 | -- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'. | ||
343 | linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double | ||
344 | linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b) | ||
345 | |||
346 | -- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'. | ||
347 | linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
348 | linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b) | ||
349 | |||
350 | |||
351 | -- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'. | ||
352 | cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double | ||
353 | cholSolveR a b = linearSolveSQAux dpotrs "cholSolveR" (fmat a) (fmat b) | ||
354 | |||
355 | -- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'. | ||
356 | cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
357 | cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b) | ||
358 | |||
359 | ----------------------------------------------------------------------------------- | ||
360 | foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM | ||
361 | foreign import ccall unsafe "linearSolveLSC_l" zgels :: TCMCMCM | ||
362 | foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM | ||
363 | foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TCMCMCM | ||
364 | |||
365 | linearSolveAux f st a b = unsafePerformIO $ do | ||
366 | r <- createMatrix ColumnMajor (max m n) nrhs | ||
367 | app3 f mat a mat b mat r st | ||
368 | return r | ||
369 | where m = rows a | ||
370 | n = cols a | ||
371 | nrhs = cols b | ||
372 | |||
373 | -- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'. | ||
374 | linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double | ||
375 | linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ | ||
376 | linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b) | ||
377 | |||
378 | -- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'. | ||
379 | linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
380 | linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ | ||
381 | linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b) | ||
382 | |||
383 | -- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. 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. | ||
384 | linearSolveSVDR :: Maybe Double -- ^ rcond | ||
385 | -> Matrix Double -- ^ coefficient matrix | ||
386 | -> Matrix Double -- ^ right hand sides (as columns) | ||
387 | -> Matrix Double -- ^ solution vectors (as columns) | ||
388 | linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ | ||
389 | linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b) | ||
390 | linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b) | ||
391 | |||
392 | -- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. 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. | ||
393 | linearSolveSVDC :: Maybe Double -- ^ rcond | ||
394 | -> Matrix (Complex Double) -- ^ coefficient matrix | ||
395 | -> Matrix (Complex Double) -- ^ right hand sides (as columns) | ||
396 | -> Matrix (Complex Double) -- ^ solution vectors (as columns) | ||
397 | linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ | ||
398 | linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b) | ||
399 | linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b) | ||
400 | |||
401 | ----------------------------------------------------------------------------------- | ||
402 | foreign import ccall unsafe "chol_l_H" zpotrf :: TCMCM | ||
403 | foreign import ccall unsafe "chol_l_S" dpotrf :: TMM | ||
404 | |||
405 | cholAux f st a = do | ||
406 | r <- createMatrix ColumnMajor n n | ||
407 | app2 f mat a mat r st | ||
408 | return r | ||
409 | where n = rows a | ||
410 | |||
411 | -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/. | ||
412 | cholH :: Matrix (Complex Double) -> Matrix (Complex Double) | ||
413 | cholH = unsafePerformIO . cholAux zpotrf "cholH" . fmat | ||
414 | |||
415 | -- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/. | ||
416 | cholS :: Matrix Double -> Matrix Double | ||
417 | cholS = unsafePerformIO . cholAux dpotrf "cholS" . fmat | ||
418 | |||
419 | -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version). | ||
420 | mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double)) | ||
421 | mbCholH = unsafePerformIO . mbCatch . cholAux zpotrf "cholH" . fmat | ||
422 | |||
423 | -- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/ ('Maybe' version). | ||
424 | mbCholS :: Matrix Double -> Maybe (Matrix Double) | ||
425 | mbCholS = unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat | ||
426 | |||
427 | ----------------------------------------------------------------------------------- | ||
428 | foreign import ccall unsafe "qr_l_R" dgeqr2 :: TMVM | ||
429 | foreign import ccall unsafe "qr_l_C" zgeqr2 :: TCMCVCM | ||
430 | |||
431 | -- | QR factorization of a real matrix, using LAPACK's /dgeqr2/. | ||
432 | qrR :: Matrix Double -> (Matrix Double, Vector Double) | ||
433 | qrR = qrAux dgeqr2 "qrR" . fmat | ||
434 | |||
435 | -- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/. | ||
436 | qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) | ||
437 | qrC = qrAux zgeqr2 "qrC" . fmat | ||
438 | |||
439 | qrAux f st a = unsafePerformIO $ do | ||
440 | r <- createMatrix ColumnMajor m n | ||
441 | tau <- createVector mn | ||
442 | app3 f mat a vec tau mat r st | ||
443 | return (r,tau) | ||
444 | where | ||
445 | m = rows a | ||
446 | n = cols a | ||
447 | mn = min m n | ||
448 | |||
449 | foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM | ||
450 | foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM | ||
451 | |||
452 | -- | build rotation from reflectors | ||
453 | qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double | ||
454 | qrgrR = qrgrAux dorgqr "qrgrR" | ||
455 | -- | build rotation from reflectors | ||
456 | qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double) | ||
457 | qrgrC = qrgrAux zungqr "qrgrC" | ||
458 | |||
459 | qrgrAux f st n (a, tau) = unsafePerformIO $ do | ||
460 | res <- createMatrix ColumnMajor (rows a) n | ||
461 | app3 f mat (fmat a) vec (subVector 0 n tau') mat res st | ||
462 | return res | ||
463 | where | ||
464 | tau' = vjoin [tau, constantD 0 n] | ||
465 | |||
466 | ----------------------------------------------------------------------------------- | ||
467 | foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM | ||
468 | foreign import ccall unsafe "hess_l_C" zgehrd :: TCMCVCM | ||
469 | |||
470 | -- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/. | ||
471 | hessR :: Matrix Double -> (Matrix Double, Vector Double) | ||
472 | hessR = hessAux dgehrd "hessR" . fmat | ||
473 | |||
474 | -- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/. | ||
475 | hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) | ||
476 | hessC = hessAux zgehrd "hessC" . fmat | ||
477 | |||
478 | hessAux f st a = unsafePerformIO $ do | ||
479 | r <- createMatrix ColumnMajor m n | ||
480 | tau <- createVector (mn-1) | ||
481 | app3 f mat a vec tau mat r st | ||
482 | return (r,tau) | ||
483 | where m = rows a | ||
484 | n = cols a | ||
485 | mn = min m n | ||
486 | |||
487 | ----------------------------------------------------------------------------------- | ||
488 | foreign import ccall unsafe "schur_l_R" dgees :: TMMM | ||
489 | foreign import ccall unsafe "schur_l_C" zgees :: TCMCMCM | ||
490 | |||
491 | -- | Schur factorization of a square real matrix, using LAPACK's /dgees/. | ||
492 | schurR :: Matrix Double -> (Matrix Double, Matrix Double) | ||
493 | schurR = schurAux dgees "schurR" . fmat | ||
494 | |||
495 | -- | Schur factorization of a square complex matrix, using LAPACK's /zgees/. | ||
496 | schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double)) | ||
497 | schurC = schurAux zgees "schurC" . fmat | ||
498 | |||
499 | schurAux f st a = unsafePerformIO $ do | ||
500 | u <- createMatrix ColumnMajor n n | ||
501 | s <- createMatrix ColumnMajor n n | ||
502 | app3 f mat a mat u mat s st | ||
503 | return (u,s) | ||
504 | where n = rows a | ||
505 | |||
506 | ----------------------------------------------------------------------------------- | ||
507 | foreign import ccall unsafe "lu_l_R" dgetrf :: TMVM | ||
508 | foreign import ccall unsafe "lu_l_C" zgetrf :: TCMVCM | ||
509 | |||
510 | -- | LU factorization of a general real matrix, using LAPACK's /dgetrf/. | ||
511 | luR :: Matrix Double -> (Matrix Double, [Int]) | ||
512 | luR = luAux dgetrf "luR" . fmat | ||
513 | |||
514 | -- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/. | ||
515 | luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int]) | ||
516 | luC = luAux zgetrf "luC" . fmat | ||
517 | |||
518 | luAux f st a = unsafePerformIO $ do | ||
519 | lu <- createMatrix ColumnMajor n m | ||
520 | piv <- createVector (min n m) | ||
521 | app3 f mat a vec piv mat lu st | ||
522 | return (lu, map (pred.round) (toList piv)) | ||
523 | where n = rows a | ||
524 | m = cols a | ||
525 | |||
526 | ----------------------------------------------------------------------------------- | ||
527 | type TW a = CInt -> PD -> a | ||
528 | type TQ a = CInt -> CInt -> PC -> a | ||
529 | |||
530 | foreign import ccall unsafe "luS_l_R" dgetrs :: TMVMM | ||
531 | foreign import ccall unsafe "luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt)))) | ||
532 | |||
533 | -- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/. | ||
534 | lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double | ||
535 | lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b) | ||
536 | |||
537 | -- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/. | ||
538 | lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
539 | lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b) | ||
540 | |||
541 | lusAux f st a piv b | ||
542 | | n1==n2 && n2==n =unsafePerformIO $ do | ||
543 | x <- createMatrix ColumnMajor n m | ||
544 | app4 f mat a vec piv' mat b mat x st | ||
545 | return x | ||
546 | | otherwise = error $ st ++ " on LU factorization of nonsquare matrix" | ||
547 | where n1 = rows a | ||
548 | n2 = cols a | ||
549 | n = rows b | ||
550 | m = cols b | ||
551 | piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double | ||
552 | |||