From 49d718705d205d62aea2762445f95735a671f305 Mon Sep 17 00:00:00 2001 From: Dominic Steinitz Date: Tue, 21 Mar 2017 17:35:43 +0000 Subject: Add tridiagonal solver and tests for it and triagonal solver. --- packages/tests/src/Numeric/LinearAlgebra/Tests.hs | 107 ++++++++++++++++++++++ 1 file changed, 107 insertions(+) (limited to 'packages/tests/src/Numeric') 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 -- cgit v1.2.3