Particle-resolved direct numerical flow solvers predominantly use a projection method to decouple the non-linear mass and momentum conservation equations. The computing performance of such solvers often decays beyond $\mathcal{O}(1000)$ cores due to the cost of solving at least one large three-dimensional pressure Poisson problem per time step. The parallelization may perform moderately well only or even poorly sometimes despite using an efficient algebraic multigrid preconditioner [38]. We present an accurate and scalable solver using a direction splitting algorithm [12] to transform all three-dimensional parabolic/elliptic problems (and in particular the elliptic pressure Poisson problem) into a sequence of three one-dimensional parabolic sub-problems, thus improving its scalability up to multiple thousands of cores. We employ this algorithm to solve mass and momentum conservation equations in flows laden with fixed non-spherical rigid bodies. We consider the presence of rigid bodies on the (uniform or non-uniform) fixed Cartesian fluid grid by modifying the diffusion and divergence stencils on the impacted grid node near the rigid body boundary. Compared to [12], we use a higher-order interpolation scheme for the velocity field to maintain a second-order stress estimation on the particle boundary, resulting in more accurate dimensionless coefficients such as drag $C_d$ and lift $C_l.$ We also correct the interpolation scheme due to the presence of any nearby particle to maintain an acceptable accuracy, making the solver robust even when particles are densely packed in a sub-region of the computational domain. We present classical validation tests involving a single or multiple (up to $\mathcal{O}(1000))$ rigid bodies and assess the robustness, accuracy and computing speed of the solver. We further show that the Direction Splitting solver is ∼5 times faster on 5120 cores than our solver [38] based on a classical projection method [5].