Consider an undirected weighted graph $G = (V,E,w)$. We study the problem of computing $(1+ε)$-approximate shortest paths for $S \times V$, for a subset $S \subseteq V$ of $|S| = n^r$ sources, for some $0 < r \le 1$. We devise a significantly improved algorithm for this problem in the entire range of parameter $r$, in both the classical centralized and the parallel (PRAM) models of computation, and in a wide range of $r$ in the distributed (Congested Clique) model. Specifically, our centralized algorithm for this problem requires time $\tilde{O}(|E| \cdot n^{o(1)} + n^{ω(r)})$, where $n^{ω(r)}$ is the time required to multiply an $n^r \times n$ matrix by an $n \times n$ one. Our PRAM algorithm has polylogarithmic time $(\log n)^{O(1/ρ)}$, and its work complexity is $\tilde{O}(|E| \cdot n^ρ+ n^{ω(r)})$, for any arbitrarily small constant $ρ>0$. In particular, for $r \le 0.313\ldots$, our centralized algorithm computes $S \times V$ $(1+ε)$-approximate shortest paths in $n^{2 + o(1)}$ time. Our PRAM polylogarithmic-time algorithm has work complexity $O(|E| \cdot n^ρ+ n^{2+o(1)})$, for any arbitrarily small constant $ρ>0$. Previously existing solutions either require centralized time/parallel work of $O(|E| \cdot |S|)$ or provide much weaker approximation guarantees. In the Congested Clique model, our algorithm solves the problem in polylogarithmic time for $|S| = n^r$ sources, for $r \le 0.655$, while previous state-of-the-art algorithms did so only for $r \le 1/2$. Moreover, it improves previous bounds for all $r > 1/2$. For unweighted graphs, the running time is improved further to $poly(\log\log n)$.