In this paper, the coordinate descent (CD) method integrated with k-means clustering is proposed to address the least-squares problem of the large-scale matrix equation $AX=B$, where the coefficient matrix $A$ is assumed to be of full column rank. Rigorous theoretical analysis shows that the iteration sequence generated by this method converges to the unique least-norm least-squares solution of the problem. The performed numerical results demonstrate that the method developed here is feasible and efficient.