$\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.