I do not fully understand when and why I should organize the vectors in a matrix as rows and when and why as columns.
As long as you use the same logical convention throughout, the choice is arbitrary, because
$$\bbox[#ffffef, 1em]{ \mathbf{M} \mathbf{v} = \left( \mathbf{v}^T \mathbf{M}^T \right)^T }$$
That is, if you use column vectors, the vector is on the right and the matrix on the left in matrix-vector multiplication; if you use row vectors, the vector is on the left and the matrix transposed. The result only differs by transpose.
Thus, it is only a matter of convention; mathematically, you can use either one.
Do remember to ensure you use the same convention throughout, however. In most practical use cases, we use math to model phenomena or operations; as a language to describe the problem and the solution we would like to find. If we change the convention midway through, we break the mapping between the model and its mathematical description. Breakage ensues.
However, it is completely okay to switch conventions between problems. In fact, it is a good idea to choose the convention based on the context, and ease of notation.
For example, we do Gaussian elimination as "row reduction", because that yields nice, compact representation using our typical mathematical notation. You can do it transposed, as "column reduction", but the notation will be messier. In fact, if you do need to do column reduction, it is easier to first take the transpose of the matrix, then do row reduction (Gaussian elimination), and transpose the result!
This is a very common issue for programmers learning how to implement 3D transforms and projection using linear algebra. For simplicity and efficiency, 3D coordinates are usually described using homogenous coordinates, with an implicit fourth component always being 1. If we use column vectors, then
$$\bbox{ \mathbf{M} = \left [ \begin{matrix} x_x & y_x & z_x & t_x \\ x_y & y_y & z_y & t_y \\ x_z & y_z & z_z & t_z \\ 0 & 0 & 0 & 1 \end{matrix} \right ]}, \quad
\bbox{ \vec{v} = \left [ \begin{matrix} x \\ y \\ z \\ 1 \end{matrix} \right ]}, \quad
\bbox{ \mathbf{M} \vec{v} = \left [ \begin{matrix} x_x x + y_x y + z_x z + t_x \\ x_y x + y_y y + z_y z + t_y \\ x_z x + y_z y + z_z z + t_z \\ 1 \end{matrix} \right ]}$$
which can be described as a transformation and translation by
$$\bbox {
\hat{e}_x = \left [ \begin{matrix} x_x \\ x_y \\ x_z \end{matrix} \right ]
}, \quad \bbox {
\hat{e}_y = \left [ \begin{matrix} y_x \\ y_y \\ y_z \end{matrix} \right ]
}, \quad \bbox {
\hat{e}_z = \left [ \begin{matrix} z_x \\ z_y \\ z_z \end{matrix} \right ]
}, \quad \bbox {
\hat{t} = \left [ \begin{matrix} t_x \\ t_y \\ t_z \end{matrix} \right ]
}$$
where $x$, $y$, and $z$ axes are transformed to $\hat{e}_x$, $\hat{e}_y$, $\hat{e}_z$, and the result translated by $\hat{t}$.
However, the above applies to column vectors. This convention is common in linear algebra books and tutorials, but for example OpenGL uses row vectors instead. Then,
$$\bbox {
\vec{v} = \left [ \begin{matrix} x & y & z & 1 \end{matrix} \right ]
}, \quad \bbox {
\mathbf{M} = \left [ \begin{matrix} x_x & x_y & x_z & 0 \\ y_x & y_y & y_z & 0 \\ z_x & z_y & z_z & 0 \\ t_x & t_y & t_z & 1 \end{matrix} \right ]
}, \quad \bbox {
\vec{v} \mathbf{M} = \left [ \begin{matrix} x_x x + y_x y + z_x z + t_x & x_y x + y_y y + z_y z + t_y & x_z x + y_z y + z_z z + t_z & 1 \end{matrix} \right ]
}$$
Note how these conventions only differ by the transpose:
$$\bbox { \begin{aligned}
\hat{e}_x &= \left [ \begin{matrix} x_x & x_y & x_z \end{matrix} \right ] \\
\hat{e}_y &= \left [ \begin{matrix} y_x & y_y & y_z \end{matrix} \right ] \\
\hat{e}_z &= \left [ \begin{matrix} z_x & z_y & z_z \end{matrix} \right ] \\
\hat{t} &= \left [ \begin{matrix} t_x & t_y & t_z \end{matrix} \right ] \\
\end{aligned} }$$
The interpretation of the vectors and the matrix stays exactly the same. Only the notation differs (by transpose, and whether vector is to the left or to the right of the matrix in matrix-vector multiplication).
A further tripping point is that the upper left 3×3 part of the transformations (the rotation/reflection part) is very often a pure rotation, and therefore orthonormal:
$$\bbox { \begin{aligned}
\hat{e}_x \cdot \hat{e}_x &= 1 \\
\hat{e}_y \cdot \hat{e}_y &= 1 \\
\hat{e}_z \cdot \hat{e}_z &= 1 \\
\hat{e}_x \cdot \hat{e}_y &= 0 \\
\hat{e}_x \cdot \hat{e}_z &= 0 \\
\hat{e}_y \cdot \hat{e}_z &= 0 \\
\end{aligned} } \quad \quad \bbox { \begin{aligned}
\hat{e}_x \times \hat{e}_y &= \hat{e}_z \\
\hat{e}_y \times \hat{e}_z &= \hat{e}_x \\
\hat{e}_z \times \hat{e}_x &= \hat{e}_y \\
\hat{e}_x \times \hat{e}_z &= -\hat{e}_y \\
\hat{e}_z \times \hat{e}_y &= -\hat{e}_x \\
\hat{e}_y \times \hat{e}_x &= -\hat{e}_z \\
\end{aligned} }$$
In that case, the inverse of that 3×3 matrix is its transpose! In other words, if we have a pure 3×3 rotation matrix $\mathbf{R}$,
$$\bbox { \mathbf{R} = \left [ \begin{matrix}
x_x & y_x & z_x \\
x_y & y_y & z_y \\
x_z & y_z & z_z \\
\end{matrix} \right ] }, \quad
\bbox { \mathbf{R}^{-1} = \mathbf{R}^T = \left [ \begin{matrix}
x_x & x_y & x_z \\
y_x & y_y & y_z \\
z_x & z_y & z_z \\
\end{matrix} \right ] }$$
Because of this, if you have a pure rotation transform, and you get the convention (of whether row or column vectors are used) wrong, you get inverse rotations of what you expected.
It is an especially annoying tripping point for programmers, because the first reaction is obviously to suspect a programming error, perhaps a sign inversion.
Even when the transform is not a pure rotation, the basis vectors tend to be orthogonal to each other (their pairwise dot products zero), in which case there is a scaling factor applied to each axis, but otherwise the same applies.
