  We present a fast direct solver for boundary integral equations on complex
surfaces in three dimensions using an extension of the recently introduced
recursive strong skeletonization scheme. For problems that are not highly
oscillatory, our algorithm computes an ${LU}$-like hierarchical factorization
of the dense system matrix, permitting application of the inverse in $\mathcal
O(n)$ time, where $n$ is the number of unknowns on the surface. The
factorization itself also scales linearly with the system size, albeit with a
somewhat larger constant. The scheme is built on a level-restricted adaptive
octree data structure, and therefore it is compatible with highly nonuniform
discretizations. Furthermore, the scheme is coupled with high-order accurate
locally-corrected Nystr\"om quadrature methods to integrate the singular and
weakly-singular Green's functions used in the integral representations. Our
method has immediate applications to a variety of problems in computational
physics. We concentrate here on studying its performance in acoustic scattering
(governed by the Helmholtz equation) at low to moderate frequencies, and
provide rigorous justification for compression of submatrices via proxy
surfaces.
