A few more updates for GPU compatibility for TensorKit #100
Conversation
There was a problem hiding this comment.
Is it only the Diagonal giving issues?
We could also just add a specialization for Diagonal in general, which we should be able to handle in a GPU friendly way:
function project_hermitian_native!(A::Diagonal, B::Diagonal, ::Val{anti}) where {anti}
if anti
diagview(A) .= imag.(diagview(B)) .* im
else
diagview(A) .= real.(diagview(B))
end
end[edit] I think we have a utility function somewhere for the imaginary part, I can't remember the name though
|
That's also fine, I wrote this in the post lunch carb coma there's probably a better way |
|
I guess the diagonal specialization is something we want anyways so that might make sense? |
Codecov Report❌ Patch coverage is
🚀 New features to boost your workflow:
|
src/implementations/svd.jl
Outdated
|
|
||
| return USVᴴ | ||
| end | ||
| svd_full!(A::Diagonal, USVᴴ, alg::GPU_SVDAlgorithm) = svd_full!(diagm(A.diag), USVᴴ, alg) |
There was a problem hiding this comment.
Is this actually being used? This is to cover the case where somehow a Diagonal ends up with a GPU_SVDAlgorithm?
There was a problem hiding this comment.
We currently don't have a DiagonalGPUAlgorithm, which is why this happens, I suppose
There was a problem hiding this comment.
I think this might also be a problem with the specific testing setup: we are manually picking the GPUAlgorithm and then calling this on a Diagonal, while the automated algorithm selection would have picked DiagonalAlgorithm instead.
I agree that we should make sure that the diagonal gpu arrays work, but I don't think we have to have specializations to make sure the GPU_SVDAlgorithm implementations work for Diagonal.
TLDR, instead of these specializations, I think I would rather just remove the tests, since this really is supposed to error?
| MatrixAlgebraKit.ishermitian_approx(A::StridedROCMatrix; kwargs...) = | ||
| @invoke MatrixAlgebraKit.ishermitian_approx(A::Any; kwargs...) |
There was a problem hiding this comment.
Is this kind of definition useful? What happens if it is removed?
There was a problem hiding this comment.
This is to avoid the StridedMatrix implementation that does scalar indexing, and instead use norm(project_hermitian(A)) <= ..., as this is the fallback definition for Any.
I updated this and simply copied over that code now, in hindsight the benefit of not repeating code is probably not outweighed by the clarity.
|
Yes, I think the max of the two tolerances is the best compromise… the
original I think was some late night coding again.
…On Wed, Nov 26, 2025 at 10:57 AM Jutho ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl
<#100 (comment)>
:
> MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; kwargs...) =
@invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...)
+function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedROCVector{T}}; atol, rtol, kwargs...) where {T <: Real}
+ return sum(abs2, A.diag) ≤ max(atol, rtol * norm(A))
Ok, maybe max(atol, rtol). Although that is not really equivalent. I
think it is natural that this cannot be satisfied if atol is zero (except
for norm(A) == 0, so it should definitely be norm(A) <= atol instead of a
strictly smaller than).
But the general condition is, for a matrix A that can always be split as
A_full = A_hermitian + A_antihermitian
Then our condition for isantihermitian_approx when atol == 0 amounts to
norm(A_hermitan) <= rtol * norm(A_full)
But in the case of a real diagonal matrix, there is no anti-hermitian part
and A_full == A_hermitian. So the condition above can never be satisfied
(for rtol < 1, whereas it is always satisfied for rtol >= 1, which isn't
really a sane choice).
—
Reply to this email directly, view it on GitHub
<#100 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AAGKJY4OG6JXW2F6HMPUEJT36V2SDAVCNFSM6AAAAACMLW45ISVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZTKMJQGEZDSMZUGE>
.
You are receiving this because you authored the thread.Message ID:
***@***.***>
|
Co-authored-by: Jutho <Jutho@users.noreply.github.com>
|
SVD now has tests, but I rather dislike the approach here and I think will implement |
|
Your PR no longer requires formatting changes. Thank you for your contribution! |
|
Apologies for the PR hijack here, I hope everything now turns green, and would be happy to merge.
|
|
I think the followup PR #106 can handle the SVD case, so I'm happy to turn on automerge if you are :) |
|
As a note for a potential (follow-up) improvement: I think the polar decompositions are now defaulting to |
Base.require_one_based_indexingreturns aBoolon some versions, notn(integer).Converting to the full
CuMatrix/ROCMatrixis inefficient/ugly but dispatching on aSubArrayof aDiagonalof aCuVector(orROCVector) is also really ugly. Not sure what the best approach is here, but this works for now (we can revisit if it's really problematic).