When a discretized Reynolds equation is to be solved iteratively at least three subjects have to be determined first. These are the iterative solution method, the size of gridwork for the numerical model, and the stopping criterion for the iterative computing. The truncation error analysis of the Reynolds equation is used to provide the stopping criterion, as well as to estimate an adequate grid size based on a required relative precision or grid convergence index. In the simulated lubrication analyses, the convergent rate of the solution is further improved by combining a simple multilevel computing, the modified Chebyshev acceleration, and multithreaded computing. The best case is obtained by using the parallel three-level red-black successive-over-relaxation (SOR) with Chebyshev acceleration. The speedups of the best case relative to the case using sequential SOR with optimal relaxation factor are around 210 and 135, respectively, for the slider and journal bearing simulations.