Euclidean distance matrix optimization with ordinal constraints (EDMOC) has found important applications in sensor network localization and molecular conformation. It can also be viewed as a matrix formulation of multidimensional scaling, which is to embed n points in a $r$-dimensional space such that the resulting distances follow the ordinal constraints. The ordinal constraints, though proved to be quite useful, may result in only zero solution when too many are added, leaving the feasibility of EDMOC as a question. In this paper, we first study the feasibility of EDMOC systematically. We show that if $rge n-2$, EDMOC always admits a nontrivial solution. Otherwise, it may have only zero solution. The latter interprets the numerical observations of crowding phenomenon. Next we overcome two obstacles in designing fast algorithms for EDMOC, i.e., the low-rankness and the potential huge number of ordinal constraints. We apply the technique developed to take the low rank constraint as the conditional positive semidefinite cone with rank cut. This leads to a majorization penalty approach. The ordinal constraints are left to the subproblem, which is exactly the weighted isotonic regression, and can be solved by the enhanced implementation of Pool Adjacent Violators Algorithm (PAVA). Extensive numerical results demonstrate {the} superior performance of the proposed approach over some state-of-the-art solvers.