function X = BAgmres(A,B,b,its,Bits)
%BAgmres solves A x = b using GMRES with "left preconditioner" B
%
% X = BAgmres(A,B,b,maxiter,its)
%
% Implements the BA-GMRES method that solves the problem
%       min_x || A x - b ||_2
% by appying the GMRES method to the problem
%       min_x || B A x - B b ||_2 ,
% where B acts as a "left 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 BA-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.
%
% 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 numnber of iteations for preconditioner')
end

% Prepare for iterations.
f = waitbar(0,'BAgmres');
waitbar(1/(its+1),f);
if isa(B,'function_handle')  % Note that we redifine b here.
    b = B(A,b,Bits);
else
    b = B*b;
end
normb  = norm(b);
Q = zeros(length(b),its+1);
Q(:,1) = b / normb;
n = size(A,2);
X = zeros(n,its);
H = zeros(its+1,its);

% Run its iterations (there is no stopping rule except breakdown).
for k = 1:its
    waitbar((k+1)/(its+1),f);
    % Arnoldi step calculations.
    z = A*Q(:,k);
    if isa(B,'function_handle')
        v = B(A,z,Bits);
    else
        v = B*z;
    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 \.
    X(:,k) = Q(:,1:k)*y;    
end

close(f)