diff options
Diffstat (limited to 'packages/base/src/Internal/Algorithms.hs')
-rw-r--r-- | packages/base/src/Internal/Algorithms.hs | 33 |
1 files changed, 20 insertions, 13 deletions
diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs index aaf6fbb..1235da3 100644 --- a/packages/base/src/Internal/Algorithms.hs +++ b/packages/base/src/Internal/Algorithms.hs | |||
@@ -29,7 +29,9 @@ import Internal.Conversion | |||
29 | import Internal.LAPACK as LAPACK | 29 | import Internal.LAPACK as LAPACK |
30 | import Internal.Numeric | 30 | import Internal.Numeric |
31 | import Data.List(foldl1') | 31 | import Data.List(foldl1') |
32 | import Data.Array | 32 | import qualified Data.Array as A |
33 | import Internal.ST | ||
34 | import Internal.Vectorized(range) | ||
33 | 35 | ||
34 | {- | Generic linear algebra functions for double precision real and complex matrices. | 36 | {- | Generic linear algebra functions for double precision real and complex matrices. |
35 | 37 | ||
@@ -578,11 +580,6 @@ eps = 2.22044604925031e-16 | |||
578 | peps :: RealFloat x => x | 580 | peps :: RealFloat x => x |
579 | peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x) | 581 | peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x) |
580 | 582 | ||
581 | |||
582 | -- | The imaginary unit: @i = 0.0 :+ 1.0@ | ||
583 | i :: Complex Double | ||
584 | i = 0:+1 | ||
585 | |||
586 | ----------------------------------------------------------------------- | 583 | ----------------------------------------------------------------------- |
587 | 584 | ||
588 | -- | The nullspace of a matrix from its precomputed SVD decomposition. | 585 | -- | The nullspace of a matrix from its precomputed SVD decomposition. |
@@ -796,13 +793,23 @@ signlp r vals = foldl f 1 (zip [0..r-1] vals) | |||
796 | where f s (a,b) | a /= b = -s | 793 | where f s (a,b) | a /= b = -s |
797 | | otherwise = s | 794 | | otherwise = s |
798 | 795 | ||
799 | swap (arr,s) (a,b) | a /= b = (arr // [(a, arr!b),(b,arr!a)],-s) | 796 | fixPerm r vals = (fromColumns $ A.elems res, sign) |
800 | | otherwise = (arr,s) | 797 | where |
801 | 798 | v = [0..r-1] | |
802 | fixPerm r vals = (fromColumns $ elems res, sign) | 799 | t = toColumns (ident r) |
803 | where v = [0..r-1] | 800 | (res,sign) = foldl swap (A.listArray (0,r-1) t, 1) (zip v vals) |
804 | s = toColumns (ident r) | 801 | swap (arr,s) (a,b) |
805 | (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals) | 802 | | a /= b = (arr A.// [(a, arr A.! b),(b,arr A.! a)],-s) |
803 | | otherwise = (arr,s) | ||
804 | |||
805 | fixPerm' :: [Int] -> Vector I | ||
806 | fixPerm' s = res $ mutable f s0 | ||
807 | where | ||
808 | s0 = reshape 1 (range (length s)) | ||
809 | res = flatten . fst | ||
810 | swap m i j = rowOper (SWAP i j AllCols) m | ||
811 | f :: (Num t, Element t) => (Int, Int) -> STMatrix s t -> ST s () -- needed because of TypeFamilies | ||
812 | f _ p = sequence_ $ zipWith (swap p) [0..] s | ||
806 | 813 | ||
807 | triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]] | 814 | triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]] |
808 | where el p q = if q-p>=h then v else 1 - v | 815 | where el p q = if q-p>=h then v else 1 - v |