Solver for generalized eigenvalue problem with multipoint constraints

We have the following generalized eigenvalue (set of) problem(s)

$$[K_R(\kappa)]\{u_R\} = \omega^2[M_R(\kappa)]\{u_R\}\quad \forall \kappa \in [\kappa_0, \kappa_1]$$

with

\begin{align} &K_R(\kappa) = T^H(\kappa) K T(\kappa)\, ,\\ &M_R(\kappa) = T^H(\kappa) M T(\kappa)\, ,\\ &u = T(\kappa) u_R\, . \end{align}

Where $$K$$ and $$M$$ are sparse and symmetric and come from a PDE, and $$T(\kappa)$$ is also sparse, non-square, complex and represents a set of multipoint constraints. Forming the matrix $$T$$ is relatively cheap compared with the other matrices since they have a fixed structure and are sparse.

If we consider a regular mesh with $$n^2$$ nodes, we have that:

\begin{align} &K\in \mathbb{R}^{n^2\times n^2}\, ,&M\in \mathbb{R}^{n^2\times n^2}\, ,\\ &u\in \mathbb{C}^{n^2}\, , & &\\ &K_R\in \mathbb{C}^{(n^2 - 2n + 1)\times (n^2 - 2n + 1)}\, ,&M_R\in \mathbb{C}^{(n^2 - 2n + 1)\times (n^2 - 2n + 1)}\, ,&\\ &u_R\in \mathbb{C}^{(n^2 - 2n + 1)}\, , & &\\ &T(\kappa) \in \mathbb{C}^{(2n -1)\times(n^2 - 2n + 1)}\, .& & \end{align}

We normally handle the problem in one of the following ways:

1. Assemble $$K$$ and $$M$$ and for each value of $$\kappa$$ form the product matrices $$K_R$$ and $$M_R$$. The main advantage of this method is that we have to assemble once, but then we loose the sparse nature of the problem.
2. Assemble $$K_R$$ and $$M_R$$ for each $$\kappa$$ value, conserving the sparse nature of the problem.

Question

Is there any way of solving the generalized eigenvalue problem

$$[T^H(\kappa) K T(\kappa)]\{u_R\} = \omega^2[T^H(\kappa) M T(\kappa)]\{u_R\}\quad \forall \kappa \in [\kappa_0, \kappa_1]\, ,$$

conserving the sparse nature of the system and assembling the matrices $$K$$ and $$M$$ only once?

• I think I'm missing something, but can't you use an iterative method for the eigenproblem, can't you work with the product of matrices as it is written? – VorKir Mar 23 at 3:47
• @VorKir, I think that you are right. I think I could just define the linear operators and pass that to the solver. – nicoguaro Mar 23 at 15:42
• @nicoguaro did it work? it seems to me that it should have. – Anton Menshov Jun 5 at 21:37
• Yes it did. I think I should write an answer about that. – nicoguaro Jun 5 at 23:25