Sample least squares problem using LAPACK routines from Intel MKL

mkl_least_squares_example.f90

! # Check out the link below for specific compile instructions for your system:
! # https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor/
!
! To compile use:
! LINK="${MKLROOT}/lib/intel64/libmkl_lapack95_lp64.a -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_lp64.a ${MKLROOT}/lib/intel64/libmkl_sequential.a ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl"
! INC="-I${MKLROOT}/include/intel64/lp64/ -I${MKLROOT}/include/"
! ifort -warn all -o example3 $INC example3.f90 $LINK 

module m_least_squares

    use f95_precision, only: wp => dp
    use lapack95, only: gelss, gels
    implicit none
    private

    public :: wp, quadratic_fit

contains

    function quadratic_fit(n,x,y,rcond,info) result(alpha)
        integer, intent(in) :: n        ! The number of points.
        real(wp), intent(in) :: x(n)    ! x-values of the points.
        real(wp), intent(in) :: y(n)    ! y-values of the points.
        real(wp), intent(in), optional :: rcond     ! Used to determine the effective rank of A. 
                                                    ! Singular values s(i) ≤ rcond*s(1) are treated as zero. 
        integer, intent(out), optional :: info      ! Zero, if routine is succesful.
        real(wp) :: alpha(3)
        
        integer :: info_, rank_
        real(wp) :: A(n,3), b(n), rcond_
        real(wp), target :: s_(3)

        rcond_ = 100*epsilon(1.0_wp)
        if (present(rcond)) rcond_ = rcond

        A(:,1) = 1.0_wp
        A(:,2) = x
        A(:,3) = x**2
        b = y

        ! SVD: https://software.intel.com/en-us/mkl-developer-reference-fortran-gelss
        call gelss(A=A,b=b,rank=rank_,rcond=rcond_,info=info_)

        ! QR: https://software.intel.com/en-us/mkl-developer-reference-fortran-gels
        ! call gels(A,b,trans="N",info=info_)

        if (present(info)) info = info_
        alpha = b(1:3)

    end function

end module

program test_least_squares

    use m_least_squares
    implicit none

    integer, parameter :: n = 4
    integer :: i
    real(8) :: x(n), y(n), a(3)
    integer :: ierr

    ! Fitting a parabola to a sine function

    x = [1.0_wp,2.0_wp,3.0_wp,4.0_wp]
    y = sin(x)

    a = quadratic_fit(n,x,y,info=ierr)

    print *, "# coefficients = ", a
    print *, "# ierr = ", ierr
    do i = 1, n
        print *, x(i), y(i), a(1) + a(2)*x(i) + a(3)*x(i)**2
    end do

end program