I had the same issue when computing "lowess" regression. I ended up using a Fortran subroutine that called the LAPACK subroutine DGELS. (It is possible to vectorize for multiple cores using OpenMP or link with a multithreaded LAPACK/BLAS. What works best is dependent on the problem.) Sturla