summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra/Algorithms.hs
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2009-04-18 18:46:57 +0000
committerAlberto Ruiz <aruiz@um.es>2009-04-18 18:46:57 +0000
commit3ccfccc82aa1cf374af1b822087fb6c4fd41b3c3 (patch)
tree7a150b0559bb770cced0d2d82e4964d8310691d1 /lib/Numeric/LinearAlgebra/Algorithms.hs
parent1b5b7b8252fddcda816722b55eae1d5f1ca9a80e (diff)
Algorithms module reorganized to fix documentation structure
Diffstat (limited to 'lib/Numeric/LinearAlgebra/Algorithms.hs')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs205
1 files changed, 124 insertions, 81 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index d978fa4..e22246f 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -3,7 +3,7 @@
3----------------------------------------------------------------------------- 3-----------------------------------------------------------------------------
4{- | 4{- |
5Module : Numeric.LinearAlgebra.Algorithms 5Module : Numeric.LinearAlgebra.Algorithms
6Copyright : (c) Alberto Ruiz 2006-7 6Copyright : (c) Alberto Ruiz 2006-9
7License : GPL-style 7License : GPL-style
8 8
9Maintainer : Alberto Ruiz (aruiz at um dot es) 9Maintainer : Alberto Ruiz (aruiz at um dot es)
@@ -20,27 +20,33 @@ imported from "Numeric.LinearAlgebra.LAPACK".
20----------------------------------------------------------------------------- 20-----------------------------------------------------------------------------
21 21
22module Numeric.LinearAlgebra.Algorithms ( 22module Numeric.LinearAlgebra.Algorithms (
23-- * Linear Systems 23-- * Supported types
24 Field(),
25-- * Products
24 multiply, dot, 26 multiply, dot,
27 outer, kronecker,
28-- * Linear Systems
25 linearSolve, 29 linearSolve,
30 luSolve,
31 linearSolveSVD,
26 inv, pinv, 32 inv, pinv,
27 pinvTol, det, rank, rcond, 33 det, rank, rcond,
28-- * Matrix factorizations 34-- * Matrix factorizations
29-- ** Singular value decomposition 35-- ** Singular value decomposition
30 svd, 36 svd,
31 full, economy, --thin, 37 full, economy, --thin,
32-- ** Eigensystems 38-- ** Eigensystems
33 eig, eigSH, 39 eig, eigSH, eigSH',
34-- ** QR 40-- ** QR
35 qr, 41 qr,
36-- ** Cholesky 42-- ** Cholesky
37 chol, 43 chol, cholSH,
38-- ** Hessenberg 44-- ** Hessenberg
39 hess, 45 hess,
40-- ** Schur 46-- ** Schur
41 schur, 47 schur,
42-- ** LU 48-- ** LU
43 lu, luPacked, luSolve, 49 lu, luPacked,
44-- * Matrix functions 50-- * Matrix functions
45 expm, 51 expm,
46 sqrtm, 52 sqrtm,
@@ -53,11 +59,10 @@ module Numeric.LinearAlgebra.Algorithms (
53-- * Misc 59-- * Misc
54 ctrans, 60 ctrans,
55 eps, i, 61 eps, i,
56 outer, kronecker,
57-- * Util 62-- * Util
58 haussholder, 63 haussholder,
59 unpackQR, unpackHess, 64 unpackQR, unpackHess,
60 Field(linearSolveSVD,eigSH',cholSH) 65 pinvTol
61) where 66) where
62 67
63 68
@@ -71,82 +76,120 @@ import Data.List(foldl1')
71import Data.Array 76import Data.Array
72 77
73 78
74
75-- | Auxiliary typeclass used to define generic computations for both real and complex matrices. 79-- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
76class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where 80class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
77 -- | Singular value decomposition using lapack's dgesvd or zgesvd. 81 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
78 svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) 82 luPacked' :: Matrix t -> (Matrix t, [Int])
79 -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. 83 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t
80 luPacked :: Matrix t -> (Matrix t, [Int]) 84 linearSolve' :: Matrix t -> Matrix t -> Matrix t
81 -- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization 85 linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
82 -- obtained by 'luPacked'. 86 eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
83 luSolve :: (Matrix t, [Int]) -> Matrix t -> Matrix t 87 eigSH'' :: Matrix t -> (Vector Double, Matrix t)
84 -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv. 88 cholSH' :: Matrix t -> Matrix t
85 -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system. 89 qr' :: Matrix t -> (Matrix t, Matrix t)
86 -- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK". 90 hess' :: Matrix t -> (Matrix t, Matrix t)
87 linearSolve :: Matrix t -> Matrix t -> Matrix t 91 schur' :: Matrix t -> (Matrix t, Matrix t)
88 linearSolveSVD :: Matrix t -> Matrix t -> Matrix t 92 ctrans' :: Matrix t -> Matrix t
89 -- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev. 93 multiply' :: Matrix t -> Matrix t -> Matrix t
90 -- 94
91 -- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ 95
92 eig :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) 96-- | Singular value decomposition using lapack's dgesvd or zgesvd.
93 -- | Similar to eigSH without checking that the input matrix is hermitian or symmetric. 97svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
94 eigSH' :: Matrix t -> (Vector Double, Matrix t) 98svd = svd'
95 -- | Similar to chol without checking that the input matrix is hermitian or symmetric. 99
96 cholSH :: Matrix t -> Matrix t 100-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
97 -- | QR factorization using lapack's dgeqr2 or zgeqr2. 101luPacked :: Field t => Matrix t -> (Matrix t, [Int])
98 -- 102luPacked = luPacked'
99 -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. 103
100 qr :: Matrix t -> (Matrix t, Matrix t) 104-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization
101 -- | Hessenberg factorization using lapack's dgehrd or zgehrd. 105-- obtained by 'luPacked'.
102 -- 106luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
103 -- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary 107luSolve = luSolve'
104 -- and h is in upper Hessenberg form. 108
105 hess :: Matrix t -> (Matrix t, Matrix t) 109-- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv.
106 -- | Schur factorization using lapack's dgees or zgees. 110-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
107 -- 111-- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".
108 -- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary 112linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
109 -- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is 113linearSolve = linearSolve'
110 -- upper triangular in 2x2 blocks. 114linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
111 -- 115linearSolveSVD = linearSolveSVD'
112 -- \"Anything that the Jordan decomposition can do, the Schur decomposition 116
113 -- can do better!\" (Van Loan) 117-- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev.
114 schur :: Matrix t -> (Matrix t, Matrix t) 118--
115 -- | Conjugate transpose. 119-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
116 ctrans :: Matrix t -> Matrix t 120eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
117 -- | Matrix product. 121eig = eig'
118 multiply :: Matrix t -> Matrix t -> Matrix t 122
123-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric.
124eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t)
125eigSH' = eigSH''
126
127-- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.
128cholSH :: Field t => Matrix t -> Matrix t
129cholSH = cholSH'
130
131-- | QR factorization using lapack's dgeqr2 or zgeqr2.
132--
133-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
134qr :: Field t => Matrix t -> (Matrix t, Matrix t)
135qr = qr'
136
137-- | Hessenberg factorization using lapack's dgehrd or zgehrd.
138--
139-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary
140-- and h is in upper Hessenberg form.
141hess :: Field t => Matrix t -> (Matrix t, Matrix t)
142hess = hess'
143
144-- | Schur factorization using lapack's dgees or zgees.
145--
146-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary
147-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
148-- upper triangular in 2x2 blocks.
149--
150-- \"Anything that the Jordan decomposition can do, the Schur decomposition
151-- can do better!\" (Van Loan)
152schur :: Field t => Matrix t -> (Matrix t, Matrix t)
153schur = schur'
154
155-- | Generic conjugate transpose.
156ctrans :: Field t => Matrix t -> Matrix t
157ctrans = ctrans'
158
159-- | Matrix product.
160multiply :: Field t => Matrix t -> Matrix t -> Matrix t
161multiply = multiply'
119 162
120 163
121instance Field Double where 164instance Field Double where
122 svd = svdR 165 svd' = svdR
123 luPacked = luR 166 luPacked' = luR
124 luSolve (l_u,perm) = lusR l_u perm 167 luSolve' (l_u,perm) = lusR l_u perm
125 linearSolve = linearSolveR -- (luSolve . luPacked) ?? 168 linearSolve' = linearSolveR -- (luSolve . luPacked) ??
126 linearSolveSVD = linearSolveSVDR Nothing 169 linearSolveSVD' = linearSolveSVDR Nothing
127 ctrans = trans 170 ctrans' = trans
128 eig = eigR 171 eig' = eigR
129 eigSH' = eigS 172 eigSH'' = eigS
130 cholSH = cholS 173 cholSH' = cholS
131 qr = unpackQR . qrR 174 qr' = unpackQR . qrR
132 hess = unpackHess hessR 175 hess' = unpackHess hessR
133 schur = schurR 176 schur' = schurR
134 multiply = multiplyR 177 multiply' = multiplyR
135 178
136instance Field (Complex Double) where 179instance Field (Complex Double) where
137 svd = svdC 180 svd' = svdC
138 luPacked = luC 181 luPacked' = luC
139 luSolve (l_u,perm) = lusC l_u perm 182 luSolve' (l_u,perm) = lusC l_u perm
140 linearSolve = linearSolveC 183 linearSolve' = linearSolveC
141 linearSolveSVD = linearSolveSVDC Nothing 184 linearSolveSVD' = linearSolveSVDC Nothing
142 ctrans = conj . trans 185 ctrans' = conj . trans
143 eig = eigC 186 eig' = eigC
144 eigSH' = eigH 187 eigSH'' = eigH
145 cholSH = cholH 188 cholSH' = cholH
146 qr = unpackQR . qrC 189 qr' = unpackQR . qrC
147 hess = unpackHess hessC 190 hess' = unpackHess hessC
148 schur = schurC 191 schur' = schurC
149 multiply = multiplyC 192 multiply' = multiplyC
150 193
151 194
152-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev. 195-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.
@@ -193,8 +236,8 @@ pinv m = linearSolveSVD m (ident (rows m))
193-- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@. 236-- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@.
194full :: Element t 237full :: Element t
195 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) 238 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)
196full svd' m = (u, d ,v) where 239full svdFun m = (u, d ,v) where
197 (u,s,v) = svd' m 240 (u,s,v) = svdFun m
198 d = diagRect s r c 241 d = diagRect s r c
199 r = rows m 242 r = rows m
200 c = cols m 243 c = cols m
@@ -204,8 +247,8 @@ full svd' m = (u, d ,v) where
204-- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. 247-- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@.
205economy :: Element t 248economy :: Element t
206 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) 249 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t)
207economy svd' m = (u', subVector 0 d s, v') where 250economy svdFun m = (u', subVector 0 d s, v') where
208 (u,s,v) = svd' m 251 (u,s,v) = svdFun m
209 sl@(g:_) = toList s 252 sl@(g:_) = toList s
210 s' = fromList . filter (>tol) $ sl 253 s' = fromList . filter (>tol) $ sl
211 t = 1 254 t = 1