+{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
-- | Classical iterative methods to solve the system Ax = b.
-module Linear.Iteration
+module Linear.Iteration (
+ gauss_seidel_iteration,
+ gauss_seidel_iterations,
+ gauss_seidel_method,
+ jacobi_iteration,
+ jacobi_iterations,
+ jacobi_method,
+ rayleigh_quotient,
+ sor_iteration,
+ sor_iterations,
+ sor_method )
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 Data.List ( find )
+import Data.Maybe ( fromJust )
+import Data.Vector.Fixed ( Arity, N1, S )
+import NumericPrelude hiding ( (*) )
+import qualified Algebra.Algebraic as Algebraic ( C )
+import qualified Algebra.Field as Field ( C )
+import qualified Algebra.RealField as RealField ( C )
+import qualified Algebra.ToRational as ToRational ( C )
+
+import Linear.Matrix (
+ Col,
+ Mat(..),
+ (!!!),
+ (*),
+ diagonal_part,
+ dot,
+ lt_part_strict,
+ transpose )
+import Linear.System ( forward_substitute )
+import Normed ( Normed(..) )
-import Linear.Matrix
-import Linear.System
-import Normed
-- | A generalized implementation for Jacobi, Gauss-Seidel, etc. All
-- that we really need to know is how to construct the matrix M, so we
-- take a function that does it as an argument.
-classical_iteration :: (Field.C a, Arity m)
- => (Mat m m a -> Mat m m a)
- -> Mat m m a
- -> Mat m N1 a
- -> Mat m N1 a
- -> Mat m N1 a
+classical_iteration :: (Eq a, Field.C a, m ~ S n, Arity n)
+ => (Mat m m a -> Mat m m a)
+ -> Mat m m a
+ -> Col m a
+ -> Col m a
+ -> Col m a
classical_iteration m_function matrix b x_current =
x_next
where
-- | Perform one iteration of successive over-relaxation.
--
-sor_iteration :: forall m a.
- (Field.C a, Arity m)
+sor_iteration :: forall m n a.
+ (Eq a, Field.C a, m ~ S n, Arity n)
=> a -- ^ Omega
-> Mat m m a -- ^ Matrix A
-> Mat m N1 a -- ^ Vector b
-- | Compute an infinite list of SOR iterations starting with the
-- vector x0.
-sor_iterations :: (Field.C a, Arity m)
+sor_iterations :: (Eq a, Field.C a, m ~ S n, Arity n)
=> a
-> Mat m m a
-> Mat m N1 a
-- | Perform one iteration of Gauss-Seidel.
-gauss_seidel_iteration :: (Field.C a, Arity m)
+gauss_seidel_iteration :: (Eq a, Field.C a, m ~ S n, Arity n)
=> Mat m m a
-> Mat m N1 a
-> Mat m N1 a
-- | Compute an infinite list of Gauss-Seidel iterations starting with
-- the vector x0.
-gauss_seidel_iterations :: (Field.C a, Arity m)
+gauss_seidel_iterations :: (Eq a, Field.C a, m ~ S n, Arity n)
=> Mat m m a
-> Mat m N1 a
-> Mat m N1 a
--
-- 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)
-- >>> jacobi_iteration m b x1
-- ((0.0),(0.25))
--
-jacobi_iteration :: (Field.C a, Arity m)
+jacobi_iteration :: (Eq a, Field.C a, m ~ S n, Arity n)
=> Mat m m a
-> Mat m N1 a
-> Mat m N1 a
-- | Compute an infinite list of Jacobi iterations starting with the
-- vector x0.
-jacobi_iterations :: (Field.C a, Arity m)
+jacobi_iterations :: (Eq a, Field.C a, m ~ S n, Arity n)
=> Mat m m a
-> Mat m N1 a
-> Mat m N1 a
--
-- 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)
--
-- 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)
--
-- 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)
--
-- Examples:
--
+-- >>> import Linear.Matrix ( Mat2, fromList, vec2d )
+--
-- >>> let m = fromList [[3,1],[1,2]] :: Mat2 Rational
-- >>> let v = vec2d (1, 1::Rational)
-- >>> rayleigh_quotient m v