From 221805e711475880d37fb7bb6d7a0b04689b1ebe Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sat, 23 Feb 2013 19:05:26 -0500 Subject: [PATCH] Add triangular functions, determinant, minor (all half-baked) to Linear.Matrix. --- src/Linear/Matrix.hs | 142 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 139 insertions(+), 3 deletions(-) diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs index 4f2e845..5bb5519 100644 --- a/src/Linear/Matrix.hs +++ b/src/Linear/Matrix.hs @@ -13,10 +13,10 @@ import Data.List (intercalate) import Data.Vector.Fixed ( Dim, + N1, Vector ) import qualified Data.Vector.Fixed as V ( - N1, and, fromList, length, @@ -239,6 +239,139 @@ cholesky m = construct r | otherwise = 0 +-- | Returns True if the given matrix is upper-triangular, and False +-- otherwise. +-- +-- Examples: +-- +-- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int +-- >>> is_upper_triangular m +-- False +-- +-- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int +-- >>> is_upper_triangular m +-- True +-- +is_upper_triangular :: (Eq a, Ring.C a, Vector w a) => Mat v w a -> Bool +is_upper_triangular m = + and $ concat results + where + results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ] + + test :: Int -> Int -> Bool + test i j + | i <= j = True + | otherwise = m !!! (i,j) == 0 + + +-- | Returns True if the given matrix is lower-triangular, and False +-- otherwise. +-- +-- Examples: +-- +-- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int +-- >>> is_lower_triangular m +-- True +-- +-- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int +-- >>> is_lower_triangular m +-- False +-- +is_lower_triangular :: (Eq a, + Ring.C a, + Vector w a, + Vector w (v a), + Vector v a) + => Mat v w a + -> Bool +is_lower_triangular = is_upper_triangular . transpose + + +-- | Returns True if the given matrix is triangular, and False +-- otherwise. +-- +-- Examples: +-- +-- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int +-- >>> is_triangular m +-- True +-- +-- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int +-- >>> is_triangular m +-- True +-- +-- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int +-- >>> is_triangular m +-- False +-- +is_triangular :: (Eq a, + Ring.C a, + Vector w a, + Vector w (v a), + Vector v a) + => Mat v w a + -> Bool +is_triangular m = is_upper_triangular m || is_lower_triangular m + + +-- | Return the (i,j)th minor of m. +-- +-- Examples: +-- +-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int +-- >>> minor m 0 0 :: Mat2 Int +-- ((5,6),(8,9)) +-- >>> minor m 1 1 :: Mat2 Int +-- ((1,3),(7,9)) +-- +minor :: (Dim v ~ S (Dim u), + Dim w ~ S (Dim z), + Vector z a, + Vector u (w a), + Vector u (z a)) + => Mat v w a + -> Int + -> Int + -> Mat u z a +minor (Mat rows) i j = m + where + rows' = delete rows i + m = Mat $ V.map ((flip delete) j) rows' + + +determinant :: (Eq a, + Ring.C a, + Vector w a, + Vector w (v a), + Vector v a, + Dim v ~ S r, + Dim w ~ S t) + => Mat v w a + -> a +determinant m + | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ] + | otherwise = undefined --determinant_recursive m + +{- +determinant_recursive :: forall v w a r c. + (Eq a, + Ring.C a, + Vector w a) + => Mat (v r) (w c) a + -> a +determinant_recursive m + | (ncols m) == 0 || (nrows m) == 0 = error "don't do that" + | (ncols m) == 1 && (nrows m) == 1 = m !!! (0,0) -- Base case + | otherwise = + sum [ (-1)^(1+(toInteger j)) NP.* (m' 1 j) NP.* (det_minor 1 j) + | j <- [0..(ncols m)-1] ] + where + m' i j = m !!! (i,j) + + det_minor :: Int -> Int -> a + det_minor i j = determinant (minor m i j) +-} + -- | Matrix multiplication. Our 'Num' instance doesn't define one, and -- we need additional restrictions on the result type anyway. -- @@ -359,6 +492,9 @@ instance (Algebraic.C a, -- >>> fixed_point g eps u0 -- ((1.0728549599342185),(1.0820591495686167)) -- +vec1d :: (a) -> Mat D1 D1 a +vec1d (x) = Mat (D1 (D1 x)) + vec2d :: (a,a) -> Mat D2 D1 a vec2d (x,y) = Mat (D2 (D1 x) (D1 y)) @@ -374,7 +510,7 @@ scalar :: a -> Mat D1 D1 a scalar x = Mat (D1 (D1 x)) dot :: (RealRing.C a, - Dim w ~ V.N1, + Dim w ~ N1, Dim v ~ S n, Vector v a, Vector w a, @@ -397,7 +533,7 @@ v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0) -- angle :: (Transcendental.C a, RealRing.C a, - Dim w ~ V.N1, + Dim w ~ N1, Dim v ~ S n, Vector w (w a), Vector v [a], -- 2.44.2