An elastoplastic damage model is established based on Mohr-Coulomb(M-C) yield criterion considering the plastic deformation mechanism and stiffness degradation of rock material in the framework of elastoplastic damage theory. The internal variable- equivalent plastic strain is used to characterize the evolution of damage variable of rock material. The M-C yield criterion is a six-pyramid in stress space, stress update exist a "singularity" problem on the six corners and ridge of a pyramid in the process of the numerical implementation. Angular point smoothing method is used to deal with this problem; but it inevitable lead to approximate results. The implicit return mapping algorithm of elastoplastic damage constitutive equation is deduced on the basis of the M-C constitutive numerical integral algorithm in principal stress space. It includes three main computational steps of elastic trial, plasticity correction and damage correction. Three conditions of flow vector returns to the main plane, ridge and sharp point are discussed in the plastic correction steps respectively. To solve the "singularity" problem from the angle of principal stress space, the elastoplastic damage constitutive solver(RDM-C) is compiled using C++ language with object-oriented programming method. The uniaxial compression test, the foundation and cavity example are used to analyze and compare the results of calculation. The results show that the established elastoplastic damage constitutive model can well describe the main mechanics and deformation characteristics of quasi-brittle material, plastic zone and damage zone change trend; the program can carry out numerical analysis for geotechnical engineering problems, so as to give certain guidance to the site construction.