Comments (9)
from mmmfe-st-dd.
The local gmres code is a matrix free adaptation of the following matlab implementation which can be found here.
function [x, e] = gmres( A, b, x, max_iterations, threshold)
n = length(A);
m = max_iterations;
% use x as the initial vector
r = b - A * x;
b_norm = norm(b);
error = norm(r) / b_norm;
% initialize the 1D vectors
sn = zeros(m, 1);
cs = zeros(m, 1);
%e1 = zeros(n, 1);
e1 = zeros(m+1, 1);
e1(1) = 1;
e = [error];
r_norm = norm(r);
Q(:,1) = r / r_norm;
beta = r_norm * e1; %Note: this is not the beta scalar in section "The method" above but the beta scalar multiplied by e1
for k = 1:m
% run arnoldi
[H(1:k+1, k) Q(:, k+1)] = arnoldi(A, Q, k);
% eliminate the last element in H ith row and update the rotation matrix
[H(1:k+1, k) cs(k) sn(k)] = apply_givens_rotation(H(1:k+1,k), cs, sn, k);
% update the residual vector
beta(k + 1) = -sn(k) * beta(k);
beta(k) = cs(k) * beta(k);
error = abs(beta(k + 1)) / b_norm;
% save the error
e = [e; error];
if (error <= threshold)
break;
end
end
% if threshold is not reached, k = m at this point (and not m+1)
% calculate the result
y = H(1:k, 1:k) \ beta(1:k);
x = x + Q(:, 1:k) * y;
end
%----------------------------------------------------%
% Arnoldi Function %
%----------------------------------------------------%
function [h, q] = arnoldi(A, Q, k)
q = A*Q(:,k); % Krylov Vector
for i = 1:k % Modified Gramm-Schmidt, keeping the Hessenberg matrix
h(i) = q' * Q(:, i);
q = q - h(i) * Q(:, i);
end
h(k + 1) = norm(q);
q = q / h(k + 1);
end
%---------------------------------------------------------------------%
% Applying Givens Rotation to H col %
%---------------------------------------------------------------------%
function [h, cs_k, sn_k] = apply_givens_rotation(h, cs, sn, k)
% apply for ith column
for i = 1:k-1
temp = cs(i) * h(i) + sn(i) * h(i + 1);
h(i+1) = -sn(i) * h(i) + cs(i) * h(i + 1);
h(i) = temp;
end
% update the next sin cos values for rotation
[cs_k sn_k] = givens_rotation(h(k), h(k + 1));
% eliminate H(i + 1, i)
h(k) = cs_k * h(k) + sn_k * h(k + 1);
h(k + 1) = 0.0;
end
%%----Calculate the Given rotation matrix----%%
function [cs, sn] = givens_rotation(v1, v2)
% if (v1 == 0)
% cs = 0;
% sn = 1;
% else
t = sqrt(v1^2 + v2^2);
% cs = abs(v1) / t;
% sn = cs * v2 / v1;
cs = v1 / t; % see http://www.netlib.org/eispack/comqr.f
sn = v2 / t;
% end
end
- It could be confusing that in local_gmres, there is division by
r_norm
instead ofb_norm
, compared to the snippet above. This is because our initial guess x_0 is zero and hence the residual norm,r_norm=norm(Ax_0-b)
, is same asb_norm
. - Regarding what is being printed in the output, you are welcome to change it and send me a pull request, I will merge the changes to the master.
Let me know if you have any other questions before closing this issue.
from mmmfe-st-dd.
Also, if there is a possibility that division of error by r_norm
will be confusing some time in the future, we could simply add a variable b_norm which will hold the place for error in first iteration. Of course this will be redundant, but will increase the readability of the algorithm.
from mmmfe-st-dd.
It's fairly clear that r_norm is actually norm(b), though you're right that it could get confusing later on. You're probably better off by erring on the side of clarity vs redundancy.
But the main issue is not the printing, but the test:
in the matlab code the test is
if (error <= threshold)
whereas the C++ code has
if (combined_error_iter/e_all_iter[0] < tolerance)
but the division by e_all_iter[0]
has already taken place (e_all_iter
is the same thing as error
, and both actually contain residual reduction )
I might still be mistaken, so I want to make sure we both agree before I make any changes.
(I also removed a comment above that was for the wrong issue)
from mmmfe-st-dd.
I now see what you are pointing to. Sorry I overlooked it the first time. There indeed are two divisions by b_norm
.
This seems be a bug, which might affect the convergence and number of iterations. Please go ahead and make the changes. Thanks.
from mmmfe-st-dd.
This issue has been resolved, using a pull request from @MichelKern .
from mmmfe-st-dd.
from mmmfe-st-dd.
I have looked at the changes carefully before merging the request. I wasn't sure I was supposed to open a review formally before merging. May be I will do that from next time. Thanks for pointing out.
from mmmfe-st-dd.
from mmmfe-st-dd.
Related Issues (4)
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from mmmfe-st-dd.