Add the rayleigh_quotient function.

diff --git a/src/Linear/Iteration.hs b/src/Linear/Iteration.hs
index 89f2f0c0eab7788a088642afd81e98525478e610..7daf3b69b5f66c536d1a2567b7896fe54c972c57 100644 (file)
--- a/src/Linear/Iteration.hs
+++ b/src/Linear/Iteration.hs
@@ -261,3 +261,28 @@ 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:
+--
+--   >>> 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)