function [X,Z] = ABgmres(A,B,b,its,Bits)
%ABgmres solves A x = b using GMRES with "right preconditioner" B.
%
% X = ABgmres(A,B,b,its,Bits)
%
% Implements the AB-GMRES method that solves the problem
%       min_x || A x - b ||_2
% by appying the GMRES method to the problem
%       min_z || A B z - b ||_2 ,  x = B z ,
% where B acts as a "right preconditiner."
%
% Input:  A    - Coefficient matrix.
%         B    - "Preconditioner," either a matrix or a function handle
%                to an iterative solver.
%         b    - Right-hand side.
%         its  - Number of AB-GMRES iterations.
%         Bits - Number of iterations in the iterative solver for the
%                preconditioner, in case B is a function handle.
%
% Output: X    - All iteration vectors x as columns in X.
%         Z    - All "internal" iteration vectors z.
%
% For simplicity we use backslash to solve the projected problem (i.e., we
% do not update a QR factorization of the Hessenberg matrix).

% Per Christian Hansen, DTU Compute, October 1, 2021.

% K. Hayami, J.-F. Yin, and T. Ito, "GMRES methods for least squares
% problems," SIAM J. Matrix Anal. Appl., 31 (2010), pp. 2400-2430.
% K. Morikuni and K. Hayami, "Convergence of inner-iteration GMRES methods
% for rank-deficient least squares problems," SIAM J. Matrix Anal. Appl.,
% 36 (2015), pp. 225--250.

% Initialization.
if nargin < 4, error('Not enough input arguments'), end
if isa(B,'function_handle') && nargin < 5
    error('Must specify number of iteations for preconditioner')
end

% Prepare for iterations.
Q = zeros(length(b),its+1);
H = zeros(its+1,its);
normb  = norm(b);
Q(:,1) = b / normb;
[m,n] = size(A);
X = zeros(n,its);
Z = zeros(m,its);

% Run its iterations (there is no stopping rule except breakdown).
f = waitbar(0,'ABgmres');
for k = 1:its
    waitbar(k/its,f);
    % Arnoldi step calculations.
    if isa(B,'function_handle')
        z = B(A,Q(:,k),Bits);
        v = A*z;
    else
        v = A*(B*Q(:,k));
    end
    for j = 1:k
        H(j,k) = Q(:,j)'* v;
        v = v  - H(j,k)*Q(:,j);
    end
    H(k+1,k) = norm(v);
    if H(k+1,k) == 0 % Breakdown so stop.
        disp('Breakdown')
        break
    end
    Q(:,k+1) = v / H(k+1,k);
    e1 = [1;zeros(k,1)];
    y = H(1:k+1,1:k) \ (normb*e1);  % For simplicity we just use Matlab's \.
    z = Q(:,1:k)*y;
    if nargout==2, Z(:,k) = z; end
    if isa(B,'function_handle')
        X(:,k) = B(A,z,Bits);
    else
        X(:,k) = B*z;
    end
end

close(f)