src/Linear/Matrix.hs

   1 {-# LANGUAGE ExistentialQuantification #-}
   2 {-# LANGUAGE FlexibleContexts #-}
   3 {-# LANGUAGE FlexibleInstances #-}
   4 {-# LANGUAGE MultiParamTypeClasses #-}
   5 {-# LANGUAGE NoMonomorphismRestriction #-}
   6 {-# LANGUAGE ScopedTypeVariables #-}
   7 {-# LANGUAGE TypeFamilies #-}
   8 {-# LANGUAGE RebindableSyntax #-}
   9
  10 -- | Boxed matrices; that is, boxed m-vectors of boxed n-vectors. We
  11 --   assume that the underlying representation is
  12 --   Data.Vector.Fixed.Boxed.Vec for simplicity. It was tried in
  13 --   generality and failed.
  14 --
  15 module Linear.Matrix
  16 where
  17
  18 import Data.List (intercalate)
  19
  20 import Data.Vector.Fixed (
  21   (!),
  22   N1,
  23   N2,
  24   N3,
  25   N4,
  26   N5,
  27   S,
  28   Z,
  29   generate,
  30   mk1,
  31   mk2,
  32   mk3,
  33   mk4,
  34   mk5
  35   )
  36 import qualified Data.Vector.Fixed as V (
  37   and,
  38   fromList,
  39   head,
  40   length,
  41   map,
  42   maximum,
  43   replicate,
  44   toList,
  45   zipWith )
  46 import Data.Vector.Fixed.Cont ( Arity, arity )
  47 import Linear.Vector ( Vec, delete, element_sum )
  48 import Normed ( Normed(..) )
  49
  50 import NumericPrelude hiding ( (*), abs )
  51 import qualified NumericPrelude as NP ( (*) )
  52 import qualified Algebra.Absolute as Absolute ( C )
  53 import Algebra.Absolute ( abs )
  54 import qualified Algebra.Additive as Additive ( C )
  55 import qualified Algebra.Algebraic as Algebraic ( C )
  56 import Algebra.Algebraic ( root )
  57 import qualified Algebra.Ring as Ring ( C )
  58 import qualified Algebra.Module as Module ( C )
  59 import qualified Algebra.RealRing as RealRing ( C )
  60 import qualified Algebra.ToRational as ToRational ( C )
  61 import qualified Algebra.Transcendental as Transcendental ( C )
  62 import qualified Prelude as P ( map )
  63
  64 data Mat m n a = (Arity m, Arity n) => Mat (Vec m (Vec n a))
  65 type Mat1 a = Mat N1 N1 a
  66 type Mat2 a = Mat N2 N2 a
  67 type Mat3 a = Mat N3 N3 a
  68 type Mat4 a = Mat N4 N4 a
  69 type Mat5 a = Mat N5 N5 a
  70
  71 instance (Eq a) => Eq (Mat m n a) where
  72   -- | Compare a row at a time.
  73   --
  74   --   Examples:
  75   --
  76   --   >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
  77   --   >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int
  78   --   >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int
  79   --   >>> m1 == m2
  80   --   True
  81   --   >>> m1 == m3
  82   --   False
  83   --
  84   (Mat rows1) == (Mat rows2) =
  85     V.and $ V.zipWith comp rows1 rows2
  86     where
  87       -- Compare a row, one column at a time.
  88       comp row1 row2 = V.and (V.zipWith (==) row1 row2)
  89
  90
  91 instance (Show a) => Show (Mat m n a) where
  92   -- | Display matrices and vectors as ordinary tuples. This is poor
  93   --   practice, but these results are primarily displayed
  94   --   interactively and convenience trumps correctness (said the guy
  95   --   who insists his vector lengths be statically checked at
  96   --   compile-time).
  97   --
  98   --   Examples:
  99   --
 100   --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 101   --   >>> show m
 102   --   ((1,2),(3,4))
 103   --
 104   show (Mat rows) =
 105     "(" ++ (intercalate "," (V.toList row_strings)) ++ ")"
 106     where
 107       row_strings = V.map show_vector rows
 108       show_vector v1 =
 109         "(" ++ (intercalate "," element_strings) ++ ")"
 110         where
 111           v1l = V.toList v1
 112           element_strings = P.map show v1l
 113
 114
 115 -- | Convert a matrix to a nested list.
 116 toList :: Mat m n a -> [[a]]
 117 toList (Mat rows) = map V.toList (V.toList rows)
 118
 119 -- | Create a matrix from a nested list.
 120 fromList :: (Arity m, Arity n) => [[a]] -> Mat m n a
 121 fromList vs = Mat (V.fromList $ map V.fromList vs)
 122
 123
 124 -- | Unsafe indexing.
 125 (!!!) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> a
 126 (!!!) m (i, j) = (row m i) ! j
 127
 128 -- | Safe indexing.
 129 (!!?) :: Mat m n a -> (Int, Int) -> Maybe a
 130 (!!?) m@(Mat rows) (i, j)
 131   | i < 0 || j < 0 = Nothing
 132   | i > V.length rows = Nothing
 133   | otherwise = if j > V.length (row m j)
 134                 then Nothing
 135                 else Just $ (row m j) ! j
 136
 137
 138 -- | The number of rows in the matrix.
 139 nrows :: forall m n a. (Arity m) => Mat m n a -> Int
 140 nrows _ = arity (undefined :: m)
 141
 142 -- | The number of columns in the first row of the
 143 --   matrix. Implementation stolen from Data.Vector.Fixed.length.
 144 ncols :: forall m n a. (Arity n) => Mat m n a -> Int
 145 ncols _ = arity (undefined :: n)
 146
 147
 148 -- | Return the @i@th row of @m@. Unsafe.
 149 row :: Mat m n a -> Int -> (Vec n a)
 150 row (Mat rows) i = rows ! i
 151
 152
 153 -- | Return the @j@th column of @m@. Unsafe.
 154 column :: Mat m n a -> Int -> (Vec m a)
 155 column (Mat rows) j =
 156   V.map (element j) rows
 157   where
 158     element = flip (!)
 159
 160
 161
 162
 163 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
 164 --   implementation.. it would be a little cleaner to use imap, but it
 165 --   doesn't seem to work.
 166 --
 167 --   TODO: Don't cheat with fromList.
 168 --
 169 --   Examples:
 170 --
 171 --   >>> let m = fromList [[1,2], [3,4]] :: Mat2 Int
 172 --   >>> transpose m
 173 --   ((1,3),(2,4))
 174 --
 175 transpose :: (Arity m, Arity n) => Mat m n a -> Mat n m a
 176 transpose m = Mat $ V.fromList column_list
 177   where
 178     column_list = [ column m i | i <- [0..(ncols m)-1] ]
 179
 180
 181 -- | Is @m@ symmetric?
 182 --
 183 --   Examples:
 184 --
 185 --   >>> let m1 = fromList [[1,2], [2,1]] :: Mat2 Int
 186 --   >>> symmetric m1
 187 --   True
 188 --
 189 --   >>> let m2 = fromList [[1,2], [3,1]] :: Mat2 Int
 190 --   >>> symmetric m2
 191 --   False
 192 --
 193 symmetric :: (Eq a, Arity m) => Mat m m a -> Bool
 194 symmetric m =
 195   m == (transpose m)
 196
 197
 198 -- | Construct a new matrix from a function @lambda@. The function
 199 --   @lambda@ should take two parameters i,j corresponding to the
 200 --   entries in the matrix. The i,j entry of the resulting matrix will
 201 --   have the value returned by lambda i j.
 202 --
 203 --   Examples:
 204 --
 205 --   >>> let lambda i j = i + j
 206 --   >>> construct lambda :: Mat3 Int
 207 --   ((0,1,2),(1,2,3),(2,3,4))
 208 --
 209 construct :: forall m n a. (Arity m, Arity n)
 210           => (Int -> Int -> a) -> Mat m n a
 211 construct lambda = Mat $ generate make_row
 212   where
 213     make_row :: Int -> Vec n a
 214     make_row i = generate (lambda i)
 215
 216
 217 -- | Create an identity matrix with the right dimensions.
 218 --
 219 --   Examples:
 220 --
 221 --   >>> identity_matrix :: Mat3 Int
 222 --   ((1,0,0),(0,1,0),(0,0,1))
 223 --   >>> identity_matrix :: Mat3 Double
 224 --   ((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))
 225 --
 226 identity_matrix :: (Arity m, Ring.C a) => Mat m m a
 227 identity_matrix =
 228   construct (\i j -> if i == j then (fromInteger 1) else (fromInteger 0))
 229
 230 -- | Given a positive-definite matrix @m@, computes the
 231 --   upper-triangular matrix @r@ with (transpose r)*r == m and all
 232 --   values on the diagonal of @r@ positive.
 233 --
 234 --   Examples:
 235 --
 236 --   >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
 237 --   >>> cholesky m1
 238 --   ((4.47213595499958,-0.22360679774997896),(0.0,4.466542286825459))
 239 --   >>> (transpose (cholesky m1)) * (cholesky m1)
 240 --   ((20.000000000000004,-1.0),(-1.0,20.0))
 241 --
 242 cholesky :: forall m n a. (Algebraic.C a, Arity m, Arity n)
 243          => (Mat m n a) -> (Mat m n a)
 244 cholesky m = construct r
 245   where
 246     r :: Int -> Int -> a
 247     r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
 248           | i < j =
 249               (((m !!! (i,j)) - sum [(r k i) NP.* (r k j) | k <- [0..i-1]]))/(r i i)
 250           | otherwise = 0
 251
 252
 253 -- | Returns True if the given matrix is upper-triangular, and False
 254 --   otherwise. The parameter @epsilon@ lets the caller choose a
 255 --   tolerance.
 256 --
 257 --   Examples:
 258 --
 259 --   >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double
 260 --   >>> is_upper_triangular m
 261 --   False
 262 --   >>> is_upper_triangular' 1e-10 m
 263 --   True
 264 --
 265 --   TODO:
 266 --
 267 --     1. Don't cheat with lists.
 268 --
 269 is_upper_triangular' :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
 270                     => a -- ^ The tolerance @epsilon@.
 271                     -> Mat m n a
 272                     -> Bool
 273 is_upper_triangular' epsilon m =
 274   and $ concat results
 275   where
 276     results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ]
 277
 278     test :: Int -> Int -> Bool
 279     test i j
 280       | i <= j = True
 281       -- use "less than or equal to" so zero is a valid epsilon
 282       | otherwise = abs (m !!! (i,j)) <= epsilon
 283
 284
 285 -- | Returns True if the given matrix is upper-triangular, and False
 286 --   otherwise. A specialized version of 'is_upper_triangular\'' with
 287 --   @epsilon = 0@.
 288 --
 289 --   Examples:
 290 --
 291 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 292 --   >>> is_upper_triangular m
 293 --   False
 294 --
 295 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 296 --   >>> is_upper_triangular m
 297 --   True
 298 --
 299 --   TODO:
 300 --
 301 --     1. The Ord constraint is too strong here, Eq would suffice.
 302 --
 303 is_upper_triangular :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
 304                     => Mat m n a -> Bool
 305 is_upper_triangular = is_upper_triangular' 0
 306
 307
 308 -- | Returns True if the given matrix is lower-triangular, and False
 309 --   otherwise. This is a specialized version of 'is_lower_triangular\''
 310 --   with @epsilon = 0@.
 311 --
 312 --   Examples:
 313 --
 314 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 315 --   >>> is_lower_triangular m
 316 --   True
 317 --
 318 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 319 --   >>> is_lower_triangular m
 320 --   False
 321 --
 322 is_lower_triangular :: (Ord a,
 323                         Ring.C a,
 324                         Absolute.C a,
 325                         Arity m,
 326                         Arity n)
 327                     => Mat m n a
 328                     -> Bool
 329 is_lower_triangular = is_upper_triangular . transpose
 330
 331
 332 -- | Returns True if the given matrix is lower-triangular, and False
 333 --   otherwise. The parameter @epsilon@ lets the caller choose a
 334 --   tolerance.
 335 --
 336 --   Examples:
 337 --
 338 --   >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
 339 --   >>> is_lower_triangular m
 340 --   False
 341 --   >>> is_lower_triangular' 1e-12 m
 342 --   True
 343 --
 344 is_lower_triangular' :: (Ord a,
 345                          Ring.C a,
 346                          Absolute.C a,
 347                          Arity m,
 348                          Arity n)
 349                     => a -- ^ The tolerance @epsilon@.
 350                     -> Mat m n a
 351                     -> Bool
 352 is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
 353
 354
 355 -- | Returns True if the given matrix is triangular, and False
 356 --   otherwise.
 357 --
 358 --   Examples:
 359 --
 360 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 361 --   >>> is_triangular m
 362 --   True
 363 --
 364 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 365 --   >>> is_triangular m
 366 --   True
 367 --
 368 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 369 --   >>> is_triangular m
 370 --   False
 371 --
 372 is_triangular :: (Ord a,
 373                   Ring.C a,
 374                   Absolute.C a,
 375                   Arity m,
 376                   Arity n)
 377               => Mat m n a
 378               -> Bool
 379 is_triangular m = is_upper_triangular m || is_lower_triangular m
 380
 381
 382 -- | Return the (i,j)th minor of m.
 383 --
 384 --   Examples:
 385 --
 386 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 387 --   >>> minor m 0 0 :: Mat2 Int
 388 --   ((5,6),(8,9))
 389 --   >>> minor m 1 1 :: Mat2 Int
 390 --   ((1,3),(7,9))
 391 --
 392 minor :: (m ~ S r,
 393           n ~ S t,
 394           Arity r,
 395           Arity t)
 396       => Mat m n a
 397       -> Int
 398       -> Int
 399       -> Mat r t a
 400 minor (Mat rows) i j = m
 401   where
 402     rows' = delete rows i
 403     m = Mat $ V.map ((flip delete) j) rows'
 404
 405
 406 class (Eq a, Ring.C a) => Determined p a where
 407   determinant :: (p a) -> a
 408
 409 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
 410   determinant (Mat rows) = (V.head . V.head) rows
 411
 412 instance (Ord a,
 413           Ring.C a,
 414           Absolute.C a,
 415           Arity n,
 416           Determined (Mat (S n) (S n)) a)
 417          => Determined (Mat (S (S n)) (S (S n))) a where
 418   -- | The recursive definition with a special-case for triangular matrices.
 419   --
 420   --   Examples:
 421   --
 422   --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 423   --   >>> determinant m
 424   --   -1
 425   --
 426   determinant m
 427     | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
 428     | otherwise = determinant_recursive
 429         where
 430           m' i j = m !!! (i,j)
 431
 432           det_minor i j = determinant (minor m i j)
 433
 434           determinant_recursive =
 435             sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (det_minor 0 j)
 436               | j <- [0..(ncols m)-1] ]
 437
 438
 439
 440 -- | Matrix multiplication.
 441 --
 442 --   Examples:
 443 --
 444 --   >>> let m1 = fromList [[1,2,3], [4,5,6]]  :: Mat N2 N3 Int
 445 --   >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat N3 N2 Int
 446 --   >>> m1 * m2
 447 --   ((22,28),(49,64))
 448 --
 449 infixl 7 *
 450 (*) :: (Ring.C a, Arity m, Arity n, Arity p)
 451         => Mat m n a
 452         -> Mat n p a
 453         -> Mat m p a
 454 (*) m1 m2 = construct lambda
 455   where
 456     lambda i j =
 457       sum [(m1 !!! (i,k)) NP.* (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]
 458
 459
 460
 461 instance (Ring.C a, Arity m, Arity n) => Additive.C (Mat m n a) where
 462
 463   (Mat rows1) + (Mat rows2) =
 464     Mat $ V.zipWith (V.zipWith (+)) rows1 rows2
 465
 466   (Mat rows1) - (Mat rows2) =
 467     Mat $ V.zipWith (V.zipWith (-)) rows1 rows2
 468
 469   zero = Mat (V.replicate $ V.replicate (fromInteger 0))
 470
 471
 472 instance (Ring.C a, Arity m, Arity n, m ~ n) => Ring.C (Mat m n a) where
 473   -- The first * is ring multiplication, the second is matrix
 474   -- multiplication.
 475   m1 * m2 = m1 * m2
 476
 477
 478 instance (Ring.C a, Arity m, Arity n) => Module.C a (Mat m n a) where
 479   -- We can multiply a matrix by a scalar of the same type as its
 480   -- elements.
 481   x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows
 482
 483
 484 instance (Algebraic.C a,
 485           ToRational.C a,
 486           Arity m)
 487          => Normed (Mat (S m) N1 a) where
 488   -- | Generic p-norms for vectors in R^n that are represented as nx1
 489   --   matrices.
 490   --
 491   --   Examples:
 492   --
 493   --   >>> let v1 = vec2d (3,4)
 494   --   >>> norm_p 1 v1
 495   --   7.0
 496   --   >>> norm_p 2 v1
 497   --   5.0
 498   --
 499   norm_p p (Mat rows) =
 500     (root p') $ sum [fromRational' (toRational x)^p' | x <- xs]
 501     where
 502       p' = toInteger p
 503       xs = concat $ V.toList $ V.map V.toList rows
 504
 505   -- | The infinity norm.
 506   --
 507   --   Examples:
 508   --
 509   --   >>> let v1 = vec3d (1,5,2)
 510   --   >>> norm_infty v1
 511   --   5
 512   --
 513   norm_infty (Mat rows) =
 514     fromRational' $ toRational $ V.maximum $ V.map V.maximum rows
 515
 516
 517 -- | Compute the Frobenius norm of a matrix. This essentially treats
 518 --   the matrix as one long vector containing all of its entries (in
 519 --   any order, it doesn't matter).
 520 --
 521 --   Examples:
 522 --
 523 --   >>> let m = fromList [[1, 2, 3],[4,5,6],[7,8,9]] :: Mat3 Double
 524 --   >>> frobenius_norm m == sqrt 285
 525 --   True
 526 --
 527 --   >>> let m = fromList [[1, -1, 1],[-1,1,-1],[1,-1,1]] :: Mat3 Double
 528 --   >>> frobenius_norm m == 3
 529 --   True
 530 --
 531 frobenius_norm :: (Algebraic.C a, Ring.C a) => Mat m n a -> a
 532 frobenius_norm (Mat rows) =
 533   sqrt $ element_sum $ V.map row_sum rows
 534   where
 535     -- | Square and add up the entries of a row.
 536     row_sum = element_sum . V.map (^2)
 537
 538
 539 -- Vector helpers. We want it to be easy to create low-dimension
 540 -- column vectors, which are nx1 matrices.
 541
 542 -- | Convenient constructor for 2D vectors.
 543 --
 544 --   Examples:
 545 --
 546 --   >>> import Roots.Simple
 547 --   >>> let fst m = m !!! (0,0)
 548 --   >>> let snd m = m !!! (1,0)
 549 --   >>> let h = 0.5 :: Double
 550 --   >>> let g1 m = 1.0 + h NP.* exp(-((fst m)^2))/(1.0 + (snd m)^2)
 551 --   >>> let g2 m = 0.5 + h NP.* atan((fst m)^2 + (snd m)^2)
 552 --   >>> let g u = vec2d ((g1 u), (g2 u))
 553 --   >>> let u0 = vec2d (1.0, 1.0)
 554 --   >>> let eps = 1/(10^9)
 555 --   >>> fixed_point g eps u0
 556 --   ((1.0728549599342185),(1.0820591495686167))
 557 --
 558 vec1d :: (a) -> Mat N1 N1 a
 559 vec1d (x) = Mat (mk1 (mk1 x))
 560
 561 vec2d :: (a,a) -> Mat N2 N1 a
 562 vec2d (x,y) = Mat (mk2 (mk1 x) (mk1 y))
 563
 564 vec3d :: (a,a,a) -> Mat N3 N1 a
 565 vec3d (x,y,z) = Mat (mk3 (mk1 x) (mk1 y) (mk1 z))
 566
 567 vec4d :: (a,a,a,a) -> Mat N4 N1 a
 568 vec4d (w,x,y,z) = Mat (mk4 (mk1 w) (mk1 x) (mk1 y) (mk1 z))
 569
 570 vec5d :: (a,a,a,a,a) -> Mat N5 N1 a
 571 vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))
 572
 573 -- Since we commandeered multiplication, we need to create 1x1
 574 -- matrices in order to multiply things.
 575 scalar :: a -> Mat N1 N1 a
 576 scalar x = Mat (mk1 (mk1 x))
 577
 578 dot :: (RealRing.C a, n ~ N1, m ~ S t, Arity t)
 579        => Mat m n a
 580        -> Mat m n a
 581        -> a
 582 v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)
 583
 584
 585 -- | The angle between @v1@ and @v2@ in Euclidean space.
 586 --
 587 --   Examples:
 588 --
 589 --   >>> let v1 = vec2d (1.0, 0.0)
 590 --   >>> let v2 = vec2d (0.0, 1.0)
 591 --   >>> angle v1 v2 == pi/2.0
 592 --   True
 593 --
 594 angle :: (Transcendental.C a,
 595           RealRing.C a,
 596           n ~ N1,
 597           m ~ S t,
 598           Arity t,
 599           ToRational.C a)
 600           => Mat m n a
 601           -> Mat m n a
 602           -> a
 603 angle v1 v2 =
 604   acos theta
 605   where
 606    theta = (recip norms) NP.* (v1 `dot` v2)
 607    norms = (norm v1) NP.* (norm v2)
 608
 609
 610 -- | Retrieve the diagonal elements of the given matrix as a \"column
 611 --   vector,\" i.e. a m-by-1 matrix. We require the matrix to be
 612 --   square to avoid ambiguity in the return type which would ideally
 613 --   have dimension min(m,n) supposing an m-by-n matrix.
 614 --
 615 --   Examples:
 616 --
 617 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 618 --   >>> diagonal m
 619 --   ((1),(5),(9))
 620 --
 621 diagonal :: (Arity m) => Mat m m a -> Mat m N1 a
 622 diagonal matrix =
 623   construct lambda
 624   where
 625     lambda i _ = matrix !!! (i,i)
 626
 627
 628 -- | Given a square @matrix@, return a new matrix of the same size
 629 --   containing only the on-diagonal entries of @matrix@. The
 630 --   off-diagonal entries are set to zero.
 631 --
 632 --   Examples:
 633 --
 634 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 635 --   >>> diagonal_part m
 636 --   ((1,0,0),(0,5,0),(0,0,9))
 637 --
 638 diagonal_part :: (Arity m, Ring.C a)
 639          => Mat m m a
 640          -> Mat m m a
 641 diagonal_part matrix =
 642   construct lambda
 643   where
 644     lambda i j = if i == j then matrix !!! (i,j) else 0
 645
 646
 647 -- | Given a square @matrix@, return a new matrix of the same size
 648 --   containing only the on-diagonal and below-diagonal entries of
 649 --   @matrix@. The above-diagonal entries are set to zero.
 650 --
 651 --   Examples:
 652 --
 653 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 654 --   >>> lt_part m
 655 --   ((1,0,0),(4,5,0),(7,8,9))
 656 --
 657 lt_part :: (Arity m, Ring.C a)
 658         => Mat m m a
 659         -> Mat m m a
 660 lt_part matrix =
 661   construct lambda
 662   where
 663     lambda i j = if i >= j then matrix !!! (i,j) else 0
 664
 665
 666 -- | Given a square @matrix@, return a new matrix of the same size
 667 --   containing only the below-diagonal entries of @matrix@. The on-
 668 --   and above-diagonal entries are set to zero.
 669 --
 670 --   Examples:
 671 --
 672 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 673 --   >>> lt_part_strict m
 674 --   ((0,0,0),(4,0,0),(7,8,0))
 675 --
 676 lt_part_strict :: (Arity m, Ring.C a)
 677         => Mat m m a
 678         -> Mat m m a
 679 lt_part_strict matrix =
 680   construct lambda
 681   where
 682     lambda i j = if i > j then matrix !!! (i,j) else 0
 683
 684
 685 -- | Given a square @matrix@, return a new matrix of the same size
 686 --   containing only the on-diagonal and above-diagonal entries of
 687 --   @matrix@. The below-diagonal entries are set to zero.
 688 --
 689 --   Examples:
 690 --
 691 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 692 --   >>> ut_part m
 693 --   ((1,2,3),(0,5,6),(0,0,9))
 694 --
 695 ut_part :: (Arity m, Ring.C a)
 696         => Mat m m a
 697         -> Mat m m a
 698 ut_part = transpose . lt_part . transpose
 699
 700
 701 -- | Given a square @matrix@, return a new matrix of the same size
 702 --   containing only the above-diagonal entries of @matrix@. The on-
 703 --   and below-diagonal entries are set to zero.
 704 --
 705 --   Examples:
 706 --
 707 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 708 --   >>> ut_part_strict m
 709 --   ((0,2,3),(0,0,6),(0,0,0))
 710 --
 711 ut_part_strict :: (Arity m, Ring.C a)
 712         => Mat m m a
 713         -> Mat m m a
 714 ut_part_strict = transpose . lt_part_strict . transpose
 715
 716
 717 -- | Compute the trace of a square matrix, the sum of the elements
 718 --   which lie on its diagonal. We require the matrix to be
 719 --   square to avoid ambiguity in the return type which would ideally
 720 --   have dimension min(m,n) supposing an m-by-n matrix.
 721 --
 722 --   Examples:
 723 --
 724 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 725 --   >>> trace m
 726 --   15
 727 --
 728 trace :: (Arity m, Ring.C a) => Mat m m a -> a
 729 trace matrix =
 730   let (Mat rows) = diagonal matrix
 731   in
 732     element_sum $ V.map V.head rows