Skip to content

Conversation

kshyatt
Copy link
Member

@kshyatt kshyatt commented Aug 26, 2025

The ReshapedArray overrides are needed to dispatch to the correct GPU algorithms. Needed to modify the type signature for the default algorithms to avoid ambiguities. Also it's nice to give some more info about dimension mismatches.

@kshyatt kshyatt requested a review from lkdvos August 26, 2025 08:29
Copy link

codecov bot commented Aug 26, 2025

Codecov Report

❌ Patch coverage is 31.03448% with 20 lines in your changes missing coverage. Please review.

Files with missing lines Patch % Lines
...MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl 11.76% 15 Missing ⚠️
src/implementations/svd.jl 44.44% 5 Missing ⚠️
Files with missing lines Coverage Δ
ext/MatrixAlgebraKitCUDAExt/yacusolver.jl 95.75% <100.00%> (ø)
src/implementations/eig.jl 98.57% <ø> (-0.03%) ⬇️
src/implementations/lq.jl 98.62% <100.00%> (ø)
src/implementations/svd.jl 91.77% <44.44%> (-2.01%) ⬇️
...MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl 48.48% <11.76%> (-24.25%) ⬇️
🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.


include("yacusolver.jl")

function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T<:StridedCuMatrix}
function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {TT<:BlasFloat, T<:StridedCuMatrix{TT}}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are probably a couple more of these somewhat complex wrapper types that can still be handled by these algorithms, how do you feel about doing something like

for MatType in [...]
@eval ...
end

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds fine to me, do we have a list of the ones we want?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Was this particular change also induced by TensorKit requirements, or simply more strictness (which I fully support)?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be equivalent to defining a new type constant

const StridedCuBLASMatrix{T} = StridedCuMatrix{T} where {T<:BlasFloat}

and then using default_xxx_algorithm(::Type{<:StridedCuBLASMatrix}; kwargs...) everywhere?

@@ -288,7 +288,7 @@ function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAl
U, S, Vᴴ = USVᴴ
if alg isa GPU_QRIteration
isempty(alg.kwargs) ||
throw(ArgumentError("GPU_QRIteration does not accept any keyword arguments"))
@warn "GPU_QRIteration does not accept any keyword arguments"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Something I've been wondering is if these kinds of warnings and checks should just be moved to the actual algorithm constructors, which might avoid having to repeat them

Copy link
Member

@Jutho Jutho Sep 3, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would need to check, but some of the LAPACK algorithm names are used for different factorisations, and might support keywords for one and not for the other, or might support a different set of keywords. We can also remove these explicit checks, just pass on alg.kwargs to the underlying LAPACK/CUSOLVER/rocSOLVER routine always, and then have these complain if the set of keyword arguments is invalid.

@@ -181,7 +181,7 @@ function lq_via_qr!(A::AbstractMatrix, L::AbstractMatrix, Q::AbstractMatrix,
qr_alg::AbstractAlgorithm)
m, n = size(A)
minmn = min(m, n)
At = adjoint!(similar(A'), A)::AbstractMatrix
At = min(m, n) > 0 ? adjoint!(similar(A'), A)::AbstractMatrix : similar(A')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this ::AbstractMatrix type assert useful?

@kshyatt
Copy link
Member Author

kshyatt commented Sep 3, 2025

Changed the format of some of the orthnull tests -- CuArray and UniformScaling and isapprox really don't play nicely together. Also added a path to use CUSOLVER_QRIteration even in the case where m < n via transpose. Happy for feedback on the implementation - should this be a separate algorithm?

Ut = similar(U')
Vᴴt = similar(Vᴴ')
if size(U) == (m, m)
_gpu_gesvd!(At, view(S, 1:minmn, 1), Vᴴt, Ut)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Didn't think long about it, but it was not immediately clear to me why this was necessary. Isn't S always of length minmn?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here I'm following what the CPU bindings have done, but I suppose we could be reusing an S over and over between differently sized arrays?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants