Abstract: High-fidelity neutron transport solutions for irregular and complex reactor cores remain a key challenge and focal point in reactor physics calculations. Traditional deterministic programs are incapable of finely resolving reactors with complex and non-standard geometries. Leveraging their superior geometric adaptability, parallel unstructured mesh solution methods are widely applied for large-scale numerical simulations of complex reactor core geometries. A three-dimensional multi-group parallel neutron transport solver on tetrahedral unstructured meshes developed on the COME coupling platform by the Virtual Reactor Coupling Analysis Laboratory at Huazhong University of Science and Technology, was presented in this study. TETRIS utilizes discrete ordinates method for angular discretization and discontinuous finite elements method for spatial discretization. Owing to the inhomogeneous distribution of materials in complex geometric reactor cores, traditional source iteration methods face significant challenges in solving neutron transport equations efficiently. To accelerate neutron transport calculations and reduce the difficulty of solving the transport equation, the TETRIS was implemented a Krylov subspace acceleration method based on a matrix-free approach on the COME platform. Moreover, due to the extensive use of unstructured meshes required for detailed modeling of nuclear reactors, which demands substantial memory and prolongs computation time, the TETRIS was implemented a parallel sweeping scheme using MPI for spatial variables, alongside MPI-based parallel resolution for energy groups, and OpenMP for angular directions to reduce both computational time and memory requirements significantly, the parallel sweeping scheme employs a pipeline sweeping algorithm, while the multi-region parallel sweeping models utilizes a block-Jacobi parallel strategy. The computational results for the PinCell model, hexagonal assembly model, C5G7 benchmark model, and KNK-Ⅱ experimental fast neutron reactor model show that TETRIS exhibits commendable computational accuracy, parallel efficiency and geometric applicability, achieving numerical accuracy comparable to that of multi-group Monte Carlo methods. For the PinCell model, the deviation in k_eff is 61 pcm, while the hexagonal geometry model shows a deviation of 86 pcm. For the KNK model under two conditions, the maximum deviation in k_eff is 39 pcm, with the maximum relative deviation in power distribution being 0.242%. In the C5G7 2D model, the deviation in k_eff is 13 pcm, and the maximum relative deviation in power distribution is 1.4%. In summary, the effective multiplication factor deviations for all four test cases are within 100 pcm, with relative power error below 1.4%. The computational time for different parallel strategies indicates that the TETRIS achieves parallel efficiency of over 50% in most cases when utilizing over 100 cores. The acceleration factor achieved by the Krylov subspace method reaches a maximum of 3.47.