各种语言/框架中最好的 FGMRES 实现是什么？特别是，是否有任何质量好的 Matlab 实现？

我指的是 GMRES 的变化，其中通过跟踪预条件向量来允许变化或非线性预条件子，正如 Saad 的论文中所讨论的，

萨德，尤瑟夫。“一种灵活的内外预处理 GMRES 算法。” SIAM 科学计算杂志 14.2 (1993): 461-469。http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.15.4854&rep=rep1&type=pdf

它看起来很容易实现，但我知道如果简单地实现 Krylov 方法有时会出现数值问题（例如，CG 中的正交性丢失）。

编辑： 这是我的简单的 MATLAB 实现，它似乎工作正常，以防阅读这个线程的任何人都需要一个可以工作但没有高度优化的 matlab 版本。

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