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