In this article, the multilevel adaptive cross approximation (MLACA) algorithm is presented to accelerate the boundary element method (BEM) for eddy current nondestructive evaluation (NDE) 3D problems involving arbitrary shapes. The Stratton-Chu formula, which does not have the low frequency breakdown issue, has been selected for modeling. The equivalent electric and magnetic surface currents are expanded with Rao-Wilton-Glisson (RWG) vector basis functions while the normal component of the magnetic field is expanded with pulse basis functions. The MLACA compresses the rank deficient matrices with the ACA and the butterfly algorithm. We improve the efficiency of MLACA by truncating the integral kernels after a certain distance and applying the multi-stage (level) algorithm adaptively based on the criteria for different operators to further decrease the memory and CPU time requirements while keeping almost the same accuracy comparing with the traditional MLACA. The proposed method is especially helpful to deal with the large solution domain issue of the BEM for eddy current problems. Numerical predictions are compared with the analytical, the semi-analytical predictions and the experimental results for 3D eddy current NDE problems of practical interest to demonstrate the robustness and efficiency of the proposed method.