The industry-standard constrained pressure residual (CPR) algorithm is often able to effectively improve the robustness behavior and the convergence speed of linear iterations for isothermal reservoir simulation. In this paper, we present and study an improved extension of CPR to the constrained pressure-temperature residual (CPTR) version for non-isothermal reservoir problems in heterogeneous porous media. In the proposed preconditioner, the corresponding approximations for the inverse of matrices are computed under a domain decomposition framework by using the restricted additive Schwarz (RAS) algorithm, to equally deal with the coupled thermal-pressure-saturation reservoir system and highly exploit the parallelism of supercomputer platforms. Moreover, we introduce and develop a family of multilevel CPTR preconditioners with suitable coarse grid corrections, to further improve the applicability of this two-stage preconditioner for large-scale computation. Numerical results for strong heterogeneous flow problems show that the new approach can dramatically improve the convergence of linear iterations, and demonstrate the superiority of CPTR over the commonly used RAS preconditioners. The parallel scalability of the non-isothermal reservoir simulator is also studied versus a supercomputer with tens of thousands of processors.