We consider the problem of minimizing the distance from a given n-dimensional vector to a set defined by constraints of the form x_i ≤ x_j. Such constraints induce a partial order of the components x_i, which can be illustrated by an acyclic directed graph. This problem is also known as the isotonic regression (IR) problem. IR has important applications in statistics, operations research and signal processing, with most of them characterized by a very large value of n. For such large-scale problems, it is of great practical importance to develop algorithms whose complexity does not rise with n too rapidly. The existing optimization-based algorithms and statistical IR algorithms have either too high computational complexity or too low accuracy of the approximation to the optimal solution they generate. We introduce a new IR algorithm, which can be viewed as a generalization of the Pool-Adjacent-Violator (PAV) algorithm from completely to partially ordered data. Our algorithm combines both low computational complexity O(n²) and high accuracy. This allows us to obtain sufficiently accurate solutions to IR problems with thousands of observations.