diff options
author | Alberto Ruiz <aruiz@um.es> | 2015-06-17 19:35:31 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2015-06-17 19:35:31 +0200 |
commit | 52009006791ee2b71530a61f4bf9e1c065c04eae (patch) | |
tree | 36c4256822d99a3abc34902a8e86150be2a0ea17 /packages/base/src/Internal/Util.hs | |
parent | 61d90ff66af8bfe53ef8cdda8dfe1e70463c213c (diff) |
improved luSolve', tests
Diffstat (limited to 'packages/base/src/Internal/Util.hs')
-rw-r--r-- | packages/base/src/Internal/Util.hs | 87 |
1 files changed, 80 insertions, 7 deletions
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index d9777ae..079663d 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs | |||
@@ -54,7 +54,7 @@ module Internal.Util( | |||
54 | -- ** 2D | 54 | -- ** 2D |
55 | corr2, conv2, separable, | 55 | corr2, conv2, separable, |
56 | block2x2,block3x3,view1,unView1,foldMatrix, | 56 | block2x2,block3x3,view1,unView1,foldMatrix, |
57 | gaussElim_1, gaussElim_2, gaussElim, luST, luSolve' | 57 | gaussElim_1, gaussElim_2, gaussElim, luST, luSolve', luSolve'', luPacked' |
58 | ) where | 58 | ) where |
59 | 59 | ||
60 | import Internal.Vector | 60 | import Internal.Vector |
@@ -64,7 +64,7 @@ import Internal.Element | |||
64 | import Internal.Container | 64 | import Internal.Container |
65 | import Internal.Vectorized | 65 | import Internal.Vectorized |
66 | import Internal.IO | 66 | import Internal.IO |
67 | import Internal.Algorithms hiding (Normed,linearSolve',luSolve') | 67 | import Internal.Algorithms hiding (Normed,linearSolve',luSolve', luPacked') |
68 | import Numeric.Matrix() | 68 | import Numeric.Matrix() |
69 | import Numeric.Vector() | 69 | import Numeric.Vector() |
70 | import Internal.Random | 70 | import Internal.Random |
@@ -686,6 +686,35 @@ luST ok (r,_) x = do | |||
686 | v <- unsafeFreezeVector p | 686 | v <- unsafeFreezeVector p |
687 | return (toList v) | 687 | return (toList v) |
688 | 688 | ||
689 | {- | Experimental implementation of 'luPacked' | ||
690 | for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'. | ||
691 | |||
692 | >>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17) | ||
693 | (5><5) | ||
694 | [ 1, 1, 2, 3, 4 | ||
695 | , 5, 7, 7, 8, 9 | ||
696 | , 10, 11, 13, 13, 14 | ||
697 | , 15, 16, 0, 2, 2 | ||
698 | , 3, 4, 5, 6, 8 ] | ||
699 | |||
700 | >>> let (l,u,p,s) = luFact $ luPacked' m | ||
701 | >>> l | ||
702 | (5><5) | ||
703 | [ 1, 0, 0, 0, 0 | ||
704 | , 6, 1, 0, 0, 0 | ||
705 | , 12, 7, 1, 0, 0 | ||
706 | , 7, 10, 7, 1, 0 | ||
707 | , 8, 2, 6, 11, 1 ] | ||
708 | >>> u | ||
709 | (5><5) | ||
710 | [ 15, 16, 0, 2, 2 | ||
711 | , 0, 13, 7, 13, 14 | ||
712 | , 0, 0, 15, 0, 11 | ||
713 | , 0, 0, 0, 15, 15 | ||
714 | , 0, 0, 0, 0, 1 ] | ||
715 | |||
716 | -} | ||
717 | luPacked' x = mutable (luST (magnit 0)) x | ||
689 | 718 | ||
690 | -------------------------------------------------------------------------------- | 719 | -------------------------------------------------------------------------------- |
691 | 720 | ||
@@ -693,35 +722,79 @@ rowRange m = [0..rows m -1] | |||
693 | 722 | ||
694 | at k = Pos (idxs[k]) | 723 | at k = Pos (idxs[k]) |
695 | 724 | ||
696 | backSust lup rhs = foldl' f (rhs?[]) (reverse ls) | 725 | backSust' lup rhs = foldl' f (rhs?[]) (reverse ls) |
697 | where | 726 | where |
698 | ls = [ (d k , u k , b k) | k <- rowRange lup ] | 727 | ls = [ (d k , u k , b k) | k <- rowRange lup ] |
699 | where | 728 | where |
700 | d k = lup ?? (at k, at k) | 729 | d k = lup ?? (at k, at k) |
701 | u k = lup ?? (at k, Drop (k+1)) | 730 | u k = lup ?? (at k, Drop (k+1)) |
702 | b k = rhs ?? (at k, All) | 731 | b k = rhs ?? (at k, All) |
703 | 732 | ||
704 | f x (d,u,b) = (b - u<>x) / d | 733 | f x (d,u,b) = (b - u<>x) / d |
705 | === | 734 | === |
706 | x | 735 | x |
707 | 736 | ||
708 | 737 | ||
709 | forwSust lup rhs = foldl' f (rhs?[]) ls | 738 | forwSust' lup rhs = foldl' f (rhs?[]) ls |
710 | where | 739 | where |
711 | ls = [ (l k , b k) | k <- rowRange lup ] | 740 | ls = [ (l k , b k) | k <- rowRange lup ] |
712 | where | 741 | where |
713 | l k = lup ?? (at k, Take k) | 742 | l k = lup ?? (at k, Take k) |
714 | b k = rhs ?? (at k, All) | 743 | b k = rhs ?? (at k, All) |
715 | 744 | ||
716 | f x (l,b) = x | 745 | f x (l,b) = x |
717 | === | 746 | === |
718 | (b - l<>x) | 747 | (b - l<>x) |
719 | 748 | ||
720 | 749 | ||
721 | luSolve' (lup,p) b = backSust lup (forwSust lup pb) | 750 | luSolve'' (lup,p) b = backSust' lup (forwSust' lup pb) |
722 | where | 751 | where |
723 | pb = b ?? (Pos (fixPerm' p), All) | 752 | pb = b ?? (Pos (fixPerm' p), All) |
724 | 753 | ||
754 | -------------------------------------------------------------------------------- | ||
755 | |||
756 | forwSust lup rhs = fst $ mutable f rhs | ||
757 | where | ||
758 | f (r,c) x = do | ||
759 | l <- unsafeThawMatrix lup | ||
760 | let go k = gemmm 1 (Slice x k 0 1 c) (-1) (Slice l k 0 1 k) (Slice x 0 0 k c) | ||
761 | mapM_ go [0..r-1] | ||
762 | |||
763 | |||
764 | backSust lup rhs = fst $ mutable f rhs | ||
765 | where | ||
766 | f (r,c) m = do | ||
767 | l <- unsafeThawMatrix lup | ||
768 | let d k = recip (lup `atIndex` (k,k)) | ||
769 | u k = Slice l k (k+1) 1 (r-1-k) | ||
770 | b k = Slice m k 0 1 c | ||
771 | x k = Slice m (k+1) 0 (r-1-k) c | ||
772 | scal k = rowOper (SCAL (d k) (Row k) AllCols) m | ||
773 | |||
774 | go k = gemmm 1 (b k) (-1) (u k) (x k) >> scal k | ||
775 | mapM_ go [r-1,r-2..0] | ||
776 | |||
777 | |||
778 | {- | Experimental implementation of 'luSolve' for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'. | ||
779 | |||
780 | >>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13) | ||
781 | (2><2) | ||
782 | [ 1, 2 | ||
783 | , 3, 5 ] | ||
784 | >>> b | ||
785 | (2><3) | ||
786 | [ 5, 1, 3 | ||
787 | , 8, 6, 3 ] | ||
788 | |||
789 | >>> luSolve' (luPacked' a) b | ||
790 | (2><3) | ||
791 | [ 4, 7, 4 | ||
792 | , 7, 10, 6 ] | ||
793 | |||
794 | -} | ||
795 | luSolve' (lup,p) b = backSust lup (forwSust lup pb) | ||
796 | where | ||
797 | pb = b ?? (Pos (fixPerm' p), All) | ||
725 | 798 | ||
726 | -------------------------------------------------------------------------------- | 799 | -------------------------------------------------------------------------------- |
727 | 800 | ||