tests/preconditioned_conjugate_gradient_method_tests.m

   1 ## Used throughout. The PCGM uses the infinity norm as the stopping
   2 ## condition, so we had better also.
   3 max_iterations = 100000;
   4 tolerance = 1e-11;
   5
   6 ## First a simple example.
   7 A = [5,1,2; ...
   8      1,6,3; ...
   9      2,3,7];
  10
  11 M = eye(3);
  12 b = [1;2;3];
  13 x0 = [1;1;1];
  14
  15 cgm  = conjugate_gradient_method(A, b, x0, tolerance, max_iterations);
  16 pcgm = preconditioned_conjugate_gradient_method(A, ...
  17                                                 M, ...
  18                                                 b, ...
  19                                                 x0, ...
  20                                                 tolerance, ...
  21                                                 max_iterations);
  22 diff = norm(cgm - pcgm, 'inf');
  23
  24 unit_test_equals("PCGM agrees with CGM when M == I", ...
  25                  true, ...
  26                  diff < 2*tolerance);
  27
  28 pcgm_simple = simple_preconditioned_cgm(A, M, b, x0, tolerance, max_iterations);
  29 diff = norm(pcgm_simple - pcgm, 'inf');
  30
  31 unit_test_equals("PCGM agrees with SimplePCGM when M == I", ...
  32                  true, ...
  33                  diff < 2*tolerance);
  34
  35 ## Needs to be symmetric!
  36 M = [0.97466, 0.24345, 0.54850; ...
  37      0.24345, 0.73251, 0.76639; ...
  38      0.54850, 0.76639, 1.47581];
  39
  40 pcgm = preconditioned_conjugate_gradient_method(A, ...
  41                                                 M, ...
  42                                                 b, ...
  43                                                 x0, ...
  44                                                 tolerance, ...
  45                                                 max_iterations);
  46 diff = norm(cgm - pcgm, 'inf');
  47
  48 unit_test_equals("PCGM agrees with CGM when M != I", ...
  49                  true, ...
  50                  diff < 2*tolerance);
  51
  52
  53 pcgm_simple = simple_preconditioned_cgm(A, M, b, x0, tolerance, max_iterations);
  54 diff = norm(pcgm_simple - pcgm, 'inf');
  55
  56 unit_test_equals("PCGM agrees with Simple PCGM when M != I", ...
  57                  true, ...
  58                  diff < 2*tolerance);
  59
  60
  61 # Test again Octave's pcg() function.
  62 for n = [ 5, 10, 25, 50, 100 ]
  63   A = random_positive_definite_matrix(5, 1000);
  64   C = random_positive_definite_matrix(5, 1000);
  65   M = C*C';
  66
  67   # Assumed by Octave's implementation when you don't supply a
  68   # preconditioner.
  69   x0 = zeros(5, 1);
  70   b  = unifrnd(-1000, 1000, 5, 1);
  71   [o_x, o_flag, o_relres, o_iter] = pcg(A, b, tolerance, max_iterations, C, C');
  72   [x, k] = preconditioned_conjugate_gradient_method(A,
  73                                                     M,
  74                                                     b,
  75                                                     x0,
  76                                                     tolerance,
  77                                                     max_iterations);
  78   diff = norm(o_x - x, 'inf');
  79   msg = sprintf("Our PCGM agrees with Octave's, n=%d.", n);
  80   unit_test_equals(msg, true, diff < 2*tolerance);
  81 end