Skip to content

Instantly share code, notes, and snippets.

@mkitti
Last active Jan 24, 2020
Embed
What would you like to do?
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment