We present a novel numerical scheme for solving the Stokes equation with variable coefficients in the unit box. Our scheme is based on a volume integral equation formulation. We employ a novel adaptive fast multipole method for volume integrals to obtain a scheme that is algorithmically optimal. To increase performance, we have integrated our code with both NVIDIA and Intel accelerators. Overall, our scheme supports non-uniform discretizations and is spectrally accurate. In our largest run we solved a problem with 20 billion unknowns, using a 14-order approximation for the velocity, on 2048 nodes of the Stampede platform at the Texas Advanced Computing Center. We achieved 0.656 PetaFLOPS for the overall code (20% efficiency) and 1 PetaFLOPS for the fast multipole part (40% efficiency). As an application example, we simulate Stokes flow in a porous medium with highly complex pore structure using a penalty formulation to enforce the no slip condition.