zgelss.jl

# Hacked toether from https://github.com/JuliaLang/julia/blob/2d5741174ce3e6a394010d2e470e4269ca54607f/stdlib/LinearAlgebra/src/lapack.jl#L1371
# http://www.netlib.org/lapack/explore-html/d6/d10/group__complex16_g_esolve_ga61e68db68886c3f80753fac87ca35a6e.html#ga61e68db68886c3f80753fac87ca35a6e

using LinearAlgebra

const liblapack = Base.liblapack_name

import LinearAlgebra.LAPACK.chklapackerror, LinearAlgebra.LAPACK.subsetrows

import ..LinearAlgebra.BLAS.@blasfunc

import ..LinearAlgebra: BlasFloat, BlasInt, LAPACKException,
    DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare

using Base: iszero, require_one_based_indexing


gelss = :zgelss_
elty = :ComplexF64
relty = :Float64

    @eval begin
        # SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
        #      $                   WORK, LWORK, RWORK, IWORK, INFO )
        # *     .. Scalar Arguments ..
        #       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
        #       DOUBLE PRECISION   RCOND
        # *     ..
        # *     .. Array Arguments ..
        #       DOUBLE PRECISION   RWORK( * ), S( * )
        #       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
        function gelss!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=-one($relty))
            require_one_based_indexing(A, B)
            chkstride1(A, B)
            m, n  = size(A)
            if size(B, 1) != m
                throw(DimensionMismatch("B has leading dimension $(size(B,1)) but needs $m"))
            end
            newB = [B; zeros($elty, max(0, n - size(B, 1)), size(B, 2))]
            s     = similar(A, $relty, min(m, n))
            rnk   = Ref{BlasInt}()
            info  = Ref{BlasInt}()
            work  = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            rwork = Vector{$relty}(undef, 5*min(m,n))
            #iwork = Vector{BlasInt}(undef, 1)
            for i = 1:2  # first call returns lwork as work[1], rwork length as rwork[1] and iwork length as iwork[1]
                ccall((@blasfunc($gelss), liblapack), Cvoid,
                      (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
                       Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$relty},
                       Ref{$relty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                       Ptr{$relty}, Ref{BlasInt}),
                      m, n, size(B,2), A,
                      max(1,stride(A,2)), newB, max(1,stride(B,2),n), s,
                      $relty(rcond), rnk, work, lwork,
                      rwork, info)
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(real(work[1]))
                    resize!(work, lwork)
                    #print(BlasInt(rwork[1]))
                    #resize!(rwork, BlasInt(rwork[1]))
                    #resize!(iwork, iwork[1])
                end
            end
            subsetrows(B, newB, n), rnk[]
        end
    end