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