+jacobi_method =
+ classical_method jacobi_iterations
+
+
+-- | Solve the system Ax = b using the Gauss-Seidel method. This will
+-- run forever if the iterations do not converge.
+--
+-- Examples:
+--
+-- >>> import Linear.Matrix ( Mat2, fromList, vec2d )
+--
+-- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
+-- >>> let x0 = vec2d (0, 0::Double)
+-- >>> let b = vec2d (1, 1::Double)
+-- >>> let epsilon = 10**(-12)
+-- >>> gauss_seidel_method m b x0 epsilon
+-- ((4.547473508864641e-13),(0.49999999999954525))
+--
+gauss_seidel_method :: (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
+gauss_seidel_method =
+ classical_method gauss_seidel_iterations
+
+
+-- | Solve the system Ax = b using the Successive Over-Relaxation
+-- (SOR) method. This will run forever if the iterations do not
+-- converge.
+--
+-- Examples:
+--
+-- >>> import Linear.Matrix ( Mat2, fromList, vec2d )
+--
+-- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
+-- >>> let x0 = vec2d (0, 0::Double)
+-- >>> let b = vec2d (1, 1::Double)
+-- >>> let epsilon = 10**(-12)
+-- >>> sor_method 1.5 m b x0 epsilon
+-- ((6.567246746413957e-13),(0.4999999999993727))
+--
+sor_method :: (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)
+ => a
+ -> Mat m m a
+ -> Mat m N1 a
+ -> Mat m N1 a
+ -> b
+ -> Mat m N1 a
+sor_method omega =
+ classical_method (sor_iterations omega)
+
+
+-- | General implementation for all classical iteration methods. For
+-- its first argument, it takes a function which generates the
+-- sequence of iterates when supplied with the remaining arguments
+-- (except for the tolerance).
+--
+classical_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 -> [Mat m N1 a])
+ -> Mat m m a
+ -> Mat m N1 a
+ -> Mat m N1 a
+ -> b
+ -> Mat m N1 a
+classical_method iterations_function matrix b x0 epsilon =