summaryrefslogtreecommitdiff
path: root/packages/base/src/Internal/Algorithms.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Internal/Algorithms.hs')
-rw-r--r--packages/base/src/Internal/Algorithms.hs71
1 files changed, 37 insertions, 34 deletions
diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index ee3ddff..c4f1a60 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -4,12 +4,6 @@
4{-# LANGUAGE UndecidableInstances #-} 4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE TypeFamilies #-} 5{-# LANGUAGE TypeFamilies #-}
6 6
7{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
8{-# LANGUAGE CPP #-}
9{-# LANGUAGE MultiParamTypeClasses #-}
10{-# LANGUAGE UndecidableInstances #-}
11{-# LANGUAGE TypeFamilies #-}
12
13----------------------------------------------------------------------------- 7-----------------------------------------------------------------------------
14{- | 8{- |
15Module : Internal.Algorithms 9Module : Internal.Algorithms
@@ -376,8 +370,8 @@ ldlPackedSH x = {-# SCC "ldlPacked" #-} LDL m p
376 (m,p) = ldlPacked' x 370 (m,p) = ldlPacked' x
377 371
378-- | Obtains the LDL decomposition of a matrix in a compact data structure suitable for 'ldlSolve'. 372-- | Obtains the LDL decomposition of a matrix in a compact data structure suitable for 'ldlSolve'.
379ldlPacked :: Field t => Her t -> LDL t 373ldlPacked :: Field t => Herm t -> LDL t
380ldlPacked (Her m) = ldlPackedSH m 374ldlPacked (Herm m) = ldlPackedSH m
381 375
382-- | Solution of a linear system (for several right hand sides) from a precomputed LDL factorization obtained by 'ldlPacked'. 376-- | Solution of a linear system (for several right hand sides) from a precomputed LDL factorization obtained by 'ldlPacked'.
383-- 377--
@@ -461,18 +455,23 @@ fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575]
4613.000 5.000 6.000 4553.000 5.000 6.000
462 456
463-} 457-}
464eigSH :: Field t => Her t -> (Vector Double, Matrix t) 458eigSH :: Field t => Herm t -> (Vector Double, Matrix t)
465eigSH (Her m) = eigSH' m 459eigSH (Herm m) = eigSH' m
466 460
467-- | Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix. 461-- | Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix.
468eigenvaluesSH :: Field t => Her t -> Vector Double 462eigenvaluesSH :: Field t => Herm t -> Vector Double
469eigenvaluesSH (Her m) = eigenvaluesSH' m 463eigenvaluesSH (Herm m) = eigenvaluesSH' m
470 464
471-------------------------------------------------------------- 465--------------------------------------------------------------
472 466
473-- | QR decomposition of a matrix in compact form. (The orthogonal matrix is not explicitly formed.) 467-- | QR decomposition of a matrix in compact form. (The orthogonal matrix is not explicitly formed.)
474data QR t = QR (Matrix t) (Vector t) 468data QR t = QR (Matrix t) (Vector t)
475 469
470instance (NFData t, Numeric t) => NFData (QR t)
471 where
472 rnf (QR m _) = rnf m
473
474
476-- | QR factorization. 475-- | QR factorization.
477-- 476--
478-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. 477-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
@@ -533,12 +532,12 @@ cholSH = cholSH'
533-- | Cholesky factorization of a positive definite hermitian or symmetric matrix. 532-- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
534-- 533--
535-- If @c = chol m@ then @c@ is upper triangular and @m == tr c \<> c@. 534-- If @c = chol m@ then @c@ is upper triangular and @m == tr c \<> c@.
536chol :: Field t => Her t -> Matrix t 535chol :: Field t => Herm t -> Matrix t
537chol (Her m) = {-# SCC "chol" #-} cholSH' m 536chol (Herm m) = {-# SCC "chol" #-} cholSH' m
538 537
539-- | Similar to 'chol', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'. 538-- | Similar to 'chol', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
540mbChol :: Field t => Her t -> Maybe (Matrix t) 539mbChol :: Field t => Herm t -> Maybe (Matrix t)
541mbChol (Her m) = {-# SCC "mbChol" #-} mbCholSH' m 540mbChol (Herm m) = {-# SCC "mbChol" #-} mbCholSH' m
542 541
543 542
544 543
@@ -977,10 +976,10 @@ relativeError norm a b = r
977-- | Generalized symmetric positive definite eigensystem Av = lBv, 976-- | Generalized symmetric positive definite eigensystem Av = lBv,
978-- for A and B symmetric, B positive definite. 977-- for A and B symmetric, B positive definite.
979geigSH :: Field t 978geigSH :: Field t
980 => Her t -- ^ A 979 => Herm t -- ^ A
981 -> Her t -- ^ B 980 -> Herm t -- ^ B
982 -> (Vector Double, Matrix t) 981 -> (Vector Double, Matrix t)
983geigSH (Her a) (Her b) = geigSH' a b 982geigSH (Herm a) (Herm b) = geigSH' a b
984 983
985geigSH' :: Field t 984geigSH' :: Field t
986 => Matrix t -- ^ A 985 => Matrix t -- ^ A
@@ -999,29 +998,33 @@ geigSH' a b = (l,v')
999 998
1000-- | A matrix that, by construction, it is known to be complex Hermitian or real symmetric. 999-- | A matrix that, by construction, it is known to be complex Hermitian or real symmetric.
1001-- 1000--
1002-- It can be created using 'sym', 'xTx', or 'trustSym', and the matrix can be extracted using 'her'. 1001-- It can be created using 'sym', 'mTm', or 'trustSym', and the matrix can be extracted using 'unSym'.
1003data Her t = Her (Matrix t) deriving Show 1002newtype Herm t = Herm (Matrix t) deriving Show
1003
1004instance (NFData t, Numeric t) => NFData (Herm t)
1005 where
1006 rnf (Herm m) = rnf m
1004 1007
1005-- | Extract the general matrix from a 'Her' structure, forgetting its symmetric or Hermitian property. 1008-- | Extract the general matrix from a 'Herm' structure, forgetting its symmetric or Hermitian property.
1006her :: Her t -> Matrix t 1009unSym :: Herm t -> Matrix t
1007her (Her x) = x 1010unSym (Herm x) = x
1008 1011
1009-- | Compute the complex Hermitian or real symmetric part of a square matrix (@(x + tr x)/2@). 1012-- | Compute the complex Hermitian or real symmetric part of a square matrix (@(x + tr x)/2@).
1010sym :: Field t => Matrix t -> Her t 1013sym :: Field t => Matrix t -> Herm t
1011sym x = Her (scale 0.5 (tr x `add` x)) 1014sym x = Herm (scale 0.5 (tr x `add` x))
1012 1015
1013-- | Compute the contraction @tr x <> x@ of a general matrix. 1016-- | Compute the contraction @tr x <> x@ of a general matrix.
1014xTx :: Numeric t => Matrix t -> Her t 1017mTm :: Numeric t => Matrix t -> Herm t
1015xTx x = Her (tr x `mXm` x) 1018mTm x = Herm (tr x `mXm` x)
1016 1019
1017instance Field t => Linear t Her where 1020instance Field t => Linear t Herm where
1018 scale x (Her m) = Her (scale x m) 1021 scale x (Herm m) = Herm (scale x m)
1019 1022
1020instance Field t => Additive (Her t) where 1023instance Field t => Additive (Herm t) where
1021 add (Her a) (Her b) = Her (a `add` b) 1024 add (Herm a) (Herm b) = Herm (a `add` b)
1022 1025
1023-- | At your own risk, declare that a matrix is complex Hermitian or real symmetric 1026-- | At your own risk, declare that a matrix is complex Hermitian or real symmetric
1024-- for usage in 'chol', 'eigSH', etc. Only a triangular part of the matrix will be used. 1027-- for usage in 'chol', 'eigSH', etc. Only a triangular part of the matrix will be used.
1025trustSym :: Matrix t -> Her t 1028trustSym :: Matrix t -> Herm t
1026trustSym x = (Her x) 1029trustSym x = (Herm x)
1027 1030