summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--examples/tests.hs1
-rw-r--r--lib/Numeric/GSL/Minimization.hs2
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs39
3 files changed, 29 insertions, 13 deletions
diff --git a/examples/tests.hs b/examples/tests.hs
index ab1b485..974ec24 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -204,6 +204,7 @@ eigTest m = complex m <> v |~| v <> diag s
204 204
205eigTestSH m = m <> v |~| v <> real (diag s) 205eigTestSH m = m <> v |~| v <> real (diag s)
206 && orthonormal v 206 && orthonormal v
207 && m |~| v <> real (diag s) <> ctrans v
207 where (s, v) = eigSH m 208 where (s, v) = eigSH m
208 209
209rank m | m |~| 0 = 0 210rank m | m |~| 0 = 0
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index 23ca51b..e93f8cb 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -150,11 +150,9 @@ minimizeConjugateGradient istep minimpar tol maxit f df xi = unsafePerformIO $ d
150 df' = (fromList . df . toList) 150 df' = (fromList . df . toList)
151 fp <- mkVecfun (iv f') 151 fp <- mkVecfun (iv f')
152 dfp <- mkVecVecfun (aux_vTov df') 152 dfp <- mkVecVecfun (aux_vTov df')
153 print "entro"
154 rawpath <- createMIO maxit (n+2) 153 rawpath <- createMIO maxit (n+2)
155 (c_minimizeConjugateGradient fp dfp istep minimpar tol maxit // vec xiv) 154 (c_minimizeConjugateGradient fp dfp istep minimpar tol maxit // vec xiv)
156 "minimizeDerivV" [xiv] 155 "minimizeDerivV" [xiv]
157 print "salgo"
158 let it = round (rawpath @@> (maxit-1,0)) 156 let it = round (rawpath @@> (maxit-1,0))
159 path = takeRows it rawpath 157 path = takeRows it rawpath
160 sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path 158 sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 5492e07..3757f33 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -9,11 +9,11 @@ Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional 9Stability : provisional
10Portability : uses ffi 10Portability : uses ffi
11 11
12A generic interface for a number of essential functions. Using it some higher level algorithms 12A generic interface for some common functions. Using it we can write higher level algorithms
13and testing properties can be written for both real and complex matrices. 13and testing properties both for real and complex matrices.
14 14
15In any case, the specific functions for particular base types can also be explicitly 15In any case, the specific functions for particular base types can also be explicitly
16imported from the LAPACK and GSL.Matrix modules. 16imported from "Numeric.LinearAlgebra.LAPACK".
17 17
18-} 18-}
19----------------------------------------------------------------------------- 19-----------------------------------------------------------------------------
@@ -29,16 +29,20 @@ module Numeric.LinearAlgebra.Algorithms (
29 full, economy, 29 full, economy,
30-- ** Eigensystems 30-- ** Eigensystems
31 eig, eigSH, 31 eig, eigSH,
32-- ** Other 32-- ** QR
33 Numeric.LinearAlgebra.Algorithms.qr, chol, 33 qr,
34-- ** Cholesky
35 chol,
34-- * Nullspace 36-- * Nullspace
35 nullspacePrec, 37 nullspacePrec,
36 nullVector, 38 nullVector,
39-- * Norms
40 Normed(..), NormType(..),
37-- * Misc 41-- * Misc
38 eps, i,
39 ctrans, 42 ctrans,
40 Normed(..), NormType(..), 43 eps, i,
41 GenMat(linearSolveSVD,lu,eigSH'), unpackQR 44-- * Util
45 GenMat(linearSolveSVD,lu,eigSH',cholSH), unpackQR, haussholder
42) where 46) where
43 47
44 48
@@ -53,15 +57,22 @@ import Data.List(foldl1')
53 57
54-- | Auxiliary typeclass used to define generic computations for both real and complex matrices. 58-- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
55class (Linear Matrix t) => GenMat t where 59class (Linear Matrix t) => GenMat t where
60 -- | Singular value decomposition using lapack's dgesvd or zgesvd.
56 svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) 61 svd :: Matrix t -> (Matrix t, Vector Double, Matrix t)
57 lu :: Matrix t -> (Matrix t, Matrix t, [Int], t) 62 lu :: Matrix t -> (Matrix t, Matrix t, [Int], t)
63 -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv and zgesv.
64 -- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".
58 linearSolve :: Matrix t -> Matrix t -> Matrix t 65 linearSolve :: Matrix t -> Matrix t -> Matrix t
59 linearSolveSVD :: Matrix t -> Matrix t -> Matrix t 66 linearSolveSVD :: Matrix t -> Matrix t -> Matrix t
67 -- | eigensystem of general square matrix using lapack's dgeev or zgeev. If @(s,v) = eig m@ then @m <> v =~= v <> diag s@
60 eig :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) 68 eig :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
69 -- | similar to eigSH without checking that the input matrix is hermitian or symmetric.
61 eigSH' :: Matrix t -> (Vector Double, Matrix t) 70 eigSH' :: Matrix t -> (Vector Double, Matrix t)
71 -- | similar to chol without checking that the input matrix is hermitian or symmetric.
62 cholSH :: Matrix t -> Matrix t 72 cholSH :: Matrix t -> Matrix t
73 -- | QR factorization using lapack's dgeqr2 or zgeqr2.
63 qr :: Matrix t -> (Matrix t, Matrix t) 74 qr :: Matrix t -> (Matrix t, Matrix t)
64 -- | conjugate transpose 75 -- | conjugate transpose.
65 ctrans :: Matrix t -> Matrix t 76 ctrans :: Matrix t -> Matrix t
66 77
67instance GenMat Double where 78instance GenMat Double where
@@ -86,12 +97,12 @@ instance GenMat (Complex Double) where
86 cholSH = cholH 97 cholSH = cholH
87 qr = unpackQR . qrC 98 qr = unpackQR . qrC
88 99
89-- | eigensystem of complex hermitian or real symmetric matrix 100-- | eigensystem of complex hermitian or real symmetric matrix using lapack's dsyev or zheev. If @(s,v) = eigSH m@ then @m =~= v <> diag s <> ctrans v@
90eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) 101eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t)
91eigSH m | m `equal` ctrans m = eigSH' m 102eigSH m | m `equal` ctrans m = eigSH' m
92 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" 103 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
93 104
94-- | Cholesky factorization of a positive definite hermitian or symmetric matrix 105-- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf.
95chol :: GenMat t => Matrix t -> Matrix t 106chol :: GenMat t => Matrix t -> Matrix t
96chol m | m `equal` ctrans m = cholSH m 107chol m | m `equal` ctrans m = cholSH m
97 | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" 108 | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
@@ -103,13 +114,16 @@ det m | square m = s * (product $ toList $ takeDiag $ u)
103 | otherwise = error "det of nonsquare matrix" 114 | otherwise = error "det of nonsquare matrix"
104 where (_,u,_,s) = lu m 115 where (_,u,_,s) = lu m
105 116
117-- | Inverse of a square matrix using lapacks' dgesv and zgesv.
106inv :: GenMat t => Matrix t -> Matrix t 118inv :: GenMat t => Matrix t -> Matrix t
107inv m | square m = m `linearSolve` ident (rows m) 119inv m | square m = m `linearSolve` ident (rows m)
108 | otherwise = error "inv of nonsquare matrix" 120 | otherwise = error "inv of nonsquare matrix"
109 121
122-- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss.
110pinv :: GenMat t => Matrix t -> Matrix t 123pinv :: GenMat t => Matrix t -> Matrix t
111pinv m = linearSolveSVD m (ident (rows m)) 124pinv m = linearSolveSVD m (ident (rows m))
112 125
126-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values: if @(u,d,v) = full svd m@ then @m =~= u \<> d \<> trans v@.
113full :: Field t 127full :: Field t
114 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) 128 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)
115full svd m = (u, d ,v) where 129full svd m = (u, d ,v) where
@@ -118,6 +132,7 @@ full svd m = (u, d ,v) where
118 r = rows m 132 r = rows m
119 c = cols m 133 c = cols m
120 134
135-- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations: if @(u,s,v) = economy svd m@ then @m =~= u \<> diag s \<> trans v@.
121economy :: Field t 136economy :: Field t
122 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) 137 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t)
123economy svd m = (u', subVector 0 d s, v') where 138economy svd m = (u', subVector 0 d s, v') where
@@ -263,9 +278,11 @@ pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
263 278
264-- many thanks, quickcheck! 279-- many thanks, quickcheck!
265 280
281haussholder :: (GenMat a) => a -> Vector a -> Matrix a
266haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) 282haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w))
267 where w = asColumn v 283 where w = asColumn v
268 284
285unpackQR :: (GenMat t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
269unpackQR (pq, tau) = (q,r) 286unpackQR (pq, tau) = (q,r)
270 where cs = toColumns pq 287 where cs = toColumns pq
271 m = rows pq 288 m = rows pq