4
$\begingroup$

What are the best FGMRES implementations in various languages/frameworks? In particular, are there any good quality Matlab implementations?

I am referring to the variation of GMRES where a changing or nonlinear preconditioner is allowed by keeping track of the preconditioned vectors as well, as discussed in Saad's paper,

Saad, Youcef. "A flexible inner-outer preconditioned GMRES algorithm." SIAM Journal on Scientific Computing 14.2 (1993): 461-469. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.15.4854&rep=rep1&type=pdf

It looks pretty simple to implement, but I know there can sometimes be numerical problems with Krylov methods if implemented naively (eg., loss of orthogonality in CG).


Edit: Here's my unsophisticated MATLAB implementation that seems to work fine, in case anyone reading this thread needs a matlab version that works but isn't highly optimized.

function [x,X] = fgmres(A,b,m,tol,maxit,P,x0)
%Flexible GMRES (FGMRES) with restart length m.
%Solves Ax = b to tolerance tol using preconditioner P and initial guess x0
%If P is a function handle, P(A) =approx= I.
%If P is a matrix, P\A =approx= I.
%Optional output X is the whole sequence of FGMRES iterates

%Written by Nick Alger 6/10/2014, released to public domain.
%Reference:
%Saad, "A flexible inner-outer preconditioned GMRES algorithm", SIAM 1993

%Parse input
N = size(b,1);

if (isa(A,'numeric'))
    if (size(A,1) == N && size(A,2) == N)
        Afct = @(u) A*u;
    else
        error(['Dimensions of A and b in Ax = b are not consistent: '...
            'A is ', num2str(size(A,1)),'-by-', num2str(size(A,2)),', and '...
            'b is ', num2str(size(b,1)),'-by-', num2str(size(b,2)),'.']);
    end
elseif (isa(A,'function_handle'))
    Afct = A;
else
    error('Forward operator needs to be a matrix or function handle.')
end

if (nargin < 3 || isempty(m))
    m = 50; %Default restart
end

if (nargin < 4 || isempty(tol))
    tol = 1e-6;    
end

if (nargin < 5 || isempty(maxit))
    maxit = min(N,200);
end

if (nargin < 6 || isempty(P))
    Pfct = @(u) u;
elseif (isa(P,'numeric'))
    Pfct = @(u) P\u;
elseif (isa(P,'function_handle'))
    Pfct = P;
else
    error('Preconditioner needs to be a matrix or function handle.')
end

if (nargin < 7 || isempty(x0))
    x0 = zeros(size(b));
end

%Run FGMRES
X = [];
nit = 1;
while (nit <= maxit)
    disp('Starting new FGMRES cycle')
    H = zeros(m+1,m);
    r0 = b - Afct(x0);
    beta = norm(r0);
    V = zeros(N, m+1);
    V(:,1) = (1/beta)*r0;
    Z = zeros(N,m);
    for j=1:m
        Z(:,j) = Pfct(V(:,j));
        w = Afct(Z(:,j));
        for i=1:j
            H(i,j) = w'*V(:,i);
            w = w - H(i,j)*V(:,i);
        end
        H(j+1, j) = norm(w);
        V(:, j+1) = (1/H(j+1, j))*w;

        e1 = zeros(j+1,1); e1(1) = 1;
        y = H(1 : j+1, 1:j)\(beta*e1);
        x = x0 + Z(:, 1:j)*y;
        if (nargout > 1)
            X = [X,x];
        end
        resnorm = norm(b-Afct(x));
        disp(['k= ', num2str(nit), ...
            ', relative residual= ', num2str(resnorm/norm(b))])
        if (resnorm < tol*norm(b))
            disp(['FGMRES converged to relative tolerance ', ...
                num2str(resnorm/norm(b)), ...
                ' at iteration ', ...
                num2str(nit)])
            return
        end

        nit = nit + 1;
        if (nit > maxit)
            break;
        end
    end
    x0 = x;
end
$\endgroup$
4
$\begingroup$

Both PETSc and Trilinos have high quality implementations. Both of these libraries also have bindings to other languages.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.