Let $A$ be a dilation matrix, an$n times n$ expansive matrix that maps a full-rank lattice $Gamma subset R^n$ into itself. Let $Lambda$ be a finite subset of$Gamma$, and for $k in Lambda$ let $c_k$ be $r times r$ complex matrices. The refinement equation corresponding to $A$,$Gamma$, $Lambda$, and $c = set{c_k}_{k in Lambda}$ is $f(x) = sum_{k in Lambda} c_k , f(Ax-k)$. A solution $f colon R^n to C^r$, if one exists, is called a refinable vector function or a vector scaling function of multiplicity $r$. In this manuscript we characterize the existence of compactly supported $L^p$ or continuous solutions of the refinement equation, in terms of the $p$-norm joint spectral radius of a finite set of finite matrices determined by the coefficients $c_k$.We obtain sufficient conditions for the $L^p$ convergence ($1 le p le infty$) of the Cascade Algorithm $f^{(i+1)}(x) = sum_{k in Lambda} c_k , f^{(i)}(Ax-k)$, and necessary conditions for the uniform convergence of the Cascade Algorithm to a continuous solution. We also characterize those compactly supported vector scaling functions which give rise to a multiresolution analysis for $L^2(R^n)$ of multiplicity $r$, and provide conditions under which there exist corresponding multiwavelets whose dilations and translations form an orthonormal basis for $L^2(R^n)$.