linear_solver Subroutine

    subroutine linear_solver(m,n,a,b,x,info)

    implicit none

    integer,intent(in)                 :: n          !! number of columns in `a`
    integer,intent(in)                 :: m          !! number of rows in `a`
    real(wp),dimension(m,n),intent(in) :: a          !! `A` matrix of the linear system
    real(wp),dimension(m),intent(in)   :: b          !! RHS of the linear system
    real(wp),dimension(n),intent(out)  :: x          !! the solution of the linear system.
    integer,intent(out)                :: info       !! output status flag (`=0` if success)

    ! LAPACK routine interfaces:
    interface
        subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
            !! See: [?gesv](https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/lapack-routines/lapack-linear-equation-routines/lapack-linear-equation-driver-routines/gesv.html)
            import :: wp
            implicit none
            integer :: info
            integer :: lda
            integer :: ldb
            integer :: n
            integer :: nrhs
            integer :: ipiv(*)
            real(wp) :: a(lda,*)
            real(wp) :: b(ldb,*)
        end subroutine dgesv
        subroutine dgels(trans,m,n,nrhs,a,lda,b,ldb,work,lwork,info)
            !! See: [?gels](https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/lapack-routines/lapack-least-squares-and-eigenvalue-problem/lapack-least-squares-eigenvalue-problem-driver/linear-least-squares-lls-problems-lapack-driver/gels.html)
            import :: wp
            implicit none
            character :: trans
            integer :: info
            integer :: lda
            integer :: ldb
            integer :: lwork
            integer :: m
            integer :: n
            integer :: nrhs
            real(wp) :: a(lda,*)
            real(wp) :: b(ldb,*)
            real(wp) :: work(*)
        end subroutine dgels
    end interface

    integer,dimension(:),allocatable    :: ipiv   !! pivot indices array
    real(wp),dimension(:,:),allocatable :: bmat   !! copy of `b` so it won't be overwritten
    real(wp),dimension(:),allocatable   :: work   !! work array for `dgels`
    real(wp),dimension(:,:),allocatable :: amat   !! copy of `a` so it won't be overwritten
    integer                             :: lwork  !! size of `work`

    allocate(amat(m,n))
    allocate(bmat(max(1,m,n),1))

    if (n==m) then  !normal inverse

        allocate(ipiv(n))

        amat = a
        bmat(1:n,1) = b
        call dgesv(n, 1, amat, n, ipiv, bmat, n, info)
        x = bmat(1:n,1)

    else

        amat = a
        bmat = zero
        bmat(1:m,1) = b
        lwork = min(m,n)+max(1,m,n)
        allocate(work(lwork))
        call dgels('N', m, n, 1, amat, m, bmat, max(1,m,n), work, lwork, info)
        x = bmat(1:n,1)

    end if

    end subroutine linear_solver