% Solution of a sequence of linear systems. If BiCGStab is convergent throughout the sequence,
% one can try this code.

% Suppose ml(n)bicgstab is used to solve a sequence of m systems A_i x = b_i
% where A_i is N-by-N. This code dynamically searches values for n in an 
% user-provided interval [n_min, n_max] so that the time per iteration is 
% as small as possible.

n_min = 2; 
n_max = 20; 
step = 3;          % step size, an integer >= 1

Q = sign(randn(N, n_max));    % a random sign-matrix
max_it = 3*N; 
tol = 10^{-7}; 
kappa = 0;      % or, say, kappa = 0.7

n = 10;       % initial value for n, an integer picked from [n_min, n_max].

walk = 1;     % walk =  1, search forward; = -1, search backward.
t1 = inf;     % solution time of the previous system.

for i = 1:m 
    x = zeros(N,1);   % initial guess for the i-th system.
    Q(:,1) = b_i - A_i*x;
    tic
    [x,err,iter,flag] = mlbicgstabt(A_i,x,b_i,Q(:,1:n),M_i,max_it,tol, kappa);
    t2 = toc / iter;   % time per iteration of the current system.
    if walk == 1
        if t2 < t1
            n = min(n + step, n_max); 
            t1 = t2;
        else
            n = max(n - step, n_min); 
            t1 = t2; 
            walk = -1;
        end
    else
        if t2 < t1
            n = max(n - step, n_min);
            t1 = t2;
        else
            n = min(n + step, n_max); 
            t1 = t2; 
            walk = 1;
        end
    end
end