X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FLinear%2FIteration.hs;h=6117770a52313c21546ed57e729854a7d1fc5f8c;hb=64e43e504a8716cb1784de5fc33d7f02e915e2ac;hp=89f2f0c0eab7788a088642afd81e98525478e610;hpb=eb2fb001d7b7a7769844b6e5c3623e7501a7ac12;p=numerical-analysis.git

diff --git a/src/Linear/Iteration.hs b/src/Linear/Iteration.hs
index 89f2f0c..6117770 100644
--- a/src/Linear/Iteration.hs
+++ b/src/Linear/Iteration.hs
@@ -1,34 +1,53 @@
+{-# 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
@@ -41,8 +60,8 @@ classical_iteration m_function matrix b x_current =
 
 -- | 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
@@ -61,7 +80,7 @@ sor_iteration omega =
 
 -- | 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
@@ -72,7 +91,7 @@ sor_iterations omega matrix b =
 
 
 -- | 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
@@ -82,7 +101,7 @@ gauss_seidel_iteration = sor_iteration one
 
 -- | 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
@@ -97,6 +116,8 @@ gauss_seidel_iterations matrix b =
 --
 --   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)
@@ -106,7 +127,7 @@ gauss_seidel_iterations matrix b =
 --   >>> 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
@@ -117,7 +138,7 @@ jacobi_iteration =
 
 -- | 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
@@ -131,6 +152,8 @@ jacobi_iterations matrix b =
 --
 --   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)
@@ -160,6 +183,8 @@ jacobi_method =
 --
 --   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)
@@ -190,6 +215,8 @@ gauss_seidel_method =
 --
 --   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)
@@ -261,3 +288,30 @@ classical_method iterations_function matrix b x0 epsilon =
 
     error_small_enough :: (Mat m N1 a, Mat m N1 a, b)-> Bool
     error_small_enough (_,_,err) = err < epsilon
+
+
+
+-- | Compute the Rayleigh quotient of  @matrix@ and @vector@.
+--
+--   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
+--   7 % 2
+--
+rayleigh_quotient :: (RealField.C a,
+                      Arity m,
+                      Arity n,
+                      m ~ S n)
+                  => (Mat m m a)
+                  -> (Mat m N1 a)
+                  -> a
+rayleigh_quotient matrix vector =
+  (vector `dot` (matrix * vector)) / (norm_squared vector)
+  where
+    -- We don't use the norm function here to avoid the algebraic
+    -- requirement on our field.
+    norm_squared v = ((transpose v) * v) !!! (0,0)