In this study, we investigate the sedimentation of rigid spheroidal particles suspensions in Stokes flow at volume fractions ranging from dilute ($\phi = 0.1$) up to dense regimes ($\phi = 0.6$) and Reynolds numbers $Re < 1$. While the sedimentation of spherical particles is well documented, showing a monotonic decrease of the average settling velocity with increasing volume fraction, the behavior of spheroidal suspensions remains sparsely studied at low Galileo number and is mostly limited to dilute conditions. Fully resolved simulations are performed using LEDDS, an open source lattice Boltzmann discrete element method (LBM-DEM) framework accelerated on GPUs. Large scale simulations involving more than $10^4$ particles allow collective dynamics to emerge and provide statistically converged quantities. Both prolate and oblate spheroids are considered, and sedimentation is studied in vertical columns as well as in fully periodic domains to access long time dynamics. Particular attention is paid to the generation of initial configurations in order to assess their influence on sedimentation regimes. The study focuses on the evolution of particle orientations, average settling velocities, and velocity fluctuations as functions of volume fraction and aspect ratio. In dilute regimes, spheroidal particles may exhibit settling velocities exceeding that of an isolated particle, whereas at higher concentrations hindrance effects are expected to dominate, similarly to spherical suspensions but with anisotropy induced deviations. The present work aims to quantify these effects in dense regimes and to identify scaling laws linking structural organization, orientation statistics, and macroscopic settling behavior. This work demonstrates how large-scale, fully resolved simulations capture the collective dynamics of spheroidal particle sedimentation in viscous flows, providing a quantitative reference for dense suspensions and a basis for future studies of polydisperse systems and collective sedimentation instabilities