summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
Diffstat (limited to 'packages')
-rw-r--r--packages/base/src/Internal/Algorithms.hs33
-rw-r--r--packages/base/src/Internal/Util.hs43
-rw-r--r--packages/base/src/Numeric/LinearAlgebra.hs3
3 files changed, 61 insertions, 18 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
29import Internal.LAPACK as LAPACK 29import Internal.LAPACK as LAPACK
30import Internal.Numeric 30import Internal.Numeric
31import Data.List(foldl1') 31import Data.List(foldl1')
32import Data.Array 32import qualified Data.Array as A
33import Internal.ST
34import 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
578peps :: RealFloat x => x 580peps :: RealFloat x => x
579peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x) 581peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x)
580 582
581
582-- | The imaginary unit: @i = 0.0 :+ 1.0@
583i :: Complex Double
584i = 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
799swap (arr,s) (a,b) | a /= b = (arr // [(a, arr!b),(b,arr!a)],-s) 796fixPerm r vals = (fromColumns $ A.elems res, sign)
800 | otherwise = (arr,s) 797 where
801 798 v = [0..r-1]
802fixPerm 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
805fixPerm' :: [Int] -> Vector I
806fixPerm' 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
807triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]] 814triang 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
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index f08f710..d9777ae 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -1,6 +1,5 @@
1{-# LANGUAGE FlexibleContexts #-} 1{-# LANGUAGE FlexibleContexts #-}
2{-# LANGUAGE FlexibleInstances #-} 2{-# LANGUAGE FlexibleInstances #-}
3{-# LANGUAGE TypeFamilies #-}
4{-# LANGUAGE MultiParamTypeClasses #-} 3{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-} 4{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE ViewPatterns #-} 5{-# LANGUAGE ViewPatterns #-}
@@ -55,7 +54,7 @@ module Internal.Util(
55 -- ** 2D 54 -- ** 2D
56 corr2, conv2, separable, 55 corr2, conv2, separable,
57 block2x2,block3x3,view1,unView1,foldMatrix, 56 block2x2,block3x3,view1,unView1,foldMatrix,
58 gaussElim_1, gaussElim_2, gaussElim, luST 57 gaussElim_1, gaussElim_2, gaussElim, luST, luSolve'
59) where 58) where
60 59
61import Internal.Vector 60import Internal.Vector
@@ -65,7 +64,7 @@ import Internal.Element
65import Internal.Container 64import Internal.Container
66import Internal.Vectorized 65import Internal.Vectorized
67import Internal.IO 66import Internal.IO
68import Internal.Algorithms hiding (i,Normed,swap,linearSolve') 67import Internal.Algorithms hiding (Normed,linearSolve',luSolve')
69import Numeric.Matrix() 68import Numeric.Matrix()
70import Numeric.Vector() 69import Numeric.Vector()
71import Internal.Random 70import Internal.Random
@@ -73,7 +72,7 @@ import Internal.Convolution
73import Control.Monad(when,forM_) 72import Control.Monad(when,forM_)
74import Text.Printf 73import Text.Printf
75import Data.List.Split(splitOn) 74import Data.List.Split(splitOn)
76import Data.List(intercalate,sortBy) 75import Data.List(intercalate,sortBy,foldl')
77import Control.Arrow((&&&)) 76import Control.Arrow((&&&))
78import Data.Complex 77import Data.Complex
79import Data.Function(on) 78import Data.Function(on)
@@ -690,6 +689,42 @@ luST ok (r,_) x = do
690 689
691-------------------------------------------------------------------------------- 690--------------------------------------------------------------------------------
692 691
692rowRange m = [0..rows m -1]
693
694at k = Pos (idxs[k])
695
696backSust lup rhs = foldl' f (rhs?[]) (reverse ls)
697 where
698 ls = [ (d k , u k , b k) | k <- rowRange lup ]
699 where
700 d k = lup ?? (at k, at k)
701 u k = lup ?? (at k, Drop (k+1))
702 b k = rhs ?? (at k, All)
703
704 f x (d,u,b) = (b - u<>x) / d
705 ===
706 x
707
708
709forwSust lup rhs = foldl' f (rhs?[]) ls
710 where
711 ls = [ (l k , b k) | k <- rowRange lup ]
712 where
713 l k = lup ?? (at k, Take k)
714 b k = rhs ?? (at k, All)
715
716 f x (l,b) = x
717 ===
718 (b - l<>x)
719
720
721luSolve' (lup,p) b = backSust lup (forwSust lup pb)
722 where
723 pb = b ?? (Pos (fixPerm' p), All)
724
725
726--------------------------------------------------------------------------------
727
693instance Testable (Matrix I) where 728instance Testable (Matrix I) where
694 checkT _ = test 729 checkT _ = test
695 730
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 2e6b8ca..e899445 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -77,6 +77,7 @@ module Numeric.LinearAlgebra (
77 linearSolveLS, 77 linearSolveLS,
78 linearSolveSVD, 78 linearSolveSVD,
79 luSolve, 79 luSolve,
80 luSolve',
80 cholSolve, 81 cholSolve,
81 cgSolve, 82 cgSolve,
82 cgSolve', 83 cgSolve',
@@ -158,7 +159,7 @@ import Numeric.Vector()
158import Internal.Matrix 159import Internal.Matrix
159import Internal.Container hiding ((<>)) 160import Internal.Container hiding ((<>))
160import Internal.Numeric hiding (mul) 161import Internal.Numeric hiding (mul)
161import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked',linearSolve') 162import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked',linearSolve',luSolve')
162import qualified Internal.Algorithms as A 163import qualified Internal.Algorithms as A
163import Internal.Util 164import Internal.Util
164import Internal.Random 165import Internal.Random