diff options
author | Alberto Ruiz <aruiz@um.es> | 2015-06-08 21:50:04 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2015-06-08 21:50:04 +0200 |
commit | 660e1eb875a3517c85285a9341e9e07e647e4e10 (patch) | |
tree | 38f28421bc2327435bb3413af5c0d1b1526ac1c0 /packages/base/src/Internal | |
parent | 7c446e77a52dc0a979c0fe570f1e7b6127d9ff42 (diff) |
gaussElim using foldMatrix
Diffstat (limited to 'packages/base/src/Internal')
-rw-r--r-- | packages/base/src/Internal/Util.hs | 84 |
1 files changed, 49 insertions, 35 deletions
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index e9cf155..1b639e1 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs | |||
@@ -53,6 +53,7 @@ module Internal.Util( | |||
53 | corr, conv, corrMin, | 53 | corr, conv, corrMin, |
54 | -- ** 2D | 54 | -- ** 2D |
55 | corr2, conv2, separable, | 55 | corr2, conv2, separable, |
56 | block2x2,block3x3,view1,unView1,foldMatrix, | ||
56 | gaussElim | 57 | gaussElim |
57 | ) where | 58 | ) where |
58 | 59 | ||
@@ -71,8 +72,7 @@ import Internal.Convolution | |||
71 | import Control.Monad(when) | 72 | import Control.Monad(when) |
72 | import Text.Printf | 73 | import Text.Printf |
73 | import Data.List.Split(splitOn) | 74 | import Data.List.Split(splitOn) |
74 | import Data.List(intercalate,sortBy) | 75 | import Data.List(intercalate,) |
75 | import Data.Function(on) | ||
76 | import Control.Arrow((&&&)) | 76 | import Control.Arrow((&&&)) |
77 | import Data.Complex | 77 | import Data.Complex |
78 | 78 | ||
@@ -498,50 +498,64 @@ dispShort maxr maxc dec m = | |||
498 | 498 | ||
499 | -------------------------------------------------------------------------------- | 499 | -------------------------------------------------------------------------------- |
500 | 500 | ||
501 | -- | generic reference implementation of gaussian elimination | 501 | -- matrix views |
502 | -- | ||
503 | -- @a <> gauss a b = b@ | ||
504 | -- | ||
505 | gaussElim | ||
506 | :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t) | ||
507 | => Matrix t -> Matrix t -> Matrix t | ||
508 | 502 | ||
509 | gaussElim x y = dropColumns (rows x) (flipud $ fromRows s2) | 503 | block2x2 r c m = [[m11,m12],[m21,m22]] |
510 | where | 504 | where |
511 | rs = toRows $ fromBlocks [[x , y]] | 505 | m11 = m ?? (Take r, Take c) |
512 | s1 = pivotDown (rows x) 0 rs | 506 | m12 = m ?? (Take r, Drop c) |
513 | s2 = pivotUp (rows x-1) (reverse s1) | 507 | m21 = m ?? (Drop r, Take c) |
508 | m22 = m ?? (Drop r, Drop c) | ||
514 | 509 | ||
515 | pivotDown t n xs | 510 | block3x3 r nr c nc m = [[m ?? (er !! i, ec !! j) | j <- [0..2] ] | i <- [0..2] ] |
516 | | t == n = [] | ||
517 | | otherwise = y : pivotDown t (n+1) ys | ||
518 | where | 511 | where |
519 | y:ys = redu (pivot n xs) | 512 | er = [ Range 0 1 (r-1), Range r 1 (r+nr-1), Drop (nr+r) ] |
513 | ec = [ Range 0 1 (c-1), Range c 1 (c+nc-1), Drop (nc+c) ] | ||
514 | |||
515 | view1 :: Numeric t => Matrix t -> Maybe (t, Vector t, Vector t, Matrix t) | ||
516 | view1 m | ||
517 | | rows m > 0 && cols m > 0 = Just (e, flatten m12, flatten m21 , m22) | ||
518 | | otherwise = Nothing | ||
519 | where | ||
520 | [[m11,m12],[m21,m22]] = block2x2 1 1 m | ||
521 | e = m11 `atIndex` (0, 0) | ||
522 | |||
523 | unView1 :: Numeric t => (t, Vector t, Vector t, Matrix t) -> Matrix t | ||
524 | unView1 (e,r,c,m) = fromBlocks [[scalar e, asRow r],[asColumn c, m]] | ||
525 | |||
526 | |||
527 | foldMatrix g f ( (f <$>) . view1 . g -> Just (e,r,c,m)) = unView1 (e, r, c, foldMatrix g f m) | ||
528 | foldMatrix _ _ m = m | ||
520 | 529 | ||
521 | pivot k = (const k &&& id) | 530 | sortRowsBy h j m = m ?? (Pos (sortIndex (h (tr m ! j))), All) |
522 | . reverse . sortBy (compare `on` (abs. (!k))) -- FIXME | ||
523 | 531 | ||
524 | redu (k,x:zs) | 532 | splitColsAt n = (takeColumns n &&& dropColumns n) |
525 | | p == 0 = error "gauss: singular!" -- FIXME | 533 | |
526 | | otherwise = u : map f zs | 534 | |
535 | down a = foldMatrix g f a | ||
536 | where | ||
537 | g = sortRowsBy (negate.abs) 0 | ||
538 | f (e,r,c,m) | ||
539 | | e /= 0 = (1, r', 0, m - outer c r') | ||
540 | | otherwise = error "singular!" | ||
527 | where | 541 | where |
528 | p = x!k | 542 | r' = r / scalar e |
529 | u = scale (recip (x!k)) x | ||
530 | f z = z - scale (z!k) u | ||
531 | redu (_,[]) = [] | ||
532 | 543 | ||
533 | 544 | ||
534 | pivotUp n xs | 545 | -- | generic reference implementation of gaussian elimination |
535 | | n == -1 = [] | 546 | -- |
536 | | otherwise = y : pivotUp (n-1) ys | 547 | -- @a <> gaussElim a b = b@ |
548 | -- | ||
549 | gaussElim | ||
550 | :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t) | ||
551 | => Matrix t -> Matrix t -> Matrix t | ||
552 | |||
553 | gaussElim a b = r | ||
537 | where | 554 | where |
538 | y:ys = redu' (n,xs) | 555 | go x y = splitColsAt (cols a) (down $ fromBlocks [[x,y]]) |
556 | (a1,b1) = go a b | ||
557 | ( _, r) = go (flipud . fliprl $ a1) (flipud . fliprl $ b1) | ||
539 | 558 | ||
540 | redu' (k,x:zs) = u : map f zs | ||
541 | where | ||
542 | u = x | ||
543 | f z = z - scale (z!k) u | ||
544 | redu' (_,[]) = [] | ||
545 | 559 | ||
546 | -------------------------------------------------------------------------------- | 560 | -------------------------------------------------------------------------------- |
547 | 561 | ||