From: Michael Orlitzky Date: Fri, 22 Mar 2013 07:48:56 +0000 (-0400) Subject: Test our (P)CGM implementation against Octave's. X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=commitdiff_plain;h=e1b71b4ca7cfa08ac76744a17a3778d4ccfaa7e2 Test our (P)CGM implementation against Octave's. --- diff --git a/tests/conjugate_gradient_method_tests.m b/tests/conjugate_gradient_method_tests.m index 59e137a..e04b2e2 100644 --- a/tests/conjugate_gradient_method_tests.m +++ b/tests/conjugate_gradient_method_tests.m @@ -14,3 +14,23 @@ diff = norm(actual - expected); unit_test_equals("CGM works on an example", ... true, ... norm(diff) < 1e-6); + + +# Let's test Octave's pcg() against our method on some easy matrices. +max_iterations = 100000; +tolerance = 1e-11; + +for n = [ 5, 10, 25, 50, 100 ] + A = random_positive_definite_matrix(5, 1000); + + # Assumed by Octave's implementation when you don't supply a + # preconditioner. + x0 = zeros(5, 1); + b = unifrnd(-1000, 1000, 5, 1); + [o_x, o_flag, o_relres, o_iter] = pcg(A, b, tolerance, max_iterations); + [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations); + + diff = norm(o_x - x); + msg = sprintf("Our CGM agrees with Octave's, n=%d.", n); + unit_test_equals(msg, true, norm(diff) < 1e-10); +end diff --git a/tests/preconditioned_conjugate_gradient_method_tests.m b/tests/preconditioned_conjugate_gradient_method_tests.m index d44ee42..2c8ff84 100644 --- a/tests/preconditioned_conjugate_gradient_method_tests.m +++ b/tests/preconditioned_conjugate_gradient_method_tests.m @@ -40,3 +40,29 @@ diff = norm(pcgm_simple - pcgm); unit_test_equals("PCGM agrees with Simple PCGM when M != I", ... true, ... norm(diff) < 1e-6); + + +# Test again Octave's pcg() function. +max_iterations = 100000; +tolerance = 1e-11; +C = random_positive_definite_matrix(5, 1000); +M = C*C'; + +for n = [ 5, 10, 25, 50, 100 ] + A = random_positive_definite_matrix(5, 1000); + + # Assumed by Octave's implementation when you don't supply a + # preconditioner. + x0 = zeros(5, 1); + b = unifrnd(-1000, 1000, 5, 1); + [o_x, o_flag, o_relres, o_iter] = pcg(A, b, tolerance, max_iterations, C, C'); + [x, k] = preconditioned_conjugate_gradient_method(A, + M, + b, + x0, + tolerance, + max_iterations); + diff = norm(o_x - x); + msg = sprintf("Our PCGM agrees with Octave's, n=%d.", n); + unit_test_equals(msg, true, norm(diff) < 1e-10); +end