Let \(A\) be a dilation matrix, an\(n \times n\) expansive matrix that maps a full-rank lattice \(\Gamma \subset \mathbf{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 = \{c_k\}_{k \in \Lambda}\) is \(f(x) = \sum_{k \in \Lambda} c_k \, f(Ax-k)\). A solution \(f \,\colon\, \mathbf{R}^n \to \mathbf{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(\mathbf{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(\mathbf{R}^n)\).