A new approach for the solution of the Holstein-Biberman equation based on the advanced matrix method is developed. It allows for the consideration of radiation trapping in arbitrary finite 3D plasma domains for the various shapes of line contours and in a wide range of optical depths. Homogeneous and inhomogeneous distributions of absorbing atoms are considered. To solve the equation, an arbitrary plasma domain is discretized on a Cartesian voxel grid. The distances between the cells which are crossed by photons are computed by means of an efficient ray traversal algorithm. The algorithm is optimized for parallel computation on a graphical processing unit (GPU). For the Lorentzian shape of emission and absorption lines, the analytical expressions (which significantly decrease the computation time) have been derived. In the high opacity limit, the matrix is transformed to the universal form with an escape factor as a multiplier. The method is validated against a previously developed matrix approach by comparing the solutions for a finite cylinder geometry. The applicability range of the old method is specified. This range is defined by the asymptotics of Lorentz line wings at high optical depths. The capability of the method is illustrated with several complex geometries which are typical for various plasma sources. The effects related to the presence of photon blocking barriers are demonstrated. The proposed method allows for the demonstration of the fundamental differences between radiation and diffusion transport processes in the plasma domains of a complex shape. The method can be integrated into multi-component collisional-radiative models.