Multiwavelets are generated from refinable function vectors by using multiresolution analysis. In this paper we investigate the approximation properties of a multivariate refinable function vector associated with a general dilation matrix in terms of both the subdivision operator and the order of sum rules satisfied by the matrix refinement mask. Based on a fact about the sum rules of biorthogonal multiwavelets, a coset by coset (CBC) algorithm is presented to construct biorthogonal multiwavelets with arbitrary order of vanishing moments. More precisely, to obtain biorthogonal multiwavelets, we have to construct primal and dual masks. Given any primal matrix mask $a$ and a general dilation matrix $M$, the proposed CBC algorithm reduces the construction of all dual masks of $a$, which satisfy the sum rules of arbitrary order, to a problem of solving a well organized system of linear equations. We prove in a constructive way that for any given primal mask $a$ with a dilation matrix $M$ and for any positive integer $k$, we can always construct a dual mask $\wt a$ of $a$ such that $\wt a$ satisfies the sum rules of order $k$. In addition, we provide a general way for the construction of Hermite interpolatory matrix masks in the univariate setting with any dilation factors. From such Hermite interpolatory masks, smooth Hermite interpolants, including the well known cubic Hermite splines as a special case, are obtained and are used to construct biorthogonal multiwavelets. As an example, a $C^3$ Hermite interpolant with support $[-3,3]$ is presented. Then we shall apply the CBC algorithm to such Hermite interpolatory masks to construct biorthogonal multiwavelets. Several examples of biorthogonal multiwavelets are provided to illustrate the general theory. In particular, a $C^1$ dual function vector with support $[-4,4]$ of the cubic Hermite splines is given.