It is messy because you have misunderstood the problem. While $q(\underline{v})$ is induced by the bilinear form $f(\underline{u}, \underline{v})=\underline{v}^TA\underline{u}$, where $A$ is your $3\times 3$ coefficient matrix, $q$ is quadratic, not bilinear, also not a linear transformation. So, what you are asked to do is to find a decomposition of the form $A = P^TDP$ (where $P$ is invertible and the diagonal of $D$ does not necessarily contain any eigenvalue of $A$), but you have confused this with an eigenvalue decomposition $A = P^{-1}DP$. Surely, as your matrix $A$ is real symmetric, you can do both by performing an orthogonal decomposition $A=Q^TDQ$ where $QQ^T=I$ and $D$ contains the eigenvalues of $A$, but this is simply not required.
In general, you can find a decomposition $A = P^TDP$ by using elementary row/column operations. This is somewhat akin to finding a row-reduced echelon form of a matrix, but here we need to perform both an elementary row operation and a corresponding elementary column operation at each step. In other words, if, in a certain step, you multiply $A$ by an elementary matrix $E$ on the left, you should also mutiply $A$ by $E^T$ on the right.
For the problem you describe, however, simple inspection plus some completing-square trick is enough. Note that $ \begin{eqnarray} &&x_1^2 + x_2^2 + 9x_3^2 + 2x_1x_2 - 6x_1x_3 - 5x_2x_3\\ &=&(x_1 + x_2 - 3x_3)^2 + x_2x_3\\ &=&(x_1 + x_2 - 3x_3)^2 + \frac14[(x_2 + x_3)^2 - (x_2 - x_3)^2]. \end{eqnarray} $ So you may take $B=\{(x_1 + x_2 - 3x_3),\ (x_2 + x_3),\ (x_2 - x_3)\}$. You may verify that $A = P^TDP$ where $ P=\begin{pmatrix} 1&1&-3\\0&1&1\\0&1&-1 \end{pmatrix}, \ D=\begin{pmatrix} 1\\&\frac14\\&&-\frac14 \end{pmatrix}. $