I'm trying to justify why the method presented by Ivan Niven (and co-authors) to solve linear diophantine equations actually works. I refer to the device that uses identity matrix, presented on page 218, of the fifth edition of the book An Introduction to The Theory of Numbers.
Given the diophantine equation ax + by = c, the table is formed:
a b c
1 0
0 1
The following three operations can be used on the first two columns:
(C1) Add an integral multiple of one of the first two columns to the other;
(C2) Exchange the first two columns;
(C3) Multiply all elements of one of the first two columns by -1.
The process terminates when a value of zero is obtained in the position initially occupied by a or b.
For example, given the equation 147x + 258y = 369, we finally arrive at:
0 3 369
86 -7
-49 4
From the first line, we take 3v = 369, that is: v = 123.
From the second line: x = 86 u - 7v, that is: x = 86u -861
of the third line: y = -49u + 4v, that is: y = -49u + 492
I ask for help.



First look at what happens just to $a$ and $b$ in this process. You will find that it is simply the Euclidean algorithm, in a sort of pedestrian manner, where instead of writing $$a=bQ+R$$ you say "let's multiply $b$ by $-1$ and then add the result $Q$ times to $a$, so as to get $a+Q(-b)=R$ instead of $a$.
The couple $a,b$ is replaced by $R,b$ and you keep going, until one of the numbers in that row is zero, that is when you reach a zero remainder, which means that the other number, $3$ in your example, is going to be $\gcd(a,b)$. If you ever have solved linear Diophantine equations with one or another method this should ring a bell.
Now what about the second and last rows?
The array that you gave is the skeleton of the following equation (which is the one you want to solve!): $$\left(\matrix{a&b\\1&0\\0&1}\right)\left(\matrix{x\\y}\right)=\left(\matrix{c\\x\\y}\right)$$
The second and last rows look stupid right now. But they are going to keep track of the change of variables that you make when you perform column operations.
If, say, you need to subtract the second column twice from the first. This corresponds to multiplying the matrix $\left(\matrix{a&b\\1&0\\0&1}\right)$ by $\left(\matrix{1&0\\-2&1}\right)$. The original equation is equivalent to $$\underbrace{\left(\matrix{a&b\\1&0\\0&1}\right)\left(\matrix{1&0\\-2&1}\right)}_{\text{Your new $3\times 2$ matrix}}\underbrace{\left(\matrix{1&0\\2&1}\right)\left(\matrix{x\\y}\right)}_{\text{Your new variable vector}\\\text{call it $\left(\matrix{u\\v}\right)$} }=\left(\matrix{c\\x\\y}\right)$$
The exact same process is going on in a recent answer that I gave on how to compute quotients of finite abelian groups.