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 (Dim, N1, Vector)
   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 import Debug.Trace (trace, traceShow)
  17
  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.
  20 --
  21 --   Examples:
  22 --
  23 --   >>> let identity = fromList [[1,0,0],[0,1,0],[0,0,1]] :: Mat3 Double
  24 --   >>> let b = vec3d (1,2,3)
  25 --   >>> forward_substitute identity b
  26 --   ((1.0),(2.0),(3.0))
  27 --   >>> (forward_substitute identity b) == b
  28 --   True
  29 --
  30 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Double
  31 --   >>> let b = vec2d (1,1)
  32 --   >>> forward_substitute m b
  33 --   ((1.0),(0.0))
  34 --
  35 forward_substitute :: forall a v w z.
  36                       (Show a, Field.C a,
  37                        Vector z a,
  38                        Vector w (z a),
  39                        Vector w a,
  40                        Dim z ~ N1,
  41                        v ~ w)
  42                    => Mat v w a
  43                    -> Mat w z a
  44                    -> Mat w z a
  45 forward_substitute m' b' = x'
  46   where
  47     x' = construct lambda
  48
  49     -- Convenient accessor for the elements of b'.
  50     b :: Int -> a
  51     b k = b' !!! (k, 0)
  52
  53     -- Convenient accessor for the elements of m'.
  54     m :: Int -> Int -> a
  55     m i j = m' !!! (i, j)
  56
  57     -- Convenient accessor for the elements of x'.
  58     x :: Int -> a
  59     x k = x' !!! (k, 0)
  60
  61     -- The second argument to lambda should always be zero here, so we
  62     -- ignore it.
  63     lambda :: Int -> Int -> a
  64     lambda 0 _ = (b 0) / (m 0 0)
  65     lambda k _ = ((b k) - sum [ (m k j) NP.* (x j) |
  66                                   j <- [0..k-1] ]) / (m k k)
  67
  68
  69 -- | Solve the system m*x = b, where m is lower-triangular. Will
  70 --   probably crash if m is non-singular. The result is the vector x.
  71 --
  72 --   Examples:
  73 --
  74 --   >>> let identity = fromList [[1,0,0],[0,1,0],[0,0,1]] :: Mat3 Double
  75 --   >>> let b = vec3d (1,2,3)
  76 --   >>> backward_substitute identity b
  77 --   ((1.0),(2.0),(3.0))
  78 --   >>> (backward_substitute identity b) == b
  79 --   True
  80 --
  81 backward_substitute :: (Show a, Field.C a,
  82                         Vector z a,
  83                         Vector v (w a),
  84                         Vector w (z a),
  85                         Vector w a,
  86                         Dim z ~ N1,
  87                         v ~ w)
  88                     => Mat v w a
  89                     -> Mat w z a
  90                     -> Mat w z a
  91 backward_substitute m b =
  92   forward_substitute (transpose m) b
  93
  94
  95 -- | Solve the linear system m*x = b where m is positive definite.
  96 {-
  97 solve_positive_definite :: Mat v w a -> Mat w z a
  98 solve_positive_definite m b = x
  99   where
 100     r = cholesky m
 101     -- First we solve r^T * y == b for y. Then let y=r*x
 102     rx = forward_substitute (transpose r) b
 103     -- Now solve r*x == b.
 104 -}