This research addresses the contact problem with friction often presenting a real challenge for efficient numerical simulations in engineering. Here, we focus on a rare methodology, flexibility method [1] based on concepts from the Boundary Element Method and thus requiring assembly and solving of dense matrices. This study explores the use of hierarchical or H-matrices to accelerate calculations within this class of methods. The paper redefines contact problem as an auxiliary problem thus enabling to consider contact tractions as external forces. Based on Hertz-Signorini-Moreau conditions for contact and eventual frictional contribution, this method integrates flexibility matrices in finite element formulation and solves this nonlinear problem in a separate loop employing integral operators inspired by boundary element methods. The main drawback of this method is that the construction of the dense flexibility matrix and its resolution could be too computationally extensive to ensure optimal performance. We suggest to overcome this drawback by an adapted use of hierarchical or H-matrices employing low-rank approximation for some blocks [2]. The method reorganizes the discretized contact surfaces hierarchically and simplifies interactions between distant elements using methods like Adaptive Cross Approximation (ACA) combined with Singular Value Decomposition (SVD). This hierarchical structure drastically reduces storage requirements and speeds up resolution by orders of magnitude. The algorithm and H-matrix construction technique will be analyzed through various examples, comparing performance with traditional methods. Potential extension to friction problems using the bi-potential method [3] is also considered.