Add tridiagonal solver and tests for it and triagonal solver.

1 files changed, 107 insertions, 0 deletions

diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 043ebf3..55a5f74 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -131,6 +131,111 @@ mbCholTest = utest "mbCholTest" (ok1 && ok2) where
     ok1 = mbChol (trustSym m1) == Nothing
     ok2 = mbChol (trustSym m2) == Just (chol $ trustSym m2)
 
+-----------------------------------------------------
+
+triTest = utest "triTest" ok1 where
+
+  a :: Matrix R
+  a = (4><4)
+    [
+       4.30,  0.00,  0.00, 0.00,
+      -3.96, -4.87,  0.00, 0.00,
+       0.40,  0.31, -8.02, 0.00,
+      -0.27,  0.07, -5.95, 0.12
+    ]
+
+  w :: Matrix R
+  w = (4><2)
+    [
+      -12.90, -21.50,
+       16.75,  14.93,
+      -17.55,   6.33,
+      -11.04,   8.09
+    ]
+
+  v :: Matrix R
+  v = triSolve Lower a w
+
+  e :: Matrix R
+  e = (4><2)
+    [
+      -3.0000, -5.0000,
+      -1.0000,  1.0000,
+       2.0000, -1.0000,
+       1.0000,  6.0000
+    ]
+
+  ok1 = (maximum $ map abs $ concat $ toLists $ e - v) <= 1e-14
+
+-----------------------------------------------------
+
+triDiagTest = utest "triDiagTest" (ok1 && ok2) where
+
+  dL, d, dU :: Vector Double
+  dL =  fromList [3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0]
+  d  =  fromList [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
+  dU =  fromList [4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0]
+
+  b :: Matrix R
+  b = (9><3)
+    [
+      1.0,   1.0,   1.0,
+      1.0,  -1.0,   2.0,
+      1.0,   1.0,   3.0,
+      1.0,  -1.0,   4.0,
+      1.0,   1.0,   5.0,
+      1.0,  -1.0,   6.0,
+      1.0,   1.0,   7.0,
+      1.0,  -1.0,   8.0,
+      1.0,   1.0,   9.0
+    ]
+
+  y :: Matrix R
+  y = (9><9)
+    [
+      1.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+      3.0, 1.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+      0.0, 3.0, 1.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+      0.0, 0.0, 3.0, 1.0, 4.0, 0.0, 0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0, 3.0, 1.0, 4.0, 0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0, 0.0, 3.0, 1.0, 4.0, 0.0, 0.0,
+      0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 1.0, 4.0, 0.0,
+      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 1.0, 4.0,
+      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 1.0
+    ]
+
+  x :: Matrix R
+  x = triDiagSolve dL d dU b
+
+  z :: Matrix C
+  z = (4><4)
+    [
+      1.0 :+ 1.0, 4.0 :+ 4.0, 0.0 :+ 0.0, 0.0 :+ 0.0,
+      3.0 :+ 3.0, 1.0 :+ 1.0, 4.0 :+ 4.0, 0.0 :+ 0.0,
+      0.0 :+ 0.0, 3.0 :+ 3.0, 1.0 :+ 1.0, 4.0 :+ 4.0,
+      0.0 :+ 0.0, 0.0 :+ 0.0, 3.0 :+ 3.0, 1.0 :+ 1.0
+    ]
+
+  zDL, zD, zDu :: Vector C
+  zDL = fromList [3.0 :+ 3.0, 3.0 :+ 3.0, 3.0 :+ 3.0]
+  zD  = fromList [1.0 :+ 1.0, 1.0 :+ 1.0, 1.0 :+ 1.0, 1.0 :+ 1.0]
+  zDu = fromList [4.0 :+ 4.0, 4.0 :+ 4.0, 4.0 :+ 4.0]
+
+  zB :: Matrix C
+  zB = (4><3)
+    [
+      1.0 :+ 1.0,   1.0  :+   1.0,  1.0 :+ (-1.0),
+      1.0 :+ 1.0, (-1.0) :+ (-1.0), 1.0 :+ (-1.0),
+      1.0 :+ 1.0,   1.0  :+   1.0,  1.0 :+ (-1.0),
+      1.0 :+ 1.0, (-1.0) :+ (-1.0), 1.0 :+ (-1.0)
+    ]
+
+  u :: Matrix C
+  u = triDiagSolve zDL zD zDu zB
+
+  ok1 = (maximum $ map abs $ concat $ toLists $ b - (y <> x)) <= 1e-15
+  ok2 = (maximum $ map magnitude $ concat $ toLists $ zB - (z <> u)) <= 1e-15
+
 ---------------------------------------------------------------------
 
 randomTestGaussian = (unSym c) :~3~: unSym (snd (meanCov dat))
@@ -715,6 +820,8 @@ runTests n = do
                        && rank ((2><3)[1,0,0,1,7*peps,0::Double]) == 2
         , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
         , mbCholTest
+        , triTest
+        , triDiagTest
         , utest "offset" offsetTest
         , normsVTest
         , normsMTest