1function [lam, V, info] = plate_fem_eigs(K, M, nwant, nb, maxit, tol)
2% PLATE_FEM_EIGS The lowest eigenpairs of K v = lambda M v by subspace iteration.
3%
4% [lam, V, info] = plate_fem_eigs(K, M, nwant, nb, maxit, tol) finds the
5% nwant smallest eigenvalues of the sparse symmetric pencil (K, M), K
6% positive semidefinite and M positive definite, with a block of nb >
7% nwant vectors. Each iteration maps the block through the shifted
8% inverse,
9%
10% Y = (K - sigma M)^-1 M X, sigma = -(2 pi 1000 Hz)^2 < 0,
11%
12% then projects K and M onto Y (Rayleigh-Ritz) and takes the Ritz vectors
13% as the next block. A negative shift keeps K - sigma M positive definite,
14% so one sparse Cholesky factorization after an approximate-minimum-degree
15% ordering serves every iteration, and the eigenvalues nearest the shift,
16% which are the lowest, converge first. The extra nb - nwant vectors
17% guard the wanted ones: eigenvalue i converges like
18% ((lambda_i - sigma) / (lambda_{nb+1} - sigma))^iteration. Degenerate
19% pairs, which a square plate has many of, are found as reliably as
20% simple eigenvalues, which is why this is a block method.
21%
22% The iteration stops when the wanted eigenvalues change by less than the
23% relative tolerance tol between iterations, or after maxit iterations;
24% info reports the count, the largest relative residual
25% |K v - lambda M v| / |K v| among the wanted pairs, and the shift. The
26% rigid-body modes sit at lambda = 0, so both measures are floored at the
27% scale of a 20 Hz mode rather than taken relative to nothing.
28% Eigenvalues ascend and the eigenvectors are M-normalized.
29%
30% numbl's own eigs is far slower on this problem; this routine is written
31% so that the whole method stays in double precision in numbl.
33ndof = size(K, 1);
34nb = min(nb, ndof);
35nwant = min(nwant, nb);
36sigma = -(2 * pi * 1000)^2;
37Ks = K - sigma * M;
38p = amd(Ks);
39R = chol(Ks(p, p));
40rng(12345);
41X = randn(ndof, nb);
42lam_prev = Inf * ones(nb, 1);
43lam_floor = (2 * pi * 20)^2;
44converged = false;
45it = 0;
46lam = zeros(nb, 1);
47while it < maxit
48 it = it + 1;
49 B = M * X;
50 Yp = R \ (R' \ B(p, :));
51 Y = zeros(ndof, nb);
52 Y(p, :) = Yp;
53 % the columns shrink like 1 / (lambda - sigma), four orders of magnitude across the
54 % block; scaling them keeps the projected mass matrix well conditioned
55 Y = Y ./ repmat(sqrt(sum(Y.^2, 1)), ndof, 1);
56 KY = K * Y;
57 MY = M * Y;
58 Kp = Y' * KY;
59 Mp = Y' * MY;
60 Kp = (Kp + Kp') / 2;
61 Mp = (Mp + Mp') / 2;
62 Rm = chol(Mp);
63 Cm = Rm' \ Kp / Rm;
64 Cm = (Cm + Cm') / 2;
65 [Vc, Dc] = eig(Cm);
66 [lam, order] = sort(diag(Dc));
67 Vc = Vc(:, order);
68 X = Y * (Rm \ Vc);
69 % The rigid-body modes sit at zero with rounding noise of order eps |K|, so
70 % they are only asked to be stable in absolute terms; the modes above the
71 % scale of a 20 Hz mode are judged relatively.
72 wanted = lam(1:nwant);
73 prev = lam_prev(1:nwant);
74 big = abs(wanted) >= lam_floor;
75 rel = abs(wanted - prev) ./ max(abs(wanted), lam_floor);
76 change = max([rel(big); 0]);
77 small_ok = all(abs(wanted(~big) - prev(~big)) < 1e-3 * lam_floor);
78 lam_prev = lam;
79 if change < tol && small_ok
80 converged = true;
81 break
82 end
83end
84V = X(:, 1:nwant);
85lam = lam(1:nwant);
86% residuals of the wanted pairs
87res = zeros(nwant, 1);
88KV = K * V;
89MV = M * V;
90for k = 1:nwant
91 r = KV(:, k) - lam(k) * MV(:, k);
92 res(k) = norm(r) / max(norm(KV(:, k)), lam_floor * norm(MV(:, k)));
93end
94info.iterations = it;
95info.converged = converged;
96info.residual = max(res);
97info.sigma = sigma;
98info.block = nb;
99end