In this paper we consider the fully discrete wavelet Galerkin scheme for the fast solution of boundary integral equations in three dimensions. It produces approximate solutions within discretization error accuracy offered by the underlying Galerkin method at a computational expense that stays proportional to the number of unknowns. We focus on implementational details of the scheme, in particular on numerical integration of relevant matrix coefficients. We illustrate the proposed algorithms by numerical results.