Consider an undirected weighted graph G = (V, E, w). We study the problem of computing (1 + ϵ)approximate shortest paths for S × V , for a subset S ⊆ V of |S| = nr sources, for some 0 < r ≤ 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 Õ(|E| · no(1) + nω(r)), where nω(r) is the time required to multiply an nr × n matrix by an n × n one. Our PRAM algorithm has polylogarithmic time (log n)O(1/ρ), and its work complexity is Õ(|E| · nρ + nω(r)), for any arbitrarily small constant ρ > 0. In particular, for r ≤ 0.313 . . ., our centralized algorithm computes S × V (1 + ϵ)-approximate shortest paths in n2+o(1) time. Our PRAM polylogarithmic-time algorithm has work complexity O(|E| · nρ + n2+o(1)), for any arbitrarily small constant ρ > 0. Previously existing solutions either require centralized time/parallel work of O(|E| · |S|) or provide much weaker approximation guarantees. In the Congested Clique model, our algorithm solves the problem in polylogarithmic time for |S| = nr sources, for r ≤ 0.655, while previous state-of-the-art algorithms did so only for r ≤ 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) for r ≤ 0.655. Previously this running time was known for r ≤ 1/2.