What "reducing a matrix" means
One issue here is what is meant by "reducing a matrix". The only purpose I can see for reducing a matrix is to decide if the matrix has some specific property. Whatever property you are interested in has to remain unchanged under the operations you do or you have not succeeded in your purpose.
If you investigate Serge Lang's book closely, you will find that he had a specific purpose for reducing the matrix,(to find the rank, discussed below) which allows him to use column operations to achieve his aim.
It is important to note that for most people, the phrase "reducing a matrix" refers specifically to finding the Reduced Row Echelon Form (also known as RREF). As the name implies, RREF is defined using the rows of the matrix:
1. The leftmost nonzero entry in any row is a 1 (called a "leading 1").
2. Leading 1's on lower rows are further to the right than leading 1's on higher rows.
3. Other entries in a column with a leading 1 are zero.
4. All rows consisting entirely of zeros are at the bottom of the matrix.
Since you are defining this matrix in terms of rows, you must do row operations to achieve it. In fact, if you do use row operations, then given a particular starting matrix (which may or may not be square, by the way), there is precisely one RREF that it will reduce to. If you do column operations you will not arrive at this same matrix.
Relationships to solutions of equations
The point others have made about finding solutions is valid. Every matrix can be used to represent the coefficients of a system of linear equations, regardless of whether you want it to or not. Many of the important properties of matrices are strongly related to this system of linear equations, and so it is a good idea to have methods that preserve these properties.
For example, the set of solutions to $A\mathbf{x} = \mathbf{0}$ is the same as the set of solutions to the (homogeneous) linear equations, and the nonzero rows of the RREF form a basis for the space perpendicular to this space of solutions. Neither of these properties are preserved by column operations, but they are preserved by row operations.
Most importantly, just to reiterate: the solution to your homogeneous linear equations is definitely not preserved by column operations, and this is an important consideration regardless of whether you want your matrix to represent linear equations or not.
Note that there are things that column operations do preserve, such as the range of the linear transformation defined by your matrix. However a mixture of row and column operations will not preserve this, nor most of the properties you care about.
The inverse of a square matrix
One property that a mixture of row and column operations does preserve is invertibility. You can see this from the idea of elementary matrices.
Doing elementary row operations corresponds to multiplying on the left by an elementary matrix. For example, the row operation of "new R2 = R2 - 3R1" is produced on a 3 by n matrix when you multiply on the left by $\begin{pmatrix} 1 & 0 & 0 \\ -3 & 1 & 0\\ 0 & 0 & 1 \end{pmatrix}$. Column operations, on the other hand, are produced when you multiply by a matrix on the right hand side. The column operation of "new C2 = C2 - 3C1" is produced on an m by 3 matrix when you multiply on the right by $\begin{pmatrix} 1 & -3 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{pmatrix}$.
The process of finding the inverse of a square matrix by augmenting an identity and doing matching row operations on both works because of these elementary matrices. If you perform row operations on $A$ to get $I$ and the elementary matrices corresponding to these operations are $E_1$, ..., $E_k$, then we have $E_k ... E_1 A = I$ and so the inverse of $A$ must be $E_k ... E_1$. But this is equal to $E_k ... E_1 I$, which is what you get when you do all those row operations on the identity!
This would work perfectly well if you did column operations, however a mixture of the two is a bit harder to understand. Suppose you did row operations with matrices $E_1$, ..., $E_k$ and column operations with matrices $F_1$, ..., $F_h$ and you produced the identity as a result. Then
$$
\begin{align}
E_k \dots E_1 A F_1 \dots F_h &= I\\
A &= (E_1)^{-1} \dots (E_k)^{-1} I (F_h)^{-1} \dots (F_1)^{-1}\\
A^{-1} &= F_1 \dots F_h I E_k \dots E_1
\end{align}
$$
This is not the result of doing the matching column operations and row operations on $I$ because the matrices are on the wrong side!
So if you can row-and-or-column reduce a matrix to the identity, then this proves that the matrix is invertible, but it is tricky to tell what the inverse actually is unless you do only row or only column operations.
Finally note that row and column operations simultaneously will allow you to calculate the determinant of a square matrix easily, and so there is a certain advantage to doing this. In particular a row/column operation of the type "new Ri = Ri + k Rj" or "new Ci = Ci + k Cj" will not change the determinant, so if you restrict yourself to those operations, you can get your matrix into a form where it is clear what the determinant is more quickly than restricting yourself to just one. However, it is worth stressing that this relates to my original comment about the purpose you are trying to achieve -- this is only really an advantage in the situation where you are finding determinants.
The rank of a matrix
Finally, this row-and-or-column reduction also does preserve the rank of your matrix. If you do this kind of reduction on any matrix, what you are guaranteed to produce eventually if you try is a matrix with some 1's down its main diagonal (but not necessarily all the way down) and zeros everywhere else. Something like these:
$$
\begin{pmatrix} 1 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0& 0\end{pmatrix} \quad \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0& 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}
$$
The number of 1's on that diagonal is the rank of your matrix. That is, it is the dimension of the column space and row space, and the dimension of the range of the linear transformation defined by your matrix, and n minus the dimension of the solution space to your system of equations. So the two matrices above have rank 3 and 2 respectively.
However, there is nothing in the final reduced form that tells you what the actual solutions to the equations are or what the bases for the spaces are. These can all be changed in ways that are not easy to follow if you do both row and column operations. If you did just column operations, then the leading-1 columns would be a basis for the range of the linear transformation, but with both row and column operations, these columns are all standard basis vectors and won't necessarily be correct. If you did just row operations, then the non-leading-1 columns would appear as the first several coordinates of the vectors in a basis for the null space of the matrix, but with both row and column operations these columns are all zero.
Summary
- Most people mean Reduced Row Echelon Form when they say "reduced matrix", which is defined by doing only row operations.
- Row operations preserve many useful properties of matrices and tell you certain information about matrices.
- Column operations preserve some useful properties of matrices and tell you certain information about matrices.
- Mixtures of row and column operations preserve only a small number of properties. In particular, a mixture will preserve the rank and the invertibility of a matrix (but will not provide bases for associated subspaces or the actual inverse itself).
Best Answer
Yes, if $\ k_1,k_2,\dots,k_n\in K\ $, with $\ \det(M)=\prod_\limits{i=1}^nk_i\ $, then there exists a sequence, $\ R_1,R_2,\dots,R_m\ $ of row operations of the specified form such that $$ R_mR_{m-1}\dots R_1M=\text{diag}\big(k_1,k_2,\dots,k_n\big)\ . $$ The key to doing this is that there exist sequences of row operations of the specified kind that will allow you to:
Here is how you can carry out the first of these operations: \begin{align} r_i'&=r_i+r_j\\ r_j'&=r_j-r_i'=-r_i\\ r_i''&=r_i'+r_j'=r_j\ . \end{align} After this sequence of operations, the new $\ i^\text{th}\ $ row, $\ r_i''=r_j\ ,$ will be the original $\ j^\text{th}\ $ row, and the new $\ j^\text{th}\ $ row, $\ r_j'={-}r_i\ ,$ will be the negative of the original $\ i^\text{th}\ $ row.
The second of the two above-mentioned operations can be performed as follows \begin{align} r_i'&=r_i+r_j\\ r_j'&=r_j+(\alpha-1)r_i'=\alpha r_j+(\alpha-1)r_i\\ r_i''&=r_i'-\alpha^{-1}r_j'=\alpha^{-1}r_i\\ r_j''&=r_j'-\alpha(1-\alpha)r_i''=\alpha r_j\ . \end{align} After this sequence of operations, the new $\ i^\text{th}\ $ row, $\ r_i''=\alpha^{-1}r_i\ $, will be the original $\ i^\text{th}\ $ row multiplied by $\ \alpha^{-1}\ $, and the new $\ j^\text{th}\ $ row, $\ r_j''=\alpha r_j\ $, will be the the original $\ j^\text{th}\ $ row multiplied by $\ \alpha\ $.
Since any matrix obtained by applying row or column operations of the specified form to $\ M\ $ must always have full rank, every such matrix must have at least one non-zero entry in every column and can never have a column which is a linear combination of other columns. Therefore, by subtracting suitable multiples of a row with a non-zero entry in the $\ j^\text{th}\ $ column from the other rows, you can reduce all the other entries in that column to zero. If the non-zero entry is not in the $\ j^\text{th}\ $ row you can get it there by applying an operation of the first kind described above. After doing this you will then have a matrix in which the only non-zero entry in the $\ j^\text{th}\ $ column is in the $\ j^\text{th}\ $ row. By applying this procedure successively to all the columns of $\ M\ $ (in any order you please) you can therefore reduce it to a matrix of the form $\ \text{diag}\big(\mu_1,\mu_2,\dots,\mu_n\big)\ $ with $\ \prod_\limits{i=1}^n\mu_i=\det(M)\ $. Note that as you apply this procedure you will never have a column whose only non-zero entries are in the same rows as those columns containing a single non-zero entry, because such a column would then be a linear combination of those latter columns.
You can now apply operations of the second kind described above to the rows of $\ \text{diag}\big(\mu_1,\mu_2,\dots,\mu_n\big)\ $ to multiply the $\ i^\text{th}\ $ row by $\ \frac{k_i}{\mu_i}\ $ and the $\ n^\text{th}\ $ row by $\ \frac{\mu_i}{k_i}\ $. The matrix you will then end up with is $$ \text{diag}\left(k_1,k_2,\dots,\frac{\prod_\limits{i=1}^n\mu_i}{\prod_\limits{i=1}^{n-1}k_i}\right)=\text{diag}\big(k_1,k_2,\dots,k_n\big)\ , $$ since $\ \prod_\limits{i=1}^n\mu_i=\det(M)\ $ and $\ \prod_\limits{i=1}^{n-1}k_i=\frac{\det(M)}{k_n}\ $.