It helps to understand the idea behind the algorithm. It uses the subtractive Euclidean algorithm to change basis, iteratively transforming $\,m\Bbb Z + n\Bbb Z\,$ into $\,d\Bbb Z + d\Bbb Z = d\Bbb Z,\,$ where $\,d = \gcd(m,n),\,$ while simultaneously keeping track of the coordinates of some given number $\,k = v_0 m + u_0 n\,$ in each transformed basis. The (loop) invariant is that the bases are equivalent (same span), i.e. $\,x_j \Bbb Z + y_j\Bbb Z = m\Bbb Z + n\Bbb Z$. and that $\,v_j x_j + u_j y_j = k = v_0 x_0 + u_0 y_0,\,$ i.e. $\,(v_j,u_j)\,$ are the coordinates of $\,k\,$ in the current basis $\,(x_j,y_j).\,$ Span persistence follows because the subtractive Euclidean algorithm changes bases via elementary transformations (i.e. determinant $=\pm1),\,$ which are invertible. The coordinate persistence follows by simple algebra, naturally viewed as persistence of the determinant under column operations, i.e. if we represent $(x,y,u,v)$ by $\begin{bmatrix} x & y\\ -u & v \end{bmatrix}$ then the Euclidean reduction step subtracts the column with least first row entry from the other column, which preserves the determinant, e.g. if $\,x < y\,$ then the reduction step simply subtracts the $\rm\color{#c00}{second}$ column from the first, as displayed below
$$\begin{align} (x,y,u,v)\ &\to\ (x\!-\!y,y, u\!+\!v,v)\,\ \ \text{has matrix form below:}\\[.5em]
\begin{bmatrix} \ \ x & \color{#c00}y\\ -u & \color{#c00}v \end{bmatrix}
\ &\to \ \begin{bmatrix} \,\ \ x-\color{#c00}y & y\\ -(u+\color{#c00}v) & v \end{bmatrix}
= \begin{bmatrix} \ \ x & y\\ -u & v \end{bmatrix}
\begin{bmatrix} \ \ 1 & 0\\ -1 & 1 \end{bmatrix} \end{align}\qquad$$
Since the reduction step amounts to multiplying by a matrix with determinant $=1$ it preserves the determinant (our loop invariant $I)$. Thus the algorithm may be viewed as column subtractions on these matrices, which amounts to augmenting the subtractive Euclidean algorithm (on rows $[x,y])$ to $\,2\times 2\,$ matrices, by appending a 2nd row $[-u,v]$ for the coordinates of the given number.
Finally, to answer your question about the invariant at termination. note that in the OP we have $\,v_0,u_0 = n,m\,$ so $\, k = 2mn,\,$ so when the algorithm terminates with $\,x = y = d\,$ we have
$$\begin{align} \overbrace{\color{#c00}d(v\!+\!u)}^{\textstyle \color{#c00}xv+\color{#c00}yu} &=\, 2mn\\[.4em] \iff \frac{v+u}2 &= \frac{mn}d = {\rm lcm}(m,n)\end{align}\qquad$$
where we used $\,d\,{\rm lcm}(m,n) = \gcd(m,n)\,{\rm lcm}(m,n) = mn\,$ by here.
Remark $ $ To learn more about the linear algebra at the heart of the matter see the Hermite - Smith normal form and closely related topics.
