function [alpha, info] = rpsas(A, epsln, opts)
% RPSAS: Subspace accelerated criss-cross type method for computing the 
%        real pseudospectral abscissa of large (sparse) matrices.
%
%   alpha = rpsas(A, epsln) returns the real pseudospectral abscissa 
%       alpha of a real square matrix A with perturbation level epsln.
%
%   [alpha, info] = rpsa(A, epsln) also returns an info structure
%       - info.succ:    false if there was an error during the iteration
%       - info.its      total number of iterations (number of svds/eigs).
%       - info.z        a vector of complex numbers 
%                               z_i = alpha_i + beta_i  (1<=i<=its+1) 
%                       that are the approximate rightmost points in the 
%                       reduced real pseudospectrum in the ith iteration.
%                       z0 is the starting (rightmost) eigenvalue.
%       - info.V        subspace used for the iterations.
%       - info.message  contains the stopping message.
%
%   [alpha, info] = rpsa(A, epsln, opts) allows specification of the 
%       parameters opts:
%       - opts.maxit (30)       maximum iterations.
%       - opts.tol (1.0E-8)     stopping (relative) increment tolerance.
%       - opts.mode ('s')       subspace expansion mode: 's'(svd) or 'e'(eig).
%       - opts.gmode ('a')      gamma parameter used for subspace expansion:
%                               gmode = 'a'(approximate), 'e'(exact).
%                               The exact option is computationally expensive. 
%       - opts.neigs (6)        number of eigenvalues for eigs fun.
%       - opts.nsvds (6)        number of singular values for eigs fun.
%       - opts.dodisp(0)        0 no output convergence history.
%                               1 display convergence history of outer loop.
%                               2 also display convergence of inner loop 
%                                   (criss-cross iterations).
%
%       
%   References: 
%       [1] Ding Lu and Bart Vandereycken, Criss-cross type algorithms 
%           for computing the real pseudospectral abscissa,
%           manuscript available at http://www.unige.ch/math/vandereycken.
%
%   Authors:
%       Ding Lu and Bart Vandereycken, University of Geneva.
%
%   Version:
%       0.1, Jun. 13, 2016.

%   POSSIBLE IMPROVEMENTS: 
%       IF n < 60, automatically run rpsa of dense problem instead
%       Handle error message of eigs, and svds in subspace expansion.

alpha = NaN;

if ~isreal(A) || size(A,1)~=size(A,2)
    info = Check_Convergence(0, [], [], [], []);
    return;
end
if nargin < 3, opts = []; end

% PARAMETER CHECK
[maxit, tol, dodisp, mode, neigs, nsvds, gmode] = Check_Opts(opts);

% INITIALIZATION 
[V, e0, msg] = rpsas_init_subspace(A, neigs);  % initial subspace
if isempty(e0) 
    info = Check_Convergence(-1, [], [], [], msg);
    return
end

n = length(A);
z = e0;

alpha = real(e0);                   % starting value (rightmost eigenvalue)
beta = imag(e0); 
opts0.tol = max([eps, 2*tol]);      % tolerance for inner loop criss-cross
                                    % ( Slightly less accurate than tol, say 2*tol
                                    %   IMPORTANT: outer loop itr benefits 
                                    %   from slightly higher tolerance compared to that of inner. )
opts0.do_acc = true;                % Apply inner loop acceleration
                                    % IMPORTANT for problem
                                    %   rdbrusselator_demo(3200), epsln=.1

if dodisp
    Display_Result(2, e0, 0);   % print first two lines
    if dodisp == 2, opts0.dodisp = true; end;  % set inner loop disp
end

% MAIN LOOP
for i = 1 : maxit
    % COMPUTE REDUCED RPSA
    opts0.alpha0 = alpha;   
    opts0.beta0 = beta; 
    [alpha, beta, gamma, info_rpsa] = rpsa_pencil(A*V, V, epsln, opts0); 
    z = [z; alpha + 1j * beta];

    if dodisp, Display_Result(2, z(end), i); end

    % CONVERGENCE CHECK 
    if i > 1 && alpha - real(z(i)) < tol * max([1, abs(alpha)])
        break; 
    end

    % SUBSPACE EXPANSION
    if gmode == 'e'     % compute exact gamma in gmode == e 
        [~, gamma] = eval_mu(A, speye(n), alpha, beta); 
    end
    %   Compute approximate maximizing singular vector(s) 
    [~, v, u] = second_smallest_sv(A, speye(n), alpha, beta, gamma, nsvds);
    v = reshape(v, n, []);

    %   Compute eigenvectors in mode 's' 
    if mode~='s', v = rpsas_expand_eigvec(A, u, v, epsln, neigs); end

    [V,~] = qr([V,v],0);
end

% RETURN INFO STRUCTURE
info = Check_Convergence(i, maxit, z, V, info_rpsa);

if dodisp, Display_Result(inf, [], [i-1,1 + (i-1)*(mode ~='s')]); end

% END
return;




% ================================================================ 
% Subfunctions
% ================================================================ 

function [maxit, tol, dodisp, mode, neigs, nsvds, gmode] = Check_Opts(opts);
if ~isempty(opts) && isfield(opts,'maxit') && ~isempty(opts.maxit) && opts.maxit > 0 % number of maximum iterations
    maxit = opts.maxit;
else
    maxit = 30;
end
if ~isempty(opts) && isfield(opts,'tol') && ~isempty(opts.tol) && opts.tol > 0 % stopping tolerance
    tol = opts.tol;
else
    tol = 1.0E-8;
end
if ~isempty(opts) && isfield(opts,'dodisp') && ~isempty(opts.dodisp)
    dodisp = opts.dodisp;
else
    dodisp = false;
end
if ~isempty(opts) && isfield(opts,'mode') && ~isempty(opts.mode) % subspace expansion mode
    mode = opts.mode;
else
    mode = 's';
end
if ~isempty(opts) && isfield(opts,'neigs') && ~isempty(opts.neigs)  % number of eigenvalues for eigs
    neigs = opts.neigs;
else
    neigs = 30;
end
if ~isempty(opts) && isfield(opts,'nsvds') && ~isempty(opts.nsvds)  % number of singular values for svds
    nsvds = opts.nsvds;
else
    nsvds = [];
end
if ~isempty(opts) && isfield(opts,'gmode') % gamma mode
    gmode = opts.gmode;
else
    gmode = 'a';
end

warning('off', 'MATLAB:normest:notconverge');

return


% -----------------------------------------------------------------

function Display_Result(mode, z, itr)
if ~isinf(mode)
    fprintf(' RPSAS Its %2i:  alpha = %+16.15e  beta = %8.7e \n', ...
        [itr, real(z), imag(z)]);
else
    fprintf(' ------------------------------------------\n');
    fprintf(' --- RPSAS svds calls %2i,  eigs calls %2i \n', itr);
    fprintf(' ------------------------------------------\n\n');
end;

return


% -----------------------------------------------------------------

function info = Check_Convergence(i, maxit, z, V, info_rpsa)

if i<0
    info.succ = false;
    info.its = 0;
    info.z = [];
    info.V = [];
    info.message = ['RPSAS subspace initilization fail: ', info_rpsa];
    return; end
info.succ = true;
info.its = i;
info.z = z; 
info.V = V;

if i == 0
    info.succ = false;
    info.message = ' Input matrix A must be a real square matrix! ';
    return;
end

if i~=maxit && info_rpsa.succ
    info.message = 'success';
elseif info_rpsa.succ
    info.succ = false;
    info.message =  ['Maximun number of iteration has achieved,'...
                     ' solution can be inaccurate.'];
else
    info.succ = false;
    info.message = ['RPSA_PENCIL does not converged: ',...
                     info_rpsa.message];
end

return
