The motivation of this work is to develop a computational treatment of eddy current simulation applicable to complex 3D geometries. To this end, an approach based on surface integral equations and involving the homogeneous full-space Green's function and dyad is developed. A Galerkin variant of the method of moments, using bilinear basis functions, is employed for discretizing the surface integral equations, leading to a linear system. Solving the linear system using a direct matrix solver is prohibitively expensive in computing time and memory resources. This difficulty is circumvented by using an iterative solver (typically GMRES) together with the Fast Multipole Method (FMM). The response of an eddy current probe is finally expressed with the Auld reciprocity theorem.