Given n pairwise distinct and arbitrarily spaced points P_i in a domain D of the x-y plane and n real numbers f_i, consider the problem of computing a bivariate function f(x,y) of class C1 in D whose values in P_i are exactly f_i, i=1,...,n, and whose first or second order partial derivatives satisfy appropriate equality and inequality constraints on a given set of p points Q_l in D.In this paper we present a method for solving the above problem, which is designed for extremely large data sets. A step of this method requires the solution of a large scale quadratic programming (QP) problem. The main purpose of this work is to analyse an iterative method for determining the solution of this QP problem: such a method is very efficient and well suited for parallel implementation on a multiprocessor system.