We present a monolithic preconditioned iterative solver for implicit discrete ordinate equations of three-dimensional grey thermal radiation transport and parallel codes are developed. A strategy of assembling linear algebraic systems is used to obtain radiation intensity in all discrete directions simultaneously. With preconditioned Krylov subspace iterative methods, the solver avoids possible mesh cycles in complex grids associated with sweep algorithms, which improves robustness and computational efficiency. First order upwind finite volume scheme is used for space discretization. Numerical experiments verify convergence rate on distorted hexahedral grids and assess performance of preconditioned iterative methods. Problems with coupled radiation and matter are simulated. Simulation results of three-dimensional crooked pipe and hohlraum problems are shown. It shows validity of the codes and flexibility of the method.