Merge pull request #41 from nidtec-una/13-add-basic-tests-for-the-alg…

…orithms

BUGFIX Arnoldi + test matrices + plane rotations encapsulation

src/modified_gram_schmidt_arnoldi.m

@@ -1,4 +1,4 @@
-function [H, V] = modified_gram_schmidt_arnoldi(A, v, m)
+function [H, V, mUpdated] = modified_gram_schmidt_arnoldi(A, v, m)
     % Modified Gram-Schmidt Arnoldi iteration
     %
     %   Description:
@@ -10,28 +10,34 @@
     %   Syntaxis:
     %   ---------
     %
-    %   [H, V] = modified_gram_schmidt_arnoldi(A, v, m)
+    %   [H, V, mUpdated] = modified_gram_schmidt_arnoldi(A, v, m)
     %
     %   Input parameters:
     %   -----------------
     %
-    %   A:  n-by-n matrix
-    %       Matrix of coefficients.
+    %   A:          n-by-n matrix
+    %               Matrix of coefficients.
     %
-    %   v:  n-by-1 vector
-    %       Normalized residual vector.
+    %   v:          n-by-1 vector
+    %               Normalized residual vector.
     %
-    %   m:  int
-    %       Restart parameter.
+    %   m:          int
+    %               Restart parameter.
     %
     %   Output parameters:
     %   ------------------
     %
-    %   H:  m+1-by-m matrix
-    %       Upper Hessenberg matrix
+    %   H:          m+1-by-m matrix
+    %               Upper Hessenberg matrix.
     %
-    %   V:  n-by-m matrix
-    %       Orthonormal basis of Krylov subspace
+    %   V:          n-by-m matrix
+    %               Orthonormal basis of Krylov subspace.
+    %
+    %   mUpdated:   int
+    %               Updated value of the restart parameter 'm'. This
+    %               parameter changes only if the last element of H is exactly
+    %               0 and the matrix is trunctated. See, in particular, the
+    %               if-else statement inside the Modified Gram-Schmidt block.
     %
     %   References:
     %   -----------
@@ -85,6 +91,7 @@
             % Slice matrices and return outputs.
             H = H(1:m + 1, 1:m);
             V = V(:, 1:m);
+            mUpdated = m;
             return
         else
             V(:, j + 1) = W(:, j) / H(j + 1, j);
@@ -93,5 +100,9 @@
     end
 
     % Slice matrix V since we are only interested in the first 'm' columns
-    V = V(1:n, 1:m);
+    V = V(:, 1:m);
+
+    % m will remain the same in this case, but we need it in the output
+    mUpdated = m;
+
 end

src/pd_gmres.m

@@ -271,7 +271,7 @@
         tic();
 
         % Calll MATLAB bult-in gmres
-        [x, flag, relres, ~, resvec] = ...
+        [x, flag, ~, ~, resvec] = ...
             gmres(A, b, [], tol, maxit, [], [], xInitial);
 
         % gmres uses a flag system. We only care wheter the solution has
@@ -310,6 +310,7 @@
 
     while flag == 0
 
+        % Control block
         if iter(size(iter, 1), :) ~= 1
             [miter] = pd_rule(m, n, mInitial, mMin, mMax, ...
                               mStep, res, iter(size(iter, 1), :), ...
@@ -319,58 +320,43 @@
         else
             m = mInitial;
         end
-
         mvec(iter(size(iter, 1), :) + 1, 1) = m;
 
         % Compute normalized residual vector
         r = b - A * xInitial;
         beta = norm(r);
         v1 = r / beta;
 
-        % Modifed Gram-Schmidt Arnoldi iteration
-        [H, V] = modified_gram_schmidt_arnoldi(A, v1, m);
-
-        g = zeros(m + 1, 1);
-        g(1, 1) = beta;
-
-        % Plane rotations (QR decompostion)
-        for j = 1:m
-            P = eye(m + 1);
-            sin = H(j + 1, j) / (sqrt(H(j + 1, j)^2 + H(j, j)^2));
-            cos = H(j, j) / (sqrt(H(j + 1, j)^2 + H(j, j)^2));
-            P(j, j) = cos;
-            P(j + 1, j + 1) = cos;
-            P(j, j + 1) = sin;
-            P(j + 1, j) = -sin;
-            H = P * H;
-            g = P * g;
-        end
-
-        R = zeros(m, m);
-        G = zeros(m, 1);
+        % Modified Gram-Schmidt Arnoldi iteration
+        [H, V, m] = modified_gram_schmidt_arnoldi(A, v1, m);
 
-        for k = 1:m
-            G(k) = g(k);
-            for i = 1:m
-                R(k, i) = H(k, i);
-            end
-        end
+        % Plane rotations
+        [HUpTri, g] = plane_rotations(H, beta);
 
-        minimizer = R \ G;
+        % Solve the least-squares problem
+        Rm = HUpTri(1:m, 1:m);
+        gm = g(1:m);
+        minimizer = Rm \ gm;
         xm = xInitial + V * minimizer;
+
+        % Update residual norm, iterations, and relative residual vector
         res(restart + 1, :) = abs(g(m + 1, 1));
         iter(restart + 1, :) = restart + 1;
-        relresvec(size(relresvec, 1) + 1, :) = abs(g(m + 1, 1) / res(1, 1));
-        % Use last component of g as residual
-        if abs(g(m + 1, 1)) / res(1, 1) < tol || size(relresvec, 1) == maxit
-            flag = 1;  % solution has converged
-            x = xm;  % solution vector
+        relresvec(size(relresvec, 1) + 1, :) = res(restart + 1, :) / res(1, 1);
+
+        % Check convergence
+        if relresvec(end) < tol || size(relresvec, 1) == maxit
+            % We reached convergence.
+            flag = 1;
+            x = xm;
         else
-            xInitial = xm;  % update and restart
+            % We have not reached convergence. Update and restart.
+            xInitial = xm;
             restart = restart + 1;
         end
+
     end
 
-    time = toc();  % record CPU time
+    time = toc();
 
 end

src/plane_rotations.m

@@ -0,0 +1,92 @@
+function [HUpTri, g] = plane_rotations(H, beta)
+    % Performs plane rotations.
+    %
+    %   Description:
+    %   ------------
+    %
+    %   Implementation of plane rotations (a.k.a Givens rotation) following [1].
+    %
+    %   Syntaxis:
+    %   ---------
+    %
+    %   [HUpTri, g] = plane_rotations(H, beta)
+    %
+    %   Input parameters:
+    %   -----------------
+    %
+    %   H:          m+1-by-m matrix
+    %               Upper Hessenberg matrix.
+    %
+    %   beta:       int
+    %               Norm of the last cycle residual vector.
+    %
+    %   Output parameters:
+    %   ------------------
+    %
+    %   HUpTri:     m+1-by-m matrix
+    %               Upper-triangular Hessenberg matrix.
+    %
+    %   g:          m+1-by-1 vector
+    %               Resulting right-hand-side vector.
+    %
+    %   Notes:
+    %   ------
+    %
+    %   HUpTri and g (with the last rows deleted) will enter into the
+    %   least-squares problem.
+    %
+    %   References:
+    %   -----------
+    %
+    %   [1] Saad, Y. (2003). Iterative methods for sparse linear systems.
+    %   Society for Industrial and Applied Mathematics.
+    %
+    %   Copyright:
+    %   ----------
+    %
+    %   This file is part of the KrySBAS MATLAB Toolbox.
+    %
+    %   Copyright 2023 CC&MA - NIDTec - FP - UNA
+    %
+    %   KrySBAS is free software: you can redistribute it and/or modify it under
+    %   the terms of the GNU General Public License as published by the Free
+    %   Software Foundation, either version 3 of the License, or (at your
+    %   option) any later version.
+    %
+    %   KrySBAS is distributed in the hope that it will be useful, but WITHOUT
+    %   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    %   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    %   for more details.
+    %
+    %   You should have received a copy of the GNU General Public License along
+    %   with this file.  If not, see <http://www.gnu.org/licenses/>.
+    %
+
+    % Infer 'm' from the size of the upper Hessenber matrix H
+    [~, m] = size(H);
+
+    % Create rhs vector g
+    g = zeros(m + 1, 1);
+    g(1, 1) = beta;
+
+    % Plane rotations
+    for j = 1:m
+        % Obtain sines and cosines
+        s = H(j + 1, j) / (sqrt(H(j + 1, j)^2 + H(j, j)^2));
+        c = H(j, j) / (sqrt(H(j + 1, j)^2 + H(j, j)^2));
+
+        % Build rotation matrix
+        P = eye(m + 1);
+        P(j, j) = c;
+        P(j + 1, j + 1) = c;
+        P(j, j + 1) = s;
+        P(j + 1, j) = -s;
+
+        % Update HUpTri and g
+        H = P * H;
+        g = P * g;
+    end
+
+    HUpTri = H;
+
+end

tests/test_pd_gmres.m

@@ -292,7 +292,7 @@ function test_outputs_unrestarted_identity_matrix()
     A = eye(3);
     b = ones(3, 1);
 
-    % Call function
+    % Call PD-GMRES
     [x, flag, relresvec, mvec, time] = pd_gmres(A, b);
 
     % Compare with expected outputs
@@ -307,19 +307,18 @@ function test_outputs_restarted_identity_matrix()
     % Test whether the correct outputs are returned using an identity
     % matrix when the restarted algorithm is called
 
+    % Setup a trivial linear system
     A = eye(3);
     b = [2; 3; 4];
+
+    % Setup PD-GMRES
     mInitial = 1;
-    mMinMax = [1; 3];
-    mStep = 1;
     tol = 1e-9;
     maxit = 10;
-    xInitial = zeros(3, 1);
-    alphaPD = [-3; 5];
 
-    % Call function
+    % Call PD-GMRES
     [x, flag, relresvec, mvec, time] = ...
-        pd_gmres(A, b, mInitial, mMinMax, mStep, tol, maxit, xInitial, alphaPD);
+        pd_gmres(A, b, mInitial, [], [], tol, maxit, [], []);
 
     % Compare with expected outputs
     assertElementsAlmostEqual(x, [2; 3; 4]);
@@ -329,3 +328,97 @@ function test_outputs_restarted_identity_matrix()
     assert(time > 0 && time < 5);
 
 end
+
+function test_embree_three_by_three_toy_example()
+    % Test Embree's 3x3 linear system from https://www.jstor.org/stable/25054403
+    % with PD-GMRES. We construct two checks, one with mInitial = 1, and
+    % one with mInitial = 2. Note that the gmres(m) algorithm with m = 2,
+    % does not converge (this was proven by Embree in his paper). This, however,
+    % does not occur with the PD-GMRES, which adaptively changes the value
+    % of m to avoid stagnation.
+
+    % Load A and b
+    load('data/embree3.mat', 'Problem');
+    A = Problem.A;
+    b = Problem.b;
+
+    % Setup PD-GMRES
+    mInitials = [1; 2];
+    tol = 1e-9;
+    maxit = 20;
+
+    % Loop over the values of mInitials, i.e., {1, 2}.
+    for i = 1:length(mInitials)
+        % Call PD-GMRES
+        [x, flag, relresvec, mvec, ~] = ...
+            pd_gmres(A, b, mInitials(i), [], [], tol, maxit, [], []);
+
+        % Compare with expected outputs
+        assertElementsAlmostEqual(x, [8; -7; 1]);
+        assertEqual(flag, 1);
+        if mInitials(i) == 1
+            relresvecExpected = [1.000000000000000
+                                 0.925820099772551
+                                 0.654653670707977
+                                 0];
+            mvecExpected = [1; 1; 1; 2];
+        else
+            relresvecExpected = [1.000000000000000
+                                 0.462910049886276
+                                 0.377189160453301
+                                 0.000000000000001];
+            mvecExpected = [2; 2; 2; 3];
+        end
+        assertElementsAlmostEqual(relresvec, relresvecExpected);
+        assertEqual(mvec, mvecExpected);
+    end
+
+end
+
+function test_sherman_one()
+    % Test pd_gmres with sherman1 matrix
+
+    % Load A and b
+    load('data/sherman1.mat', 'Problem');
+    A = Problem.A;
+    b = Problem.b;
+
+    % Setup PD-GMRES
+    mInitial = 30;
+    tol = 1e-9;
+    mStep = 3;
+    maxit = 1000;
+
+    % Call PD-GMRES
+    [~, flag, ~, mvec, ~] = ...
+            pd_gmres(A, b, mInitial, [], mStep, tol, maxit, [], []);
+
+    % We check if it has converged and the total sum of outer iterations
+    assertEqual(flag, 1);
+    assertEqual(sum(mvec), 955);
+
+end
+
+function test_sherman_four()
+    % Test pd_gmres with sherman4 matrix
+
+    % Load A and b
+    load('data/sherman4.mat', 'Problem');
+    A = Problem.A;
+    b = Problem.b;
+
+    % Setup PD-GMRES
+    mInitial = 30;
+    tol = 1e-9;
+    mStep = 3;
+    maxit = 1000;
+
+    % Call PD-GMRES
+    [~, flag, ~, mvec, ~] = ...
+            pd_gmres(A, b, mInitial, [], mStep, tol, maxit, [], []);
+
+    % We check if it has converged and the total sum of outer iterations
+    assertEqual(flag, 1);
+    assertEqual(sum(mvec), 440);
+
+end