1 {-# LANGUAGE NoImplicitPrelude #-}
2 {-# LANGUAGE ScopedTypeVariables #-}
3 {-# LANGUAGE TypeFamilies #-}
5 -- | Classical iterative methods to solve the system Ax = b.
7 module Linear.Iteration (
8 gauss_seidel_iteration,
9 gauss_seidel_iterations,
20 import Data.List ( find )
21 import Data.Maybe ( fromJust )
22 import Data.Vector.Fixed ( Arity, N1, S )
23 import NumericPrelude hiding ( (*) )
24 import qualified Algebra.Algebraic as Algebraic ( C )
25 import qualified Algebra.Field as Field ( C )
26 import qualified Algebra.RealField as RealField ( C )
27 import qualified Algebra.ToRational as ToRational ( C )
29 import Linear.Matrix (
37 import Linear.System ( forward_substitute )
38 import Normed ( Normed(..) )
41 -- | A generalized implementation for Jacobi, Gauss-Seidel, etc. All
42 -- that we really need to know is how to construct the matrix M, so we
43 -- take a function that does it as an argument.
44 classical_iteration :: (Field.C a, Arity m)
45 => (Mat m m a -> Mat m m a)
50 classical_iteration m_function matrix b x_current =
53 big_m = m_function matrix
54 big_n = big_m - matrix
55 rhs = big_n*x_current + b
56 -- TODO: Should be solve below! M might not be lower-triangular.
57 x_next = forward_substitute big_m rhs
60 -- | Perform one iteration of successive over-relaxation.
62 sor_iteration :: forall m a.
65 -> Mat m m a -- ^ Matrix A
66 -> Mat m N1 a -- ^ Vector b
67 -> Mat m N1 a -- ^ Vector x_current
68 -> Mat m N1 a -- ^ Output vector x_next
70 classical_iteration m_function
72 m_function :: Mat m m a -> Mat m m a
74 let diag = (recip omega) *> (diagonal_part matrix)
75 lt = lt_part_strict matrix
80 -- | Compute an infinite list of SOR iterations starting with the
82 sor_iterations :: (Field.C a, Arity m)
88 sor_iterations omega matrix b =
89 iterate (sor_iteration omega matrix b)
92 -- | Perform one iteration of Gauss-Seidel.
93 gauss_seidel_iteration :: (Field.C a, Arity m)
98 gauss_seidel_iteration = sor_iteration one
101 -- | Compute an infinite list of Gauss-Seidel iterations starting with
103 gauss_seidel_iterations :: (Field.C a, Arity m)
108 gauss_seidel_iterations matrix b =
109 iterate (gauss_seidel_iteration matrix b)
112 -- | Perform one Jacobi iteration,
114 -- x1 = M^(-1) * (N*x0 + b)
118 -- >>> import Linear.Matrix ( Mat2, fromList, vec2d )
120 -- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
121 -- >>> let x0 = vec2d (0, 0::Double)
122 -- >>> let b = vec2d (1, 1::Double)
123 -- >>> jacobi_iteration m b x0
125 -- >>> let x1 = jacobi_iteration m b x0
126 -- >>> jacobi_iteration m b x1
129 jacobi_iteration :: (Field.C a, Arity m)
135 classical_iteration diagonal_part
138 -- | Compute an infinite list of Jacobi iterations starting with the
140 jacobi_iterations :: (Field.C a, Arity m)
145 jacobi_iterations matrix b =
146 iterate (jacobi_iteration matrix b)
149 -- | Solve the system Ax = b using the Jacobi method. This will run
150 -- forever if the iterations do not converge.
154 -- >>> import Linear.Matrix ( Mat2, fromList, vec2d )
156 -- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
157 -- >>> let x0 = vec2d (0, 0::Double)
158 -- >>> let b = vec2d (1, 1::Double)
159 -- >>> let epsilon = 10**(-6)
160 -- >>> jacobi_method m b x0 epsilon
161 -- ((0.0),(0.4999995231628418))
163 jacobi_method :: (RealField.C a,
164 Algebraic.C a, -- Normed instance
165 ToRational.C a, -- Normed instance
169 Arity n, -- Normed instance
177 classical_method jacobi_iterations
180 -- | Solve the system Ax = b using the Gauss-Seidel method. This will
181 -- run forever if the iterations do not converge.
185 -- >>> import Linear.Matrix ( Mat2, fromList, vec2d )
187 -- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
188 -- >>> let x0 = vec2d (0, 0::Double)
189 -- >>> let b = vec2d (1, 1::Double)
190 -- >>> let epsilon = 10**(-12)
191 -- >>> gauss_seidel_method m b x0 epsilon
192 -- ((4.547473508864641e-13),(0.49999999999954525))
194 gauss_seidel_method :: (RealField.C a,
195 Algebraic.C a, -- Normed instance
196 ToRational.C a, -- Normed instance
200 Arity n, -- Normed instance
207 gauss_seidel_method =
208 classical_method gauss_seidel_iterations
211 -- | Solve the system Ax = b using the Successive Over-Relaxation
212 -- (SOR) method. This will run forever if the iterations do not
217 -- >>> import Linear.Matrix ( Mat2, fromList, vec2d )
219 -- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
220 -- >>> let x0 = vec2d (0, 0::Double)
221 -- >>> let b = vec2d (1, 1::Double)
222 -- >>> let epsilon = 10**(-12)
223 -- >>> sor_method 1.5 m b x0 epsilon
224 -- ((6.567246746413957e-13),(0.4999999999993727))
226 sor_method :: (RealField.C a,
227 Algebraic.C a, -- Normed instance
228 ToRational.C a, -- Normed instance
232 Arity n, -- Normed instance
241 classical_method (sor_iterations omega)
244 -- | General implementation for all classical iteration methods. For
245 -- its first argument, it takes a function which generates the
246 -- sequence of iterates when supplied with the remaining arguments
247 -- (except for the tolerance).
249 classical_method :: forall m n a b.
251 Algebraic.C a, -- Normed instance
252 ToRational.C a, -- Normed instance
256 Arity n, -- Normed instance
258 => (Mat m m a -> Mat m N1 a -> Mat m N1 a -> [Mat m N1 a])
264 classical_method iterations_function matrix b x0 epsilon =
265 -- fromJust is "safe," because the list is infinite. If the
266 -- algorithm doesn't converge, 'find' will search forever and never
268 fst' $ fromJust $ find error_small_enough diff_pairs
270 x_n = iterations_function matrix b x0
272 pairs :: [(Mat m N1 a, Mat m N1 a)]
273 pairs = zip (tail x_n) x_n
275 append_diff :: (Mat m N1 a, Mat m N1 a)
276 -> (Mat m N1 a, Mat m N1 a, b)
277 append_diff (cur,prev) =
280 diff = norm (cur - prev)
282 diff_pairs :: [(Mat m N1 a, Mat m N1 a, b)]
283 diff_pairs = map append_diff pairs
285 fst' :: (c, d, e) -> c
288 error_small_enough :: (Mat m N1 a, Mat m N1 a, b)-> Bool
289 error_small_enough (_,_,err) = err < epsilon
293 -- | Compute the Rayleigh quotient of @matrix@ and @vector@.
297 -- >>> import Linear.Matrix ( Mat2, fromList, vec2d )
299 -- >>> let m = fromList [[3,1],[1,2]] :: Mat2 Rational
300 -- >>> let v = vec2d (1, 1::Rational)
301 -- >>> rayleigh_quotient m v
304 rayleigh_quotient :: (RealField.C a,
311 rayleigh_quotient matrix vector =
312 (vector `dot` (matrix * vector)) / (norm_squared vector)
314 -- We don't use the norm function here to avoid the algebraic
315 -- requirement on our field.
316 norm_squared v = ((transpose v) * v) !!! (0,0)