Completing the square when the pivot are zeros. I watched the 18.06 MIT course on linear Algebra online. It shows a nice method to complete the square using Gauss elimination.
For example if I consider the quadratic form \begin{equation}x^2+y^2+4xy=(x,y)^T\cdot  \underbrace{\begin{pmatrix} 1&2\\2 & 1 \end{pmatrix}}_{A} \cdot(x,y) \end{equation}
I can complete the square by finding the $LU$ decomposition of the matrix $A$:
\begin{equation}A=\begin{pmatrix} 1&2\\2 & 1 \end{pmatrix}=\underbrace{\begin{pmatrix} 1&0\\2 & 1 \end{pmatrix}}_{L} \cdot \underbrace{\begin{pmatrix} 1&2\\0 & -3 \end{pmatrix}}_U \end{equation}
and then I complete the square by defining new variables given by the columns of $L$ and using as multiplier the pivot of $U$ such that the solution is  \begin{equation}x^2+y^2+4xy=(x+2y)^2-3y^2 \end{equation}
I like very much this method as it is easy to do without thinking too much. However I am not able to complete the square with this method when the pivots are zero. For example if I take the following symmetric matrix \begin{equation}B=\begin{pmatrix} 0&1\\1 & 0 \end{pmatrix} \end{equation} I was lucky and was able for this simple example to find that
\begin{equation}(x,y)^T\cdot  B \cdot(x,y) =2xy=\frac{1}{2} (x+y)^2-\frac{1}{2}(x-y)^2\end{equation} but I cannot find it from the Gauss elimination ($PLU$ factorization). It is really what I would like to be able to do so that I can handle matrices that are more complex. Someone knows how to do?

A more comprehensive approach is to diagonalize your matrix.  This is possible as long as it is symmetric.  $B$ is clearly symmetric, and its diagonalisation is given by
$$B = \begin{bmatrix}
-1 & 1\\
1 & 1
\end{bmatrix}
\begin{bmatrix}
-1 & 0\\
0 & 1
\end{bmatrix}
\begin{bmatrix}
-1/2 & 1/2\\
1/2 & 1/2
\end{bmatrix}.$$
The first matrix gives the eigenvectors.  Notice that the vectors $[1,0]^T$ and $[0,1]^T$ represent $x$ and $y$, respectively.
Thus, the eigenvectors are $[-1,1]^T := -x+y$ and $[1,1]^T := x+y$.
Thus,
$$B[x,y] = \begin{bmatrix}
-x + y,x+y
\end{bmatrix}
\begin{bmatrix}
-1 & 0\\
0 & 1
\end{bmatrix}
\begin{bmatrix}
-1/2x + 1/2y\\
1/2x + 1/2y
\end{bmatrix}.$$
$$=[-(y-x),(y+x)]\begin{bmatrix}
-1/2x+1/2y\\
1/2x+1/2y]\end{bmatrix} = -\frac{1}{2}(x-y)^2+\frac{1}{2}(x+y)^2.$$
The issue with the LU approach is that you are requiring one of the vectors in parenthesis to have a zero element (upper element of $L$ is zero), which is not always going to work.
