High quality flexible GMRES (FGMRES) implementation

Question

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

Answer 1

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