Computing a solution of the Laplace-Eigenvalueproblem with Neumann-b.c. Good day!  
I was considering the Laplace-Eigenvalueproblem with Neumann b.c., i.e. find $u \in H^1(\Omega) \setminus \{0\}$ and $\lambda \in \mathbb{R}$, such that:
\begin{eqnarray}
 -\Delta u  \ &=& \ \lambda u \ \ \mathrm{in} \ \Omega,\\
 \frac{\partial u}{\partial \overrightarrow{n}} \ &=& \ 0 \ \ \ \ \mathrm{on} \ \partial \Omega.
\end{eqnarray}
(the following is just a summary for Finite-Element-Methods for solving this problem, if you're already familiar with it, the question starts below)
Or in variational formulation in the proper Testspace for Neumann-problems $\displaystyle V:= \{v \in H^1(\Omega) : \int_{\Omega}v(x) \mathrm{dx} = 0 \}$, we want to find $u \in V \setminus \{0\}$ and $\lambda \in \mathbb{R}$, such that:
\begin{equation}
   \langle \nabla u, \nabla v \rangle_{L_2(\Omega)} \ = \ \lambda \langle u, v \rangle_{L_2(\Omega)} \ \ \forall v \in V
\end{equation}
Now suppose $\Omega \subset \mathbb{R}^2$ is some bounded, polygonal domain. Making the usual Galerkin-Ansatz by chosing linear Finite Triangle-Elements, we only solve this on the finite-dimensional subspace $S^0_1 \subset H^1(\Omega)$ (all functions that are global continuous and element-wise linear polynomials). Taking the nodal (lagrange) basis $[\vartheta_1, \dots, \vartheta_N]$ for this space on the N nodes of the mesh, this leads us to the matrix eigenvalue problem:
\begin{equation}
 K_h \ \underline{u} \ = \ \lambda_h \ M_h \ \underline{u}
\end{equation}
With ,,stiffiness matrix'' 
\begin{equation*}
 K_h := \left[ \langle \nabla \vartheta_j, \nabla \vartheta_i \rangle_{L_2(\Omega)} \right]^{N}_{i,j=1} \in \mathbb{R}^{N,N}
\end{equation*}
and ,,mass matrix''
\begin{equation*}
 M_h := \left[ \langle \vartheta_j, \vartheta_i \rangle_{L_2(\Omega)} \right]^{N}_{i,j=1} \in \mathbb{R}^{N,N}.
\end{equation*}
$\textbf{Question:}$ 
Now if i want to solve the discrete problem, i first have to bring in the conformal ,,$\displaystyle \int_{\Omega}u(x) \mathrm{dx} = 0$'' into the FE-space, but what exactly is the proper way to do that? 
If i express that condition in something like $\displaystyle \sum_{k=1}^{N} \mu_k u_k = 0$ with $\underline{u} = [u_k]^N_{k=1}$ the values on node k, i could add a row to matrices $K_h$ and $M_h$, but this would destroy the symmetry of the matrices and i don't even know if that eigenvalue problem for rectangular matrices would even make sense anymore. 
I would be grateful for help!
 A: OK, 
i came to the conclusion that you don't need to enter the average-condition ,,$\displaystyle \int_{\Omega}u(x) \mathrm{dx} = 0$'' at all. 
The variational formulation on the space V is equivalent in solving the variational formulation on the whole space $H^1(\Omega)$, i.e. find $u \in H^1(\Omega) \setminus \{0\}$ and $\lambda \in \mathbb{R}$, such that:
\begin{equation}
 \langle \nabla u, \nabla v \rangle_{L_2(\Omega)} \ = \ \lambda \langle u, v \rangle_{L_2(\Omega)} \ \ \forall v \in H^1(\Omega)
\end{equation}
By chosing $v \equiv 1 \in H^1(\Omega)$, the above equation reads:
\begin{align}
 \langle \nabla u, \nabla 1 \rangle_{L_2(\Omega)} \ &= \ \lambda \langle u, 1 \rangle_{L_2(\Omega)} \\
 \Leftrightarrow \hspace{13ex} 0 \ &= \ \lambda \int_{\Omega} u(x) \mathrm{dx}
\end{align}
For all $\lambda > 0$, this gives us already that the eigenfunctions u satisfy the average condition $\Rightarrow u \in V$. The solutions for $\lambda = 0$ of the Laplace-Eigenvalueproblem are uninteresting to compute, since we already know that the Eigenspace to $\lambda = 0$ consists just of the constant functions on $\Omega$. 
$M_h$ should be symmetric and positiv definit, since $\langle \cdot, \cdot \rangle_{L_2(\Omega)}$ is a scalar product on $L_2(\Omega) \supset S^0_1$. So we could just solve the standard symmetric eigenvalueproblem of the form:
\begin{equation}
 L^{-1} K_h L^{-T} \ \underline{u} \ = \ \lambda_h \  \ \underline{u}
\end{equation}
With $M_h = LL^T$ Cholesky-decomposition, so $L^{-1} K_h L^{-T}$ is symmetric and similar to $M_h^{-1} K_h$ ($\rightarrow$ same eigenvalues, eigenvectors of the original problem can be get by the transformation $u \to L^Tu$).
