Compute OLS with Fortran and Lapack

ols.f90

function ols(y, x, n, k) result (beta) 

    implicit none

    external DGELS

    integer, intent(in) :: n, k
    real(dp), allocatable, intent(in) :: y(:), x(:, :)
    integer :: info, lwork
    real(dp) :: beta(k)
    real(dp), allocatable :: work(:)
    allocate(work(100 * n * k))
    lwork = 100 * n * k

    call DGELS('No transpose', n, k, 1, x, n, y, n, work, lwork, info)
    beta = y(1:k)
    deallocate(work)

end function ols