X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FLinear%2FMatrix.hs;h=5bb5519664c5fa4028d284775493bb7da4f5598e;hb=221805e711475880d37fb7bb6d7a0b04689b1ebe;hp=ef0e9b6ff3ccba65c310b81126ea915d543e0411;hpb=2d0ecad8695e129443e53311d6494b2465f1a672;p=numerical-analysis.git diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs index ef0e9b6..5bb5519 100644 --- a/src/Linear/Matrix.hs +++ b/src/Linear/Matrix.hs @@ -13,23 +13,19 @@ import Data.List (intercalate) import Data.Vector.Fixed ( Dim, + N1, Vector ) import qualified Data.Vector.Fixed as V ( - Fun(..), - N1, and, - eq, - foldl, fromList, length, map, - maximum, replicate, toList, zipWith ) -import Data.Vector.Fixed.Internal (Arity, arity, S, Dim) +import Data.Vector.Fixed.Internal (Arity, arity, S) import Linear.Vector import Normed @@ -37,13 +33,9 @@ import NumericPrelude hiding ((*), abs) import qualified NumericPrelude as NP ((*)) import qualified Algebra.Algebraic as Algebraic import Algebra.Algebraic (root) -import qualified Algebra.Absolute as Absolute import qualified Algebra.Additive as Additive import qualified Algebra.Ring as Ring -import Algebra.Absolute (abs) -import qualified Algebra.Field as Field import qualified Algebra.Module as Module -import qualified Algebra.RealField as RealField import qualified Algebra.RealRing as RealRing import qualified Algebra.ToRational as ToRational import qualified Algebra.Transcendental as Transcendental @@ -247,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. -- @@ -294,7 +419,7 @@ instance (Ring.C a, => Ring.C (Mat v w a) where -- The first * is ring multiplication, the second is matrix -- multiplication. - m1 * m2 = m1 * m1 + m1 * m2 = m1 * m2 instance (Ring.C a, @@ -359,14 +484,17 @@ instance (Algebraic.C a, -- -- >>> import Roots.Simple -- >>> let h = 0.5 :: Double --- >>> let g1 (Mat (D2 (D1 x) (D1 y))) = 1.0 + h*exp(-(x^2))/(1.0 + y^2) --- >>> let g2 (Mat (D2 (D1 x) (D1 y))) = 0.5 + h*atan(x^2 + y^2) +-- >>> let g1 (Mat (D2 (D1 x) (D1 y))) = 1.0 + h NP.* exp(-(x^2))/(1.0 + y^2) +-- >>> let g2 (Mat (D2 (D1 x) (D1 y))) = 0.5 + h NP.* atan(x^2 + y^2) -- >>> let g u = vec2d ((g1 u), (g2 u)) -- >>> let u0 = vec2d (1.0, 1.0) -- >>> let eps = 1/(10^9) -- >>> fixed_point g eps u0 --- (1.0728549599342185,1.0820591495686167) +-- ((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)) @@ -382,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, @@ -405,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],