In iterative methods, are matrix decompositions considered useful for implementation?

When we study an iterative method from textbooks, for example, see the Gauss-Seidel Method, the given matrix is decomposed with suitable splittings. In the example, $$A = L+U$$. So we can proceed with that decomposition if we solve the exercise through pen and paper.

I don't understand, why when we implement that method on a computer, that decomposition is not considered anymore?

For example, on this page- implementation with Matlab, it is not used the decomposition introduced before. Why? Why we have to struggle to find a suitable decomposition that will not be used within the code?

In Gauss-Seidel, you are using this $$A=L+U$$ splitting implicitly. So, you never form $$L$$, $$U$$, or $$L^{-1}$$ explicitly. Which is extremely good, since forming an additional matrix (not even talking about a calculation of an explicit inverse) is a huge burden.

Instead, since $$L$$ is lower triangular, you can change the explicit inverse of $$L$$, by performing a row-by-row forward subsitution.

Notice, that solving for (via the forward substitution when $$L$$ is lower-triangular):

$$$$L\mathbf{y}=\mathbf{b} \label{1} \tag{1}$$$$

is theoretically the same as performing a matrix-vector product: $$$$\mathbf{y}=L^{-1}\mathbf{b} \label{2} \tag{2}$$$$ However, numerically you always want $$\eqref{1}$$, not $$\eqref{2}$$. The reasons are simple: computation of the matrix inverse is numerically unstable and has a huge cost. (see Q1 (especially this answer), Q2 for some additional details).

With that in mind, take a look at the expression from Wikipedia article on Gauss-Seidel. The iteration: $$\mathbf{x}^{(k+1)}=L^{-1}(\mathbf{b}-U\mathbf{x}^{(k)})$$ is totally equivalent to a for-loop for $$i$$ $$$$x^{(k+1)}_i=\frac{1}{a_{ii}}\left(b_i-\sum\limits_{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}\right)-\frac{1}{a_{ii}}\left(\sum\limits_{j=i+1}^{n}a_{ij}x_j^{(k)}\right) \label{3} \tag{3}$$$$

Here, $$a_{ij}$$ are the entries of the original matrix $$A=L+U$$, $$k$$ is the iteration number, $$n$$ is the dimension of the square matrix $$A$$. By looking at $$\eqref{3}$$, you still see those $$L$$ and $$U$$; however, they are present only implicitly, and all operations are done with the entries of the original matrix $$A$$ stored exactly as it was before.

side-note:

Matrix decomposition usually implies factorization, so $$A=L+U$$ is not really (strictly) one according to the definition.

$$A=L+U$$ splitting (used in this question) has nothing to do with a more famous $$A=LU$$ LU decomposition.