summaryrefslogtreecommitdiff
path: root/packages/base/src
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-08 12:18:56 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-08 12:18:56 +0200
commit561a6c0e21bb77c21114ccbbd86d3af5ddb5a3f1 (patch)
treed49b67d75d63938229f2d5cbed5c49d06dc02bcf /packages/base/src
parent5992d92357cfd911c8f2e9f5faaa4fd8a323fd9a (diff)
Conversion, LAPACK -> base
Diffstat (limited to 'packages/base/src')
-rw-r--r--packages/base/src/Data/Packed/Development.hs2
-rw-r--r--packages/base/src/Data/Packed/Vector.hs13
-rw-r--r--packages/base/src/Numeric/Conversion.hs91
-rw-r--r--packages/base/src/Numeric/LAPACK.hs552
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
30import Data.Packed.Internal 30import 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
19module Data.Packed.Vector ( 19module 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
29import Data.Packed.Internal.Vector 29import Data.Packed.Internal.Vector
30import Data.Packed.Internal.Matrix
30import Foreign.Storable 31import Foreign.Storable
31 32
32------------------------------------------------------------------- 33-------------------------------------------------------------------
@@ -94,3 +95,13 @@ unzipVector = unzipVectorWith id
94join :: Storable t => [Vector t] -> Vector t 95join :: Storable t => [Vector t] -> Vector t
95join = vjoin 96join = vjoin
96 97
98{- | creates a vector with a given number of equal components:
99
100@> constant 2 7
1017 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
102-}
103constant :: Element a => a -> Int -> Vector a
104-- constant x n = runSTVector (newVector x n)
105constant = 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
22module Numeric.Conversion (
23 Complexable(..), RealElement,
24 module Data.Complex
25) where
26
27import Data.Packed.Internal.Vector
28import Data.Packed.Internal.Matrix
29import Data.Complex
30import Control.Arrow((***))
31
32-------------------------------------------------------------------
33
34-- | Supported single-double precision type pairs
35class (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
39instance Precision Float Double where
40 double2FloatG = double2FloatV
41 float2DoubleG = float2DoubleV
42
43instance Precision (Complex Float) (Complex Double) where
44 double2FloatG = asComplex . double2FloatV . asReal
45 float2DoubleG = asComplex . float2DoubleV . asReal
46
47-- | Supported real types
48class (Element t, Element (Complex t), RealFloat t
49-- , RealOf t ~ t, RealOf (Complex t) ~ t
50 )
51 => RealElement t
52
53instance RealElement Double
54instance RealElement Float
55
56
57-- | Structures that may contain complex numbers
58class 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
66instance 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
75toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a)
76toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]
77
78-- | the inverse of 'toComplex'
79fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a)
80fromComplexV z = (r,i) where
81 [r,i] = toColumns $ reshape 2 $ asReal z
82
83
84instance 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
14module 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
43import Data.Packed.Development
44import Data.Packed
45import Data.Packed.Internal
46import Numeric.Conversion
47
48import Foreign.Ptr(nullPtr)
49import Foreign.C.Types
50import Control.Monad(when)
51import System.IO.Unsafe(unsafePerformIO)
52
53-----------------------------------------------------------------------------------
54
55foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM
56foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM
57foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM
58foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM
59
60isT Matrix{order = ColumnMajor} = 0
61isT Matrix{order = RowMajor} = 1
62
63tt x@Matrix{order = ColumnMajor} = x
64tt x@Matrix{order = RowMajor} = trans x
65
66multiplyAux 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/.
74multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
75multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b
76
77-- | Matrix product based on BLAS's /zgemm/.
78multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
79multiplyC a b = multiplyAux zgemmc "zgemmc" a b
80
81-- | Matrix product based on BLAS's /sgemm/.
82multiplyF :: Matrix Float -> Matrix Float -> Matrix Float
83multiplyF a b = multiplyAux sgemmc "sgemmc" a b
84
85-- | Matrix product based on BLAS's /cgemm/.
86multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float)
87multiplyQ a b = multiplyAux cgemmc "cgemmc" a b
88
89-----------------------------------------------------------------------------
90foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM
91foreign import ccall unsafe "svd_l_C" zgesvd :: TCMCMVCM
92foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TMMVM
93foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TCMCMVCM
94
95-- | Full SVD of a real matrix using LAPACK's /dgesvd/.
96svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
97svdR = svdAux dgesvd "svdR" . fmat
98
99-- | Full SVD of a real matrix using LAPACK's /dgesdd/.
100svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
101svdRd = svdAux dgesdd "svdRdd" . fmat
102
103-- | Full SVD of a complex matrix using LAPACK's /zgesvd/.
104svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
105svdC = svdAux zgesvd "svdC" . fmat
106
107-- | Full SVD of a complex matrix using LAPACK's /zgesdd/.
108svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
109svdCd = svdAux zgesdd "svdCdd" . fmat
110
111svdAux 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\'.
122thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
123thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat
124
125-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'.
126thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
127thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat
128
129-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'.
130thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
131thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat
132
133-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'.
134thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
135thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat
136
137thinSVDAux 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\'.
149svR :: Matrix Double -> Vector Double
150svR = svAux dgesvd "svR" . fmat
151
152-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'.
153svC :: Matrix (Complex Double) -> Vector Double
154svC = svAux zgesvd "svC" . fmat
155
156-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'.
157svRd :: Matrix Double -> Vector Double
158svRd = svAux dgesdd "svRd" . fmat
159
160-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'.
161svCd :: Matrix (Complex Double) -> Vector Double
162svCd = svAux zgesdd "svCd" . fmat
163
164svAux 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\'.
175rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
176rightSVR = 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\'.
179rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
180rightSVC = rightSVAux zgesvd "rightSVC" . fmat
181
182rightSVAux 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\'.
194leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
195leftSVR = 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\'.
198leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
199leftSVC = leftSVAux zgesvd "leftSVC" . fmat
200
201leftSVAux 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
213foreign import ccall unsafe "eig_l_R" dgeev :: TMMCVM
214foreign import ccall unsafe "eig_l_C" zgeev :: TCMCMCVCM
215foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> TMVM
216foreign import ccall unsafe "eig_l_H" zheev :: CInt -> TCMVCM
217
218eigAux 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.
229eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
230eigC = eigAux zgeev "eigC" . fmat
231
232eigOnlyAux 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.
241eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
242eigOnlyC = 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.
246eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
247eigR 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
253eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
254eigRaux 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
262fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s))
263 where r = dim s
264
265fixeig [] _ = []
266fixeig [_] [v] = [comp' v]
267fixeig ((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)
270fixeig _ _ = 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.
275eigOnlyR :: Matrix Double -> Vector (Complex Double)
276eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat
277
278
279-----------------------------------------------------------------------------
280
281eigSHAux 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).
291eigS :: Matrix Double -> (Vector Double, Matrix Double)
292eigS m = (s', fliprl v)
293 where (s,v) = eigS' (fmat m)
294 s' = fromList . reverse . toList $ s
295
296-- | 'eigS' in ascending order
297eigS' :: Matrix Double -> (Vector Double, Matrix Double)
298eigS' = 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).
303eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
304eigH m = (s', fliprl v)
305 where (s,v) = eigH' (fmat m)
306 s' = fromList . reverse . toList $ s
307
308-- | 'eigH' in ascending order
309eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
310eigH' = 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.
315eigOnlyS :: Matrix Double -> Vector Double
316eigOnlyS = 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.
320eigOnlyH :: Matrix (Complex Double) -> Vector Double
321eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat
322
323vrev = flatten . flipud . reshape 1
324
325-----------------------------------------------------------------------------
326foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM
327foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM
328foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM
329foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM
330
331linearSolveSQAux 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'.
343linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
344linearSolveR 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'.
347linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
348linearSolveC 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'.
352cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
353cholSolveR 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'.
356cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
357cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b)
358
359-----------------------------------------------------------------------------------
360foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM
361foreign import ccall unsafe "linearSolveLSC_l" zgels :: TCMCMCM
362foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM
363foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TCMCMCM
364
365linearSolveAux 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'.
374linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
375linearSolveLSR 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'.
379linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
380linearSolveLSC 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.
384linearSolveSVDR :: Maybe Double -- ^ rcond
385 -> Matrix Double -- ^ coefficient matrix
386 -> Matrix Double -- ^ right hand sides (as columns)
387 -> Matrix Double -- ^ solution vectors (as columns)
388linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
389 linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b)
390linearSolveSVDR 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.
393linearSolveSVDC :: 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)
397linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
398 linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b)
399linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b)
400
401-----------------------------------------------------------------------------------
402foreign import ccall unsafe "chol_l_H" zpotrf :: TCMCM
403foreign import ccall unsafe "chol_l_S" dpotrf :: TMM
404
405cholAux 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/.
412cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
413cholH = unsafePerformIO . cholAux zpotrf "cholH" . fmat
414
415-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/.
416cholS :: Matrix Double -> Matrix Double
417cholS = unsafePerformIO . cholAux dpotrf "cholS" . fmat
418
419-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version).
420mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
421mbCholH = unsafePerformIO . mbCatch . cholAux zpotrf "cholH" . fmat
422
423-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/ ('Maybe' version).
424mbCholS :: Matrix Double -> Maybe (Matrix Double)
425mbCholS = unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat
426
427-----------------------------------------------------------------------------------
428foreign import ccall unsafe "qr_l_R" dgeqr2 :: TMVM
429foreign import ccall unsafe "qr_l_C" zgeqr2 :: TCMCVCM
430
431-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/.
432qrR :: Matrix Double -> (Matrix Double, Vector Double)
433qrR = qrAux dgeqr2 "qrR" . fmat
434
435-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/.
436qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
437qrC = qrAux zgeqr2 "qrC" . fmat
438
439qrAux 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
449foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM
450foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM
451
452-- | build rotation from reflectors
453qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double
454qrgrR = qrgrAux dorgqr "qrgrR"
455-- | build rotation from reflectors
456qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double)
457qrgrC = qrgrAux zungqr "qrgrC"
458
459qrgrAux 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-----------------------------------------------------------------------------------
467foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM
468foreign import ccall unsafe "hess_l_C" zgehrd :: TCMCVCM
469
470-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/.
471hessR :: Matrix Double -> (Matrix Double, Vector Double)
472hessR = hessAux dgehrd "hessR" . fmat
473
474-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/.
475hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
476hessC = hessAux zgehrd "hessC" . fmat
477
478hessAux 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-----------------------------------------------------------------------------------
488foreign import ccall unsafe "schur_l_R" dgees :: TMMM
489foreign import ccall unsafe "schur_l_C" zgees :: TCMCMCM
490
491-- | Schur factorization of a square real matrix, using LAPACK's /dgees/.
492schurR :: Matrix Double -> (Matrix Double, Matrix Double)
493schurR = schurAux dgees "schurR" . fmat
494
495-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/.
496schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
497schurC = schurAux zgees "schurC" . fmat
498
499schurAux 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-----------------------------------------------------------------------------------
507foreign import ccall unsafe "lu_l_R" dgetrf :: TMVM
508foreign import ccall unsafe "lu_l_C" zgetrf :: TCMVCM
509
510-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/.
511luR :: Matrix Double -> (Matrix Double, [Int])
512luR = luAux dgetrf "luR" . fmat
513
514-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/.
515luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
516luC = luAux zgetrf "luC" . fmat
517
518luAux 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-----------------------------------------------------------------------------------
527type TW a = CInt -> PD -> a
528type TQ a = CInt -> CInt -> PC -> a
529
530foreign import ccall unsafe "luS_l_R" dgetrs :: TMVMM
531foreign 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/.
534lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
535lusR 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/.
538lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
539lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b)
540
541lusAux 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