The problem can be often formulated as a weighted p-median problem. Real instances of the problem are characterized by big numbers of possible service center locations, which can take the value of several hundreds or thousands. The optimal solution can be obtained by the universal IP solvers only for smaller instances of the problem. The universal IP solvers are very time-consuming and often fail when solving a large instance. Our approach to the problem is based on the Erlenkotter procedure for solving of the uncapacitated facility location problem and on the Lagrangean relaxation of the constraint which limits number of the located center. The suggested approach finds the optimal solution in most of the studied instances. The quality and the feasibility of the resulting solutions of the suggested approach depend on the setting of the Lagrangean multiplier. A suitable value of the multiplier can be obtained by a bisection algorithm. The resulting multiplier cannot guarantee an optimal solution, but provides a near-to-optimal solution and a lower bound. If our approach does not obtain the optimal solution, then a heuristic improves the near-to-optimal solution. The resulting solution of our approach and the optimal solution obtained by the universal IP solver XPRESS-IVE are compared in the computational time and the quality of solutions.