From 2731e4bbaca816eae92b4093f326eb197540192f Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Tue, 13 Oct 2026 05:07:14 -0400 Subject: [PATCH 1/5] Make gaugefixing for svd_full and eigh GPU friendly --- src/common/gauge.jl | 36 +++++++++++++++--------------------- 1 file changed, 15 insertions(+), 21 deletions(-) diff --git a/src/common/gauge.jl b/src/common/gauge.jl index 4c756b087..580d78742 100644 --- a/src/common/gauge.jl +++ b/src/common/gauge.jl @@ -36,32 +36,26 @@ function gaugefix!(::typeof(lq_householder!), L, Q, Ld) end function gaugefix!(::Union{typeof(eig_full!), typeof(eigh_full!), typeof(gen_eig_full!)}, V::AbstractMatrix) - for j in axes(V, 2) - v = view(V, :, j) - s = sign(_argmaxabs(v)) - @inbounds v .*= conj(s) - end + signs = reduce(_largest, V; dims = 1, init = zero(eltype(V))) + @. signs = sign(signs) + @. V = V * conj(signs) return V end function gaugefix!(::typeof(svd_full!), U, Vᴴ) m, n = size(U, 2), size(Vᴴ, 1) - for j in 1:max(m, n) - if j <= min(m, n) - u = view(U, :, j) - v = view(Vᴴ, j, :) - s = sign(_argmaxabs(u)) - u .*= conj(s) - v .*= s - elseif j <= m - u = view(U, :, j) - s = sign(_argmaxabs(u)) - u .*= conj(s) - else - v = view(Vᴴ, j, :) - s = sign(_argmaxabs(v)) - v .*= conj(s) - end + k = min(m, n) + signs = reduce(_largest, U; dims = 1, init = zero(eltype(U))) + @. signs = sign(signs) + @. U = U * conj(signs) + Vᴴk = view(Vᴴ, axes(Vᴴ, 1)[1:k], :) + signs_t = transpose(view(signs, :, axes(signs, 2)[1:k])) + @. Vᴴk = signs_t * Vᴴk + if n > m + Vᴴr = view(Vᴴ, axes(Vᴴ, 1)[(k + 1):n], :) + signs_r = reduce(_largest, Vᴴr; dims = 2, init = zero(eltype(Vᴴ))) + @. signs_r = sign(signs_r) + @. Vᴴr = Vᴴr * conj(signs_r) end return (U, Vᴴ) end From 54471069f3609dc8be36a3afeb9664b0e5423247 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Tue, 13 Oct 2026 11:49:11 +0200 Subject: [PATCH 2/5] Update src/common/gauge.jl Co-authored-by: Jutho --- src/common/gauge.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/common/gauge.jl b/src/common/gauge.jl index 580d78742..505ef65ce 100644 --- a/src/common/gauge.jl +++ b/src/common/gauge.jl @@ -55,7 +55,7 @@ function gaugefix!(::typeof(svd_full!), U, Vᴴ) Vᴴr = view(Vᴴ, axes(Vᴴ, 1)[(k + 1):n], :) signs_r = reduce(_largest, Vᴴr; dims = 2, init = zero(eltype(Vᴴ))) @. signs_r = sign(signs_r) - @. Vᴴr = Vᴴr * conj(signs_r) + @. Vᴴr = conj(signs_r) * Vᴴr end return (U, Vᴴ) end From d855454786d3ad5d21763ecfe0cfc7af621c1438 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Tue, 13 Oct 2026 14:18:07 +0200 Subject: [PATCH 3/5] Update gauge.jl Add a comment reminding us to investigate further sign related optimizations --- src/common/gauge.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/common/gauge.jl b/src/common/gauge.jl index 505ef65ce..2889c4e52 100644 --- a/src/common/gauge.jl +++ b/src/common/gauge.jl @@ -42,6 +42,9 @@ function gaugefix!(::Union{typeof(eig_full!), typeof(eigh_full!), typeof(gen_eig return V end +# TODO: can do the multiplication @. Vᴴk = signs_t * Vᴴk in one fell swoop +# earlier, or compute signs in the n > m case entirely with reference to V? +# This might require modifying forward-mode AD rules. function gaugefix!(::typeof(svd_full!), U, Vᴴ) m, n = size(U, 2), size(Vᴴ, 1) k = min(m, n) From 19fe0c86994a6fdc1e26009254cf71b3418f0bc1 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Tue, 13 Oct 2026 09:35:38 -0400 Subject: [PATCH 4/5] Seed rng for V in eig test --- test/testsuite/decompositions/eig.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/testsuite/decompositions/eig.jl b/test/testsuite/decompositions/eig.jl index 6ddfef5d6..64af4eba7 100644 --- a/test/testsuite/decompositions/eig.jl +++ b/test/testsuite/decompositions/eig.jl @@ -117,7 +117,7 @@ function test_eig_trunc( m4 = 4 atol = sqrt(eps(real(eltype(T)))) - V = T <: Diagonal ? I : randn!(similar(A, eltype(T), m4, m4)) + V = T <: Diagonal ? I : copyto!(similar(A, eltype(T), m4, m4), randn(rng, eltype(T), m4, m4)) Ddiag = similar(A, real(eltype(T)), m4) copyto!(Ddiag, real(eltype(T))[0.9, 0.3, 0.1, 0.01]) D = Diagonal(Ddiag) @@ -158,7 +158,7 @@ function test_eig_trunc_algs( m4 = 4 atol = sqrt(eps(real(eltype(T)))) - V = T <: Diagonal ? I : randn!(similar(A, eltype(T), m4, m4)) + V = T <: Diagonal ? I : copyto!(similar(A, eltype(T), m4, m4), randn(rng, eltype(T), m4, m4)) Ddiag = similar(A, real(eltype(T)), m4) copyto!(Ddiag, real(eltype(T))[0.9, 0.3, 0.1, 0.01]) D = Diagonal(Ddiag) From eae633c9ce56e363134cba62d8230dae175031cd Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Tue, 13 Oct 2026 10:11:46 -0400 Subject: [PATCH 5/5] Fixup eig(h) --- test/testsuite/decompositions/eig.jl | 4 ++-- test/testsuite/decompositions/eigh.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/testsuite/decompositions/eig.jl b/test/testsuite/decompositions/eig.jl index 64af4eba7..78d5550e4 100644 --- a/test/testsuite/decompositions/eig.jl +++ b/test/testsuite/decompositions/eig.jl @@ -117,7 +117,7 @@ function test_eig_trunc( m4 = 4 atol = sqrt(eps(real(eltype(T)))) - V = T <: Diagonal ? I : copyto!(similar(A, eltype(T), m4, m4), randn(rng, eltype(T), m4, m4)) + V = instantiate_matrix(T, (m4, m4)) Ddiag = similar(A, real(eltype(T)), m4) copyto!(Ddiag, real(eltype(T))[0.9, 0.3, 0.1, 0.01]) D = Diagonal(Ddiag) @@ -158,7 +158,7 @@ function test_eig_trunc_algs( m4 = 4 atol = sqrt(eps(real(eltype(T)))) - V = T <: Diagonal ? I : copyto!(similar(A, eltype(T), m4, m4), randn(rng, eltype(T), m4, m4)) + V = instantiate_matrix(T, (m4, m4)) Ddiag = similar(A, real(eltype(T)), m4) copyto!(Ddiag, real(eltype(T))[0.9, 0.3, 0.1, 0.01]) D = Diagonal(Ddiag) diff --git a/test/testsuite/decompositions/eigh.jl b/test/testsuite/decompositions/eigh.jl index 087ee7341..e1654b974 100644 --- a/test/testsuite/decompositions/eigh.jl +++ b/test/testsuite/decompositions/eigh.jl @@ -122,7 +122,7 @@ function test_eigh_trunc( @testset "specify truncation algorithm" begin atol = sqrt(eps(real(eltype(T)))) m4 = 4 - smallA = randn!(similar(A, (m4, m4))) + smallA = instantiate_matrix(T, (m4, m4)) V = T <: Diagonal ? I : qr_compact(smallA)[1] Ddiag = similar(A, real(eltype(T)), m4) copyto!(Ddiag, real(eltype(T))[0.9, 0.3, 0.1, 0.01]) @@ -161,7 +161,7 @@ function test_eigh_trunc_algs( # truncrank atol = sqrt(eps(real(eltype(T)))) m4 = 4 - smallA = randn!(similar(A, (m4, m4))) + smallA = instantiate_matrix(T, (m4, m4)) V = T <: Diagonal ? I : qr_compact(smallA)[1] Ddiag = similar(A, real(eltype(T)), m4) copyto!(Ddiag, real(eltype(T))[0.9, 0.3, 0.1, 0.01])