In this paper, we consider solving a class of matrix inequality constrained matrix least squares problems of the formwhere · is the Frobenius norm, matrices A i ∈ R l×m , B i ∈ R n×s (i = 1, . . . , t), C ∈ R l×s , E ∈ R p×m , F ∈ R n×q and L , U ∈ R p×q are given. An inexact version of alternating direction method (ADM) with truly implementable inexactness criteria is proposed for solving this problem and its several reduced versions which are applicable in image restoration. Numerical experiments are performed to illustrate the feasibility and efficiency of the proposed algorithm, including when the algorithm is tested with randomly generated data and on some image restoration problems. Comparisons with some existing methods (with necessary modifications) are also given.