  Consider an undirected weighted graph $G = (V,E,w)$. We study the problem of
computing $(1+\epsilon)$-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^{\omega(r)})$, where
$n^{\omega(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/\rho)}$, and its work complexity is $\tilde{O}(|E| \cdot n^\rho +
n^{\omega(r)})$, for any arbitrarily small constant $\rho >0$.
  In particular, for $r \le 0.313\ldots$, our centralized algorithm computes $S
\times V$ $(1+\epsilon)$-approximate shortest paths in $n^{2 + o(1)}$ time. Our
PRAM polylogarithmic-time algorithm has work complexity $O(|E| \cdot n^\rho +
n^{2+o(1)})$, for any arbitrarily small constant $\rho >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)$.
