+-- >>> let f1 = [0.0418]
+-- >>> let f2 = [0.0805]
+-- >>> let f3 = [0.1007]
+-- >>> let f4 = [-0.0045]
+-- >>> let f5 = [-0.0332]
+-- >>> let f6 = [-0.0054]
+-- >>> let f7 = [-0.0267]
+-- >>> let big_F = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col N7 Double
+-- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
+-- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
+-- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
+-- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
+-- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
+-- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
+-- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
+-- >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
+-- >>> let r = cholesky big_K
+-- >>> let rt = transpose r
+-- >>> let e1 = [0.0170647785413895] :: [Double]
+-- >>> let e2 = [0.0338] :: [Double]
+-- >>> let e3 = [0.07408] :: [Double]
+-- >>> let e4 = [-0.00183711730708738] :: [Double]
+-- >>> let e5 = [-0.0135538432434003] :: [Double]
+-- >>> let e6 = [-0.00220454076850486] :: [Double]
+-- >>> let e7 = [-0.00689391035624920] :: [Double]
+-- >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Col N7 Double
+-- >>> let actual = forward_substitute rt big_F
+-- >>> frobenius_norm (actual - expected) < 1e-10
+-- True
+--
+forward_substitute :: forall a m. (Eq a, Field.C a, Arity m)
+ => Mat (S m) (S m) a
+ -> Col (S m) a
+ -> Col (S m) a
+forward_substitute matrix b
+ | not (is_lower_triangular matrix) =
+ error "forward substitution on non-lower-triangular matrix"
+ | otherwise = ifoldl2 f zero pairs
+ where
+ -- Pairs (m_ii, b_i) needed at each step.
+ pairs :: Col (S m) (a,a)
+ pairs = zip2 (diagonal matrix) b