Large-scale fluid flow in porous media demands intense computations and occurs in the most diverse applications, including groundwater flow and oil recovery. This article presents novel computational strategies applied to reservoir geomechanics. Advances are proposed for the efficient assembly of finite element matrices and the solution of linear systems using highly vectorized code in MATLAB. In the CPU version, element matrix assembly is performed using conventional vectorization procedures, based on two strategies: the explicit matrices, and the multidimensional products. Further assembly of the global sparse matrix is achieved using the native sparse function. For the GPU version, computation of the complete set of element matrices is performed with the same strategies as the CPU approach, using gpuArray structures and the native CUDA support provided by MATLAB Parallel Computing Toolbox. Solution of the resulting linear system in CPU and GPU versions is obtained with two strategies using a one-way approach: the native conjugate gradient solver (pcg), and the one provided by the Eigen library. A broad discussion is presented in a dedicated benchmark, where the different strategies using CPU and GPU are compared in processing time and memory requirements. These analyses present significant speedups over serial codes. Finally, recommendations are provided.