diff options
Diffstat (limited to 'packages/base/src/Internal')
-rw-r--r-- | packages/base/src/Internal/Algorithms.hs | 71 | ||||
-rw-r--r-- | packages/base/src/Internal/Util.hs | 2 |
2 files changed, 38 insertions, 35 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 | {- | |
15 | Module : Internal.Algorithms | 9 | Module : 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'. |
379 | ldlPacked :: Field t => Her t -> LDL t | 373 | ldlPacked :: Field t => Herm t -> LDL t |
380 | ldlPacked (Her m) = ldlPackedSH m | 374 | ldlPacked (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] | |||
461 | 3.000 5.000 6.000 | 455 | 3.000 5.000 6.000 |
462 | 456 | ||
463 | -} | 457 | -} |
464 | eigSH :: Field t => Her t -> (Vector Double, Matrix t) | 458 | eigSH :: Field t => Herm t -> (Vector Double, Matrix t) |
465 | eigSH (Her m) = eigSH' m | 459 | eigSH (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. |
468 | eigenvaluesSH :: Field t => Her t -> Vector Double | 462 | eigenvaluesSH :: Field t => Herm t -> Vector Double |
469 | eigenvaluesSH (Her m) = eigenvaluesSH' m | 463 | eigenvaluesSH (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.) |
474 | data QR t = QR (Matrix t) (Vector t) | 468 | data QR t = QR (Matrix t) (Vector t) |
475 | 469 | ||
470 | instance (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@. |
536 | chol :: Field t => Her t -> Matrix t | 535 | chol :: Field t => Herm t -> Matrix t |
537 | chol (Her m) = {-# SCC "chol" #-} cholSH' m | 536 | chol (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'. |
540 | mbChol :: Field t => Her t -> Maybe (Matrix t) | 539 | mbChol :: Field t => Herm t -> Maybe (Matrix t) |
541 | mbChol (Her m) = {-# SCC "mbChol" #-} mbCholSH' m | 540 | mbChol (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. |
979 | geigSH :: Field t | 978 | geigSH :: 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) |
983 | geigSH (Her a) (Her b) = geigSH' a b | 982 | geigSH (Herm a) (Herm b) = geigSH' a b |
984 | 983 | ||
985 | geigSH' :: Field t | 984 | geigSH' :: 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'. |
1003 | data Her t = Her (Matrix t) deriving Show | 1002 | newtype Herm t = Herm (Matrix t) deriving Show |
1003 | |||
1004 | instance (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. |
1006 | her :: Her t -> Matrix t | 1009 | unSym :: Herm t -> Matrix t |
1007 | her (Her x) = x | 1010 | unSym (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@). |
1010 | sym :: Field t => Matrix t -> Her t | 1013 | sym :: Field t => Matrix t -> Herm t |
1011 | sym x = Her (scale 0.5 (tr x `add` x)) | 1014 | sym 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. |
1014 | xTx :: Numeric t => Matrix t -> Her t | 1017 | mTm :: Numeric t => Matrix t -> Herm t |
1015 | xTx x = Her (tr x `mXm` x) | 1018 | mTm x = Herm (tr x `mXm` x) |
1016 | 1019 | ||
1017 | instance Field t => Linear t Her where | 1020 | instance 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 | ||
1020 | instance Field t => Additive (Her t) where | 1023 | instance 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. |
1025 | trustSym :: Matrix t -> Her t | 1028 | trustSym :: Matrix t -> Herm t |
1026 | trustSym x = (Her x) | 1029 | trustSym x = (Herm x) |
1027 | 1030 | ||
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index 36b7855..cf42961 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs | |||
@@ -462,7 +462,7 @@ null1 :: Matrix R -> Vector R | |||
462 | null1 = last . toColumns . snd . rightSV | 462 | null1 = last . toColumns . snd . rightSV |
463 | 463 | ||
464 | -- | solution of overconstrained homogeneous symmetric linear system | 464 | -- | solution of overconstrained homogeneous symmetric linear system |
465 | null1sym :: Her R -> Vector R | 465 | null1sym :: Herm R -> Vector R |
466 | null1sym = last . toColumns . snd . eigSH | 466 | null1sym = last . toColumns . snd . eigSH |
467 | 467 | ||
468 | -------------------------------------------------------------------------------- | 468 | -------------------------------------------------------------------------------- |