The algorithm and implementation details of a new parallel CFD-DEM solver that uses both CPU and GPU resources for computations are presented in this study. The new solver supports complex geometries, various shapes of mesh, converges at very high fluid velocities and executes with available computational resources of a simple desktop computer. At first, verification studies were performed for the solver. Among three parts of the solver, the DEM solver and coupling solver results were verified by simple problems for which analytical solutions exist. Then, the solver results were evaluated in three cases of gas-solid flows with various geometries (from simple to complex) and various numbers of particles. For a large system that contains 870 k particles, it took around 6 h (with two CPU cores) to complete each second of simulation while in a smaller system with 47 k particles, it took only 30 min (with one CPU core). Results obtained by the new solver were compared with experimental data. The solver can properly predict the frequency of bubble passage and intensity/share of macrostructures in the bubbling regime, while underestimates intensity/share of meso-structures. The hydrodynamics of a spouted bed with draft tube was well-captured in terms of the fountain shape and its maximum height, particle concentration profile in different parts and time-averaged upward velocity of particles. The solver was also used to investigate the particle dynamics in a Wurster bed. Two velocity profiles were tested. In the first velocity profile, an unstable fountain with slugging conveying in the draft tube was observed while in the second velocity profile, a stable fountain with dispersed conveying in the draft tube. This was found to be attributed to the formation of the bubbles in the annulus region. (C) 2016 Elsevier B.V. All rights reserved.