Before revision:  We devise a significantly improved algorithm for this problem in the entire range of parameter r, in both the classical centralized and in the parallel (PRAM) models of computation  . Specifically, our centralized algorithm for this problem requires time 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.
After revision:  We devise a significantly improved algorithm for this problem in the entire range of parameter r, in both the classical centralized and in 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 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.