Smoothness and symmetry are two important properties of a refinable function. It is known that the Sobolev smoothness exponent of a refinable function can be estimated by computing the spectral radius of certain finite matrix which is generated from a mask. However, the increase of dimension and the support of a mask tremendously increases the size of the matrix and therefore make the computation very expensive. In this paper, we shall present a simple algorithm to efficiently numerically compute the smoothness exponent of a symmetric refinable function with a general dilation matrix. By taking into account of symmetry of a refinable function, our algorithm greatly reduces the size of the matrix and enables us to numerically compute the Sobolev smoothness exponents of a large class of symmetric refinable functions. Step by step numerically stable algorithms and details of the numerical implementation are given. To illustrate our results by performing some numerical experiments, we construct a family of dyadic interpolatory masks in any dimension and we compute the smoothness exponents of their refinable functions in dimension three. Several examples will also be presented for computing smoothness exponents of symmetric refinable functions on the quincunx lattice and on the hexagonal lattice.