A univariate compactly supported refinable function $f$ can always be factored into $B_kast f$, with $B_k$ the B-spline of order $k$, $f$ a compactly supported distribution, and $k$ the approximation orders provided by the underlying shift-invariant space $S(f)$. Factorizations of univariate refinable {it vectors} $ff$ were also studied and utilized in the literature. One of the by-products of this article is a rigorous analysis of that factorization notion, including, possibly, the first precise definition of that process. The main goal of this article is the introduction of a special factorization algorithm of refinable vectors that generalizes the scalar case as closely (and unexpectedly) as possible: the original vector $ff$ is shown to be `almost' in the form $B_kastfff$, with $fff$ still compactly supported and refinable, and $k$ the approximation order of $S(ff)$: `almost' in the sense that $ff$ and $B_kastfff$ differ at most in one entry. The algorithm guarantees $fff$ to retain the possible favorable properties of $ff$, such as the stability of the shifts of $ff$ and/or the polynomiality of the mask symbol. At the same time, the theory and the algorithm are derived under relatively mild conditions and, in particular, apply to $ff$ whose shifts are not stable, as well as to refinable vectors which are not compactly supported. The usefulness of this specific factorization for the study of the smoothness of FSI wavelets (known also as `multiwavelets' and `multiple wavelets') is explained. The analysis invokes in an essential way the theory of finitely generated shift-invariant (FSI) spaces, and, in particular, the tool of {it superfunction theory}