Unitarily diagonalize $P=V\Lambda V^{*}$, where the adjoint $V^{*}$ is just $V^T$ if your matrix is real. Having done this, do not store $V^{*}$ separately.
The matrix $\Lambda$ is to have the eigenvalues of $P$ along its main diagonal and is to be null elsewhere. The columns of $V$ are to be mutually orthogonal unit vectors, whereby $V^{*}V = I$.
The reason this works is that it lends $V$ the useful property that $V^{-1}=V^{*}.$ Therefore, because $\Lambda$ is also trivial to invert, you never ask your software to invoke the Gauss-Jordan algorithm that is desymmetrifying your matrix. Don't let your software decide for itself how to invert anything, but rather work around the problem as here explained by telling your software what the inverses of the various elements are.
EXPLANATION
A matrix represents a linear operation and, if the matrix is invertible, then it has certain properties. Moreover, if the matrix is symmetric, it has other properties. One cannot use a symmetric matrix effectively without exploiting the properties that belong to all symmetric matrices. If you try to treat a symmetric matrix like any old matrix, then, numerically, weird things are likely to happen.
When I said symmetric, I really meant self-adjoint, which is a jargon word but here is an example:
$$\left[\begin{array}{cc}2 & 4+i3 \\ 4-i3 & 7\end{array}\right].$$
So, as you see, a self-adjoint matrix is symmetric with complex conjugation. If the matrix is real, then "symmetric" and "self-adjoint" mean the same.
Another word for "self-adjoint" is Hermitian.
The symbol for the transpose is $P^T$ as you know. The symbol for the adjoint is $P^{*}=\text{conjugate}\{P^{T}\}$, but for real matrices $P^{*}=P^T$, so you can read the two symbols as though they were the same.
Now, most matrices are diagonalizable. Significantly, all self-adjoint (real symmetric) matrices are diagonalizable. This means that you can factor the self-adjoint matrix in a special form that resembles
$$P =
\left[\begin{array}{cc}0.8 & 0.6 \\ -0.6 & 0.8\end{array}\right]
\left[\begin{array}{cc}1 & 0 \\ 0 & 2\end{array}\right]
\left[\begin{array}{cc}0.8 & -0.6 \\ 0.6 & 0.8\end{array}\right].
$$
For convenience, we give names to the factors, abbreviating the last line as
$$
P = V\Lambda V^{*},
$$
the inverse of which conveniently is
$$
P^{-1} = V\Lambda^{-1} V^{*}.
$$
Writing out the factors,
$$
P^{-1} =
\left[\begin{array}{cc}0.8 & 0.6 \\ -0.6 & 0.8\end{array}\right]
\left[\begin{array}{cc}1 & 0 \\ 0 & 0.5\end{array}\right]
\left[\begin{array}{cc}0.8 & -0.6 \\ 0.6 & 0.8\end{array}\right].
$$
So, as matrices go, that's pretty easy.
Interestingly and usefully,
$$
P^n = V\Lambda^n V^{*},
$$
wherein, if you think about it,
$$
\Lambda^n = \left[\begin{array}{cc}1^n & 0 \\ 0 & 2^n\end{array}\right],
$$
so the computer really cannot very well mess that up.
My advice to you is this: don't let the computer directly handle $P$. Instead, make the computer deal in $\Lambda$ and $V$.
One last question is, how do you diagonalize? The answer is that your software almost certainly has a built-in function to diagonalize, but if you wish to learn how to do it, yourself, look up eigenvalue, eigenvector and/or diagonalization in your linear-algebra textbook.
Good luck.
