  This paper presents a general method for applying hierarchical matrix
skeletonization factorizations to the numerical solution of boundary integral
equations with possibly rank-deficient integral operators. Rank-deficient
operators arise in boundary integral approaches to elliptic partial
differential equations with multiple boundary components, such as in the case
of multiple vesicles in a viscous fluid flow. Our generalized skeletonization
factorization retains the locality property afforded by the "proxy point
method," and allows for a parallelized implementation where different
processors work on different parts of the boundary simultaneously. Further,
when the boundary undergoes local geometric perturbations (such as movement of
an interior hole), the factorization can be recomputed efficiently with respect
to the number of modified discretization nodes. We present an application that
leverages a parallel implementation of skeletonization with updates in a shape
optimization regime.
