You basically want to solve for the homogeneous solutions of the system
$$
\textbf{M}\, \ddot{\textbf{q}} + \textbf{D}\, \dot{\textbf{q}} + \textbf{K}\, \textbf{q} = \textbf{0}.
$$
It can be shown that each solution will be of the form
$$
\textbf{q}_i(t) = \textbf{u}_i\,e^{\mu_i\,t},
$$
where $\textbf{u}_i$ is a vector of the same dimension as $\textbf{q}$ (usually normalized to a length of one) and $\lambda_i$ is a scalar. In general any homogeneous solution can be constructed by taking a linear combination these solution (I will ignore for now that you can also have Jordan blocks of bigger then one by one, which will increase the possible solutions).
Substituting the solution into the homogeneous equation yields
$$
\textbf{M}\, \textbf{u}_i\,\mu_i^2\,\,e^{\mu_i\,t} + \textbf{D}\, \textbf{u}_i\,\mu_i\,\,e^{\mu_i\,t} + \textbf{K}\, \textbf{u}_i\,e^{\mu_i\,t} = \left(\textbf{M}\, \mu_i^2 + \textbf{D}\, \mu_i + \textbf{K}\right) \textbf{u}_i\,e^{\mu_i\,t} = \textbf{0}.
$$
Since $e^{\mu_i\,t}$ is always non-zero (for finite time) it can be ignored since the remaining terms on the left hand side of the equation still has equal to zero in order to satisfy the right hand side. Doing so already starts to make the equation look more like a generalized eigenvalue problem
$$
\textbf{A}\, \textbf{v} = \lambda\, \textbf{B}\, \textbf{v},
$$
however $\mu_i$ appears both linearly and quadratically. In the special case when $\textbf{B} = 0$ then the generalized eigenvalue problem can be formulated using $\textbf{A} = \textbf{K}$, $\textbf{B} = \textbf{M}$, $\textbf{v} = \textbf{u}_i$ and $\lambda = -\mu_i^2$. Another special case is when $\textbf{M}$, $\textbf{D}$ and $\textbf{K}$ are all simultaneously diagonalizable by pre and post multiplying them by the same matrices. Namely in that case the differential equation can be written as a system of decoupled second order differential equations, like a single mass-spring-damper system. One way such transformation could be constructed when $\textbf{M}$, $\textbf{D}$ and $\textbf{K}$ are symmetric is
$$
\mathbf{U} = \begin{bmatrix}\mathbf{u}_1 & \mathbf{u}_2 & \cdots & \mathbf{u}_n\end{bmatrix},
$$
where $\mathbf{u}_i$ are the eigenvectors of the eigenvalue problem with $\textbf{B} = 0$. The matrices with which you pre and post multiply would then be $\mathbf{U}^\top$ and $\mathbf{U}$ respectively. Often if the system is not simultaneously diagonalizable and the damping is small, then assuming it is can be a good approximation (so setting the non-diagonal elements of $\mathbf{U}^\top \mathbf{B}\, \mathbf{U}$ to zero). A specific form for which this is the case is called Rayleigh damping, when $\mathbf{B} = \alpha\,\mathbf{M} + \beta\,\mathbf{K}$.
The advantage of all the previous mentioned methods is that the generalized eigenvalue problem you have to solve are of the same dimension as $\mathbf{q}$. However if none of cases mentioned previously apply to your system you would have to solve the generalized eigenvalue problem of twice the dimension of $\mathbf{q}$. There are multiple ways of formulating this. A few of them are
$$
\begin{bmatrix}
\mathbf{I} & \mathbf{0} \\
\mathbf{0} & \mathbf{M}
\end{bmatrix}
\begin{bmatrix}
\dot{\mathbf{q}} \\ \ddot{\mathbf{q}}
\end{bmatrix} =
\begin{bmatrix}
\mathbf{0} & \mathbf{I} \\
-\mathbf{K} & -\mathbf{D}
\end{bmatrix}
\begin{bmatrix}
\mathbf{q} \\ \dot{\mathbf{q}}
\end{bmatrix},
$$
$$
\begin{bmatrix}
\mathbf{I} & \mathbf{0} \\
\mathbf{0} & \mathbf{I}
\end{bmatrix}
\begin{bmatrix}
\dot{\mathbf{q}} \\ \ddot{\mathbf{q}}
\end{bmatrix} =
\begin{bmatrix}
\mathbf{0} & \mathbf{I} \\
-\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{D}
\end{bmatrix}
\begin{bmatrix}
\mathbf{q} \\ \dot{\mathbf{q}}
\end{bmatrix},
$$
$$
\begin{bmatrix}
\mathbf{D} & \mathbf{M} \\
\mathbf{M} & \mathbf{0}
\end{bmatrix}
\begin{bmatrix}
\dot{\mathbf{q}} \\ \ddot{\mathbf{q}}
\end{bmatrix} =
\begin{bmatrix}
-\mathbf{K} & \mathbf{0} \\
\mathbf{0} & \mathbf{M}
\end{bmatrix}
\begin{bmatrix}
\mathbf{q} \\ \dot{\mathbf{q}}
\end{bmatrix}.
$$
For this extended state space it can easily be seen that the solutions will be of the form
$$
\begin{bmatrix}
\mathbf{q}_i(t) \\ \dot{\mathbf{q}}_i(t)
\end{bmatrix} =
\begin{bmatrix}
\mathbf{u}_i \\ \mu_i\,\mathbf{u}_i
\end{bmatrix} e^{\mu_i\,t}.
$$
Substituting this into the three equations above and factoring out $e^{\mu_i\,t}$ yields
$$
\begin{bmatrix}
\mathbf{I} & \mathbf{0} \\
\mathbf{0} & \mathbf{M}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}_i \\ \mu_i\,\mathbf{u}_i
\end{bmatrix} \mu_i =
\begin{bmatrix}
\mathbf{0} & \mathbf{I} \\
-\mathbf{K} & -\mathbf{D}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}_i \\ \mu_i\,\mathbf{u}_i
\end{bmatrix},
$$
$$
\begin{bmatrix}
\mathbf{I} & \mathbf{0} \\
\mathbf{0} & \mathbf{I}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}_i \\ \mu_i\,\mathbf{u}_i
\end{bmatrix} \mu_i =
\begin{bmatrix}
\mathbf{0} & \mathbf{I} \\
-\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{D}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}_i \\ \mu_i\,\mathbf{u}_i
\end{bmatrix},
$$
$$
\begin{bmatrix}
\mathbf{D} & \mathbf{M} \\
\mathbf{M} & \mathbf{0}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}_i \\ \mu_i\,\mathbf{u}_i
\end{bmatrix} \mu_i =
\begin{bmatrix}
-\mathbf{K} & \mathbf{0} \\
\mathbf{0} & \mathbf{M}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}_i \\ \mu_i\,\mathbf{u}_i
\end{bmatrix}.
$$
These are all just generalized eigenvalue problems. I am not an expert on numerically solving these problems, however I have few things that can be noted. Namely the second problem formulation requires a matrix inverse, which might not exist, but does reduce it to a normal eigenvalue problem (not sure if this requires any less computation time in general). The third problem formulation will yield symmetric matrices if $\textbf{M}$, $\textbf{D}$ and $\textbf{K}$ are symmetric as well, which might have some attractive properties when doing numerical calculations. So depending on your system and solver a different formulation might give more accurate results.
