Add the Linear.Iteration module.

diff --git a/src/Linear/Iteration.hs b/src/Linear/Iteration.hs
new file mode 100644 (file)
index 0000000..31c7c73
--- /dev/null
+++ b/src/Linear/Iteration.hs
@@ -0,0 +1,114 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeFamilies #-}
+
+-- | Classical iterative methods to solve the system Ax = b.
+
+module Linear.Iteration
+where
+
+import Data.List (find)
+import Data.Maybe (fromJust)
+import Data.Vector.Fixed (Arity, N1, S)
+import NumericPrelude hiding ((*))
+import qualified Algebra.Algebraic as Algebraic
+import qualified Algebra.Field as Field
+import qualified Algebra.RealField as RealField
+import qualified Algebra.ToRational as ToRational
+import qualified Prelude as P
+
+import Linear.Matrix
+import Linear.System
+import Normed
+
+-- | Perform one Jacobi iteration,
+--
+--     x1 = M^(-1) * (N*x0 + b)
+--
+--   Examples:
+--
+--   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
+--   >>> let x0 = vec2d (0, 0::Double)
+--   >>> let b  = vec2d (1, 1::Double)
+--   >>> jacobi_iteration m b x0
+--   ((0.25),(0.5))
+--   >>> let x1 = jacobi_iteration m b x0
+--   >>> jacobi_iteration m b x1
+--   ((0.0),(0.25))
+--
+jacobi_iteration :: (Field.C a, Arity m)
+                 => Mat m m  a
+                 -> Mat m N1 a
+                 -> Mat m N1 a
+                 -> Mat m N1 a
+jacobi_iteration matrix b x_current =
+  x_next
+  where
+    big_m = diagonal matrix
+    big_n = big_m - matrix
+    rhs = big_n*x_current + b
+    x_next = forward_substitute big_m rhs
+
+
+-- | Compute an infinite list of Jacobi iterations starting with the
+--   vector x0.
+jacobi_iterations :: (Field.C a, Arity m)
+                  => Mat m m  a
+                  -> Mat m N1 a
+                  -> Mat m N1 a
+                  -> [Mat m N1 a]
+jacobi_iterations matrix b =
+  iterate (jacobi_iteration matrix b)
+
+
+-- | Solve the system Ax = b using the Jacobi method. This will run
+--   forever if the iterations do not converge.
+--
+--   Examples:
+--
+--   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
+--   >>> let x0 = vec2d (0, 0::Double)
+--   >>> let b  = vec2d (1, 1::Double)
+--   >>> let epsilon = 10**(-6)
+--   >>> jacobi_method m b x0 epsilon
+--   ((0.0),(0.4999995231628418))
+--
+jacobi_method :: forall m n a b.
+                 (RealField.C a,
+                  Algebraic.C a, -- Normed instance
+                  ToRational.C a, -- Normed instance
+                  Algebraic.C b,
+                  RealField.C b,
+                  Arity m,
+                  Arity n, -- Normed instance
+                  m ~ S n)
+              => Mat m m  a
+              -> Mat m N1 a
+              -> Mat m N1 a
+              -> b
+              -> Mat m N1 a
+jacobi_method matrix b x0 epsilon =
+  -- fromJust is "safe," because the list is infinite. If the
+  -- algorithm doesn't converge, 'find' will search forever and never
+  -- return Nothing.
+  fst' $ fromJust $ find error_small_enough diff_pairs
+  where
+    x_n = jacobi_iterations matrix b x0
+
+    pairs :: [(Mat m N1 a, Mat m N1 a)]
+    pairs = zip (tail x_n) x_n
+
+    append_diff :: (Mat m N1 a, Mat m N1 a)
+                -> (Mat m N1 a, Mat m N1 a, b)
+    append_diff (cur,prev) =
+      (cur,prev,diff)
+      where
+        diff = norm (cur - prev)
+
+    diff_pairs :: [(Mat m N1 a, Mat m N1 a, b)]
+    diff_pairs = map append_diff pairs
+
+    fst' :: (c, d, e) -> c
+    fst' (x,_,_) = x
+
+    error_small_enough :: (Mat m N1 a, Mat m N1 a, b)-> Bool
+    error_small_enough (_,_,err) = err < epsilon