summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs116
-rw-r--r--lib/Numeric/LinearAlgebra/Instances.hs4
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK.hs225
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c159
-rw-r--r--lib/Numeric/LinearAlgebra/Tests.hs21
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/Instances.hs8
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/Properties.hs78
7 files changed, 456 insertions, 155 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index e0cc0a0..1544d7c 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -10,10 +10,9 @@ Maintainer : Alberto Ruiz (aruiz at um dot es)
10Stability : provisional 10Stability : provisional
11Portability : uses ffi 11Portability : uses ffi
12 12
13A generic interface for some common functions. Using it we can write higher level algorithms 13Generic interface for the most common functions. Using it we can write higher level algorithms and testing properties for both real and complex matrices.
14and testing properties both for real and complex matrices.
15 14
16In any case, the specific functions for particular base types can also be explicitly 15Specific functions for particular base types can also be explicitly
17imported from "Numeric.LinearAlgebra.LAPACK". 16imported from "Numeric.LinearAlgebra.LAPACK".
18 17
19-} 18-}
@@ -34,7 +33,11 @@ module Numeric.LinearAlgebra.Algorithms (
34-- * Matrix factorizations 33-- * Matrix factorizations
35-- ** Singular value decomposition 34-- ** Singular value decomposition
36 svd, 35 svd,
37 full, economy, --thin, 36 fullSVD,
37 thinSVD,
38 compactSVD,
39 sv,
40 leftSV, rightSV,
38-- ** Eigensystems 41-- ** Eigensystems
39 eig, eigSH, eigSH', 42 eig, eigSH, eigSH',
40-- ** QR 43-- ** QR
@@ -54,6 +57,7 @@ module Numeric.LinearAlgebra.Algorithms (
54-- * Nullspace 57-- * Nullspace
55 nullspacePrec, 58 nullspacePrec,
56 nullVector, 59 nullVector,
60 nullspaceSVD,
57-- * Norms 61-- * Norms
58 Normed(..), NormType(..), 62 Normed(..), NormType(..),
59-- * Misc 63-- * Misc
@@ -63,8 +67,8 @@ module Numeric.LinearAlgebra.Algorithms (
63 haussholder, 67 haussholder,
64 unpackQR, unpackHess, 68 unpackQR, unpackHess,
65 pinvTol, 69 pinvTol,
66 rankSVD, ranksv, 70 ranksv,
67 nullspaceSVD 71 full, economy
68) where 72) where
69 73
70 74
@@ -80,6 +84,8 @@ import Data.Array
80-- | Auxiliary typeclass used to define generic computations for both real and complex matrices. 84-- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
81class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where 85class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
82 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) 86 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
87 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
88 sv' :: Matrix t -> Vector Double
83 luPacked' :: Matrix t -> (Matrix t, [Int]) 89 luPacked' :: Matrix t -> (Matrix t, [Int])
84 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t 90 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t
85 linearSolve' :: Matrix t -> Matrix t -> Matrix t 91 linearSolve' :: Matrix t -> Matrix t -> Matrix t
@@ -94,10 +100,18 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
94 multiply' :: Matrix t -> Matrix t -> Matrix t 100 multiply' :: Matrix t -> Matrix t -> Matrix t
95 101
96 102
97-- | Singular value decomposition using lapack's dgesvd or zgesvd. 103-- | Full singular value decomposition using lapack's dgesdd or zgesdd.
98svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) 104svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
99svd = svd' 105svd = svd'
100 106
107-- | Thin singular value decomposition using lapack's dgesvd or zgesvd.
108thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
109thinSVD = thinSVD'
110
111-- | Singular values using lapack's dgesvd or zgesvd.
112sv :: Field t => Matrix t -> Vector Double
113sv = sv'
114
101-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. 115-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
102luPacked :: Field t => Matrix t -> (Matrix t, [Int]) 116luPacked :: Field t => Matrix t -> (Matrix t, [Int])
103luPacked = luPacked' 117luPacked = luPacked'
@@ -163,7 +177,9 @@ multiply = multiply'
163 177
164 178
165instance Field Double where 179instance Field Double where
166 svd' = svdR 180 svd' = svdRd
181 thinSVD' = thinSVDRd
182 sv' = svR
167 luPacked' = luR 183 luPacked' = luR
168 luSolve' (l_u,perm) = lusR l_u perm 184 luSolve' (l_u,perm) = lusR l_u perm
169 linearSolve' = linearSolveR -- (luSolve . luPacked) ?? 185 linearSolve' = linearSolveR -- (luSolve . luPacked) ??
@@ -178,7 +194,9 @@ instance Field Double where
178 multiply' = multiplyR 194 multiply' = multiplyR
179 195
180instance Field (Complex Double) where 196instance Field (Complex Double) where
181 svd' = svdC 197 svd' = svdCd
198 thinSVD' = thinSVDCd
199 sv' = svC
182 luPacked' = luC 200 luPacked' = luC
183 luSolve' (l_u,perm) = lusC l_u perm 201 luSolve' (l_u,perm) = lusC l_u perm
184 linearSolve' = linearSolveC 202 linearSolve' = linearSolveC
@@ -232,36 +250,61 @@ inv m | square m = m `linearSolve` ident (rows m)
232pinv :: Field t => Matrix t -> Matrix t 250pinv :: Field t => Matrix t -> Matrix t
233pinv m = linearSolveSVD m (ident (rows m)) 251pinv m = linearSolveSVD m (ident (rows m))
234 252
235-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. 253
236-- 254{-# DEPRECATED full "use fullSVD instead" #-}
237-- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@.
238full :: Element t
239 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)
240full svdFun m = (u, d ,v) where 255full svdFun m = (u, d ,v) where
241 (u,s,v) = svdFun m 256 (u,s,v) = svdFun m
242 d = diagRect s r c 257 d = diagRect s r c
243 r = rows m 258 r = rows m
244 c = cols m 259 c = cols m
245 260
246-- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations. 261-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
247-- 262--
248-- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. 263-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
249economy :: Element t 264fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
250 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) 265fullSVD m = (u,d,v) where
266 (u,s,v) = svd m
267 d = diagRect s r c
268 r = rows m
269 c = cols m
270
271
272{-# DEPRECATED economy "use compactSVD instead" #-}
251economy svdFun m = (u', subVector 0 d s, v') where 273economy svdFun m = (u', subVector 0 d s, v') where
252 r@(u,s,v) = svdFun m 274 (u,s,v) = svdFun m
253 d = rankSVD (1*eps) m r `max` 1 275 d = rankSVD (1*eps) m s `max` 1
276 u' = takeColumns d u
277 v' = takeColumns d v
278
279-- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations.
280--
281-- If @(u,s,v) = compactSVD m@ then @m == u \<> diag s \<> trans v@.
282compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
283compactSVD m = (u', subVector 0 d s, v') where
284 (u,s,v) = thinSVD m
285 d = rankSVD (1*eps) m s `max` 1
254 u' = takeColumns d u 286 u' = takeColumns d u
255 v' = takeColumns d v 287 v' = takeColumns d v
256 288
289vertical m = rows m >= cols m
290
291-- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd.
292rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
293rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
294 | otherwise = let (_,s,v) = svd m in (s,v)
295
296-- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd.
297leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
298leftSV m | vertical m = let (u,s,_) = svd m in (u,s)
299 | otherwise = let (u,s,_) = thinSVD m in (u,s)
257 300
258-- | Numeric rank of a matrix from the SVD decomposition. 301-- | Numeric rank of a matrix from the SVD decomposition.
259rankSVD :: Element t 302rankSVD :: Element t
260 => Double -- ^ numeric zero (e.g. 1*'eps') 303 => Double -- ^ numeric zero (e.g. 1*'eps')
261 -> Matrix t -- ^ input matrix m 304 -> Matrix t -- ^ input matrix m
262 -> (Matrix t, Vector Double, Matrix t) -- ^ 'svd' of m 305 -> Vector Double -- ^ 'sv' of m
263 -> Int -- ^ rank of m 306 -> Int -- ^ rank of m
264rankSVD teps m (_,s,_) = ranksv teps (max (rows m) (cols m)) (toList s) 307rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s)
265 308
266-- | Numeric rank of a matrix from its singular values. 309-- | Numeric rank of a matrix from its singular values.
267ranksv :: Double -- ^ numeric zero (e.g. 1*'eps') 310ranksv :: Double -- ^ numeric zero (e.g. 1*'eps')
@@ -316,12 +359,12 @@ pnormCV PNorm1 = norm1 . mapVector magnitude
316pnormCV Infinity = vectorMax . mapVector magnitude 359pnormCV Infinity = vectorMax . mapVector magnitude
317--pnormCV _ = error "pnormCV not yet defined" 360--pnormCV _ = error "pnormCV not yet defined"
318 361
319pnormRM PNorm2 m = head (toList s) where (_,s,_) = svdR m 362pnormRM PNorm2 m = sv m @> 0
320pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m 363pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m
321pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m) 364pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m)
322--pnormRM _ _ = error "p norm not yet defined" 365--pnormRM _ _ = error "p norm not yet defined"
323 366
324pnormCM PNorm2 m = head (toList s) where (_,s,_) = svdC m 367pnormCM PNorm2 m = sv m @> 0
325pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m 368pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m
326pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m) 369pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m)
327--pnormCM _ _ = error "p norm not yet defined" 370--pnormCM _ _ = error "p norm not yet defined"
@@ -354,9 +397,9 @@ nullspaceSVD :: Field t
354 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'), 397 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
355 -- or Right \"theoretical\" matrix rank. 398 -- or Right \"theoretical\" matrix rank.
356 -> Matrix t -- ^ input matrix m 399 -> Matrix t -- ^ input matrix m
357 -> (Matrix t, Vector Double, Matrix t) -- ^ 'svd' of m 400 -> (Vector Double, Matrix t) -- ^ 'rightSV' of m
358 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace 401 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
359nullspaceSVD hint a z@(_,_,v) = vs where 402nullspaceSVD hint a (s,v) = vs where
360 r = rows a 403 r = rows a
361 c = cols a 404 c = cols a
362 tol = case hint of 405 tol = case hint of
@@ -364,9 +407,8 @@ nullspaceSVD hint a z@(_,_,v) = vs where
364 _ -> eps 407 _ -> eps
365 k = case hint of 408 k = case hint of
366 Right t -> t 409 Right t -> t
367 _ -> rankSVD tol a z 410 _ -> rankSVD tol a s
368 nsol = c - k 411 vs = drop k $ toRows $ ctrans v
369 vs = drop k (toColumns v)
370 412
371 413
372-- | The nullspace of a matrix. See also 'nullspaceSVD'. 414-- | The nullspace of a matrix. See also 'nullspaceSVD'.
@@ -374,10 +416,7 @@ nullspacePrec :: Field t
374 => Double -- ^ relative tolerance in 'eps' units 416 => Double -- ^ relative tolerance in 'eps' units
375 -> Matrix t -- ^ input matrix 417 -> Matrix t -- ^ input matrix
376 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace 418 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
377nullspacePrec t m = ns where 419nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m)
378 r@(_,_,v) = svd m
379 k = rankSVD (t*eps) m r
380 ns = drop k $ toRows $ ctrans v
381 420
382-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance 1*'eps'. 421-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance 1*'eps'.
383nullVector :: Field t => Matrix t -> Vector t 422nullVector :: Field t => Matrix t -> Vector t
@@ -405,7 +444,7 @@ multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').
405-} 444-}
406--pinvTol :: Double -> Matrix Double -> Matrix Double 445--pinvTol :: Double -> Matrix Double -> Matrix Double
407pinvTol t m = v' `mXm` diag s' `mXm` trans u' where 446pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
408 (u,s,v) = svdR m 447 (u,s,v) = thinSVDRd m
409 sl@(g:_) = toList s 448 sl@(g:_) = toList s
410 s' = fromList . map rec $ sl 449 s' = fromList . map rec $ sl
411 rec x = if x < g*tol then 1 else 1/x 450 rec x = if x < g*tol then 1 else 1/x
@@ -460,15 +499,14 @@ uH (pq, tau) = (p,h)
460 499
461-------------------------------------------------------------------------- 500--------------------------------------------------------------------------
462 501
463-- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD. 502-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
464rcond :: Field t => Matrix t -> Double 503rcond :: Field t => Matrix t -> Double
465rcond m = last s / head s 504rcond m = last s / head s
466 where (_,s',_) = svd m 505 where s = toList (sv m)
467 s = toList s'
468 506
469-- | Number of linearly independent rows or columns. 507-- | Number of linearly independent rows or columns.
470rank :: Field t => Matrix t -> Int 508rank :: Field t => Matrix t -> Int
471rank m = rankSVD eps m (svd m) 509rank m = rankSVD eps m (sv m)
472 510
473{- 511{-
474expm' m = case diagonalize (complex m) of 512expm' m = case diagonalize (complex m) of
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
index 04bbd68..a96cf15 100644
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -190,8 +190,8 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati
190 190
191--------------------------------------------------------------- 191---------------------------------------------------------------
192 192
193instance (Storable a) => Monoid (Vector a) where 193instance (Storable a, Num (Vector a)) => Monoid (Vector a) where
194 mempty = V { dim = 0, fptr = undefined } 194 mempty = 0 { dim = 0 }
195 mappend a b = mconcat [a,b] 195 mappend a b = mconcat [a,b]
196 mconcat = j . filter ((>0).dim) 196 mconcat = j . filter ((>0).dim)
197 where j [] = mempty 197 where j [] = mempty
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 30a3d3b..2eae9b7 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -1,4 +1,3 @@
1{-# OPTIONS_GHC #-}
2----------------------------------------------------------------------------- 1-----------------------------------------------------------------------------
3-- | 2-- |
4-- Module : Numeric.LinearAlgebra.LAPACK 3-- Module : Numeric.LinearAlgebra.LAPACK
@@ -9,21 +8,34 @@
9-- Stability : provisional 8-- Stability : provisional
10-- Portability : portable (uses FFI) 9-- Portability : portable (uses FFI)
11-- 10--
12-- Wrappers for a few LAPACK functions (<http://www.netlib.org/lapack>). 11-- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>).
13-- 12--
14----------------------------------------------------------------------------- 13-----------------------------------------------------------------------------
15 14
16module Numeric.LinearAlgebra.LAPACK ( 15module Numeric.LinearAlgebra.LAPACK (
16 -- * Matrix product
17 multiplyR, multiplyC, 17 multiplyR, multiplyC,
18 svdR, svdRdd, svdC, 18 -- * Linear systems
19 eigC, eigR, eigS, eigH, eigS', eigH',
20 linearSolveR, linearSolveC, 19 linearSolveR, linearSolveC,
20 lusR, lusC,
21 linearSolveLSR, linearSolveLSC, 21 linearSolveLSR, linearSolveLSC,
22 linearSolveSVDR, linearSolveSVDC, 22 linearSolveSVDR, linearSolveSVDC,
23 luR, luC, lusR, lusC, 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 -- * LU
31 luR, luC,
32 -- * Cholesky
24 cholS, cholH, 33 cholS, cholH,
34 -- * QR
25 qrR, qrC, 35 qrR, qrC,
36 -- * Hessenberg
26 hessR, hessC, 37 hessR, hessC,
38 -- * Schur
27 schurR, schurC 39 schurR, schurC
28) where 40) where
29 41
@@ -61,28 +73,27 @@ multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Compl
61multiplyC a b = multiplyAux zgemmc "zgemmc" a b 73multiplyC a b = multiplyAux zgemmc "zgemmc" a b
62 74
63----------------------------------------------------------------------------- 75-----------------------------------------------------------------------------
64foreign import ccall "LAPACK/lapack-aux.h svd_l_R" dgesvd :: TMMVM 76foreign import ccall "svd_l_R" dgesvd :: TMMVM
65foreign import ccall "LAPACK/lapack-aux.h svd_l_C" zgesvd :: TCMCMVCM 77foreign import ccall "svd_l_C" zgesvd :: TCMCMVCM
66foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd" dgesdd :: TMMVM 78foreign import ccall "svd_l_Rdd" dgesdd :: TMMVM
79foreign import ccall "svd_l_Cdd" zgesdd :: TCMCMVCM
67 80
68-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix. 81-- | Full SVD of a real matrix using LAPACK's /dgesvd/.
69--
70-- @(u,s,v)=full svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
71svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) 82svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
72svdR = svdAux dgesvd "svdR" . fmat 83svdR = svdAux dgesvd "svdR" . fmat
73 84
74-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix. 85-- | Full SVD of a real matrix using LAPACK's /dgesdd/.
75-- 86svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
76-- @(u,s,v)=full svdRdd m@ so that @m=u \<\> s \<\> 'trans' v@. 87svdRd = svdAux dgesdd "svdRdd" . fmat
77svdRdd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
78svdRdd = svdAux dgesdd "svdRdd" . fmat
79 88
80-- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix. 89-- | Full SVD of a complex matrix using LAPACK's /zgesvd/.
81--
82-- @(u,s,v)=full svdC m@ so that @m=u \<\> comp s \<\> 'trans' v@.
83svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) 90svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
84svdC = svdAux zgesvd "svdC" . fmat 91svdC = svdAux zgesvd "svdC" . fmat
85 92
93-- | Full SVD of a complex matrix using LAPACK's /zgesdd/.
94svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
95svdCd = svdAux zgesdd "svdCdd" . fmat
96
86svdAux f st x = unsafePerformIO $ do 97svdAux f st x = unsafePerformIO $ do
87 u <- createMatrix ColumnMajor r r 98 u <- createMatrix ColumnMajor r r
88 s <- createVector (min r c) 99 s <- createVector (min r c)
@@ -92,7 +103,104 @@ svdAux f st x = unsafePerformIO $ do
92 where r = rows x 103 where r = rows x
93 c = cols x 104 c = cols x
94 105
106
107-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'.
108thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
109thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat
110
111-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'.
112thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
113thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat
114
115-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'.
116thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
117thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat
118
119-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'.
120thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
121thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat
122
123thinSVDAux f st x = unsafePerformIO $ do
124 u <- createMatrix ColumnMajor r q
125 s <- createVector q
126 v <- createMatrix ColumnMajor q c
127 app4 f mat x mat u vec s mat v st
128 return (u,s,trans v)
129 where r = rows x
130 c = cols x
131 q = min r c
132
133
134-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'.
135svR :: Matrix Double -> Vector Double
136svR = svAux dgesvd "svR" . fmat
137
138-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'.
139svC :: Matrix (Complex Double) -> Vector Double
140svC = svAux zgesvd "svC" . fmat
141
142-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'.
143svRd :: Matrix Double -> Vector Double
144svRd = svAux dgesdd "svRd" . fmat
145
146-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'.
147svCd :: Matrix (Complex Double) -> Vector Double
148svCd = svAux zgesdd "svCd" . fmat
149
150svAux f st x = unsafePerformIO $ do
151 s <- createVector q
152 app2 g mat x vec s st
153 return s
154 where r = rows x
155 c = cols x
156 q = min r c
157 g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr
158
159
160-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'.
161rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
162rightSVR = rightSVAux dgesvd "rightSVR" . fmat
163
164-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'.
165rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
166rightSVC = rightSVAux zgesvd "rightSVC" . fmat
167
168rightSVAux f st x = unsafePerformIO $ do
169 s <- createVector q
170 v <- createMatrix ColumnMajor c c
171 app3 g mat x vec s mat v st
172 return (s,trans v)
173 where r = rows x
174 c = cols x
175 q = min r c
176 g ra ca pa = f ra ca pa 0 0 nullPtr
177
178
179-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'.
180leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
181leftSVR = leftSVAux dgesvd "leftSVR" . fmat
182
183-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'.
184leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
185leftSVC = leftSVAux zgesvd "leftSVC" . fmat
186
187leftSVAux f st x = unsafePerformIO $ do
188 u <- createMatrix ColumnMajor r r
189 s <- createVector q
190 app3 g mat x mat u vec s st
191 return (u,s)
192 where r = rows x
193 c = cols x
194 q = min r c
195 g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr
196
95----------------------------------------------------------------------------- 197-----------------------------------------------------------------------------
198
199foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
200foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM
201foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM
202foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
203
96eigAux f st m 204eigAux f st m
97 | r == 1 = (fromList [flatten m `at` 0], singleton 1) 205 | r == 1 = (fromList [flatten m `at` 0], singleton 1)
98 | otherwise = unsafePerformIO $ do 206 | otherwise = unsafePerformIO $ do
@@ -104,28 +212,13 @@ eigAux f st m
104 where r = rows m 212 where r = rows m
105 213
106 214
107foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM 215-- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/.
108foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM 216-- The eigenvectors are the columns of v. The eigenvalues are not sorted.
109foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM
110foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
111
112-- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix:
113--
114-- if @(l,v)=eigC m@ then @m \<\> v = v \<\> diag l@.
115--
116-- The eigenvectors are the columns of v.
117-- The eigenvalues are not sorted.
118eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) 217eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
119eigC = eigAux zgeev "eigC" . fmat 218eigC = eigAux zgeev "eigC" . fmat
120 219
121----------------------------------------------------------------------------- 220-- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/.
122 221-- The eigenvectors are the columns of v. The eigenvalues are not sorted.
123-- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix:
124--
125-- if @(l,v)=eigR m@ then @m \<\> v = v \<\> diag l@.
126--
127-- The eigenvectors are the columns of v.
128-- The eigenvalues are not sorted.
129eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) 222eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
130eigR m = (s', v'') 223eigR m = (s', v'')
131 where (s,v) = eigRaux (fmat m) 224 where (s,v) = eigRaux (fmat m)
@@ -155,17 +248,16 @@ fixeig _ _ = error "fixeig with impossible inputs"
155 248
156----------------------------------------------------------------------------- 249-----------------------------------------------------------------------------
157 250
158-- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix: 251-- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/.
159--
160-- if @(l,v)=eigSl m@ then @m \<\> v = v \<\> diag l@.
161--
162-- The eigenvectors are the columns of v. 252-- The eigenvectors are the columns of v.
163-- The eigenvalues are sorted in descending order (use eigS' for ascending order). 253-- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order).
164eigS :: Matrix Double -> (Vector Double, Matrix Double) 254eigS :: Matrix Double -> (Vector Double, Matrix Double)
165eigS m = (s', fliprl v) 255eigS m = (s', fliprl v)
166 where (s,v) = eigS' (fmat m) 256 where (s,v) = eigS' (fmat m)
167 s' = fromList . reverse . toList $ s 257 s' = fromList . reverse . toList $ s
168 258
259-- | 'eigS' in ascending order
260eigS' :: Matrix Double -> (Vector Double, Matrix Double)
169eigS' m 261eigS' m
170 | r == 1 = (fromList [flatten m `at` 0], singleton 1) 262 | r == 1 = (fromList [flatten m `at` 0], singleton 1)
171 | otherwise = unsafePerformIO $ do 263 | otherwise = unsafePerformIO $ do
@@ -177,17 +269,16 @@ eigS' m
177 269
178----------------------------------------------------------------------------- 270-----------------------------------------------------------------------------
179 271
180-- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix: 272-- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/.
181--
182-- if @(l,v)=eigH m@ then @m \<\> s v = v \<\> diag l@.
183--
184-- The eigenvectors are the columns of v. 273-- The eigenvectors are the columns of v.
185-- The eigenvalues are sorted in descending order (use eigH' for ascending order). 274-- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order).
186eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) 275eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
187eigH m = (s', fliprl v) 276eigH m = (s', fliprl v)
188 where (s,v) = eigH' (fmat m) 277 where (s,v) = eigH' (fmat m)
189 s' = fromList . reverse . toList $ s 278 s' = fromList . reverse . toList $ s
190 279
280-- | 'eigH' in ascending order
281eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
191eigH' m 282eigH' m
192 | r == 1 = (fromList [realPart (flatten m `at` 0)], singleton 1) 283 | r == 1 = (fromList [realPart (flatten m `at` 0)], singleton 1)
193 | otherwise = unsafePerformIO $ do 284 | otherwise = unsafePerformIO $ do
@@ -212,11 +303,11 @@ linearSolveSQAux f st a b
212 r = rows b 303 r = rows b
213 c = cols b 304 c = cols b
214 305
215-- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition. 306-- | 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'.
216linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double 307linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
217linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b) 308linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b)
218 309
219-- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition. 310-- | 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'.
220linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) 311linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
221linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b) 312linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b)
222 313
@@ -234,17 +325,17 @@ linearSolveAux f st a b = unsafePerformIO $ do
234 n = cols a 325 n = cols a
235 nrhs = cols b 326 nrhs = cols b
236 327
237-- | Wrapper for LAPACK's /dgels/, which obtains the least squared error solution of an overconstrained real linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDR'. 328-- | 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'.
238linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double 329linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
239linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ 330linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $
240 linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b) 331 linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b)
241 332
242-- | Wrapper for LAPACK's /zgels/, which obtains the least squared error solution of an overconstrained complex linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDC'. 333-- | 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'.
243linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) 334linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
244linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ 335linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $
245 linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b) 336 linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b)
246 337
247-- | Wrapper for LAPACK's /dgelss/, which obtains the minimum norm solution to a real linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. 338-- | 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.
248linearSolveSVDR :: Maybe Double -- ^ rcond 339linearSolveSVDR :: Maybe Double -- ^ rcond
249 -> Matrix Double -- ^ coefficient matrix 340 -> Matrix Double -- ^ coefficient matrix
250 -> Matrix Double -- ^ right hand sides (as columns) 341 -> Matrix Double -- ^ right hand sides (as columns)
@@ -253,7 +344,7 @@ linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
253 linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b) 344 linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b)
254linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b) 345linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b)
255 346
256-- | Wrapper for LAPACK's /zgelss/, which obtains the minimum norm solution to a complex linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. 347-- | 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.
257linearSolveSVDC :: Maybe Double -- ^ rcond 348linearSolveSVDC :: Maybe Double -- ^ rcond
258 -> Matrix (Complex Double) -- ^ coefficient matrix 349 -> Matrix (Complex Double) -- ^ coefficient matrix
259 -> Matrix (Complex Double) -- ^ right hand sides (as columns) 350 -> Matrix (Complex Double) -- ^ right hand sides (as columns)
@@ -266,13 +357,11 @@ linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b)
266foreign import ccall "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM 357foreign import ccall "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM
267foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM 358foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM
268 359
269-- | Wrapper for LAPACK's /zpotrf/, which computes the Cholesky factorization of a 360-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
270-- complex Hermitian positive definite matrix.
271cholH :: Matrix (Complex Double) -> Matrix (Complex Double) 361cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
272cholH = cholAux zpotrf "cholH" . fmat 362cholH = cholAux zpotrf "cholH" . fmat
273 363
274-- | Wrapper for LAPACK's /dpotrf/, which computes the Cholesky factorization of a 364-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/.
275-- real symmetric positive definite matrix.
276cholS :: Matrix Double -> Matrix Double 365cholS :: Matrix Double -> Matrix Double
277cholS = cholAux dpotrf "cholS" . fmat 366cholS = cholAux dpotrf "cholS" . fmat
278 367
@@ -286,11 +375,11 @@ cholAux f st a = unsafePerformIO $ do
286foreign import ccall "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM 375foreign import ccall "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM
287foreign import ccall "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM 376foreign import ccall "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM
288 377
289-- | Wrapper for LAPACK's /dgeqr2/, which computes a QR factorization of a real matrix. 378-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/.
290qrR :: Matrix Double -> (Matrix Double, Vector Double) 379qrR :: Matrix Double -> (Matrix Double, Vector Double)
291qrR = qrAux dgeqr2 "qrR" . fmat 380qrR = qrAux dgeqr2 "qrR" . fmat
292 381
293-- | Wrapper for LAPACK's /zgeqr2/, which computes a QR factorization of a complex matrix. 382-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/.
294qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) 383qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
295qrC = qrAux zgeqr2 "qrC" . fmat 384qrC = qrAux zgeqr2 "qrC" . fmat
296 385
@@ -307,11 +396,11 @@ qrAux f st a = unsafePerformIO $ do
307foreign import ccall "LAPACK/lapack-aux.h hess_l_R" dgehrd :: TMVM 396foreign import ccall "LAPACK/lapack-aux.h hess_l_R" dgehrd :: TMVM
308foreign import ccall "LAPACK/lapack-aux.h hess_l_C" zgehrd :: TCMCVCM 397foreign import ccall "LAPACK/lapack-aux.h hess_l_C" zgehrd :: TCMCVCM
309 398
310-- | Wrapper for LAPACK's /dgehrd/, which computes a Hessenberg factorization of a square real matrix. 399-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/.
311hessR :: Matrix Double -> (Matrix Double, Vector Double) 400hessR :: Matrix Double -> (Matrix Double, Vector Double)
312hessR = hessAux dgehrd "hessR" . fmat 401hessR = hessAux dgehrd "hessR" . fmat
313 402
314-- | Wrapper for LAPACK's /zgehrd/, which computes a Hessenberg factorization of a square complex matrix. 403-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/.
315hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) 404hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
316hessC = hessAux zgehrd "hessC" . fmat 405hessC = hessAux zgehrd "hessC" . fmat
317 406
@@ -328,11 +417,11 @@ hessAux f st a = unsafePerformIO $ do
328foreign import ccall "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM 417foreign import ccall "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM
329foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM 418foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM
330 419
331-- | Wrapper for LAPACK's /dgees/, which computes a Schur factorization of a square real matrix. 420-- | Schur factorization of a square real matrix, using LAPACK's /dgees/.
332schurR :: Matrix Double -> (Matrix Double, Matrix Double) 421schurR :: Matrix Double -> (Matrix Double, Matrix Double)
333schurR = schurAux dgees "schurR" . fmat 422schurR = schurAux dgees "schurR" . fmat
334 423
335-- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix. 424-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/.
336schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double)) 425schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
337schurC = schurAux zgees "schurC" . fmat 426schurC = schurAux zgees "schurC" . fmat
338 427
@@ -347,11 +436,11 @@ schurAux f st a = unsafePerformIO $ do
347foreign import ccall "LAPACK/lapack-aux.h lu_l_R" dgetrf :: TMVM 436foreign import ccall "LAPACK/lapack-aux.h lu_l_R" dgetrf :: TMVM
348foreign import ccall "LAPACK/lapack-aux.h lu_l_C" zgetrf :: TCMVCM 437foreign import ccall "LAPACK/lapack-aux.h lu_l_C" zgetrf :: TCMVCM
349 438
350-- | Wrapper for LAPACK's /dgetrf/, which computes a LU factorization of a general real matrix. 439-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/.
351luR :: Matrix Double -> (Matrix Double, [Int]) 440luR :: Matrix Double -> (Matrix Double, [Int])
352luR = luAux dgetrf "luR" . fmat 441luR = luAux dgetrf "luR" . fmat
353 442
354-- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix. 443-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/.
355luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int]) 444luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
356luC = luAux zgetrf "luC" . fmat 445luC = luAux zgetrf "luC" . fmat
357 446
@@ -370,11 +459,11 @@ type TQ a = CInt -> CInt -> PC -> a
370foreign import ccall "LAPACK/lapack-aux.h luS_l_R" dgetrs :: TMVMM 459foreign import ccall "LAPACK/lapack-aux.h luS_l_R" dgetrs :: TMVMM
371foreign import ccall "LAPACK/lapack-aux.h luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt)))) 460foreign import ccall "LAPACK/lapack-aux.h luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt))))
372 461
373-- | Wrapper for LAPACK's /dgetrs/, which solves a general real linear system (for several right-hand sides) from a precomputed LU decomposition. 462-- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/.
374lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double 463lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
375lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b) 464lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b)
376 465
377-- | Wrapper for LAPACK's /zgetrs/, which solves a general real linear system (for several right-hand sides) from a precomputed LU decomposition. 466-- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/.
378lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) 467lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
379lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b) 468lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b)
380 469
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index f18bbed..db94cc6 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -48,7 +48,27 @@ int svd_l_R(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
48 integer m = ar; 48 integer m = ar;
49 integer n = ac; 49 integer n = ac;
50 integer q = MIN(m,n); 50 integer q = MIN(m,n);
51 REQUIRES(ur==m && uc==m && sn==q && vr==n && vc==n,BAD_SIZE); 51 REQUIRES(sn==q,BAD_SIZE);
52 REQUIRES(up==NULL || ur==m && (uc==m || uc==q),BAD_SIZE);
53 char* jobu = "A";
54 if (up==NULL) {
55 jobu = "N";
56 } else {
57 if (uc==q) {
58 jobu = "S";
59 }
60 }
61 REQUIRES(vp==NULL || vc==n && (vr==n || vr==q),BAD_SIZE);
62 char* jobvt = "A";
63 integer ldvt = n;
64 if (vp==NULL) {
65 jobvt = "N";
66 } else {
67 if (vr==q) {
68 jobvt = "S";
69 ldvt = q;
70 }
71 }
52 DEBUGMSG("svd_l_R"); 72 DEBUGMSG("svd_l_R");
53 double *B = (double*)malloc(m*n*sizeof(double)); 73 double *B = (double*)malloc(m*n*sizeof(double));
54 CHECK(!B,MEM); 74 CHECK(!B,MEM);
@@ -57,25 +77,21 @@ int svd_l_R(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
57 integer res; 77 integer res;
58 // ask for optimal lwork 78 // ask for optimal lwork
59 double ans; 79 double ans;
60 //printf("ask zgesvd\n"); 80 dgesvd_ (jobu,jobvt,
61 char* job = "A";
62 dgesvd_ (job,job,
63 &m,&n,B,&m, 81 &m,&n,B,&m,
64 sp, 82 sp,
65 up,&m, 83 up,&m,
66 vp,&n, 84 vp,&ldvt,
67 &ans, &lwork, 85 &ans, &lwork,
68 &res); 86 &res);
69 lwork = ceil(ans); 87 lwork = ceil(ans);
70 //printf("ans = %d\n",lwork);
71 double * work = (double*)malloc(lwork*sizeof(double)); 88 double * work = (double*)malloc(lwork*sizeof(double));
72 CHECK(!work,MEM); 89 CHECK(!work,MEM);
73 //printf("dgesdd\n"); 90 dgesvd_ (jobu,jobvt,
74 dgesvd_ (job,job,
75 &m,&n,B,&m, 91 &m,&n,B,&m,
76 sp, 92 sp,
77 up,&m, 93 up,&m,
78 vp,&n, 94 vp,&ldvt,
79 work, &lwork, 95 work, &lwork,
80 &res); 96 &res);
81 CHECK(res,res); 97 CHECK(res,res);
@@ -90,25 +106,36 @@ int svd_l_Rdd(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
90 integer m = ar; 106 integer m = ar;
91 integer n = ac; 107 integer n = ac;
92 integer q = MIN(m,n); 108 integer q = MIN(m,n);
93 REQUIRES(ur==m && uc==m && sn==q && vr==n && vc==n,BAD_SIZE); 109 REQUIRES(sn==q,BAD_SIZE);
110 REQUIRES(up == NULL && vp == NULL
111 || ur==m && vc==n
112 && (uc == q && vr == q
113 || uc == m && vc==n),BAD_SIZE);
114 char* jobz = "A";
115 integer ldvt = n;
116 if (up==NULL) {
117 jobz = "N";
118 } else {
119 if (uc==q && vr == q) {
120 jobz = "S";
121 ldvt = q;
122 }
123 }
94 DEBUGMSG("svd_l_Rdd"); 124 DEBUGMSG("svd_l_Rdd");
95 double *B = (double*)malloc(m*n*sizeof(double)); 125 double *B = (double*)malloc(m*n*sizeof(double));
96 CHECK(!B,MEM); 126 CHECK(!B,MEM);
97 memcpy(B,ap,m*n*sizeof(double)); 127 memcpy(B,ap,m*n*sizeof(double));
98 integer* iwk = (integer*) malloc(8*q*sizeof(int)); 128 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
99 CHECK(!iwk,MEM); 129 CHECK(!iwk,MEM);
100 integer lwk = -1; 130 integer lwk = -1;
101 integer res; 131 integer res;
102 // ask for optimal lwk 132 // ask for optimal lwk
103 double ans; 133 double ans;
104 //printf("ask dgesdd\n"); 134 dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res);
105 dgesdd_ ("A",&m,&n,B,&m,sp,up,&m,vp,&n,&ans,&lwk,iwk,&res); 135 lwk = ans;
106 lwk = 2*ceil(ans); // ????? otherwise 50x100 rejects lwk
107 //printf("lwk = %d\n",lwk);
108 double * workv = (double*)malloc(lwk*sizeof(double)); 136 double * workv = (double*)malloc(lwk*sizeof(double));
109 CHECK(!workv,MEM); 137 CHECK(!workv,MEM);
110 //printf("dgesdd\n"); 138 dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res);
111 dgesdd_ ("A",&m,&n,B,&m,sp,up,&m,vp,&n,workv,&lwk,iwk,&res);
112 CHECK(res,res); 139 CHECK(res,res);
113 free(iwk); 140 free(iwk);
114 free(workv); 141 free(workv);
@@ -120,17 +147,36 @@ int svd_l_Rdd(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
120 147
121// not in clapack.h 148// not in clapack.h
122 149
123int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n, 150int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
124 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, 151 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
125 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, 152 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
126 integer *lwork, doublereal *rwork, integer *info); 153 integer *lwork, doublereal *rwork, integer *info);
127 154
128int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) { 155int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
129 integer m = ar; 156 integer m = ar;
130 integer n = ac; 157 integer n = ac;
131 integer q = MIN(m,n); 158 integer q = MIN(m,n);
132 REQUIRES(ur==m && uc==m && sn==q && vr==n && vc==n,BAD_SIZE); 159 REQUIRES(sn==q,BAD_SIZE);
133 DEBUGMSG("svd_l_C"); 160 REQUIRES(up==NULL || ur==m && (uc==m || uc==q),BAD_SIZE);
161 char* jobu = "A";
162 if (up==NULL) {
163 jobu = "N";
164 } else {
165 if (uc==q) {
166 jobu = "S";
167 }
168 }
169 REQUIRES(vp==NULL || vc==n && (vr==n || vr==q),BAD_SIZE);
170 char* jobvt = "A";
171 integer ldvt = n;
172 if (vp==NULL) {
173 jobvt = "N";
174 } else {
175 if (vr==q) {
176 jobvt = "S";
177 ldvt = q;
178 }
179 }DEBUGMSG("svd_l_C");
134 double *B = (double*)malloc(2*m*n*sizeof(double)); 180 double *B = (double*)malloc(2*m*n*sizeof(double));
135 CHECK(!B,MEM); 181 CHECK(!B,MEM);
136 memcpy(B,ap,m*n*2*sizeof(double)); 182 memcpy(B,ap,m*n*2*sizeof(double));
@@ -141,26 +187,22 @@ int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
141 integer res; 187 integer res;
142 // ask for optimal lwork 188 // ask for optimal lwork
143 doublecomplex ans; 189 doublecomplex ans;
144 //printf("ask zgesvd\n"); 190 zgesvd_ (jobu,jobvt,
145 char* job = "A";
146 zgesvd_ (job,job,
147 &m,&n,(doublecomplex*)B,&m, 191 &m,&n,(doublecomplex*)B,&m,
148 sp, 192 sp,
149 (doublecomplex*)up,&m, 193 (doublecomplex*)up,&m,
150 (doublecomplex*)vp,&n, 194 (doublecomplex*)vp,&ldvt,
151 &ans, &lwork, 195 &ans, &lwork,
152 rwork, 196 rwork,
153 &res); 197 &res);
154 lwork = ceil(ans.r); 198 lwork = ceil(ans.r);
155 //printf("ans = %d\n",lwork);
156 doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double)); 199 doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double));
157 CHECK(!work,MEM); 200 CHECK(!work,MEM);
158 //printf("zgesvd\n"); 201 zgesvd_ (jobu,jobvt,
159 zgesvd_ (job,job,
160 &m,&n,(doublecomplex*)B,&m, 202 &m,&n,(doublecomplex*)B,&m,
161 sp, 203 sp,
162 (doublecomplex*)up,&m, 204 (doublecomplex*)up,&m,
163 (doublecomplex*)vp,&n, 205 (doublecomplex*)vp,&ldvt,
164 work, &lwork, 206 work, &lwork,
165 rwork, 207 rwork,
166 &res); 208 &res);
@@ -171,7 +213,64 @@ int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
171 OK 213 OK
172} 214}
173 215
216int zgesdd_ (char *jobz, integer *m, integer *n,
217 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
218 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
219 integer *lwork, doublereal *rwork, integer* iwork, integer *info);
174 220
221int svd_l_Cdd(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
222 //printf("entro\n");
223 integer m = ar;
224 integer n = ac;
225 integer q = MIN(m,n);
226 REQUIRES(sn==q,BAD_SIZE);
227 REQUIRES(up == NULL && vp == NULL
228 || ur==m && vc==n
229 && (uc == q && vr == q
230 || uc == m && vc==n),BAD_SIZE);
231 char* jobz = "A";
232 integer ldvt = n;
233 if (up==NULL) {
234 jobz = "N";
235 } else {
236 if (uc==q && vr == q) {
237 jobz = "S";
238 ldvt = q;
239 }
240 }
241 DEBUGMSG("svd_l_Cdd");
242 doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
243 CHECK(!B,MEM);
244 memcpy(B,ap,m*n*sizeof(doublecomplex));
245 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
246 CHECK(!iwk,MEM);
247 int lrwk;
248 if (0 && *jobz == 'N') {
249 lrwk = 5*q; // does not work, crash at free below
250 } else {
251 lrwk = 5*q*q + 7*q;
252 }
253 double *rwk = (double*)malloc(lrwk*sizeof(double));;
254 CHECK(!rwk,MEM);
255 //printf("%s %ld %d\n",jobz,q,lrwk);
256 integer lwk = -1;
257 integer res;
258 // ask for optimal lwk
259 doublecomplex ans;
260 zgesdd_ (jobz,&m,&n,B,&m,sp,(doublecomplex*)up,&m,(doublecomplex*)vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
261 lwk = ans.r;
262 //printf("lwk = %ld\n",lwk);
263 doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));
264 CHECK(!workv,MEM);
265 zgesdd_ (jobz,&m,&n,B,&m,sp,(doublecomplex*)up,&m,(doublecomplex*)vp,&ldvt,workv,&lwk,rwk,iwk,&res);
266 //printf("res = %ld\n",res);
267 CHECK(res,res);
268 free(workv); // printf("freed workv\n");
269 free(rwk); // printf("freed rwk\n");
270 free(iwk); // printf("freed iwk\n");
271 free(B); // printf("freed B, salgo\n");
272 OK
273}
175 274
176//////////////////// general complex eigensystem //////////// 275//////////////////// general complex eigensystem ////////////
177 276
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index e591285..efc8c40 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -22,6 +22,7 @@ module Numeric.LinearAlgebra.Tests(
22) where 22) where
23 23
24import Numeric.LinearAlgebra 24import Numeric.LinearAlgebra
25import Numeric.LinearAlgebra.LAPACK
25import Numeric.LinearAlgebra.Tests.Instances 26import Numeric.LinearAlgebra.Tests.Instances
26import Numeric.LinearAlgebra.Tests.Properties 27import Numeric.LinearAlgebra.Tests.Properties
27import Test.HUnit hiding ((~:),test,Testable) 28import Test.HUnit hiding ((~:),test,Testable)
@@ -186,8 +187,24 @@ runTests n = do
186 putStrLn "------ svd" 187 putStrLn "------ svd"
187 test (svdProp1 . rM) 188 test (svdProp1 . rM)
188 test (svdProp1 . cM) 189 test (svdProp1 . cM)
189 test (svdProp2 . rM) 190 test (svdProp1a svdR)
190 test (svdProp2 . cM) 191 test (svdProp1a svdC)
192 test (svdProp1a svdRd)
193 test (svdProp1a svdCd)
194 test (svdProp2 thinSVDR)
195 test (svdProp2 thinSVDC)
196 test (svdProp2 thinSVDRd)
197 test (svdProp2 thinSVDCd)
198 test (svdProp3 . rM)
199 test (svdProp3 . cM)
200 test (svdProp4 . rM)
201 test (svdProp4 . cM)
202 test (svdProp5a)
203 test (svdProp5b)
204 test (svdProp6a)
205 test (svdProp6b)
206 test (svdProp7 . rM)
207 test (svdProp7 . cM)
191 putStrLn "------ eig" 208 putStrLn "------ eig"
192 test (eigSHProp . rHer) 209 test (eigSHProp . rHer)
193 test (eigSHProp . cHer) 210 test (eigSHProp . cHer)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
index 9b18513..4995e39 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -143,8 +143,8 @@ instance (Field a, Arbitrary a) => Arbitrary (WC a) where
143 r = rows m 143 r = rows m
144 c = cols m 144 c = cols m
145 n = min r c 145 n = min r c
146 sv <- replicateM n (choose (1,100)) 146 sv' <- replicateM n (choose (1,100))
147 let s = diagRect (fromList sv) r c 147 let s = diagRect (fromList sv') r c
148 return $ WC (u <> real s <> trans v) 148 return $ WC (u <> real s <> trans v)
149 149
150#if MIN_VERSION_QuickCheck(2,0,0) 150#if MIN_VERSION_QuickCheck(2,0,0)
@@ -160,8 +160,8 @@ instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where
160 Sq m <- arbitrary 160 Sq m <- arbitrary
161 let (u,_,v) = svd m 161 let (u,_,v) = svd m
162 n = rows m 162 n = rows m
163 sv <- replicateM n (choose (1,100)) 163 sv' <- replicateM n (choose (1,100))
164 let s = diag (fromList sv) 164 let s = diag (fromList sv')
165 return $ SqWC (u <> real s <> trans v) 165 return $ SqWC (u <> real s <> trans v)
166 166
167#if MIN_VERSION_QuickCheck(2,0,0) 167#if MIN_VERSION_QuickCheck(2,0,0)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index d4c2770..d4dff34 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -29,7 +29,8 @@ module Numeric.LinearAlgebra.Tests.Properties (
29 pinvProp, 29 pinvProp,
30 detProp, 30 detProp,
31 nullspaceProp, 31 nullspaceProp,
32 svdProp1, svdProp2, 32 svdProp1, svdProp1a, svdProp2, svdProp3, svdProp4,
33 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
33 eigProp, eigSHProp, 34 eigProp, eigSHProp,
34 qrProp, 35 qrProp,
35 hessProp, 36 hessProp,
@@ -41,10 +42,12 @@ module Numeric.LinearAlgebra.Tests.Properties (
41) where 42) where
42 43
43import Numeric.LinearAlgebra 44import Numeric.LinearAlgebra
45import Numeric.LinearAlgebra.LAPACK
46import Debug.Trace
44#include "quickCheckCompat.h" 47#include "quickCheckCompat.h"
45-- import Debug.Trace
46 48
47-- debug x = trace (show x) x 49
50debug x = trace (show x) x
48 51
49-- relative error 52-- relative error
50dist :: (Normed t, Num t) => t -> t -> Double 53dist :: (Normed t, Num t) => t -> t -> Double
@@ -71,7 +74,10 @@ a :~n~: b = dist a b < 10^^(-n)
71 74
72square m = rows m == cols m 75square m = rows m == cols m
73 76
74unitary m = square m && m <> ctrans m |~| ident (rows m) 77-- orthonormal columns
78orthonormal m = ctrans m <> m |~| ident (cols m)
79
80unitary m = square m && orthonormal m
75 81
76hermitian m = square m && m |~| ctrans m 82hermitian m = square m && m |~| ctrans m
77 83
@@ -119,12 +125,64 @@ nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c))
119 r = rows m 125 r = rows m
120 c = cols m - rank m 126 c = cols m - rank m
121 127
122svdProp1 m = u <> real d <> trans v |~| m 128------------------------------------------------------------------
123 && unitary u && unitary v 129
124 where (u,d,v) = full svd m 130-- fullSVD
125 131svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v
126svdProp2 m = (m |~| 0) `trivial` ((m |~| 0) || u <> real (diag s) <> trans v |~| m) 132 where (u,d,v) = fullSVD m
127 where (u,s,v) = economy svd m 133
134svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where
135 (u,s,v) = svdfun m
136 d = diagRect s (rows m) (cols m)
137
138-- thinSVD
139svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m)
140 where (u,s,v) = thinSVDfun m
141
142-- compactSVD
143svdProp3 m = (m |~| u <> real (diag s) <> trans v
144 && orthonormal u && orthonormal v)
145 where (u,s,v) = compactSVD m
146
147svdProp4 m' = m |~| u <> real (diag s) <> trans v
148 && orthonormal u && orthonormal v
149 && (dim s == r || r == 0 && dim s == 1)
150 where (u,s,v) = compactSVD m
151 m = m' <-> m'
152 r = rank m'
153
154svdProp5a m = and (map (s1|~|) [s2,s3,s4,s5,s6]) where
155 s1 = svR m
156 s2 = svRd m
157 (_,s3,_) = svdR m
158 (_,s4,_) = svdRd m
159 (_,s5,_) = thinSVDR m
160 (_,s6,_) = thinSVDRd m
161
162svdProp5b m = and (map (s1|~|) [s2,s3,s4,s5,s6]) where
163 s1 = svC m
164 s2 = svCd m
165 (_,s3,_) = svdC m
166 (_,s4,_) = svdCd m
167 (_,s5,_) = thinSVDC m
168 (_,s6,_) = thinSVDCd m
169
170svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
171 where (u,s,v) = svdR m
172 (s',v') = rightSVR m
173 (u',s'') = leftSVR m
174
175svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
176 where (u,s,v) = svdC m
177 (s',v') = rightSVC m
178 (u',s'') = leftSVC m
179
180svdProp7 m = s |~| s' && u |~| u' && v |~| v'
181 where (u,s,v) = svd m
182 (s',v') = rightSV m
183 (u',s'') = leftSV m
184
185------------------------------------------------------------------
128 186
129eigProp m = complex m <> v |~| v <> diag s 187eigProp m = complex m <> v |~| v <> diag s
130 where (s, v) = eig m 188 where (s, v) = eig m