Although I'm quite late, I think I've got the solution acording to your algorithm.
As far as I know the algorithm is known as "symmetric Gaussian elimination". It's basically the normal Gaussian elimination for finding the inverse of a matrix, where you record the row operations on the right in the identity matrix. But for every row operation you do the exact same operation for the columns without recording it in the identity matrix and you stop as soon as you have gotten a diagonal matrix.
So here you go (rows in Latin, columns in Arabic):
$$ \left[
\begin{array}{cccc|cccc}
1&0&0&0 & 1&0&0&0\\
0&0&1&0 & 0&1&0&0\\
0&1&0&0 & 0&0&1&0\\
0&0&0&1 & 0&0&0&1\\
\end{array}
\right]
\quad (\;A\;\;|\;\;\mathbb{1}\;) $$
$$
\left[
\begin{array}{cccc|cccc}
1&0&0&0 & 1&0&0&0\\
0&1&1&0 & 0&1&1&0\\
0&1&0&0 & 0&0&1&0\\
0&0&0&1 & 0&0&0&1\\
\end{array}
\right]
\quad \text{II = II+III}\\
\left[
\begin{array}{cccc|cccc}
1&0&0&0 & 1&0&0&0\\
0&2&1&0 & 0&1&1&0\\
0&1&0&0 & 0&0&1&0\\
0&0&0&1 & 0&0&0&1\\
\end{array}
\right]
\quad \text{2 = 2 + 3}\\ $$
$$
\left[
\begin{array}{cccc|cccc}
1&0&0&0 & 1&0&0&0\\
0&2&1&0 & 0&1&1&0\\
0&3&1&0 & 0&1&2&0\\
0&0&0&1 & 0&0&0&1\\
\end{array}
\right]
\quad \text{III = III+II}\\
\left[
\begin{array}{cccc|cccc}
1&0&0&0 & 1&0&0&0\\
0&2&3&0 & 0&1&1&0\\
0&3&4&0 & 0&1&2&0\\
0&0&0&1 & 0&0&0&1\\
\end{array}
\right]
\quad \text{3 = 3 + 2}\\ $$
$$
\left[
\begin{array}{cccc|cccc}
1&0&0&0 & 1&0&0&0\\
0&2&3&0 & 0&1&1&0\\
0&0&-1&0 & 0&-1&1&0\\
0&0&0&1 & 0&0&0&1\\
\end{array}
\right]
\quad \text{III = 2*III - 3*II}\\
\left[
\begin{array}{cccc|cccc}
1&0&0&0 & 1&0&0&0\\
0&2&0&0 & 0&1&1&0\\
0&0&-2&0 & 0&-1&1&0\\
0&0&0&1 & 0&0&0&1\\
\end{array}
\right]
\quad \text{3 = 2*3 - 3*2}\\ $$
$$
\left[
\begin{array}{cccc|cccc}
1&0&0&0 & 1&0&0&0\\
0&2&0&0 & 0&1&1&0\\
0&0&-2&0 & 0&-1&1&0\\
0&0&0&1 & 0&0&0&1\\
\end{array}
\right]
\quad (\;D\;\;|\;\;P\;) $$
Please let me know if I made any errors in the calculation.