src/Linear/System.hs

   1 {-# LANGUAGE RebindableSyntax #-}
   2 {-# LANGUAGE ScopedTypeVariables #-}
   3 {-# LANGUAGE TypeFamilies #-}
   4
   5 module Linear.System
   6 where
   7
   8 import Data.Vector.Fixed (Arity, N1)
   9
  10 import Linear.Matrix
  11
  12 import NumericPrelude hiding ((*), abs)
  13 import qualified NumericPrelude as NP ((*))
  14 import qualified Algebra.Field as Field
  15
  16
  17 -- | Solve the system m' * x = b', where m' is upper-triangular. Will
  18 --   probably crash if m' is non-singular. The result is the vector x.
  19 --
  20 --   Examples:
  21 --
  22 --   >>> let identity = fromList [[1,0,0],[0,1,0],[0,0,1]] :: Mat3 Double
  23 --   >>> let b = vec3d (1, 2, 3::Double)
  24 --   >>> forward_substitute identity b
  25 --   ((1.0),(2.0),(3.0))
  26 --   >>> (forward_substitute identity b) == b
  27 --   True
  28 --
  29 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Double
  30 --   >>> let b = vec2d (1, 1::Double)
  31 --   >>> forward_substitute m b
  32 --   ((1.0),(0.0))
  33 --
  34 --   >>> let m = fromList [[4,0],[0,2]] :: Mat2 Double
  35 --   >>> let b = vec2d (2, 1.5 :: Double)
  36 --   >>> forward_substitute m b
  37 --   ((0.5),(0.75))
  38 --
  39 forward_substitute :: forall a m. (Field.C a, Arity m)
  40                    => Mat m m a
  41                    -> Mat m N1 a
  42                    -> Mat m N1 a
  43 forward_substitute m' b' = x'
  44   where
  45     x' = construct lambda
  46
  47     -- Convenient accessor for the elements of b'.
  48     b :: Int -> a
  49     b k = b' !!! (k, 0)
  50
  51     -- Convenient accessor for the elements of m'.
  52     m :: Int -> Int -> a
  53     m i j = m' !!! (i, j)
  54
  55     -- Convenient accessor for the elements of x'.
  56     x :: Int -> a
  57     x k = x' !!! (k, 0)
  58
  59     -- The second argument to lambda should always be zero here, so we
  60     -- ignore it.
  61     lambda :: Int -> Int -> a
  62     lambda 0 _ = (b 0) / (m 0 0)
  63     lambda k _ = ((b k) - sum [ (m k j) NP.* (x j) |
  64                                   j <- [0..k-1] ]) / (m k k)
  65
  66
  67 -- | Solve the system m*x = b, where m is lower-triangular. Will
  68 --   probably crash if m is non-singular. The result is the vector x.
  69 --
  70 --   Examples:
  71 --
  72 --   >>> let identity = fromList [[1,0,0],[0,1,0],[0,0,1]] :: Mat3 Double
  73 --   >>> let b = vec3d (1, 2, 3::Double)
  74 --   >>> backward_substitute identity b
  75 --   ((1.0),(2.0),(3.0))
  76 --   >>> (backward_substitute identity b) == b
  77 --   True
  78 --
  79 backward_substitute :: (Field.C a, Arity m)
  80                     => Mat m m a
  81                     -> Mat m N1 a
  82                     -> Mat m N1 a
  83 backward_substitute m =
  84   forward_substitute (transpose m)
  85
  86
  87 -- | Solve the linear system m*x = b where m is positive definite.
  88 {-
  89 solve_positive_definite :: Mat v w a -> Mat w z a
  90 solve_positive_definite m b = x
  91   where
  92     r = cholesky m
  93     -- First we solve r^T * y == b for y. Then let y=r*x
  94     rx = forward_substitute (transpose r) b
  95     -- Now solve r*x == b.
  96 -}