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| = n^r 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| ⋅ n^{o(1)} + n^{ω(r)}), where n^{ω(r)} is the time required to multiply an n^r × 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 n^{2 + o(1)} time. Our PRAM polylogarithmic-time algorithm has work complexity O(|E| ⋅ n^ρ + n^{2+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| = n^r 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.