Algorithm 886: Padua2D-Lagrange Interpolation at Padua Points on Bivariate Domains

We present a stable and efficient Fortran implementation of polynomial interpolation at the Padua points on the square [−1, 1]^2 These points are unisolvent and their Lebesgue constant has . minimal order of growth (log square of the degree). The algorithm is based on the representation of the Lagrange interpolation formula in a suitable orthogonal basis, and takes advantage of a new matrix formulation together with the machine-specific optimized BLAS subroutine DGEMM for the matrix-matrix product. Extension to interpolation on rectangles, triangles and ellipses is also described.