We provide a numerical package for the computation of a d-variate near G-optimal polynomial regression design of degree m on a finite design space X ⊂ R d , by few iterations of a basic multiplicative algorithm followed by Tchakaloff-like compression of the discrete measure keeping the reached G-efficiency, via an accelerated version of the Lawson-Hanson algorithm for Non-Negative Least Squares (NNLS) problems. This package can solve on a personal computer large-scale problems where c a r d ( X ) × dim ( P 2 m d ) is up to 10 8 – 10 9 , being dim ( P 2 m d ) = 2 m + d d = 2 m + d 2 m . Several numerical tests are presented on complex shapes in d = 3 and on hypercubes in d > 3 .