summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--lib/Numeric/LinearAlgebra.hs4
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs205
2 files changed, 127 insertions, 82 deletions
diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs
index cc2c1d3..b8fd01e 100644
--- a/lib/Numeric/LinearAlgebra.hs
+++ b/lib/Numeric/LinearAlgebra.hs
@@ -1,7 +1,7 @@
1----------------------------------------------------------------------------- 1-----------------------------------------------------------------------------
2{- | 2{- |
3Module : Numeric.LinearAlgebra 3Module : Numeric.LinearAlgebra
4Copyright : (c) Alberto Ruiz 2006-7 4Copyright : (c) Alberto Ruiz 2006-9
5License : GPL-style 5License : GPL-style
6 6
7Maintainer : Alberto Ruiz (aruiz at um dot es) 7Maintainer : Alberto Ruiz (aruiz at um dot es)
@@ -10,6 +10,8 @@ Portability : uses ffi
10 10
11Basic matrix computations implemented by BLAS, LAPACK and GSL. 11Basic matrix computations implemented by BLAS, LAPACK and GSL.
12 12
13This is module reexports the most comon functions (including "Numeric.LinearAlgebra.Instances").
14
13-} 15-}
14----------------------------------------------------------------------------- 16-----------------------------------------------------------------------------
15module Numeric.LinearAlgebra ( 17module Numeric.LinearAlgebra (
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