Skip to content

Use FastLapackInterface.jl? #300

@nathanaelbosch

Description

@nathanaelbosch

https://github.com/DynareJulia/FastLapackInterface.jl
We're already doing something similar in

"""
getupperright!(A)
Get the upper right part of a matrix `A` without allocating.
This function is mostly there to make accessing `R` from a `QR` object more efficient, as
`qr!(A).R` allocates since it does not use views. This function is not exported and is only
ever used internally by `triangularize!`(@ref).
"""
function getupperright!(A)
m, n = size(A)
return UpperTriangular(triu!(A[1:min(m, n), 1:n]))
end
"""
triangularize!(A; [cachemat])
Compute `qr(A).R` in the most efficient and allocation-free way possible.
The fallback implementation essentially computes `qr!(A).R`. But if `A` is of type
`StridedMatrix{<:LinearAlgebra.BlasFloat}`, we can make things more efficient by calling
LAPACK directly and using the preallocated cache `cachemat`.
"""
triangularize!(A; cachemat)
function triangularize!(A; cachemat=nothing)
QR = qr!(A)
return getupperright!(getfield(QR, :factors))
end
function triangularize!(A::StridedMatrix{<:LinearAlgebra.BlasFloat}; cachemat)
D = size(A, 2)
BLOCKSIZE = 36
R, _ = LinearAlgebra.LAPACK.geqrt!(A, @view cachemat[1:min(BLOCKSIZE, D), :])
return getupperright!(R)
end

but not quite. FastLapackInterface.jl might be the cleaner way to do this.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions