  Geometric matching is an important topic in computational geometry and has
been extensively studied over decades. In this paper, we study a
geometric-matching problem, known as geometric many-to-many matching. In this
problem, the input is a set $S$ of $n$ colored points in $\mathbb{R}^d$, which
implicitly defines a graph $G = (S,E(S))$ where $E(S) = \{(p,q): p,q \in S
\text{ have different colors}\}$, and the goal is to compute a minimum-cost
subset $E^* \subseteq E(S)$ of edges that cover all points in $S$. Here the
cost of $E^*$ is the sum of the costs of all edges in $E^*$, where the cost of
a single edge $e$ is the Euclidean distance (or more generally, the
$L_p$-distance) between the two endpoints of $e$. Our main result is a
$(1+\varepsilon)$-approximation algorithm with an optimal running time
$O_\varepsilon(n \log n)$ for geometric many-to-many matching in any fixed
dimension, which works under any $L_p$-norm. This is the first near-linear
approximation scheme for the problem in any $d \geq 2$. Prior to this work,
only the bipartite case of geometric many-to-many matching was considered in
$\mathbb{R}^1$ and $\mathbb{R}^2$, and the best known approximation scheme in
$\mathbb{R}^2$ takes $O_\varepsilon(n^{1.5} \cdot \mathsf{poly}(\log n))$ time.
