1 {-# LANGUAGE RebindableSyntax #-}
2 {-# LANGUAGE ScopedTypeVariables #-}
3 {-# LANGUAGE TypeFamilies #-}
10 import Data.Vector.Fixed ( Arity, N1 )
11 import NumericPrelude hiding ( (*), abs )
12 import qualified NumericPrelude as NP ( (*) )
13 import qualified Algebra.Field as Field ( C )
15 import Linear.Matrix ( Mat(..), (!!!), construct, transpose )
18 -- | Solve the system m' * x = b', where m' is upper-triangular. Will
19 -- probably crash if m' is non-singular. The result is the vector x.
23 -- >>> import Linear.Matrix ( Mat2, Mat3, fromList, vec2d, vec3d )
25 -- >>> let identity = fromList [[1,0,0],[0,1,0],[0,0,1]] :: Mat3 Double
26 -- >>> let b = vec3d (1, 2, 3::Double)
27 -- >>> forward_substitute identity b
28 -- ((1.0),(2.0),(3.0))
29 -- >>> (forward_substitute identity b) == b
32 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Double
33 -- >>> let b = vec2d (1, 1::Double)
34 -- >>> forward_substitute m b
37 -- >>> let m = fromList [[4,0],[0,2]] :: Mat2 Double
38 -- >>> let b = vec2d (2, 1.5 :: Double)
39 -- >>> forward_substitute m b
42 forward_substitute :: forall a m. (Field.C a, Arity m)
46 forward_substitute m' b' = x'
50 -- Convenient accessor for the elements of b'.
54 -- Convenient accessor for the elements of m'.
58 -- Convenient accessor for the elements of x'.
62 -- The second argument to lambda should always be zero here, so we
64 lambda :: Int -> Int -> a
65 lambda 0 _ = (b 0) / (m 0 0)
66 lambda k _ = ((b k) - sum [ (m k j) NP.* (x j) |
67 j <- [0..k-1] ]) / (m k k)
70 -- | Solve the system m*x = b, where m is lower-triangular. Will
71 -- probably crash if m is non-singular. The result is the vector x.
75 -- >>> import Linear.Matrix ( Mat3, fromList, vec3d )
77 -- >>> let identity = fromList [[1,0,0],[0,1,0],[0,0,1]] :: Mat3 Double
78 -- >>> let b = vec3d (1, 2, 3::Double)
79 -- >>> backward_substitute identity b
80 -- ((1.0),(2.0),(3.0))
81 -- >>> (backward_substitute identity b) == b
84 backward_substitute :: (Field.C a, Arity m)
88 backward_substitute m =
89 forward_substitute (transpose m)
92 -- | Solve the linear system m*x = b where m is positive definite.
94 solve_positive_definite :: Mat v w a -> Mat w z a
95 solve_positive_definite m b = x
98 -- First we solve r^T * y == b for y. Then let y=r*x
99 rx = forward_substitute (transpose r) b
100 -- Now solve r*x == b.