An　elastic　domain　with　multiple　yield　surfaces　brings　about　huge　troubles　to　numerical　computations　due　to　singularities,　while　non-associative　plasticity　leads　further　to　the　open　issue　of　whether　constitutive　integration　is　well　posed.　This　study　reduces　the　constitutive　integration　of　non-associative　plasticity　with　multiple　yield　surfaces　to　a　mixed　complementarity　problem,　represented　by　MiCP,　a　special　case　of　finite-dimensional　variational　inequalities.　By　means　of　the　projection-　contraction　algorithm　for　finite-dimensional　variational　inequalities　and　the　idea　of　the　Gauss-Seidel　method,　a　new　projection-contraction　algorithm　for　MiCP　is　designed　and　denoted　by　GSPC.　Applying　the　monotonicity　of　the　mapping　of　the　MiCP,　GSPC　is　proved　convergent　theoretically　for　associative　plasticity.　For　non-associative　plasticity,　the　sufficient　condition　for　GSPC　to　be　convergent　is　also　established　if　the　tension　part　of　the　Mohr-Coulomb　elastic　domain　is　cut　off.　Typical　examples　are　designed　to　illustrate　GSPC　is　highly　efficient,　accurate　and　stable,　some　of　which　cannot　be　solved　by　the　conventional　return-mapping　method.　In　all　the　cases,　the　computational　efficiency　of　GSPC　is　apparently　above　the　mapping　return　method.　(C)　2019　Elsevier　B.Y.　All　rights　reserved.