The boundary integral equation method (BIEM) is applied to a wide spectrum of numerical simulations. The method of hierarchical matrices (H-matrices) combined with the adaptive cross approximation (ACA) is an efficient approximation technique that can be applied to the dense matrices arising in the integral equation method. We have been developing the HACApK library for the H-matrices with ACA. To conduct large-sized simulations, we adopt new techniques including improvement and parallelization of H-matrices with ACA. These techniques are examined through practical numerical experiments for electric field analyses and earthquake cycle simulations.