summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2007-10-31 12:47:19 +0000
committerAlberto Ruiz <aruiz@um.es>2007-10-31 12:47:19 +0000
commitbf838323545fe0878382f8f4d41b0f36714afa43 (patch)
tree296198dfcea483ffa933018471a01482066d96b0 /lib
parentf637161ac988979b35ab7254f753a67df8ec812a (diff)
added general expm described in Golub and Van Loan
Diffstat (limited to 'lib')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs49
1 files changed, 43 insertions, 6 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 52f9b6f..794ef69 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -26,7 +26,7 @@ module Numeric.LinearAlgebra.Algorithms (
26-- * Matrix factorizations 26-- * Matrix factorizations
27-- ** Singular value decomposition 27-- ** Singular value decomposition
28 svd, 28 svd,
29 full, economy, 29 full, economy, --thin,
30-- ** Eigensystems 30-- ** Eigensystems
31 eig, eigSH, 31 eig, eigSH,
32-- ** QR 32-- ** QR
@@ -101,8 +101,6 @@ class (Normed (Matrix t), Linear Matrix t) => GenMat t where
101 schur :: Matrix t -> (Matrix t, Matrix t) 101 schur :: Matrix t -> (Matrix t, Matrix t)
102 -- | Conjugate transpose. 102 -- | Conjugate transpose.
103 ctrans :: Matrix t -> Matrix t 103 ctrans :: Matrix t -> Matrix t
104 -- | Matrix exponential (currently implemented only for diagonalizable matrices).
105 expm :: Matrix t -> Matrix t
106 104
107 105
108instance GenMat Double where 106instance GenMat Double where
@@ -117,7 +115,6 @@ instance GenMat Double where
117 qr = GSL.unpackQR . qrR 115 qr = GSL.unpackQR . qrR
118 hess = unpackHess hessR 116 hess = unpackHess hessR
119 schur = schurR 117 schur = schurR
120 expm = fst . fromComplex . expm' . comp
121 118
122instance GenMat (Complex Double) where 119instance GenMat (Complex Double) where
123 svd = svdC 120 svd = svdC
@@ -131,7 +128,6 @@ instance GenMat (Complex Double) where
131 qr = unpackQR . qrC 128 qr = unpackQR . qrC
132 hess = unpackHess hessC 129 hess = unpackHess hessC
133 schur = schurC 130 schur = schurC
134 expm = expm'
135 131
136-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev. 132-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.
137-- 133--
@@ -245,7 +241,11 @@ pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` const
245--pnormCM _ _ = error "p norm not yet defined" 241--pnormCM _ _ = error "p norm not yet defined"
246 242
247-- | Objects which have a p-norm. 243-- | Objects which have a p-norm.
248-- Using it you can define convenient shortcuts: @norm2 = pnorm PNorm2@, @frobenius = norm2 . flatten@, etc. 244-- Using it you can define convenient shortcuts:
245--
246-- @norm2 x = pnorm PNorm2 x@
247--
248-- @frobenius m = norm2 . flatten $ m@
249class Normed t where 249class Normed t where
250 pnorm :: NormType -> t -> Double 250 pnorm :: NormType -> t -> Double
251 251
@@ -385,3 +385,40 @@ matFunc :: GenMat t => (Complex Double -> Complex Double) -> Matrix t -> Matrix
385matFunc f m = case diagonalize (complex m) of 385matFunc f m = case diagonalize (complex m) of
386 Just (l,v) -> v `mXm` diag (liftVector f l) `mXm` inv v 386 Just (l,v) -> v `mXm` diag (liftVector f l) `mXm` inv v
387 Nothing -> error "Sorry, matFunc requieres a diagonalizable matrix" 387 Nothing -> error "Sorry, matFunc requieres a diagonalizable matrix"
388
389--------------------------------------------------------------
390
391golubeps :: Integer -> Integer -> Double
392golubeps p q = a * fromIntegral b / fromIntegral c where
393 a = 2^^(3-p-q)
394 b = fact p * fact q
395 c = fact (p+q) * fact (p+q+1)
396 fact n = product [1..n]
397
398epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]]
399
400geps delta = head [ k | (k,g) <- epslist, g<delta]
401
402expGolub m = iterate msq f !! j
403 where j = max 0 $ floor $ log2 $ pnorm Infinity m
404 log2 x = log x / log 2
405 a = m */ fromIntegral (2^j)
406 q = geps eps -- 7 steps
407 eye = ident (rows m)
408 work (k,c,x,n,d) = (k',c',x',n',d')
409 where k' = k+1
410 c' = c * fromIntegral (q-k+1) / fromIntegral ((2*q-k+1)*k)
411 x' = a <> x
412 n' = n `add` (c' `scale` x')
413 d' = d `add` (((-1)^k * c') `scale` x')
414 (_,_,_,n,d) = iterate work (1,1,eye,eye,eye) !! q
415 f = linearSolve d n
416 msq m = m <> m
417 (<>) = multiply
418 v */ x = scale (recip x) v
419
420{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
421 based on a scaled Pade approximation.
422-}
423expm :: GenMat t => Matrix t -> Matrix t
424expm = expGolub