+
+
+-- | Unsafely set the (i,j) element of the given matrix.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> set_idx m (1,1) 17
+-- ((1,2,3),(4,17,6),(7,8,9))
+--
+set_idx :: forall m n a.
+ (Arity m, Arity n)
+ => Mat m n a
+ -> (Int, Int)
+ -> a
+ -> Mat m n a
+set_idx matrix (i,j) newval =
+ imap2 updater matrix
+ where
+ updater :: Int -> Int -> a -> a
+ updater k l existing =
+ if k == i && l == j
+ then newval
+ else existing
+
+
+-- | Compute the i,jth cofactor of the given @matrix@. This simply
+-- premultiplues the i,jth minor by (-1)^(i+j).
+cofactor :: (Arity m, Determined (Mat m m) a)
+ => Mat (S m) (S m) a
+ -> Int
+ -> Int
+ -> a
+cofactor matrix i j =
+ (-1)^(toInteger i + toInteger j) NP.* (minor matrix i j)
+
+
+-- | Compute the inverse of a matrix using cofactor expansion
+-- (generalized Cramer's rule).
+--
+-- Examples:
+--
+-- >>> let m1 = fromList [[37,22],[17,54]] :: Mat2 Double
+-- >>> let e1 = [54/1624, -22/1624] :: [Double]
+-- >>> let e2 = [-17/1624, 37/1624] :: [Double]
+-- >>> let expected = fromList [e1, e2] :: Mat2 Double
+-- >>> let actual = inverse m1
+-- >>> frobenius_norm (actual - expected) < 1e-12
+-- True
+--
+inverse :: (Arity m,
+ Determined (Mat (S m) (S m)) a,
+ Determined (Mat m m) a,
+ Field.C a)
+ => Mat (S m) (S m) a
+ -> Mat (S m) (S m) a
+inverse matrix =
+ (1 / (determinant matrix)) *> (transpose $ construct lambda)
+ where
+ lambda i j = cofactor matrix i j
+
+
+
+-- | Retrieve the rows of a matrix as a column matrix. If the given
+-- matrix is m-by-n, the result would be an m-by-1 column whose
+-- entries are 1-by-n row matrices.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+-- >>> (rows2 m) !!! (0,0)
+-- ((1,2))
+-- >>> (rows2 m) !!! (1,0)
+-- ((3,4))
+--
+rows2 :: (Arity m, Arity n)
+ => Mat m n a
+ -> Col m (Row n a)
+rows2 (Mat rows) =
+ Mat $ V.map (mk1. Mat . mk1) rows
+
+
+
+-- | Sum the elements of a matrix.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,-1],[3,4]] :: Mat2 Int
+-- >>> element_sum2 m
+-- 7
+--
+element_sum2 :: (Arity m, Arity n, Additive.C a) => Mat m n a -> a
+element_sum2 = foldl2 (+) zero