Skip to content

Conversation

@kshyatt
Copy link
Member

@kshyatt kshyatt commented Nov 17, 2025

Base.require_one_based_indexing returns a Bool on some versions, not n (integer).

Converting to the full CuMatrix/ROCMatrix is inefficient/ugly but dispatching on a SubArray of a Diagonal of a CuVector (or ROCVector) 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).

@kshyatt kshyatt requested a review from lkdvos November 17, 2025 18:31
Copy link
Member

@lkdvos lkdvos left a comment

Choose a reason for hiding this comment

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

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

@kshyatt
Copy link
Member Author

kshyatt commented Nov 17, 2025

That's also fine, I wrote this in the post lunch carb coma there's probably a better way

@lkdvos
Copy link
Member

lkdvos commented Nov 17, 2025

I guess the diagonal specialization is something we want anyways so that might make sense?

@codecov
Copy link

codecov bot commented Nov 17, 2025

Codecov Report

❌ Patch coverage is 66.66667% with 2 lines in your changes missing coverage. Please review.

Files with missing lines Patch % Lines
...ixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl 0.00% 1 Missing ⚠️
...MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl 0.00% 1 Missing ⚠️
Files with missing lines Coverage Δ
src/implementations/projections.jl 96.20% <100.00%> (+0.09%) ⬆️
...ixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl 80.88% <0.00%> (-1.21%) ⬇️
...MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl 56.66% <0.00%> (-0.97%) ⬇️
🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

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))
Copy link
Member

Choose a reason for hiding this comment

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

This is a strange condition, as the left hand side is norm(A)^2, so this is saying norm(A)^2 < max(atol, rtol * norm(A))

This amounts to norm(A)^2 < atol or norm(A) < rtol, both of which look off.

Because a real diagonal matrix equals its hermitian project and has zero non-hermitian projection, I think the condition simply needs to be norm(A) < atol to be compatible with how approximate antithermiticity is checked elsewhere.

Copy link
Member Author

Choose a reason for hiding this comment

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

Fixed!

Copy link
Member Author

Choose a reason for hiding this comment

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

Actually, doesn't this break things in the case that atol is 0, but rtol is not? (As CI shows)

Copy link
Member

Choose a reason for hiding this comment

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

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).


return USVᴴ
end
svd_full!(A::Diagonal, USVᴴ, alg::GPU_SVDAlgorithm) = svd_full!(diagm(A.diag), USVᴴ, alg)
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 actually being used? This is to cover the case where somehow a Diagonal ends up with a GPU_SVDAlgorithm?

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes

Copy link
Member Author

Choose a reason for hiding this comment

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

We currently don't have a DiagonalGPUAlgorithm, which is why this happens, I suppose

Comment on lines 135 to 136
MatrixAlgebraKit.ishermitian_approx(A::StridedROCMatrix; kwargs...) =
@invoke MatrixAlgebraKit.ishermitian_approx(A::Any; kwargs...)
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 kind of definition useful? What happens if it is removed?

@kshyatt
Copy link
Member Author

kshyatt commented Nov 26, 2025 via email

@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))
return norm(A) max(atol, rtol)
Copy link
Member

Choose a reason for hiding this comment

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

I do think just having atol in the right hand side is what we want to be compatible with what is happening for non-GPU Diagonal arrays, or for the equivalent representation of the same matrix as dense CuMatrix or ROCMatrix. In all of these cases, the general implementation will essentially amount to norm(A) <= atol when A only has nonzero entries on the diagonal, which are furthermore purely real.

Copy link
Member

Choose a reason for hiding this comment

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

So if that causes test failures, I would argue the tests are incorrect.

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.

4 participants