The objective of this paper is to introduce a general procedure for deriving interpolatory surface subdivision schemes with “symmetric subdivision templates” (SSTs) for regular vertices. While the precise definition of “symmetry” will be clarified in the paper, the property of SSTs is instrumental to facilitate application of the standard procedure for finding symmetric weights for taking weighted averages to accommodate extraordinary (or irregular) vertices in surface subdivisions, a topic to be studied in a continuation paper. By allowing the use of matrices as weights, the SSTs introduced in this paper may be constructed to overcome the size barrier limited to scalar-valued interpolatory subdivision templates, and thus avoiding the unnecessary surface oscillation artifacts. On the other hand, while the old vertices in a (scalar) interpolatory subdivision scheme do not require a subdivision template, we will see that this is not the case for the matrix-valued setting. Here, we employ the same definition of interpolation subdivisions as in the usual scalar consideration, simply by requiring the old vertices to be stationary in the definition of matrix-valued interpolatory subdivisions. Hence, there would be another complication when the templates are extended to accommodate extraordinary vertices if the template sizes are not small. In this paper, we show that even for C 2 interpolatory subdivisions, only one “ring” is sufficient in general, for both old and new vertices. For example, for 1-to-4 split C 2 interpolatory surface subdivisions, we obtain matrix-valued symmetric interpolatory subdivision templates (SISTs) for both triangular and quadrilateral meshes with sizes that agree with those of the Loop and Catmull-Clark schemes, respectively. Matrixvalued SISTs of similar sizes are also constructed for C 2 interpolatory √ 3 and √ 2 subdivision schemes in this paper. In addition to small template sizes, an obvious feature of matrix-valued weights is the flexibility for introducing shape-control parameters. Another significance is that, in contrast to the usual scalar setting, matrix-valued SISTs can be formulated in terms of the coefficient sequence of some vector refinement equation of interpolating bivariate C 2 splines with small support. For example, by modifying the spline function vectors introduced in our previous work [3, 6], C 2 symmetric interpolatory subdivision schemes associated with refinement equations of C 2 cubic and quartic splines on the 6-directional and 4-directional meshes, respectively, are also constructed in this paper.