Abstract: A 3D neutron transport computational code, STORK, has been developed based on a small-scale multi-GPU computing platform, utilizing the coupled approach of the two-dimensional full-core layer-by-layer transport calculation by the method of characteristics (MOC) and the 3D pin-by-pin simplified P3 (SP3) calculation. In this code, firstly, the core was layered according to the axial characteristics and the two-dimensional multi-group (69-group) transport equation was solved by MOC method (with fully reflective boundary conditions in the axial direction) for each axial layer. Secondly, utilizing the results from 2D MOC calculations, based on the equivalent homogenization theory and the super-homogenization (SPH) technology, the heterogenous cells were homogenized, which produced the few-group homogenous cross sections as well as SPH factors. Finally, the 3D whole-core pin-by-pin SP3 calculation was carried out to obtain cell flux and power distribution. Moreover, the constructive solid geometry (CSG) was applied to enhance the complex geometric modeling capability in STORK. A combination of the enhanced neutron flow method and the equivalence theory was used to perform resonance calculations and a pre-produced table of resonance interference factors was adopted to handle the resonance interference effects. During 2D transport calculation, a two-level unstructured coarse mesh finite difference method was applied to accelerate the convergence of the MOC calculation. In the 3D pin-by-pin calculation, the 3D SP3 equations were solved by the transverse integration technique and the nodal expansion method with group transverse-integrated neutron fluxes approximated by the parabola expansion in the radial direction and by semi-analytical expansion in the axial direction. In terms of code development, a hybrid programming of CUDA, C++ and Python was adopted, and all the computational modules were developed based on CUDA/C++ with a large number of performance optimizations, so that 2D MOC calculations at each layer of the core could be carried out on multiple GPUs at the same time. To maximize computational efficiency, the computationally-intensive modules in STORK, including MOC calculation, CMFD, resonance calculation, burnup calculation, and SP3 calculation modules, were executed on the GPU. The validation of the SRORK code through the C5G7 3D Rodded problem and VERA benchmark problems demonstrates its high computational accuracy, with a radial assembly power error of less than 1%. However, due to the code's direct utilization of the energy spectrum of the adjacent layers' active regions for the axial reflector and the lack of consideration for neutron leakage from neighboring axial layers, significant discrepancies in axial power occur near the reflector and in fuel layers containing spacer grids, but they remain below 3%. More importantly, developed based on the CPU/GPU heterogeneous system, the code exhibits significant advantages in terms of computational efficiency and cost compared to similar neutron transport softwares.