Add big_K code to FEM.R1.

[numerical-analysis.git] / src / Linear / Matrix.hs

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index 3d4daab3dc1f7be81277238b52253c5df54ae61f..119d77089461dfae968320c79bd1375c31d46b41 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -37,6 +37,7 @@ import qualified Data.Vector.Fixed as V (
   fromList,
   head,
   ifoldl,
+  ifoldr,
   imap,
   map,
   maximum,
@@ -662,11 +663,16 @@ vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))
 scalar :: a -> Mat1 a
 scalar x = Mat (mk1 (mk1 x))
 
-dot :: (RealRing.C a, m ~ S t, Arity t)
-       => Col m a
-       -> Col m a
+-- Get the scalar value out of a 1x1 matrix.
+unscalar :: Mat1 a -> a
+unscalar (Mat rows) = V.head $ V.head rows
+
+
+dot :: (Ring.C a, Arity m)
+       => Col (S m) a
+       -> Col (S m) a
        -> a
-v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)
+v1 `dot` v2 = unscalar $ ((transpose v1) * v2)
 
 
 -- | The angle between @v1@ and @v2@ in Euclidean space.
@@ -906,7 +912,8 @@ map2 f (Mat rows) =
 
 
 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
---   (of the row/column) to the accumulation function.
+--   (of the row/column) to the accumulation function. The fold occurs
+--   from top-left to bottom-right.
 --
 --   Examples:
 --
@@ -933,6 +940,36 @@ ifoldl2 f initial (Mat rows) =
     row_function rowinit idx r = V.ifoldl (g idx) rowinit r
 
 
+-- | Fold over the entire matrix passing the coordinates @i@ and @j@
+--   (of the row/column) to the accumulation function. The fold occurs
+--   from bottom-right to top-left.
+--
+--   The order of the arguments in the supplied function are different
+--   from those in V.ifoldr; we keep them similar to ifoldl2.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+--   >>> ifoldr2 (\i j cur _ -> cur + i + j) 0 m
+--   18
+--
+ifoldr2 :: forall a b m n.
+           (Int -> Int -> b -> a -> b)
+        -> b
+        -> Mat m n a
+        -> b
+ifoldr2 f initial (Mat rows) =
+  V.ifoldr row_function initial rows
+  where
+    -- | Swap the order of arguments in @f@ so that it agrees with the
+    --   @f@ passed to ifoldl2.
+    g :: Int -> Int -> a -> b -> b
+    g w x y z = f w x z y
+
+    row_function :: Int -> Vec n a -> b -> b
+    row_function idx r rowinit = V.ifoldr (g idx) rowinit r
+
+
 -- | Map a function over a matrix of any dimensions, passing the
 --   coordinates @i@ and @j@ to the function @f@.
 --
@@ -969,3 +1006,25 @@ reverse2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a
 reverse2 (Mat rows) = Mat $ V.reverse $ V.map V.reverse rows
 
 
+-- | 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