summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra/Algorithms.hs
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2008-10-27 13:03:41 +0000
committerAlberto Ruiz <aruiz@um.es>2008-10-27 13:03:41 +0000
commitedf12982f21c56c21bfc21eb2b2fcbc406838130 (patch)
tree4f9463ceccc49dc5b9dfdf77b16dccef9bc8d3e5 /lib/Numeric/LinearAlgebra/Algorithms.hs
parentd8639b28ec9e83b54b45c987508d270d5469451c (diff)
added luSolve
Diffstat (limited to 'lib/Numeric/LinearAlgebra/Algorithms.hs')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs19
1 files changed, 13 insertions, 6 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index e2fec9d..75f4ba3 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -40,7 +40,7 @@ module Numeric.LinearAlgebra.Algorithms (
40-- ** Schur 40-- ** Schur
41 schur, 41 schur,
42-- ** LU 42-- ** LU
43 lu, 43 lu, luPacked, luSolve,
44-- * Matrix functions 44-- * Matrix functions
45 expm, 45 expm,
46 sqrtm, 46 sqrtm,
@@ -77,8 +77,13 @@ import Foreign.C.Types
77class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where 77class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
78 -- | Singular value decomposition using lapack's dgesvd or zgesvd. 78 -- | Singular value decomposition using lapack's dgesvd or zgesvd.
79 svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) 79 svd :: Matrix t -> (Matrix t, Vector Double, Matrix t)
80 -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
80 luPacked :: Matrix t -> (Matrix t, [Int]) 81 luPacked :: Matrix t -> (Matrix t, [Int])
81 -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv and zgesv. 82 -- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization
83 -- obtained by 'luPacked'.
84 luSolve :: (Matrix t, [Int]) -> Matrix t -> Matrix t
85 -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv.
86 -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
82 -- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK". 87 -- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".
83 linearSolve :: Matrix t -> Matrix t -> Matrix t 88 linearSolve :: Matrix t -> Matrix t -> Matrix t
84 linearSolveSVD :: Matrix t -> Matrix t -> Matrix t 89 linearSolveSVD :: Matrix t -> Matrix t -> Matrix t
@@ -110,13 +115,15 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
110 schur :: Matrix t -> (Matrix t, Matrix t) 115 schur :: Matrix t -> (Matrix t, Matrix t)
111 -- | Conjugate transpose. 116 -- | Conjugate transpose.
112 ctrans :: Matrix t -> Matrix t 117 ctrans :: Matrix t -> Matrix t
118 -- | Matrix product.
113 multiply :: Matrix t -> Matrix t -> Matrix t 119 multiply :: Matrix t -> Matrix t -> Matrix t
114 120
115 121
116instance Field Double where 122instance Field Double where
117 svd = svdR 123 svd = svdR
118 luPacked = luR 124 luPacked = luR
119 linearSolve = linearSolveR 125 luSolve (l_u,perm) = lusR l_u perm
126 linearSolve = linearSolveR -- (luSolve . luPacked) ??
120 linearSolveSVD = linearSolveSVDR Nothing 127 linearSolveSVD = linearSolveSVDR Nothing
121 ctrans = trans 128 ctrans = trans
122 eig = eigR 129 eig = eigR
@@ -130,6 +137,7 @@ instance Field Double where
130instance Field (Complex Double) where 137instance Field (Complex Double) where
131 svd = svdC 138 svd = svdC
132 luPacked = luC 139 luPacked = luC
140 luSolve (l_u,perm) = lusC l_u perm
133 linearSolve = linearSolveC 141 linearSolve = linearSolveC
134 linearSolveSVD = linearSolveSVDC Nothing 142 linearSolveSVD = linearSolveSVDC Nothing
135 ctrans = conj . trans 143 ctrans = conj . trans
@@ -165,7 +173,7 @@ det m | square m = s * (product $ toList $ takeDiag $ lup)
165 where (lup,perm) = luPacked m 173 where (lup,perm) = luPacked m
166 s = signlp (rows m) perm 174 s = signlp (rows m) perm
167 175
168-- | LU factorization of a general matrix using lapack's dgetrf or zgetrf. 176-- | Explicit LU factorization of a general matrix using lapack's dgetrf or zgetrf.
169-- 177--
170-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular, 178-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
171-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation. 179-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
@@ -513,7 +521,7 @@ luFact (l_u,perm) | r <= c = (l ,u ,p, s)
513 521
514-------------------------------------------------- 522--------------------------------------------------
515 523
516-- | euclidean inner product 524-- | Euclidean inner product.
517dot :: (Field t) => Vector t -> Vector t -> t 525dot :: (Field t) => Vector t -> Vector t -> t
518dot u v = multiply r c @@> (0,0) 526dot u v = multiply r c @@> (0,0)
519 where r = asRow u 527 where r = asRow u
@@ -629,5 +637,4 @@ multiplyC3 x y = unsafePerformIO $ multiply3 zgemm "cblas_zgemm" (fmat x) (fmat
629 free palpha 637 free palpha
630 free pbeta 638 free pbeta
631 return s 639 return s
632 -- if toLists s== toLists s then return s else error $ "HORROR " ++ (show (toLists s))
633 | otherwise = error $ st ++ " (matrix product) of nonconformant matrices" 640 | otherwise = error $ st ++ " (matrix product) of nonconformant matrices"