Various approximations combining coupled-cluster (CC) and many-body perturbation theories have been derived and implemented into the parallel execution programs that take into account the spin, spatial (real Abelian), and permutation symmetries and that are applicable to closed- and open-shell molecules. The implemented models range from the CCSD(T), CCSD[T], CCSD(2)(T), CCSD(2)(TQ), and CCSDT(2)(Q) methods to the completely renormalized (CR) CCSD(T) and CCSD[T] approaches, where CCSD (CCSDT) stands for the CC method with connected single and double (single, double, and triple) cluster operators, and subscripted or parenthesized 2, T, and Q indicate the perturbation order or the excitation ranks of the cluster operators included in the corrections. The derivation and computer implementation have been automated by the algebraic and symbolic manipulation program TENSOR CONTRACTION ENGINE (TCE). The TCE-synthesized subroutines generate the tensors with the highest excitation rank in a blockwise manner so that they need not be stored in their entirety, while enabling the efficient reuse of other precalculated intermediate tensors defined by prioritizing the memory optimization as well as operation minimization. Consequently, the overall storage requirements for the corrections due to connected triple and quadruple cluster operators scale as O(n(4)) and O(n(6)), respectively (n being a measure of the system size). For systems with modest multireference character of their wave functions, we found that the order of accuracy is CCSD