Skip to content

Commit

Permalink
Merge pull request #72 from nidtec-una/70-encapsulation-of-harmonic-r…
Browse files Browse the repository at this point in the history 
…itz-vectors-method-gusespinola

70 encapsulation of harmonic ritz vectors method
  • Loading branch information
@gusespinola
gusespinola authored Sep 11, 2024
2 parents b636bbe + 98acba6 commit 47b8a77
Show file tree
Hide file tree
Showing 4 changed files with 157 additions and 100 deletions.
106 changes: 8 additions & 98 deletions src/gmres_e.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,24 +233,27 @@
iter(1, :) = restart;
beta = norm(r0);
v1 = r0 / beta;
s = m; % s = m + k

tic(); % start measuring CPU time

% Modified Gram-Schmidt Arnoldi iteration
[H, V, ~] = modified_gram_schmidt_arnoldi(A, v1, s);
% This is the first run. Since we don't have
% harmonic Ritz vectors yet, we use GMRES(m + k).
s = m + k;
[H, V, sUpdated] = modified_gram_schmidt_arnoldi(A, v1, s);

% Plane rotations
[HUpTri, g] = plane_rotations(H, beta);

% Solve the least-squares problem
s = sUpdated;
Rs = HUpTri(1:s, 1:s);
gs = g(1:s);
minimizer = Rs \ gs;
xm = xInitial + V * minimizer;

% Update residual norm, iterations, and relative residual vector
res(restart + 1, :) = abs(g(m + 1, 1));
res(restart + 1, :) = abs(g(s + 1, 1));
iter(restart + 1, :) = restart + 1;
relresvec(restart + 1, :) = res(restart + 1, :) / res(1, 1);

Expand All @@ -267,53 +270,7 @@
% Eigenvalue problem setup, from [1], p. 1161, step 5
Fold = H(1:s, 1:s)';
G = Rs' * Rs;
opts.tol = tol;
opts.v0 = ones(size(Fold, 1), 1);
E = zeros(s, k);
D = zeros(k, 1);

% Compute the harmonic Ritz vectors
[E2, D2] = eigs(Fold, G, k, 'LM', opts);
for p = 1:k
D(p, 1) = abs(D2(p, p));
end
[~, I] = sort(D, 1);
for q = 1:k
E(:, q) = E2(:, I(q, 1));
end
dy0 = V * E; % yi = Q * gi???

dy = zeros(n, 1);

% If dy0 has complex components, we separate each eigenvector
% into real and complex parts and treat as two distinct vectors
if isreal(dy0) == 0
ij = 1;
jj = 0;
while size(dy, 2) <= k && ij <= k
if isreal(dy0(:, ij)) == 0 && norm(real(dy0(:, ij))) > 0
dy(:, jj + 1) = real(dy0(:, ij));
jj = size(dy, 2);
if ij <= k
dy(:, jj + 1) = abs(imag(dy0(:, ij)) * sqrt(1));
jj = size(dy, 2);
if ij < k
if dy0(:, ij) == conj(dy0(:, ij + 1))
ij = ij + 2;
else
ij = ij + 1;
end
end
end
else
dy(:, jj + 1) = dy0(:, ij);
ij = ij + 1;
jj = size(dy, 2);
end
end
else
dy = dy0;
end
dy = harmonic_ritz_vectors(Fold, G, k, V, tol);

% Update and restart.
restart = restart + 1;
Expand Down Expand Up @@ -366,59 +323,12 @@
return

elseif restart < maxit

% We have not reached convergence.
% Eigenvalue problem setup, from [1], p. 1161, step 5
W = V(1:n, 1:s);
Fold = W' * A' * W;
G = Rs' * Rs;
opts.tol = tol;
opts.v0 = ones(size(Fold, 1), 1);
E = zeros(s, k);
D = zeros(k, 1);

% Compute the harmonic Ritz vectors
[E2, D2] = eigs(Fold, G, k, 'LM', opts);
for p = 1:k
D(p, 1) = abs(D2(p, p));
end
[~, I] = sort(D, 1);
for q = 1:k
E(:, q) = E2(:, I(q, 1));
end
dy0 = V * E;

dy = [];

% If dy0 has complex components, we separate each eigenvector
% into real and complex parts and treat as two distinct vectors
if isreal(dy0) == 0
ij = 1;
jj = 0;
while size(dy, 2) <= k && ij <= k
if isreal(dy0(:, ij)) == 0 && norm(real(dy0(:, ij))) > 0
dy(:, jj + 1) = real(dy0(:, ij));
jj = size(dy, 2);
if ij <= k
dy(:, jj + 1) = abs(imag(dy0(:, ij)) * sqrt(1));
jj = size(dy, 2);
if ij < k
if dy0(:, ij) == conj(dy0(:, ij + 1))
ij = ij + 2;
else
ij = ij + 1;
end
end
end
else
dy(:, jj + 1) = dy0(:, ij);
ij = ij + 1;
jj = size(dy, 2);
end
end
else
dy = dy0;
end
dy = harmonic_ritz_vectors(Fold, G, k, V, tol);
end

% Update and restart.
Expand Down
124 changes: 124 additions & 0 deletions src/harmonic_ritz_vectors.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
function dy = harmonic_ritz_vectors(F, G, k, V, tol)
% Harmonic Ritz Vectors function
%
% This function is a modified implementation of the Rayleigh-Ritz
% method that finds good approximations to the smallest eigenvalues
% and its associated eigenvectors, and appends these ones to the next
% search subspaces. Please refer to [1] for more information about
% this technique.
%
% Signature:
% ----------
%
% dy = harmonic_ritz_vectors(F, G, k, V, tol)
%
%
% Input Parameters:
% -----------------
%
% F: s-by-s maxtrix
% Matrix F from the generalized eigenvalue problem, as
% employed in step 5, p. 1161 of [1].
%
% G: s-by-1 vector
% Matrix G from the generalized eigenvalue problem, as
% employed in step 5, p. 1161 of [1].
%
% k: int
% Number of eigenvectors corresponding to a few of the
% smallest eigenvalues in magnitude. According to [1],
% "even just a few eigenvectors can make a big difference
% if the matrix has both small and large eigenvalues".
%
% V: n-by-s matrix
% A matrix from the modified Gram-Schmidt Arnoldi method,
% whose column vectors span the search subspace.
%
% tol: float
% Convergence tolerance for the built-in 'eigs' algorithm.
%
%
% Output parameters:
% ------------------
%
% dy: n-by-k matrix
% Matrix whose column vectors are the approximate
% eigenvectors of the matrix A.
%
% References:
% -----------
%
% [1] Morgan, R. B. (1995). A restarted GMRES method augmented with
% eigenvectors. SIAM Journal on Matrix Analysis and Applications,
% 16(4), 1154-1171.
%
% 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/>.
%
opts.tol = tol;
s = size(F, 1);
opts.v0 = ones(s, 1);
E = zeros(s, k);
D = zeros(k, 1);

% Compute the approximate eigenvectors
[E2, D2] = eigs(F, G, k, 'LM', opts);
for p = 1:k
D(p, 1) = abs(D2(p, p));
end
[~, I] = sort(D, 1);
for q = 1:k
E(:, q) = E2(:, I(q, 1));
end
dy0 = V * E; % Implements yi = Q * gi, step 5, p. 1161 of [1]

dy = [];

% If dy0 has complex components, its complex conjugate is also
% an approximate eigenvector, hence we separate each eigenvector
% into its real and complex parts and treat as two distinct vectors
if isreal(dy0) == 0
ij = 1;
jj = 0;
while size(dy, 2) <= k && ij <= k
if isreal(dy0(:, ij)) == 0 && norm(real(dy0(:, ij))) > 0
dy(:, jj + 1) = real(dy0(:, ij));
jj = size(dy, 2);
if ij <= k
dy(:, jj + 1) = abs(imag(dy0(:, ij)) * sqrt(1));
jj = size(dy, 2);
if ij < k
if dy0(:, ij) == conj(dy0(:, ij + 1))
ij = ij + 2;
else
ij = ij + 1;
end
end
end
else
dy(:, jj + 1) = dy0(:, ij);
ij = ij + 1;
jj = size(dy, 2);
end
end
else
dy = dy0;
end
end
1 change: 0 additions & 1 deletion test_report.log
View file

This file was deleted.

26 changes: 25 additions & 1 deletion tests/test_gmres_e.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -235,6 +235,30 @@ function test_outputs_restarted_identity_matrix() % Linear system # 2

end

function test_issue_68() % Linear system # 2
% Test whether the bugfix resolves the issue #68 or not

% Setup a trivial linear system
A = eye(3);
b = [2; 3; 4];

% Setup GMRES-E
m = 2;
k = 1;
tol = 1e-9;
maxit = 100;

% Call GMRES-E
[x, flag, relresvec, time] = gmres_e(A, b, m, k, tol, maxit, []);

% Compare with expected outputs
assertElementsAlmostEqual(x, [2; 3; 4]);
assert(flag == 1);
assertElementsAlmostEqual(relresvec, [1; 0]);
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

Expand Down Expand Up @@ -277,7 +301,7 @@ function test_sherman_one()

% We check if it has converged and the total sum of outer iterations
assertEqual(flag, 1);
assertEqual(size(relresvec, 1), 61);
assertEqual(size(relresvec, 1), 60);
assert(time > 0 && time < 100);
end

Expand Down

0 comments on commit 47b8a77

Please sign in to comment.