In　this　paper,　a　general　algorithm　is　presented　for　the　integration　of　multisurface　plasticity　models.　The　algorithm　combines　the　return　mapping　with　the　technique　of　the　sequential　linear　complementarity　problem　(SLCP)　to　update　stress.　The　set　of　active　constraints　is　naturally　identified　by　the　classical　Lemke＇s　algorithm.　With　the　assumptions　of　isotropic　linear　elasticity,　perfect　plasticity,　and　the　associated　flow　rule,　all　details　are　provided　in　matrix　notations　to　facilitate　computer　implementation.　The　extension　to　hardening/softening　multisurface　plasticity　models　is　also　presented.　The　application　of　the　algorithm　is　demonstrated　via　simulation　of　three　types　of　geotechnical　problems　in　2D　and　3D.　Both　linear　and　nonlinear　multisurface　plasticity　models,　e.g.　Mohr-Coulomb　and　generalized　Hoek-Brown　yield　criteria,　are　examined　within　the　framework　of　the　proposed　algorithm.　The　numerical　results　are　in　good　agreement　with　the　analytic　solutions.　Moreover,　the　accuracy　of　the　proposed　stress　integration　procedure　is　investigated　through　3D　isoerror　map.　The　convergence　using　the　consistent　tangent　matrix　at　the　global　level　is　examined　by　a　one-increment　example　consisting　of　one　element.　(c)　2021　Elsevier　Inc.　All　rights　reserved.