Method 2 is just an algorithm to find a diagonal matrix which is congruent to the first, so maybe you can call it matrix congruence method.
Another method along the lines of matrices is to do orthogonal diagonalization, but method 2 would be faster...
You can do all the row operations first: but be careful then you have $P^TA$, so the safest then is to take the resulting matrix on the right, in your augmented matrix (which is $P^T$), transpose it and multiply on the right which gives you $P^TAP$ - I don't think you will gain too much when you do this - it is much simpler to just alternate row and column operations.
Completing squares I get, with $\;x,y,z\;$ instead of $\;x_i\;$ :
$$3x^2-2xy+2y^2-2yz+3z^2=3\left(x-\frac13y\right)^2+\frac43y^2+3\left(z-\frac13y\right)^2\implies$$
$$\implies \begin{cases}I\;\;\;\;\,a=x-\frac13y\\II\;\;\;b=y\\III\;c=z-\frac13y\end{cases}\;\;\;
\implies\;\;\begin{cases}x=a+\frac13b\\y=b\\z=c+\frac13b\end{cases}$$
So the base part is $\;\left(1,\,\frac13,\,0\right)\;,\;(0,1,0)\;,\;\left(0,\,\frac13,\,1\right)\;$ , and the quadratic can be expressed as
$$\color{red}{q(a,b,c)=3a^2+\frac43b^2+3c^2}$$
Observe that the eigenvalues of $\;A\;$ are
$$\det(xI-A)=\begin{vmatrix}x-3&1&0\\1&x-2&1\\0&1&x-3\end{vmatrix}=$$
$$=(x-2)(x-3)^2-2(x-3)=(x-3)(x-4)(x-1)$$
with eigenvectors:
$$\begin{align*}
&\bullet\;\lambda=1:\;\;\begin{cases}-2x+y\;\;\;\;\;\;\;=0\\\;\;\;\;\;\;\;\;\;\;y-2z=0\end{cases}\implies2x=2z=y\implies&\begin{pmatrix}1\\2\\1\end{pmatrix}\\{}
&\bullet\;\lambda=3:\;\;\begin{cases}\;\;\;\;\;\;y\;\;\;\,\;\;=0\\x+y+z=0\end{cases}\implies x=-z\,,\,y=0\implies&\begin{pmatrix}1\\0\\\!-1\end{pmatrix}\\{}
&\bullet\;\lambda=4:\;\;\begin{cases}x+y\;\;\;\;\;\;=0\\\;\;\;\;\;\,y+z=0\end{cases}\implies x= z=-y\implies&\begin{pmatrix}\;1\\\!-1\\\;1\end{pmatrix}\end{align*}$$
Now orthonormalize the above three vectors, form the matrix $\;P\;$ whose columns are the new vectors and check that
$$P^{-1}AP=\begin{pmatrix}1&0&0\\0&3&0\\0&0&4\end{pmatrix}$$
Best Answer
This is known as orthogonal diagonalisation, and the process is outlined on Wikipedia.
First write $q(x, y)$ as
$$q(x, y) = [x\quad y]A\left[\begin{array}\ x\\ y\end{array}\right]$$
where
$$A = \left[\begin{array}\ 0 & 3\\ 3 & 0\end{array}\right].$$
Note that $A$ is symmetric. As such, it is orthogonally diagonalisable. That is, there is an invertible matrix $P$ satisfying $P^T = P^{-1}$ such that $P^TAP$ is a diagonal matrix. To find $P$, you need to find an orthonormal basis of eigenvectors; these form the columns of $P$. The eigenvectors of $A$ are $\lambda_1 = 3$ and $\lambda_2 = -3$ with eigenvectors $v_1 = [1\quad 1]^T$ and $v_2 = [-1\quad 1]^T$ respectively. As eigenvectors from distinct eigenvalues are orthogonal, normalising $v_1$ and $v_2$, we obtain an orthonormal basis: $\left\{\frac{1}{\sqrt{2}}v_1, \frac{1}{\sqrt{2}}v_2\right\}$. Therefore
$$P = \frac{1}{\sqrt{2}}\left[\begin{array}{cc} 1 & -1\\ 1 & 1\end{array}\right].$$
Note that
$$P^TAP = \frac{1}{\sqrt{2}}\left[\begin{array}{cc} 1 & 1\\ -1 & 1\end{array}\right]\left[\begin{array}\ 0 & 3\\ 3 & 0\end{array}\right]\frac{1}{\sqrt{2}}\left[\begin{array}{cc} 1 & -1\\ 1 & 1\end{array}\right] = \left[\begin{array}{cc} 3 & 0\\ 0 & -3\end{array}\right]$$
which is diagonal.
Now make the variable substitution
$$\left[\begin{array}\ \hat{x}\\ \hat{y}\end{array}\right] = P^{-1}\left[\begin{array}\ x\\ y\end{array}\right].$$
Then we have
\begin{align*} q(x, y) &= 6xy\\ &= [x\quad y]A\left[\begin{array}\ x\\ y\end{array}\right]\\ &= \left[\begin{array}\ x\\ y\end{array}\right]^TA\left[\begin{array}\ x\\ y\end{array}\right]\\ &= \left(P\left[\begin{array}\ \hat{x}\\ \hat{y}\end{array}\right]\right)^TA\left(P\left[\begin{array}\ \hat{x}\\ \hat{y}\end{array}\right]\right)\\ &= [\hat{x}\quad \hat{y}]P^TAP\left[\begin{array}\ \hat{x}\\ \hat{y}\end{array}\right]\\ &= [\hat{x}\quad \hat{y}]\left[\begin{array}{cc} 3 & 0\\ 0 & -3\end{array}\right]\left[\begin{array}\ \hat{x}\\ \hat{y}\end{array}\right]\\ &= 3\hat{x}^2 - 3\hat{y}^2. \end{align*}
Geometrically, $P$ corresponds to a rotation of the coordinate system by $\frac{\pi}{4}$. To see this, note that $P$ is the rotation matrix for $\theta = \frac{\pi}{4}$.