Dense linear systems in large-scale kernel approximation in machine learning, discretization of boundary integral equations in mathematical physics, and low-rank approximation of Schur complements in large sparse matrix factorization often employ a multilevel structure of low-rank off-diagonal blocks. To solve such systems efficiently, we present a GPU-based parallel recursive skeletonization solver that utilizes batched dense linear algebra to achieve linear time complexity, proportionally to the matrix size, under fixed ranks, while allowing for tunable precision