The paper is organized as follows. In the next section, we go briefly through the setting of Petrov-Galerkin discretization of boundary integral equations and introduce the notations necessary for the sequel. Then, in Section 3, we investigate the computation of the surface integrals arising when computing the matrix entries of the linear system. First, we have to provide some analytical properties of the kernel functions that arise, which then are used in the design of efficient cubature strategies. (Note that in more than one dimension the term quadrature is replaced by cubature). In Section 4, we explain the panel clustering algorithm. After having developed the basic principle, we formulate the algorithm and present an error analysis. Estimates of the asymptotic complexity of the algorithm follow and remarks on their relevance for practical problem sizes are given. Section 6 is devoted to the design of modern BEM software in C++. It is explained which data structures are well suited to making the different algorithms efficient and flexible with respect to various kinds of integral equations, discretization schemes such as, e.g., collocation or Galerkin BEMs, and various geometries and cubature techniques.