function [V, Lam, Res, info] = pal(M, C, K, sigma, neig, m, opts)
%   [V, Lam, Res] = pal(M, C, K, sigma, neig, m) computes neig eigenvalues 
%   close to sigma, for the quadratic eigenvalue problem of low rank
%   damping Q(e) * v = (e^2 * M + e * C + K) * x = 0 with C = E * F
%   being a low rank matrix.  m specifies the Pade approximation order.
%    
%   On output, Lam is a vector of the neig computed eigenvalues, with each 
%   column of V being the corresponding (unitary) eigenvector. 
%   Res is a vector of the relative backward error of the eigenpairs 
%   |Q(e) * v|_2 / (|e|^2 * |M|_1 + |e| * |C|_1 + |K|_1).
%
%   NOTE THAT the potential `pole' eigenvalues of the REP will be 
%   automatically removed, so the number of eigenvalues returned <= neig.
%    
%   opts specify options
%      - opts.E, opts.F:    user provide rank revealing factorization 
%                           C = E * F.
%      - opts.disp[{0} | 1 | 2]: diagnostic information display level.
%       
%   IMPORTANT NOTES:
%   The damping matrix C should be extremely sparse, otherwise 
%   a rank revealing decomposition C = E * F should be provided
%   via opts.E, and opts.F.
%    
%   [V, E, Res, info] = pal(M, C, K, sigma, k, m) also return the
%   info structure with 
%       - info.rank         rank(C)
%       - info.E, F         rank revealing factorization C = E * F
%       - info.zeta         scaling factors 
%       - info.mu           REP eigenvalues mu corresponding to Lam
%       - info.poles        REP poles 
%       - info.poles_Lam:   eigenvalues removed as potential poles, with
%         info.poles_mu     corresponding mu
%         info.poles_V      and eigenvector V.
%
%
%   REFERENCES:
%
%   [1]  Ding Lu, Xin Huang, Zhaojun Bai, and Yangfeng Su, 
%        A Padé approximate linearization algorithm for solving the 
%        quadratic eigenvalue problem with low-rank damping,
%        Int. J. Numer. Methods Eng., 2015. 103(11): 840--858. 
%
%   AUTHORS:
%       Ding Lu, Fudan Univ.
%       Modified: 22-04-2016.
% 

tolp = 1e-3; % tolerance in pole checking |mu - poles| < tol |mu| 
n = length(C);

if sigma == 0
    warning('zero shift is not allowed');
    V = []; Lam = []; Res = []; info = [];
    return
end

if nargin < 7, opts = []; end; 

% Set info display level
if isfield(opts, 'disp')
   dodisp = opts.disp;
else
   dodisp = 0;
end

% Rank-revealing decomposition of damping matrix C = E * F
tic;
if isfield(opts, 'E')
    E = opts.E;  F = opts.F;
    l = size(E,2);
    % compute sigma splitting factor 
    if isequal(E, F)    % symmetric case
        ne = 1; 
        nf = 1; 
    else                % non-symmetric
        [Q, R] = qr(E, 0);
        E = Q;	ne = 1; % norm of E
        F = F * R';	nf = norm(F); % norm of F
    end
else
    [E, F, l, ne, nf] = sparse_lufac(C);
end;
t0 = toc;

% Split sigma = sigma1 * sigma2
sigma1 = sqrt(sigma * nf / ne);
sigma2 = sqrt(sigma * ne / nf);

% Matrix norms 
nm = norm(M,1);
nc = norm(C,1);
nk = norm(K,1);

% Scaling factors
z1 = abs(sigma)^2 * nm;
z2 = 2 * m * abs(sigma) * nc;
zeta = 1 / max([z1, z2, nk]);
sz = sqrt(zeta);

% Pade approximation parameters
d = (2 * m + 1);
gamma = 2 / d  * sin([1:m]' * pi / d ).^2;  
xi = cos([1:m]' * pi / d).^2;
a = (gamma./xi).^.5;
D = - diag(xi);
poles = - 1 ./ xi;

% Approximate linearization Mat-Vec operator
Qi = sigma^2 * M + sigma * C + K;
tic; [LL1 UU1 P11] = lu(Qi, 0.1); t1 = toc;
invqi = @(x) UU1\(LL1\(P11*x));

P12 = speye(n);
X1 = (sigma1 / sz) * kron(speye(l), a.' * D);
X2 = kron(speye(l), D); 
X3 = sigma2 * sz * kron(speye(l), a);

afun01 = @(x) sigma^2 * M * x(1:n)+ E * X1 * x(n+1:end);
afun1 = @(x) - invqi( afun01(x) );
afun2 = @(x, y) [x; X2 * y - X3 * (F.' * x )];
afun = @(x) afun2(afun1(x), x(n+1:end));

% Compute eigenvalues by eigs
Lsize = n + l * m;
opts0.isreal = false;
opts0.disp = dodisp;
tic; [V1 MU] = eigs(afun, Lsize, neig, 'lm', opts0);
t2 = toc;

% Recover eigenpairs
mu = 1./diag(MU);
[mu ind] = sort(mu);
Lam = sigma * sqrt(mu + 1);
V = V1(1:n, ind);  

% Potential pole removal
mu_minus_poles = mu * ones(1,m) - ones(neig,1) * poles'; % _{i,j} = mu_i - pole_j
distmp = min(abs(mu_minus_poles'), [], 1)'; % distmp_i = min_j |mu_i - pole_j|
idx = find(distmp < abs(mu) * tolp);
if ~isempty(idx)  % remove potential poles
    info.poles_Lam = Lam(idx);
    Lam(idx) = [];
    info.poles_V = V(:,idx);
    V(:,idx) = [];
else
    info.poles_Lam = []; 
    info.poles_V = []; 
end

% Compute relative backward errors 
nneig = length(Lam);
Res = zeros(nneig, 1);  
for j = 1:nneig
    e = Lam(j);
    v = V(:,j);
    v = v / norm(v);
    Res(j)  = norm( e^2 * M * v + e * C * v + K * v ) /...
        (abs(e)^2 * nm + abs(e) * nc + nk);
    V(:, j) = v;
end

% Timing statistic display 
if dodisp
    disp(['PAL timing statistcs:']); 
    disp(['  -- Rank revealing factorization ', num2str(t0),' s']);
    disp(['  -- LU factorization ', num2str(t1),' s']);
    disp(['  -- EIGS ', num2str(t2),' s']);
end

% Information structure
info.rank = l; info.E = E; info.F = F;
info.mu = mu; info.poles = poles;
info.zeta = zeta;

% END
return;



%--------------------INTERNAL FUNCTIONS-------------------%
function [L, U, r, l, u] = sparse_lufac(C)
% Compute rank-revealing factorization LU = C.
% l and u are 2-norm of L and U.
n = size(C,1);
[idx, idy, ~] = find(C); % Extract nnz matrix with 
			 % row index idx, and column index idy
idx = unique(idx); % return the value of idx and idy with no repetition
idy = unique(idy);
CC = full( C(idx, idy) ); % extract the non-zero submatrix
[Q, S, P] = svd(CC); % SVD for submatrix
s = diag(S); % estimate the numerical rank by 
	     % r = sum( s > eps*max(s)*max(size(EE)) ) 
	     % where s =diag(S)
s1 = S(1,1);
r = sum(s > max(size(CC)) * eps(s1));
if ishermitian(C) % symmetric LDLT 
   LL = Q(:,1:r) * sqrt(S(1:r,1:r));
   UU = P(:,1:r) * sqrt(S(1:r,1:r));
   L = sparse(n,r);
   L(idx,:) = LL;
   U = L;
   l = sqrt(s1);
   u = l;
else % QR factorization 
   UU = P(:,1:r) * S(1:r,1:r);
   L = sparse(n,r); 
   L(idx,:) = Q(:,1:r);
   U = sparse(n,r); 
   U(idy,:) = UU;
   l = 1;
   u = s1;
end

% END
return
