In this paper we study a multi-server queueing system with retrials and an infinite orbit. The arrival of primary customers is described by a batch Markovian arrival process (BMAP), the service times have a phase type (PH) distribution. Previously, in the literature, such a system was mainly considered under the strict assumption that the intervals between the repeated attempts from the orbit have an exponential distribution. Only a few publications deal with retrial queueing systems with non-exponential inter-retrial times. These publications assume either the rate of retrials is constant regardless of the number of customers in the orbit or this rate is constant when the number of orbital customers exceeds a certain threshold. Such assumptions essentially simplify the mathematical analysis of the system but do not reflect the nature of the majority of the real-life retrial processes. The main feature of the model under study is that we consider the classical retrial strategy under which the retrial rate is proportional to the number of orbital customers. However, in this case, the assumption of the non-exponential distribution of inter-retrial times leads to insurmountable computational difficulties. To overcome these difficulties, we suppose that inter-retrial times have a phase-type distribution if the number of customers in the orbit is less or equal to some non-negative integer (threshold) and have an exponential distribution in the contrary case. By appropriately choosing the threshold, one can obtain a sufficiently accurate approximation of the system with PH distribution of inter-retrial times. Thus, the model under study takes into account the realistic nature of retrial process and, at the same time, does not resort to restrictions such as constant retrial rate or to rough truncation methods often applied for analysis of retrial queueing systems with an infinite orbit. We describe the behavior of the system by a multi-dimensional Markov chain, derive the stability condition, calculate the steady state distribution and the main performance indicators of the system. We make sure numerically that there is a reasonable value of a threshold under which our model can be served as a good approximation of the BMAP/PH/N queueing system with PH distribution of inter-retrial times. We also numerically compare the system under consideration with the corresponding queueing system having exponentially distributed inter-retrial times and see that the latter is a poor approximation of the system with $PH$ distribution of inter-retrial times. We present a number of illustrative numerical examples to analyze the behavior of the system performance indicators depending on the system parameters, the variance of inter-retrial times and the correlation in the input flow.