Image reconstruction in diffuse optical tomography (DOT) is, in general, posed as a model-based, nonlinear optimization problem, which requires repeated use of the three-dimensional (3D) forward and inverse solvers. To cope with the computation and storage problem for some applications, such as breast tumor diagnosis, it is preferable to develop a subdomain-based parallel computation scheme. In this study, we propose a two-level image reconstruction scheme for 3D DOT, which combines the Schwarz-type domain-decomposition (DD)-based forward calculation and the matrix-decomposition (MD)-based inversion. In the forward calculation, the solution to the diffusion equation is initially obtained using a whole-domain finite difference method at a coarse grid, and then updated with a parallel DD scheme at a fine grid. The inversion procedure starts with the wavelet-decomposition-based reconstruction at a coarse grid, and then follows with a Levenberg–Marquardt least-squares solution at a fine grid, where an MD strategy is adopted for the relevant linear inversion. It is demonstrated that the combination of the DD-based forward solver and MD-based inversion allows for coarse-grain parallel implementation of both the forward and inverse issues and effectively reduces computation and storage loads for the large-scale problem. Also, both numerical simulations and phantom experiments show that MD-based linear inversion is superior to the row-fashioned algebraic reconstruction technique.