  We introduce a fast direct solver for variable-coefficient elliptic partial
differential equations on surfaces based on the hierarchical Poincar\'e-Steklov
method. The method takes as input an unstructured, high-order quadrilateral
mesh of a surface and discretizes surface differential operators on each
element using a high-order spectral collocation scheme. Elemental solution
operators and Dirichlet-to-Neumann maps tangent to the surface are precomputed
and merged in a pairwise fashion to yield a hierarchy of solution operators
that may be applied in $\mathcal{O}(N \log N)$ operations for a mesh with $N$
elements. The resulting fast direct solver may be used to accelerate high-order
implicit time-stepping schemes, as the precomputed operators can be reused for
fast elliptic solves on surfaces. On a standard laptop, precomputation for a
12th-order surface mesh with over 1 million degrees of freedom takes 17
seconds, while subsequent solves take only 0.25 seconds. We apply the method to
a range of problems on both smooth surfaces and surfaces with sharp corners and
edges, including the static Laplace-Beltrami problem, the Hodge decomposition
of a tangential vector field, and some time-dependent nonlinear
reaction-diffusion systems.
