Consider the generalized eigenproblem:
$$ [K]\{u\} = \omega^2 [M]\{u\} $$
One can transform the generalized eigenproblem into the standard eigenproblem:
$$ [D] = [K]^{-1}[M] $$
$$ [D]\{u\} =\lambda\{u\} \qquad \text{where} \qquad \lambda=\frac{1}{\omega^2} $$
And then, to compute the largest $\lambda$ (or the lowest $\omega^2$), $\lambda_1$, one can use the Power Iteration method, iterating the following expression (also obtaining the first eigenvector, $\{u\}$):
$$ \{u\}_{k+1} = [D]\{u\}_k \qquad \text{and} \qquad \frac{||\{u\}_{k+1}||}{||\{u\}_k||} \rightarrow \lambda_1 \qquad \text{as} \qquad k \rightarrow \infty $$
As answered here, to find the next largest $\lambda$, $\lambda_2$, one can replace $[D]$ by:
$$ [B] = [D] - \lambda_1 \{u\}\{u\}^T $$
And reiterate:
$$ \{v\}_{k+1} = [B]\{v\}_k \qquad \text{and} \qquad \frac{||\{v\}_{k+1}||}{||\{v\}_k||} \rightarrow \lambda_2 \qquad \text{as} \qquad k \rightarrow \infty $$
And so on for the next largest eigenvalue and associated eigenvector.
However
The presented method consists in first calculating $[D]$. When dealing with large sparse systems, where $[M]$ and $[K]$ are stored in a CSR format, calculating $[D]$ is just not an option.
As shown in this PDF, one can divide the Power Iteration method in two steps:
$$ \{x\}_k = [M]\{u\}_k $$
$$ [K]\{u\}_{k+1} = \{x\}_k $$
This works very well as one only needs to factor $[K]$ once, and only needs to solve a linear system instead of calculating $[D]$.
My question is
With the second, more efficient method, how does one calculate the second and following eigevalues and eigenvectors?
