  A fast method for solving boundary integral equations with the generalized
Neumann kernel and the adjoint generalized Neumann kernel is presented. The
method is based on discretizing the integral equations by the Nystr\"om method
with the trapezoidal rule to obtain $(m+1)n\times(m+1)n$ linear systems where
$m+1$ is the multiplicity of the multiply connected domain and $n$ is the
number of nodes in the discretization of each boundary component. The obtained
linear systems are solved by the generalized minimal residual (GMRES) method.
Each iteration of the GMRES method requires a matrix-vector product which can
be computed using the Fast Multipole Method (FMM). The complexity of the
presented method is $O((m+1)n\ln n)$ for the integral equation with the
generalized Neumann kernel and $O((m+1)n)$ for the integral equation with the
adjoint generalized Neumann kernel. The presented numerical results illustrate
that the presented method gives accurate results even for domains with high
connectivity, domains with piecewise smooth boundaries, and domains with close
boundaries.
