Eigensolver for Black-box matrix $\DeclareMathOperator{\diag}{diag}$
Consider the generalized eigenproblem $A\mathbf{x}=\lambda B\mathbf{x}$. When solving PDEs numerically (specifically, I am interested on finding the Dirichlet eigenmodes of the 2D laplacian in orthogonal coordinates, i.e. $\nabla^2 u=\lambda u$) this problem takes the form
$$ A_1X + XA_2^\text{T} = \lambda B \circ X$$
where $X\in\mathbb{R}^{m\times n}$ is our eigenvector (the eigenfunction $u$ sampled at the meshpoints), $A_1\in\mathbb{R}^{m\times m}$ and $A_2\in\mathbb{R}^{n\times n}$ are differentiation matrices, $B\in\mathbb{R}^{m\times n}$ is the Jacobian determinant (sampled at the meshpoints) and the symbol $\circ$ denotes  element-wise multiplication. I am only interested on the eigenmodes corresponding to the smallest eigenvalues.
With the help of the Kronecker product, the problem can be vectorized as follows:
$$ (I \otimes A_1 + A_2 \otimes I)\mathbf{x} = \lambda \diag(B)\mathbf{x}$$
where $\mathbf{x}\in\mathbb{R}^{mn}$ is a vector formed by stacking the columns of $X$ and $\diag{(B)}$ is a diagonal matrix formed with the columns of $B$.
The problem with this is that the resulting matrices are unnecessarily sparse and have size of $mn\times mn$, and when I try to solve it for modestly large sizes of $m$ and $n$, say 200, my computer runs out of memory. Clearly that is not the way.
I understand that modern eigenvalue algorithms are highly sophisticated, and that coding one by myself is not the best option. So I am looking for an iterative solver that only depends on the computation of the matrix-vector products with $A$ and $B$ and that allows the user to provide them as a black-box.
 A: The ARPACK library can be used to find few eigenvalues and eigenvectors. This library can solve the standard eigenvalue problem $Dx = \lambda x$ or generalized eigenvalue problem $Dx = \lambda E x$ if $E$ is a symmetric positive definite matrix.
Iterative eigensolvers may not converge or convergence can be very slow if smallest eigenvalues must be found. In this case shift-invert method should be used, i.e. find largest eigenvalues of $(D-\sigma E)^{-1}Ex = \lambda x$, where $\sigma$ is a small number such that $D-\sigma E$ is invertible. For given problem efficient eigensolver can be constructed only for $\sigma = 0$. In order to find smallest eigenvalues one needs to form $D^{-1}Ey$ for any vector $y$, i.e. solve $Dx = Ey$.
Thus, we need an efficient solver for the Sylverster equation:
$$A_1X + XA_2^T = B\circ Y \tag{1}$$
for any matrix $Y$. This equation must be solved many times with different $Y$. This can be done efficiently if one form the Schur factorization of $A_1$ and $A_2$. Then only solving linear equations with (almost) triangular matrices will be required. Required algorithms can be found in Lapack.
The operator $T = I\otimes A_1 + A_2\otimes I$ can be singular or nearly singular. The shift-invert method for finding the largest eigenvalue is backward stable even if $T$ is nearly singular. Thus, if $T$ is singular, then one can perturb slighly the matrix $T$ to make it invertible. This is equivaluent to adding  small perturbations to triangular factors of the Schur factorization of $A_1$ and $A_2$.
