The main method for calculating the gravity field involves discretizing the density sources into a stack of rectangular prisms with a regular grid distribution. The analytical formulation of the gravity anomaly for a right-angled rectangular prism is affected by depth, with the kernel function decaying as depth increases. In addition, the efficiency of the computation and the storage requirements often pose challenges. We present a fast computational method for three-dimensional gravity forward modelling of subsurface space using rectilinear grid and apply the Block–Toeplitz Toeplitz–Block method to the rectilinear grid. The size of the upright rectangles increases with depth to offset the effect of depth. We assume that the observation points are distributed on a homogeneous grid, and the kernel sensitivity matrices exhibit a Block–Toeplitz Toeplitz–Block structure, which is symmetric. For rectilinear dissections of subsurface space in MATLAB, the logarithmic interval size is used. The rectilinear mesh can offset the effect of depth to some degree allowing gravity anomalies to decrease more quickly. For the test of a single model, the gravity anomalies decrease faster and more rapidly in the case of the rectilinear grid compared to the uniform grid. In addition to this, we performed simulations on more complex models and demonstrated that using the Block–Toeplitz Toeplitz–Block method on this basis greatly improves the computational efficiency.