Matrix algorithms for modeling acoustic waves in piezoelectric multilayers are presented. All the algorithms considered are capable of resolving the numerical instability of transfer matrix at high frequency-thickness product. The formulation of basic matrices for the algorithm building blocks, and the development of recursive algorithms for the stack matrices, are systematically presented for both the conventional scattering and impedance matrices as well as the more recent hybrid matrix. Many variants of the algorithms are discussed, along with their respective usefulness and deficiency. Comparisons are made in their computational efficiency and numerical stability. For unconditional stability throughout large and small thicknesses, both scattering and hybrid matrix algorithms are applicable. For most efficiency, the algorithms that synergize both scattering and hybrid or impedance submatrices are superior, using surface matrix approach. Other aspects of algorithms such as formula conciseness, physical insight, versatility in incorporating boundary conditions, etc., are also noted.