关于灵活的 GMRES (fgmres)，我们知道它是右预条件 gmres 的变体。matlab中的健壮命令gmres如下：
>> help gmres
gmres Generalized Minimum Residual Method.
X = gmres(A,B) attempts to solve the system of linear equations A*X = B
for X. The N-by-N coefficient matrix A must be square and the right
hand side column vector B must have length N. This uses the unrestarted
method with MIN(N,10) total iterations.
我们可以看到matlab命令gmres可以支持左右预处理gmres。如何实现fgmres使用matlab的gmres.m？
这是我的简单示例，左右预处理器成功但 fgmres 失败：
clc;clear;
n = 21;
A = rand(n);
b = sum(A,2);
tol = 1e-7;
maxit = n;
M = diag(diag(A));
x_true = A\b;% exact solution
restart = n;
% left precondition
x1 = gmres(A,b,restart,tol,maxit,M);
norm(x_true-x1)
% right precondition
x2 = gmres(@(x)A*(M\x),b,restart,tol,maxit);
norm(x_true-M\x2)
% fgmres
Mfun=@(x) minres(M,x);
x3 = gmres(@(x)A*Mfun(x),b,restart,tol,maxit);
norm(x_true-Mfun(x3))
编辑：
我写了一个fgmres.m但是当它发生故障时，它没有得到正确的解决方案，即当故障发生时，迭代步数为外循环为3，内循环为1（重启= 30），即总迭代步数为61 , 但近似解是. 能给我一些帮助吗，非常感谢。无需修改任何代码即可在matlab中运行。我的matlab是2018b，8GB内存。
clc;clear;close all;
restart = 30;
maxit = 100;
tol = 1e-6;
%%
fprintf('----------------------- fgmres with inexact inner solves -----------\n');
mu = 1;q =64;
fprintf('------------------Grid = %4d, mu = %6.4f----------\n',q,mu);
fprintf('flag\t\t|\t\titer\t\t|\t\tcputime\t\t|\t\trelres\t\t|\t\t|x-x_m|_2\n');
alpha = mu;
%% generate the saddle point matrix : bigA*x = rhs
h = 1/(1+q);
n = 2*q^2;m = q^2;
N = m+n;
I = speye(q);
T = spdiags(ones(q,1).*[-1 2 -1],[-1 0 1],q,q)*mu/h^2;
F = spdiags(ones(q,1).*[-1 1 0],[-1 0 1],q,q)/h;
B = [kron(I,F);kron(F,I)]';
A = kron(I,T)+kron(T,I);
A = blkdiag(A,A);
bigA = [A, B';-B,sparse(m,m)];
x_true = ones(N,1);
rhs = bigA*x_true;x0 = zeros(N,1);
fprintf('------------------------ my fgmres --------------------\n');
%% Hss
tic;
M = @(x)hss_precd_inexact(alpha,A,B,x);% a function handle returns M_j\x
[x,flag,relres,iter,resvec]=myfgmres_right(bigA,rhs,restart,tol,maxit,M);
t=toc;
iter = (iter(1)-1)*restart+iter(2);
err = norm(x_true-x);
fprintf('%4d%19d%25.4f%20.4e%22.4e\n',flag,iter,t,relres,err);
%% the defined preconditioner which uses iterative method to solve the sub system
function z = hss_precd_inexact(alpha,A,B,r)
% HSS peconditioner for saddle point using iterative method for solving
% inner sub-linear systems
% 20191228
% P_hss = [alpha*In+A O ] [alpha*In B']
% [ O alpha*Im] * [-B alpha*Im]
[m,n]=size(B);
In = speye(n);
% Im = speye(m);
r1 = r(1:n,1);
r2 = r(n+1:end,1);
% L_A = ichol(alpha*In+A);
% L_B = chol(alpha*Im+1/alpha*(B*B'),'low');
[w1,~] = pcg(@afun1,r1);
w2 = 1/alpha*r2;
temp = 1/alpha*B*w1+w2;
t1 = 1/alpha*w1;
[t2,~] = pcg(@afun2,temp);
z1 = t1-1/alpha*B'*t2;
z2 = t2;
z = [z1;z2];
%% handle returns A*x
function y = afun1(x)
y =alpha*x+A*x;
end
function y = afun2(x)
y =alpha*x+1/alpha*(B*(B'*x));
end
end
%% my fgmres.m
function [x,flag,relres,iter,resvec] = myfgmres_right(A,b,restart,tol,maxit,M,x0)
% myfgmres.m generalized minimal residual to solve : A*x= b using right
% preconditioner i.e., A*inv(M) *u = b, u=M*x
% input
% A any real nonsingular matrix or function handle
% returns A*x
% b real right hand side
% restart the maximum of iteration (means dimension of Krylov)
% tol tolerance
% maxit outer iteration steps
% x0 initialized guess vector (default is zero vector)
% M right preconditioner: matrix or function handle
% returns M\x
% output
% x approxiamte solution: x_k
% flag 0 = converge, 1=unconverge
% relres relative residual
% iter the iteration steps
% resvec ||r_k||_2, r_k=b-A*x_k, res(1)=norm(b-A*x0)
%-------------------
% Initialization
%-------------------
% size of the problem
n = size(b,1);
if nargin==7
% do nothing
elseif nargin==6
x0 = zeros(n,1);
elseif nargin ==5
M=[];x0 = zeros(n,1);
elseif nargin ==4
maxit=n;M=[];x0 = zeros(n,1);
elseif nargin ==3
tol = 1e-6;maxit=n;M=[];x0 = zeros(n,1);
elseif nargin == 2
restart = 10;tol = 1e-6;maxit=n;M=[];x0 = zeros(n,1);
else
error('Input variables are not enough!!!!!!!!!!!!!!!!');
end
%% restart number
if isempty(restart)% full gmres
restart = maxit;
maxit=1;
elseif restart ~= 0
restart = min(restart, n);
restart = min(restart,maxit);
elseif restart == 0
error('restart ==0 is wrong!!!!!!');
else
error('restart number is wrong!!!!!!');
end
%% initialization space
m = restart;
% n = length(A);
H = zeros(m +1,m );% the upper hessenberg matrix H (m+1,m)*****
c = zeros(m,1);% the givens transformation parameters: G1,...Gm
s = zeros(m,1);
resvec = zeros(maxit*m+1,1);% preallocate the maximum space of residual norm
flag = 1;% unconverge
Z = zeros(n,m);
V = zeros(n,m+1);
%% prepare to iteration
x = x0;
% initial residual
r = b-afun(x);% r0=b-A*x0
% r = mfun(M,r);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% left precondition: M\r
resvec(1) = norm(r);% initial residual
total_iter = 0;% total iteration steps
for out = 1:maxit
r = b-afun(x);% r0=b-A*x0
% r = mfun(M,r);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% left precondition: M\r
beta = norm(r);
e1 = zeros(m+1,1);e1(1) = 1;% e1
g = beta * e1;% beta*e1
% V = zeros(n,maxit);% orthonormal basis V = [v1,v2,...v_m]
V(:,1) = r/beta;% % v1
%% begin iteration
for j = 1: m
total_iter = total_iter+1;
Z(:,j) = mfun(V(:,j)); % right precondition
w = afun(Z(:,j));% right precondition
% modified Gram-Schmidt
for i = 1:j
H(i,j) = w.'*V(:,i);% h_ij
w = w - H(i,j) * V(:,i);% w_j = w_j - ...
end
H(j+1,j) = norm(w);% ||w||_2
%% lucky breakdown
if H(j+1,j) < eps
fprintf('lucky breakdown!!!!!!!!!!!\n');
flag = 0;
% apply the first j-1 givens to the last column of H_{j+1}_{j}
for k = 1:j-1
temp = c(k)*H(k,j)+s(k)*H(k+1,j);
H(k+1,j) = -s(k)*H(k,j)+c(k)*H(k+1,j);
H(k,j) = temp;
end
% apply the givens to the last 2 elements of H(:,j)
[s(j), c(j),r] = mygivens(H(j,j), H(j+1,j));
H(j,j) = r;
H(j+1,j) = 0;
% apply givens to the last 2 elements of g= beta*e1
% g(j:j+1,1) = [c(j) s(j);-s(j) c(j)] * [g(j);0]; %20191210
%----------------- 20191227
g(j+1) = -s(j)*g(j);
g(j) = c(j)*g(j);
%----------------- 20191227
resvec(total_iter+1) = abs(g(j+1)); % obtain norm(r_k)
relres = resvec(total_iter+1)/resvec(1);% ||r_k||/||r0||
break;
end
%% generate a new orthonomal basis
V(:,j+1) = w/H(j+1,j);% v_{j+1}
% apply the first j-1 givens to the last column of H_{j+1}_{j}
for k = 1:j-1
temp = c(k)*H(k,j)+s(k)*H(k+1,j);
H(k+1,j) = -s(k)*H(k,j)+c(k)*H(k+1,j);
H(k,j) = temp;
end
% apply the givens to the last 2 elements of H(:,j)
[s(j), c(j),r] = mygivens(H(j,j), H(j+1,j));
H(j,j) = r;
H(j+1,j) = 0;
% apply givens to the last 2 elements of g= beta*e1
% g(j:j+1,1) = [c(j) s(j);-s(j) c(j)] * [g(j);0]; %20191210
%----------------- 20191227
g(j+1) = -s(j)*g(j);
g(j) = c(j)*g(j);
%----------------- 20191227
resvec(total_iter+1) = abs(g(j+1)); % obtain norm(r_k)
relres = resvec(total_iter+1)/resvec(1);% ||r_k||/||r0||
% check convergence
if relres < tol
flag = 0;
break;
end
end% end of inner iteration
%% update the new iterate
y = H(1:j,1:j)\g(1:j);
% x = x + V(:,1:j)*y;
x = x+Z(:,1:j)*y;%--------------- right precondition
if flag==0
break;
end
end% end of outer iteration
iter = [out, j];
resvec = resvec(1:total_iter+1);
% end of gmres
%% children function
%% givens transformation
function [s,c,r] = mygivens(a,b)
% function Givens transformation: make sure r >= 0
% [c s] *[a] =[r]
% -s c] [b] =[0]
% written by Sun,Zhen-Wei on 2019.6.20
if ( a==0 && b==0 )
c=1;s=0;r=0;
return;
end
if ( a==0 && b~=0 )
c = 0;
s = sign(b);
r = abs(b);
return;
end
if ( a~=0 && b==0 )
c = sign(a);
s = 0;
r = abs(a);
return;
end
%% case for a~=0 and b~=0
if abs(b) > abs(a)
tau = a/b;
s = sign(b)/sqrt(1+tau^2);
c = s*tau;
else
tau = b/a;
c = sign(a)/sqrt(1+tau^2);
s = c*tau;
end
r = sqrt(a^2+b^2);
end
%% function handle returns A*x
function y = afun(x)
if isa(A,'double')
y = A*x;
elseif isa(A,'function_handle')
y = A(x);
else
error('------- A is neither a matrix or a function hanlde');
end
end
%% preconditioner: returns M\x
function z = mfun(x)
if isempty(M)
z = x;
elseif isa(M,'double')
z = M\x;
elseif isa(M,'function_handle')
z = M(x);
else
error('----------- Precnoditioner is neither a matrix or function handle');
end
end
end
```