# How to implement flexible gmres in matlab?

About the flexible GMRES (fgmres), we know that it is a variant of right preconditioned gmres. And the robust command gmres in matlab as follows:

>> 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.


We can see that the matlab command gmres can support left and right preconditioned gmres. How to implement the fgmres using matlab's gmres.m?

Here is my simple example, the left and right preconditioner are successful but the fgmres fails:

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))


edit:

I write a fgmres.m but when it occurs breakdown, it failed to get the correct solution i.e., when the breakdown occurs, the iteration steps is 3 for outer loop, 1 for inner loop (restart = 30), i.e., the total iteration steps is 61, but the approximate solution is $$x_{61} = NaN$$. can you give me some help, thanks very much. It can run in matlab without modifying any code. My matlab is 2018b, 8GB memory.

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]     =
%  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
$$$$

• I'm voting to close this question as off-topic because it's a matlab specific question that most likely should be on a matlab forum. – EMP Dec 28 '19 at 23:39
• sorry, Please do not close this. I have tried to ask on a matlab forum, but nobody care this question, and the computational science is very related to the matlab implementation, so I ask in here, and expecte some useful answers. Thanks very much. – sunshine Dec 29 '19 at 9:10
• According to the on-topic guide, this question can be interesting and can be categorized as on-topic for this community. The worst that could happen, it will be unanswered. – Anton Menshov Dec 29 '19 at 11:26
• @AntonMenshov, thanks Prof. Anton, I have written a fgmres.m function, but when I use it to solve a saddle point, it failed, I have dubugged the matlab code, but I have not found the correct way to modify it, can you take a look at my code? below is my code. Thanks, in my example, it occurs the so called 'lucky breakdown' in Gram-Schmidt process.but the approximate is not correct. – sunshine Dec 29 '19 at 13:09
• @AntonMenshov, Thanks Prof. Anton, I have written it in the edit in my question, can you take a look at it, thanks very much. – sunshine Dec 30 '19 at 2:13

First of all, MATLAB's gmres assumes that the preconditioner you use is linear. This is important! Actually it is the main difference between FGMRES and GMRES. Right preconditioned GMRES and FGMRES are exactly the same if you use a linear preconditioner, however, FGMRES allows the use of non-linear preconditioners.
What do I mean by a non-linear preconditioner? With MINRES, i.e. Mfun=@(x) minres(M,x), it is not guaranteed that Mfun(x+av) = Mfun(x)+aMfun(v), hence, MINRES is a non-linear operator. So MATLAB's gmres will fail. And for this reason, I think your question belongs to scicomp StackExchange. Believe me, industry experts fall into this trap.
One quick fix is to set MINRES tolerance to be less than GMRES tolerance, for example tol_minres = 1e-10, tol_gmres = 1e-7. In that case, MINRES would become a linear operator to the tolerance of GMRES, i.e Mfun(x+av) = Mfun(x)+aMfun(v)+E where |E|<1e-7. So from the perspective of GMRES Mfun is a linear operator as E` will be neglected. Note that, this is handwavy, no one did the analysis for this as far as I know.