function [lam, V, info] = plate_fem_eigs(K, M, nwant, nb, maxit, tol)
% PLATE_FEM_EIGS  The lowest eigenpairs of K v = lambda M v by subspace iteration.
%
%   [lam, V, info] = plate_fem_eigs(K, M, nwant, nb, maxit, tol) finds the
%   nwant smallest eigenvalues of the sparse symmetric pencil (K, M), K
%   positive semidefinite and M positive definite, with a block of nb >
%   nwant vectors.  Each iteration maps the block through the shifted
%   inverse,
%
%     Y = (K - sigma M)^-1 M X,   sigma = -(2 pi 1000 Hz)^2 < 0,
%
%   then projects K and M onto Y (Rayleigh-Ritz) and takes the Ritz vectors
%   as the next block.  A negative shift keeps K - sigma M positive definite,
%   so one sparse Cholesky factorization after an approximate-minimum-degree
%   ordering serves every iteration, and the eigenvalues nearest the shift,
%   which are the lowest, converge first.  The extra nb - nwant vectors
%   guard the wanted ones: eigenvalue i converges like
%   ((lambda_i - sigma) / (lambda_{nb+1} - sigma))^iteration.  Degenerate
%   pairs, which a square plate has many of, are found as reliably as
%   simple eigenvalues, which is why this is a block method.
%
%   The iteration stops when the wanted eigenvalues change by less than the
%   relative tolerance tol between iterations, or after maxit iterations;
%   info reports the count, the largest relative residual
%   |K v - lambda M v| / |K v| among the wanted pairs, and the shift.  The
%   rigid-body modes sit at lambda = 0, so both measures are floored at the
%   scale of a 20 Hz mode rather than taken relative to nothing.
%   Eigenvalues ascend and the eigenvectors are M-normalized.
%
%   numbl's own eigs is far slower on this problem; this routine is written
%   so that the whole method stays in double precision in numbl.

ndof = size(K, 1);
nb = min(nb, ndof);
nwant = min(nwant, nb);
sigma = -(2 * pi * 1000)^2;
Ks = K - sigma * M;
p = amd(Ks);
R = chol(Ks(p, p));
rng(12345);
X = randn(ndof, nb);
lam_prev = Inf * ones(nb, 1);
lam_floor = (2 * pi * 20)^2;
converged = false;
it = 0;
lam = zeros(nb, 1);
while it < maxit
    it = it + 1;
    B = M * X;
    Yp = R \ (R' \ B(p, :));
    Y = zeros(ndof, nb);
    Y(p, :) = Yp;
    % the columns shrink like 1 / (lambda - sigma), four orders of magnitude across the
    % block; scaling them keeps the projected mass matrix well conditioned
    Y = Y ./ repmat(sqrt(sum(Y.^2, 1)), ndof, 1);
    KY = K * Y;
    MY = M * Y;
    Kp = Y' * KY;
    Mp = Y' * MY;
    Kp = (Kp + Kp') / 2;
    Mp = (Mp + Mp') / 2;
    Rm = chol(Mp);
    Cm = Rm' \ Kp / Rm;
    Cm = (Cm + Cm') / 2;
    [Vc, Dc] = eig(Cm);
    [lam, order] = sort(diag(Dc));
    Vc = Vc(:, order);
    X = Y * (Rm \ Vc);
    % The rigid-body modes sit at zero with rounding noise of order eps |K|, so
    % they are only asked to be stable in absolute terms; the modes above the
    % scale of a 20 Hz mode are judged relatively.
    wanted = lam(1:nwant);
    prev = lam_prev(1:nwant);
    big = abs(wanted) >= lam_floor;
    rel = abs(wanted - prev) ./ max(abs(wanted), lam_floor);
    change = max([rel(big); 0]);
    small_ok = all(abs(wanted(~big) - prev(~big)) < 1e-3 * lam_floor);
    lam_prev = lam;
    if change < tol && small_ok
        converged = true;
        break
    end
end
V = X(:, 1:nwant);
lam = lam(1:nwant);
% residuals of the wanted pairs
res = zeros(nwant, 1);
KV = K * V;
MV = M * V;
for k = 1:nwant
    r = KV(:, k) - lam(k) * MV(:, k);
    res(k) = norm(r) / max(norm(KV(:, k)), lam_floor * norm(MV(:, k)));
end
info.iterations = it;
info.converged = converged;
info.residual = max(res);
info.sigma = sigma;
info.block = nb;
end