From 93f37f5793e256cbe6efc3b233aff5d70e8d742b Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Wed, 29 Oct 2025 14:02:05 +0100 Subject: [PATCH 01/25] add GenericExt Add extension using GenericLinearAlgebra and GenericSchur to be able to deal with `BigFloat`s --- Project.toml | 7 +- ext/MatrixAlgebraKitGenericExt.jl | 230 ++++++++++++++++++++++++++++++ src/MatrixAlgebraKit.jl | 1 + src/interface/decompositions.jl | 11 ++ test/bigfloats/eig.jl | 117 +++++++++++++++ test/bigfloats/eigh.jl | 119 ++++++++++++++++ test/bigfloats/lq.jl | 150 +++++++++++++++++++ test/bigfloats/qr.jl | 135 ++++++++++++++++++ test/bigfloats/svd.jl | 210 +++++++++++++++++++++++++++ test/runtests.jl | 14 ++ 10 files changed, 993 insertions(+), 1 deletion(-) create mode 100644 ext/MatrixAlgebraKitGenericExt.jl create mode 100644 test/bigfloats/eig.jl create mode 100644 test/bigfloats/eigh.jl create mode 100644 test/bigfloats/lq.jl create mode 100644 test/bigfloats/qr.jl create mode 100644 test/bigfloats/svd.jl diff --git a/Project.toml b/Project.toml index be0bee32..aa1f3cf1 100644 --- a/Project.toml +++ b/Project.toml @@ -10,11 +10,14 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" +GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e" [extensions] MatrixAlgebraKitChainRulesCoreExt = "ChainRulesCore" MatrixAlgebraKitAMDGPUExt = "AMDGPU" MatrixAlgebraKitCUDAExt = "CUDA" +MatrixAlgebraKitGenericExt = ["GenericLinearAlgebra", "GenericSchur"] [compat] AMDGPU = "2" @@ -22,6 +25,8 @@ Aqua = "0.6, 0.7, 0.8" ChainRulesCore = "1" ChainRulesTestUtils = "1" CUDA = "5" +GenericLinearAlgebra = "0.3.19" +GenericSchur = "0.5.6" JET = "0.9, 0.10" LinearAlgebra = "1" SafeTestsets = "0.1" @@ -42,4 +47,4 @@ TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Aqua", "JET", "SafeTestsets", "Test", "TestExtras","ChainRulesCore", "ChainRulesTestUtils", "StableRNGs", "Zygote", "CUDA", "AMDGPU"] +test = ["Aqua", "JET", "SafeTestsets", "Test", "TestExtras", "ChainRulesCore", "ChainRulesTestUtils", "StableRNGs", "Zygote", "CUDA", "AMDGPU", "GenericLinearAlgebra", "GenericSchur"] diff --git a/ext/MatrixAlgebraKitGenericExt.jl b/ext/MatrixAlgebraKitGenericExt.jl new file mode 100644 index 00000000..b95741ea --- /dev/null +++ b/ext/MatrixAlgebraKitGenericExt.jl @@ -0,0 +1,230 @@ +module MatrixAlgebraKitGenericExt + +using MatrixAlgebraKit +using MatrixAlgebraKit: LAPACK_SVDAlgorithm, LAPACK_EigAlgorithm, LAPACK_EighAlgorithm, LAPACK_QRIteration +using MatrixAlgebraKit: uppertriangular! +using MatrixAlgebraKit: @algdef, Algorithm, check_input +using MatrixAlgebraKit: sign_safe +using GenericLinearAlgebra +using GenericSchur +using LinearAlgebra: I, Diagonal + +function MatrixAlgebraKit.default_svd_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} + return BigFloat_svd_QRIteration() +end + +function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix{T}, USVᴴ, alg::BigFloat_svd_QRIteration) where {T <: Union{BigFloat, Complex{BigFloat}}} + check_input(svd_compact!, A, USVᴴ, alg) + U, S, V = GenericLinearAlgebra.svd(A) + return U, Diagonal(S), V' # conjugation to account for difference in convention +end + +function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::BigFloat_svd_QRIteration)::Tuple{Matrix{T}, Matrix{BigFloat}, Matrix{T}} where {T <: Union{BigFloat, Complex{BigFloat}}} + check_input(svd_full!, A, USVᴴ, alg) + U, S, Vᴴ = USVᴴ + m, n = size(A) + minmn = min(m, n) + if minmn == 0 + MatrixAlgebraKit.one!(U) + MatrixAlgebraKit.zero!(S) + MatrixAlgebraKit.one!(Vᴴ) + return USVᴴ + end + S̃ = zeros(eltype(S), size(S)) + U_compact, S_compact, V_compact = GenericLinearAlgebra.svd(A) + S̃[1:minmn, 1:minmn] .= Diagonal(S_compact) + Ũ = _gram_schmidt(U_compact) + Ṽ = _gram_schmidt(V_compact) + + copyto!(U, Ũ) + copyto!(S, S̃) + copyto!(Vᴴ, Ṽ') + + return MatrixAlgebraKit.gaugefix!(svd_full!, U, S, Vᴴ, m, n) +end + +function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::BigFloat_svd_QRIteration) where {T <: Union{BigFloat, Complex{BigFloat}}} + check_input(svd_vals!, A, S, alg) + S̃ = GenericLinearAlgebra.svdvals!(A) + copyto!(S, S̃) + return S +end + +function MatrixAlgebraKit.default_eig_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} + return BigFloat_eig_Francis(; kwargs...) +end + +function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eig_Francis)::Tuple{Diagonal{Complex{BigFloat}}, Matrix{Complex{BigFloat}}} where {T <: Union{BigFloat, Complex{BigFloat}}} + D, V = DV + D̃, Ṽ = GenericSchur.eigen!(A) + copyto!(D, Diagonal(D̃)) + copyto!(V, Ṽ) + return D, V +end + +function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::BigFloat_eig_Francis)::Vector{Complex{BigFloat}} where {T <: Union{BigFloat, Complex{BigFloat}}} + check_input(eig_vals!, A, D, alg) + eigval = GenericSchur.eigvals!(A) + return eigval +end + + +function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} + return BigFloat_eigh_Francis(; kwargs...) +end + +function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eigh_Francis)::Tuple{Diagonal{BigFloat}, Matrix{T}} where {T <: Union{BigFloat, Complex{BigFloat}}} + check_input(eigh_full!, A, DV, alg) + eigval, eigvec = GenericLinearAlgebra.eigen(A; sortby = λ -> real(λ)) + return Diagonal(eigval), eigvec +end + +function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::BigFloat_eigh_Francis)::Vector{BigFloat} where {T <: Union{BigFloat, Complex{BigFloat}}} + check_input(eigh_vals!, A, D, alg) + D = GenericLinearAlgebra.eigvals(A; sortby = λ -> real(λ)) + return real.(D) +end + +function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} + return BigFloat_QR_Householder(; kwargs...) +end + +function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::BigFloat_QR_Householder) + Q, R = QR + m, n = size(A) + minmn = min(m, n) + computeR = length(R) > 0 + + Q_zero = zeros(eltype(Q), (m, minmn)) + R_zero = zeros(eltype(R), (minmn, n)) + Q_compact, R_compact = _bigfloat_householder_qr!(A, Q_zero, R_zero; alg.kwargs...) + copyto!(Q, _gram_schmidt(Q_compact[:, 1:min(m, n)])) + if computeR + R̃ = zeros(eltype(R), m, n) + R̃[1:minmn, 1:n] .= R_compact + copyto!(R, R̃) + end + return Q, R +end + +function MatrixAlgebraKit.lq_full!(A::AbstractMatrix, LQ, alg::BigFloat_LQ_Householder) + L, Q = LQ + m, n = size(A) + minmn = min(m, n) + computeL = length(L) > 0 + + L_zero = zeros(eltype(L), (m, minmn)) + Q_zero = zeros(eltype(Q), (minmn, n)) + L_compact, Q_compact = _bigfloat_householder_lq!(A, L_zero, Q_zero; alg.kwargs...) + copyto!(Q, _gram_schmidt(Q_compact'[:, 1:min(m, n)])') + if computeL + L̃ = zeros(eltype(L), m, n) + L̃[1:m, 1:minmn] .= L_compact + copyto!(L, L̃) + end + return L, Q +end + +function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::BigFloat_QR_Householder) + check_input(qr_compact!, A, QR, alg) + Q, R = QR + Q, R = _bigfloat_householder_qr!(A, Q, R; alg.kwargs...) + return Q, R +end + +function _bigfloat_householder_qr!(A::AbstractMatrix{T}, Q, R; positive = false, blocksize = 1, pivoted = false) where {T <: Union{BigFloat, Complex{BigFloat}}} + pivoted && throw(ArgumentError("Only pivoted = false implemented for BigFloats.")) + (blocksize == 1) || throw(ArgumentError("Only blocksize = 1 implemented for BigFloats.")) + + m, n = size(A) + k = min(m, n) + computeR = length(R) > 0 + Q̃, R̃ = GenericLinearAlgebra.qr(A) + Q̃ = convert(Array, Q̃) + if positive + @inbounds for j in 1:k + s = sign_safe(R̃[j, j]) + @simd for i in 1:m + Q̃[i, j] *= s + end + end + end + copyto!(Q, Q̃) + if computeR + if positive + @inbounds for j in n:-1:1 + @simd for i in 1:min(k, j) + R̃[i, j] = R̃[i, j] * conj(sign_safe(R̃[i, i])) + end + end + end + copyto!(R, R̃) + end + return Q, R +end + +function _gram_schmidt(Q_compact) + m, minmn = size(Q_compact) + if minmn >= m + return Q_compact + end + Q = zeros(eltype(Q_compact), (m, m)) + Q[:, 1:minmn] .= Q_compact + for j in (minmn + 1):m + v = rand(eltype(Q), m) + for i in 1:(j - 1) + r = sum([v[k] * conj(Q[k, i])] for k in 1:size(v)[1])[1] + v .= v .- r * Q[:, i] + end + Q[:, j] = v ./ MatrixAlgebraKit.norm(v) + end + return Q +end + +function MatrixAlgebraKit.default_lq_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} + return BigFloat_LQ_Householder(; kwargs...) +end + +function MatrixAlgebraKit.lq_compact!(A::AbstractMatrix, LQ, alg::BigFloat_LQ_Householder) + check_input(lq_compact!, A, LQ, alg) + L, Q = LQ + L, Q = _bigfloat_householder_lq!(A, L, Q; alg.kwargs...) + return L, Q +end + + +function _bigfloat_householder_lq!(A::AbstractMatrix{T}, L, Q; positive = false, blocksize = 1, pivoted = false) where {T <: Union{BigFloat, Complex{BigFloat}}} + pivoted && throw(ArgumentError("Only pivoted = false implemented for BigFloats.")) + (blocksize == 1) || throw(ArgumentError("Only blocksize = 1 implemented for BigFloats.")) + + m, n = size(A) + k = min(m, n) + computeL = length(L) > 0 + + Q̃, R̃ = GenericLinearAlgebra.qr(A') + Q̃ = convert(Array, Q̃) + + if positive + @inbounds for j in 1:k + s = sign_safe(R̃[j, j]) + @simd for i in 1:n + Q̃[i, j] *= s + end + end + end + copyto!(Q, Q̃') + if computeL + if positive + @inbounds for j in m:-1:1 + for i in 1:min(k, j) + R̃[i, j] = R̃[i, j] * conj(sign_safe(R̃[i, i])) + end + end + end + copyto!(L, R̃') + + end + return L, Q +end + +end diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index 5211476f..a5480a65 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -33,6 +33,7 @@ export left_orth!, right_orth!, left_null!, right_null! export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert, LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations, LAPACK_DivideAndConquer, LAPACK_Jacobi +export BigFloat_QR_Householder, BigFloat_LQ_Householder, BigFloat_eig_Francis, BigFloat_eigh_Francis, BigFloat_svd_QRIteration export LQViaTransposedQR export PolarViaSVD, PolarNewton export DiagonalAlgorithm diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index bdda6612..7a1b020c 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -33,6 +33,8 @@ elements of `L` are non-negative. @algdef LAPACK_HouseholderLQ # TODO: +@algdef BigFloat_QR_Householder +@algdef BigFloat_LQ_Householder @algdef LAPACK_HouseholderQL @algdef LAPACK_HouseholderRQ @@ -56,6 +58,9 @@ eigenvalue decomposition of a matrix. const LAPACK_EigAlgorithm = Union{LAPACK_Simple, LAPACK_Expert} +# TODO: +@algdef BigFloat_eig_Francis + # Hermitian Eigenvalue Decomposition # ---------------------------------- """ @@ -100,6 +105,9 @@ const LAPACK_EighAlgorithm = Union{ LAPACK_MultipleRelativelyRobustRepresentations, } +# TODO: +@algdef BigFloat_eigh_Francis + # Singular Value Decomposition # ---------------------------- """ @@ -117,6 +125,9 @@ const LAPACK_SVDAlgorithm = Union{ LAPACK_Jacobi, } +# TODO: +@algdef BigFloat_svd_QRIteration + # ========================= # Polar decompositions # ========================= diff --git a/test/bigfloats/eig.jl b/test/bigfloats/eig.jl new file mode 100644 index 00000000..ac244b4a --- /dev/null +++ b/test/bigfloats/eig.jl @@ -0,0 +1,117 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using LinearAlgebra: Diagonal +using MatrixAlgebraKit: TruncatedAlgorithm, diagview, norm +using GenericLinearAlgebra +using GenericSchur + +const eltypes = (BigFloat, Complex{BigFloat}) + +@testset "eig_full! for T = $T" for T in eltypes + rng = StableRNG(123) + m = 24 + alg = BigFloat_eig_Francis() + A = randn(rng, T, m, m) + Tc = complex(T) + + D, V = @constinferred eig_full(A; alg = ($alg)) + @test eltype(D) == eltype(V) == Tc + @test A * V ≈ V * D + + alg′ = @constinferred MatrixAlgebraKit.select_algorithm(eig_full!, A, $alg) + + Ac = similar(A) + D2, V2 = @constinferred eig_full!(copy!(Ac, A), (D, V), alg′) + @test D2 ≈ D + @test V2 ≈ V + @test A * V ≈ V * D + + Dc = @constinferred eig_vals(A, alg′) + @test eltype(Dc) == Tc + @test D ≈ Diagonal(Dc) +end + +@testset "eig_trunc! for T = $T" for T in eltypes + rng = StableRNG(123) + m = 6 + alg = BigFloat_eig_Francis() + A = randn(rng, T, m, m) + A *= A' # TODO: deal with eigenvalue ordering etc + # eigenvalues are sorted by ascending real component... + D₀ = sort!(eig_vals(A); by = abs, rev = true) + rmin = findfirst(i -> abs(D₀[end - i]) != abs(D₀[end - i - 1]), 1:(m - 2)) + r = length(D₀) - rmin + atol = sqrt(eps(real(T))) + + D1, V1, ϵ1 = @constinferred eig_trunc(A; alg, trunc = truncrank(r)) + D1base, V1base = @constinferred eig_full(A; alg) + + @test length(diagview(D1)) == r + @test A * V1 ≈ V1 * D1 + @test ϵ1 ≈ norm(view(D₀, (r + 1):m)) atol = atol + + s = 1 + sqrt(eps(real(T))) + trunc = trunctol(; atol = s * abs(D₀[r + 1])) + D2, V2, ϵ2 = @constinferred eig_trunc(A; alg, trunc) + @test length(diagview(D2)) == r + @test A * V2 ≈ V2 * D2 + @test ϵ2 ≈ norm(view(D₀, (r + 1):m)) atol = atol + + s = 1 - sqrt(eps(real(T))) + trunc = truncerror(; atol = s * norm(@view(D₀[r:end]), 1), p = 1) + D3, V3, ϵ3 = @constinferred eig_trunc(A; alg, trunc) + @test length(diagview(D3)) == r + @test A * V3 ≈ V3 * D3 + @test ϵ3 ≈ norm(view(D₀, (r + 1):m)) atol = atol + + # trunctol keeps order, truncrank might not + # test for same subspace + @test V1 * ((V1' * V1) \ (V1' * V2)) ≈ V2 + @test V2 * ((V2' * V2) \ (V2' * V1)) ≈ V1 + @test V1 * ((V1' * V1) \ (V1' * V3)) ≈ V3 + @test V3 * ((V3' * V3) \ (V3' * V1)) ≈ V1 +end + +@testset "eig_trunc! specify truncation algorithm T = $T" for T in eltypes + rng = StableRNG(123) + m = 4 + atol = sqrt(eps(real(T))) + V = randn(rng, T, m, m) + D = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) + A = V * D * inv(V) + alg = TruncatedAlgorithm(BigFloat_eig_Francis(), truncrank(2)) + D2, V2, ϵ2 = @constinferred eig_trunc(A; alg) + @test diagview(D2) ≈ diagview(D)[1:2] + @test ϵ2 ≈ norm(diagview(D)[3:4]) atol = atol + @test_throws ArgumentError eig_trunc(A; alg, trunc = (; maxrank = 2)) + + alg = TruncatedAlgorithm(BigFloat_eig_Francis(), truncerror(; atol = 0.2, p = 1)) + D3, V3, ϵ3 = @constinferred eig_trunc(A; alg) + @test diagview(D3) ≈ diagview(D)[1:2] + @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol +end + +@testset "eig for Diagonal{$T}" for T in eltypes + rng = StableRNG(123) + m = 24 + Ad = randn(rng, T, m) + A = Diagonal(Ad) + atol = sqrt(eps(real(T))) + + D, V = @constinferred eig_full(A) + @test D isa Diagonal{T} && size(D) == size(A) + @test V isa Diagonal{T} && size(V) == size(A) + @test A * V ≈ V * D + + D2 = @constinferred eig_vals(A) + @test D2 isa AbstractVector{T} && length(D2) == m + @test diagview(D) ≈ D2 + + A2 = Diagonal(T[0.9, 0.3, 0.1, 0.01]) + alg = TruncatedAlgorithm(DiagonalAlgorithm(), truncrank(2)) + D2, V2, ϵ2 = @constinferred eig_trunc(A2; alg) + @test diagview(D2) ≈ diagview(A2)[1:2] + @test ϵ2 ≈ norm(diagview(A2)[3:4]) atol = atol +end diff --git a/test/bigfloats/eigh.jl b/test/bigfloats/eigh.jl new file mode 100644 index 00000000..8042df20 --- /dev/null +++ b/test/bigfloats/eigh.jl @@ -0,0 +1,119 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using LinearAlgebra: LinearAlgebra, Diagonal, I +using MatrixAlgebraKit: TruncatedAlgorithm, diagview, norm +using GenericLinearAlgebra +using GenericSchur + +const eltypes = (BigFloat, Complex{BigFloat}) + +@testset "eigh_full! for T = $T" for T in eltypes + rng = StableRNG(123) + m = 54 + alg = BigFloat_eigh_Francis() + + A = randn(rng, T, m, m) + A = (A + A') / 2 + + D, V = @constinferred eigh_full(A; alg) + @test A * V ≈ V * D + @test isunitary(V) + @test all(isreal, D) + + D2, V2 = eigh_full!(copy(A), (D, V), alg) + @test D2 ≈ D + @test V2 ≈ V + + D3 = @constinferred eigh_vals(A, alg) + @test D ≈ Diagonal(D3) +end + +@testset "eigh_trunc! for T = $T" for T in eltypes + rng = StableRNG(123) + m = 54 + alg = BigFloat_eigh_Francis() + A = randn(rng, T, m, m) + A = A * A' + A = (A + A') / 2 + Ac = similar(A) + D₀ = reverse(eigh_vals(A)) + + r = m - 2 + s = 1 + sqrt(eps(real(T))) + atol = sqrt(eps(real(T))) + + D1, V1, ϵ1 = @constinferred eigh_trunc(A; alg, trunc = truncrank(r)) + Dfull, Vfull = eigh_full(A; alg) + @test length(diagview(D1)) == r + @test isisometric(V1) + @test A * V1 ≈ V1 * D1 + @test LinearAlgebra.opnorm(A - V1 * D1 * V1') ≈ D₀[r + 1] + @test ϵ1 ≈ norm(view(D₀, (r + 1):m)) atol = atol + + trunc = trunctol(; atol = s * D₀[r + 1]) + D2, V2, ϵ2 = @constinferred eigh_trunc(A; alg, trunc) + @test length(diagview(D2)) == r + @test isisometric(V2) + @test A * V2 ≈ V2 * D2 + @test ϵ2 ≈ norm(view(D₀, (r + 1):m)) atol = atol + + s = 1 - sqrt(eps(real(T))) + trunc = truncerror(; atol = s * norm(@view(D₀[r:end]), 1), p = 1) + D3, V3, ϵ3 = @constinferred eigh_trunc(A; alg, trunc) + @test length(diagview(D3)) == r + @test A * V3 ≈ V3 * D3 + @test ϵ3 ≈ norm(view(D₀, (r + 1):m)) atol = atol + + # test for same subspace + @test V1 * (V1' * V2) ≈ V2 + @test V2 * (V2' * V1) ≈ V1 + @test V1 * (V1' * V3) ≈ V3 + @test V3 * (V3' * V1) ≈ V1 +end + +@testset "eigh_trunc! specify truncation algorithm T = $T" for T in eltypes + rng = StableRNG(123) + m = 4 + atol = sqrt(eps(real(T))) + V = qr_compact(randn(rng, T, m, m))[1] + D = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) + A = V * D * V' + A = (A + A') / 2 + alg = TruncatedAlgorithm(BigFloat_eigh_Francis(), truncrank(2)) + D2, V2, ϵ2 = @constinferred eigh_trunc(A; alg) + @test diagview(D2) ≈ diagview(D)[1:2] + @test_throws ArgumentError eigh_trunc(A; alg, trunc = (; maxrank = 2)) + @test ϵ2 ≈ norm(diagview(D)[3:4]) atol = atol + + alg = TruncatedAlgorithm(BigFloat_eigh_Francis(), truncerror(; atol = 0.2)) + D3, V3, ϵ3 = @constinferred eigh_trunc(A; alg) + @test diagview(D3) ≈ diagview(D)[1:2] + @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol +end + +@testset "eigh for Diagonal{$T}" for T in eltypes + rng = StableRNG(123) + m = 54 + Ad = randn(rng, T, m) + Ad .+= conj.(Ad) + A = Diagonal(Ad) + atol = sqrt(eps(real(T))) + + D, V = @constinferred eigh_full(A) + @test D isa Diagonal{real(T)} && size(D) == size(A) + @test V isa Diagonal{T} && size(V) == size(A) + @test A * V ≈ V * D + + D2 = @constinferred eigh_vals(A) + @test D2 isa AbstractVector{real(T)} && length(D2) == m + @test diagview(D) ≈ D2 + + A2 = Diagonal(T[0.9, 0.3, 0.1, 0.01]) + alg = TruncatedAlgorithm(BigFloat_eigh_Francis(), truncrank(2)) + D2, V2, ϵ2 = @constinferred eigh_trunc(A2; alg) + @test diagview(D2) ≈ diagview(A2)[1:2] + @test ϵ2 ≈ norm(diagview(A2)[3:4]) atol = atol + +end diff --git a/test/bigfloats/lq.jl b/test/bigfloats/lq.jl new file mode 100644 index 00000000..c0e4faba --- /dev/null +++ b/test/bigfloats/lq.jl @@ -0,0 +1,150 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using LinearAlgebra: diag, I, Diagonal +using GenericLinearAlgebra +using GenericSchur + +eltypes = (BigFloat, Complex{BigFloat}) + +@testset "qr_compact! for T = $T" for T in eltypes + + rng = StableRNG(123) + m = 54 + for n in (37, m, 63) + minmn = min(m, n) + m = 54 + A = randn(rng, T, m, n) + L, Q = @constinferred lq_compact(A) + @test L isa Matrix{T} && size(L) == (m, minmn) + @test Q isa Matrix{T} && size(Q) == (minmn, n) + @test L * Q ≈ A + @test isisometric(Q; side = :right) + + Ac = similar(A) + L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q)) + @test L2 === L + @test Q2 === Q + + noL = similar(A, 0, minmn) + Q2 = similar(Q) + lq_compact!(copy!(Ac, A), (noL, Q2)) + @test Q == Q2 + + # Transposed QR algorithm + qr_alg = BigFloat_QR_Householder() + lq_alg = LQViaTransposedQR(qr_alg) + L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q), lq_alg) + @test L2 === L + @test Q2 === Q + noL = similar(A, 0, minmn) + Q2 = similar(Q) + lq_compact!(copy!(Ac, A), (noL, Q2), lq_alg) + @test Q == Q2 + + @test_throws ArgumentError lq_compact(A; blocksize = 2) + @test_throws ArgumentError lq_compact(A; pivoted = true) + + # positive + lq_compact!(copy!(Ac, A), (L, Q); positive = true) + @test L * Q ≈ A + @test isisometric(Q; side = :right) + @test all(>=(zero(real(T))), real(diag(L))) + lq_compact!(copy!(Ac, A), (noL, Q2); positive = true) + @test Q == Q2 + end +end + +@testset "lq_full! for T = $T" for T in eltypes + rng = StableRNG(123) + m = 54 + for n in (37, m, 63) + minmn = min(m, n) + A = randn(rng, T, m, n) + L, Q = lq_full(A) + @test L isa Matrix{T} && size(L) == (m, n) + @test Q isa Matrix{T} && size(Q) == (n, n) + @test L * Q ≈ A + @test isunitary(Q) + + Ac = similar(A) + L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q)) + @test L2 === L + @test Q2 === Q + @test L * Q ≈ A + @test isunitary(Q) + + noL = similar(A, 0, n) + Q2 = similar(Q) + lq_full!(copy!(Ac, A), (noL, Q2)) + @test Q[1:minmn, n] ≈ Q2[1:minmn, n] + + # Transposed QR algorithm + qr_alg = BigFloat_QR_Householder() + lq_alg = LQViaTransposedQR(qr_alg) + L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q), lq_alg) + @test L2 === L + @test Q2 === Q + @test L * Q ≈ A + @test Q * Q' ≈ I + noL = similar(A, 0, n) + Q2 = similar(Q) + lq_full!(copy!(Ac, A), (noL, Q2), lq_alg) + @test Q[1:minmn, n] ≈ Q2[1:minmn, n] + + # Argument errors for unsupported options + @test_throws ArgumentError lq_full(A; blocksize = 2) + @test_throws ArgumentError lq_full(A; pivoted = true) + + # positive + lq_full!(copy!(Ac, A), (L, Q); positive = true) + @test L * Q ≈ A + @test isunitary(Q) + @test all(>=(zero(real(T))), real(diag(L))) + lq_full!(copy!(Ac, A), (noL, Q2); positive = true) + @test Q[1:minmn, n] ≈ Q2[1:minmn, n] + + qr_alg = BigFloat_QR_Householder(; positive = true) + lq_alg = LQViaTransposedQR(qr_alg) + lq_full!(copy!(Ac, A), (L, Q), lq_alg) + @test L * Q ≈ A + @test Q * Q' ≈ I + @test all(>=(zero(real(T))), real(diag(L))) + lq_full!(copy!(Ac, A), (noL, Q2), lq_alg) + @test Q[1:minmn, n] ≈ Q2[1:minmn, n] + + # positive and blocksize 1 + lq_full!(copy!(Ac, A), (L, Q); positive = true, blocksize = 1) + @test L * Q ≈ A + @test isunitary(Q) + @test all(>=(zero(real(T))), real(diag(L))) + lq_full!(copy!(Ac, A), (noL, Q2); positive = true, blocksize = 1) + @test Q[1:minmn, n] ≈ Q2[1:minmn, n] + end +end + +@testset "lq_compact for Diagonal{$T}" for T in eltypes + rng = StableRNG(123) + atol = eps(real(T))^(3 / 4) + for m in (54, 0) + Ad = randn(rng, T, m) + A = Diagonal(Ad) + + # compact + L, Q = @constinferred lq_compact(A) + @test Q isa Diagonal{T} && size(Q) == (m, m) + @test L isa Diagonal{T} && size(L) == (m, m) + @test L * Q ≈ A + @test isunitary(Q) + + # compact and positive + Lp, Qp = @constinferred lq_compact(A; positive = true) + @test Qp isa Diagonal{T} && size(Qp) == (m, m) + @test Lp isa Diagonal{T} && size(Lp) == (m, m) + @test Lp * Qp ≈ A + @test isunitary(Qp) + @test all(≥(zero(real(T))), real(diag(Lp))) && + all(≈(zero(real(T)); atol), imag(diag(Lp))) + end +end diff --git a/test/bigfloats/qr.jl b/test/bigfloats/qr.jl new file mode 100644 index 00000000..32c8e375 --- /dev/null +++ b/test/bigfloats/qr.jl @@ -0,0 +1,135 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using LinearAlgebra: diag, I, Diagonal +using GenericLinearAlgebra +using GenericSchur + +eltypes = (BigFloat, Complex{BigFloat}) + +@testset "qr_compact! for T = $T" for T in eltypes + + rng = StableRNG(123) + m = 54 + for n in (37, m, 63) + minmn = min(m, n) + m = 54 + A = randn(rng, T, m, n) + Q, R = @constinferred qr_compact(A) + @test Q isa Matrix{T} && size(Q) == (m, minmn) + @test R isa Matrix{T} && size(R) == (minmn, n) + @test Q * R ≈ A + + Ac = similar(A) + Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R)) + @test Q2 === Q + @test R2 === R + + Q2 = similar(Q) + noR = similar(A, minmn, 0) + qr_compact!(copy!(Ac, A), (Q2, noR)) + @test Q == Q2 + + @test_throws ArgumentError qr_compact(A; blocksize = 2) + @test_throws ArgumentError qr_compact(A; pivoted = true) + + # positive + qr_compact!(copy!(Ac, A), (Q, R); positive = true) + @test Q * R ≈ A + @test isisometric(Q) + @test all(>=(zero(real(T))), real(diag(R))) + qr_compact!(copy!(Ac, A), (Q2, noR); positive = true) + @test Q == Q2 + end +end + +@testset "qr_full! for T = $T" for T in eltypes + rng = StableRNG(123) + m = 54 + for n in (37, m, 63) + minmn = min(m, n) + A = randn(rng, T, m, n) + Q, R = qr_full(A) + @test Q isa Matrix{T} && size(Q) == (m, m) + @test R isa Matrix{T} && size(R) == (m, n) + Qc, Rc = qr_compact(A) + @test Q * R ≈ A + @test isunitary(Q) + + Ac = similar(A) + Q2 = similar(Q) + noR = similar(A, m, 0) + Q2, R2 = @constinferred qr_full!(copy!(Ac, A), (Q, R)) + @test Q2 === Q + @test R2 === R + @test Q * R ≈ A + @test isunitary(Q) + qr_full!(copy!(Ac, A), (Q2, noR)) + @test Q == Q2 + + # unblocked algorithm + qr_full!(copy!(Ac, A), (Q, R); blocksize = 1) + @test Q * R ≈ A + @test isunitary(Q) + qr_full!(copy!(Ac, A), (Q2, noR); blocksize = 1) + @test Q == Q2 + if n == m + qr_full!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q + @test Q ≈ Q2 + end + + # Argument errors for unsupported options + @test_throws ArgumentError qr_full(A; blocksize = 2) + @test_throws ArgumentError qr_compact(A; pivoted = true) + + # positive + qr_full!(copy!(Ac, A), (Q, R); positive = true) + @test Q * R ≈ A + @test isunitary(Q) + @test all(>=(zero(real(T))), real(diag(R))) + qr_full!(copy!(Ac, A), (Q2, noR); positive = true) + @test Q == Q2 + # positive and blocksize 1 + qr_full!(copy!(Ac, A), (Q, R); positive = true, blocksize = 1) + @test Q * R ≈ A + @test isunitary(Q) + @test all(>=(zero(real(T))), real(diag(R))) + qr_full!(copy!(Ac, A), (Q2, noR); positive = true, blocksize = 1) + @test Q == Q2 + if n <= m + # the following test tries to find the diagonal element (in order to test positivity) + # before the column permutation. This only works if all columns have a diagonal + # element + for j in 1:n + i = findlast(!iszero, view(R, :, j)) + @test real(R[i, j]) >= zero(real(T)) + end + end + end +end + +@testset "qr_compact for Diagonal{$T}" for T in eltypes + rng = StableRNG(123) + atol = eps(real(T))^(3 / 4) + for m in (54, 0) + Ad = randn(rng, T, m) + A = Diagonal(Ad) + + # compact + Q, R = @constinferred qr_compact(A) + @test Q isa Diagonal{T} && size(Q) == (m, m) + @test R isa Diagonal{T} && size(R) == (m, m) + @test Q * R ≈ A + @test isunitary(Q) + + # compact and positive + Qp, Rp = @constinferred qr_compact(A; positive = true) + @test Qp isa Diagonal{T} && size(Qp) == (m, m) + @test Rp isa Diagonal{T} && size(Rp) == (m, m) + @test Qp * Rp ≈ A + @test isunitary(Qp) + @test all(≥(zero(real(T))), real(diag(Rp))) && + all(≈(zero(real(T)); atol), imag(diag(Rp))) + end +end diff --git a/test/bigfloats/svd.jl b/test/bigfloats/svd.jl new file mode 100644 index 00000000..5e616d37 --- /dev/null +++ b/test/bigfloats/svd.jl @@ -0,0 +1,210 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using LinearAlgebra: LinearAlgebra, Diagonal, I, isposdef, norm +using MatrixAlgebraKit: TruncatedAlgorithm, diagview, isisometric +using GenericLinearAlgebra +using GenericSchur + +eltypes = (BigFloat, Complex{BigFloat}) + +@testset "svd_compact! for T = $T" for T in eltypes + rng = StableRNG(123) + m = 54 + @testset "size ($m, $n)" for n in (37, m, 63, 0) + k = min(m, n) + alg = BigFloat_svd_QRIteration() + minmn = min(m, n) + A = randn(rng, T, m, n) + + if VERSION < v"1.11" + # This is type unstable on older versions of Julia. + U, S, Vᴴ = svd_compact(A; alg) + else + U, S, Vᴴ = @constinferred svd_compact(A; alg = ($alg)) + end + @test U isa Matrix{T} && size(U) == (m, minmn) + @test S isa Diagonal{real(T)} && size(S) == (minmn, minmn) + @test Vᴴ isa Matrix{T} && size(Vᴴ) == (minmn, n) + @test U * S * Vᴴ ≈ A + @test isisometric(U) + @test isisometric(Vᴴ; side = :right) + @test isposdef(S) + + Ac = similar(A) + Sc = similar(A, real(T), min(m, n)) + alg′ = @constinferred MatrixAlgebraKit.select_algorithm(svd_compact!, A, $alg) + U2, S2, V2ᴴ = @constinferred svd_compact!(copy!(Ac, A), (U, S, Vᴴ), alg′) + @test U2 ≈ U + @test S2 ≈ S + @test V2ᴴ ≈ Vᴴ + @test U * S * Vᴴ ≈ A + @test isisometric(U) + @test isisometric(Vᴴ; side = :right) + @test isposdef(S) + + Sd = @constinferred svd_vals(A, alg′) + @test S ≈ Diagonal(Sd) + end +end + +@testset "svd_full! for T = $T" for T in eltypes + rng = StableRNG(123) + m = 54 + @testset "size ($m, $n)" for n in (37, m, 63, 0) + alg = BigFloat_svd_QRIteration() + A = randn(rng, T, m, n) + U, S, Vᴴ = svd_full(A; alg) + @test U isa Matrix{T} && size(U) == (m, m) + @test S isa Matrix{real(T)} && size(S) == (m, n) + @test Vᴴ isa Matrix{T} && size(Vᴴ) == (n, n) + @test U * S * Vᴴ ≈ A + @test isunitary(U) + @test isunitary(Vᴴ) + @test all(isposdef, diagview(S)) + + Ac = similar(A) + U2, S2, V2ᴴ = @constinferred svd_full!(copy!(Ac, A), (U, S, Vᴴ), alg) + @test U2 === U + @test S2 === S + @test V2ᴴ === Vᴴ + @test U * S * Vᴴ ≈ A + @test isunitary(U) + @test isunitary(Vᴴ) + @test all(isposdef, diagview(S)) + + Sc = similar(A, real(T), min(m, n)) + Sc2 = svd_vals!(copy!(Ac, A), Sc, alg) + @test Sc === Sc2 + @test diagview(S) ≈ Sc + end + @testset "size (0, 0)" begin + @testset "algorithm $alg" for alg in + (LAPACK_DivideAndConquer(), LAPACK_QRIteration()) + A = randn(rng, T, 0, 0) + U, S, Vᴴ = svd_full(A; alg) + @test U isa Matrix{T} && size(U) == (0, 0) + @test S isa Matrix{real(T)} && size(S) == (0, 0) + @test Vᴴ isa Matrix{T} && size(Vᴴ) == (0, 0) + @test U * S * Vᴴ ≈ A + @test isunitary(U) + @test isunitary(Vᴴ) + @test all(isposdef, diagview(S)) + end + end +end + +@testset "svd_trunc! for T = $T" for T in eltypes + rng = StableRNG(123) + m = 54 + atol = sqrt(eps(real(T))) + alg = BigFloat_svd_QRIteration() + + @testset "size ($m, $n)" for n in (37, m, 63) + n > m && alg isa LAPACK_Jacobi && continue # not supported + A = randn(rng, T, m, n) + S₀ = svd_vals(A) + minmn = min(m, n) + r = minmn - 2 + + U1, S1, V1ᴴ, ϵ1 = @constinferred svd_trunc(A; alg, trunc = truncrank(r)) + @test length(diagview(S1)) == r + @test diagview(S1) ≈ S₀[1:r] + @test LinearAlgebra.opnorm(A - U1 * S1 * V1ᴴ) ≈ S₀[r + 1] + # Test truncation error + @test ϵ1 ≈ norm(view(S₀, (r + 1):minmn)) atol = atol + + s = 1 + sqrt(eps(real(T))) + trunc = trunctol(; atol = s * S₀[r + 1]) + + U2, S2, V2ᴴ, ϵ2 = @constinferred svd_trunc(A; alg, trunc) + @test length(diagview(S2)) == r + @test U1 ≈ U2 + @test S1 ≈ S2 + @test V1ᴴ ≈ V2ᴴ + @test ϵ2 ≈ norm(view(S₀, (r + 1):minmn)) atol = atol + + trunc = truncerror(; atol = s * norm(@view(S₀[(r + 1):end]))) + U3, S3, V3ᴴ, ϵ3 = @constinferred svd_trunc(A; alg, trunc) + @test length(diagview(S3)) == r + @test U1 ≈ U3 + @test S1 ≈ S3 + @test V1ᴴ ≈ V3ᴴ + @test ϵ3 ≈ norm(view(S₀, (r + 1):minmn)) atol = atol + end +end + +@testset "svd_trunc! mix maxrank and tol for T = $T" for T in eltypes + rng = StableRNG(123) + alg = BigFloat_svd_QRIteration() + m = 4 + U = qr_compact(randn(rng, T, m, m))[1] + S = Diagonal(T[0.9, 0.3, 0.1, 0.01]) + Vᴴ = qr_compact(randn(rng, T, m, m))[1] + A = U * S * Vᴴ + + for trunc_fun in ( + (rtol, maxrank) -> (; rtol, maxrank), + (rtol, maxrank) -> truncrank(maxrank) & trunctol(; rtol), + ) + U1, S1, V1ᴴ, ϵ1 = svd_trunc(A; alg, trunc = trunc_fun(0.2, 1)) + @test length(diagview(S1)) == 1 + @test diagview(S1) ≈ diagview(S)[1:1] + + U2, S2, V2ᴴ, ϵ2 = svd_trunc(A; alg, trunc = trunc_fun(0.2, 3)) + @test length(diagview(S2)) == 2 + @test diagview(S2) ≈ diagview(S)[1:2] + end +end + +@testset "svd_trunc! specify truncation algorithm T = $T" for T in eltypes + rng = StableRNG(123) + atol = sqrt(eps(real(T))) + m = 4 + U = qr_compact(randn(rng, T, m, m))[1] + S = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) + Vᴴ = qr_compact(randn(rng, T, m, m))[1] + A = U * S * Vᴴ + alg = TruncatedAlgorithm(BigFloat_svd_QRIteration(), trunctol(; atol = 0.2)) + U2, S2, V2ᴴ, ϵ2 = @constinferred svd_trunc(A; alg) + @test diagview(S2) ≈ diagview(S)[1:2] + @test ϵ2 ≈ norm(diagview(S)[3:4]) atol = atol + @test_throws ArgumentError svd_trunc(A; alg, trunc = (; maxrank = 2)) +end + +@testset "svd for Diagonal{$T}" for T in eltypes + rng = StableRNG(123) + atol = sqrt(eps(real(T))) + for m in (54, 0) + Ad = randn(T, m) + A = Diagonal(Ad) + + U, S, Vᴴ = @constinferred svd_compact(A) + @test U isa AbstractMatrix{T} && size(U) == size(A) + @test Vᴴ isa AbstractMatrix{T} && size(Vᴴ) == size(A) + @test S isa Diagonal{real(T)} && size(S) == size(A) + @test isunitary(U) + @test isunitary(Vᴴ) + @test all(≥(0), diagview(S)) + @test A ≈ U * S * Vᴴ + + U, S, Vᴴ = @constinferred svd_full(A) + @test U isa AbstractMatrix{T} && size(U) == size(A) + @test Vᴴ isa AbstractMatrix{T} && size(Vᴴ) == size(A) + @test S isa Diagonal{real(T)} && size(S) == size(A) + @test isunitary(U) + @test isunitary(Vᴴ) + @test all(≥(0), diagview(S)) + @test A ≈ U * S * Vᴴ + + S2 = @constinferred svd_vals(A) + @test S2 isa AbstractVector{real(T)} && length(S2) == m + @test S2 ≈ diagview(S) + + alg = TruncatedAlgorithm(DiagonalAlgorithm(), truncrank(2)) + U3, S3, Vᴴ3, ϵ3 = @constinferred svd_trunc(A; alg) + @test diagview(S3) ≈ S2[1:min(m, 2)] + @test ϵ3 ≈ norm(S2[(min(m, 2) + 1):m]) atol = atol + end +end diff --git a/test/runtests.jl b/test/runtests.jl index af4996ed..8b50996f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -107,3 +107,17 @@ if AMDGPU.functional() include("amd/orthnull.jl") end end + +@safetestset "QR / LQ Decomposition" begin + include("bigfloats/qr.jl") + include("bigfloats/lq.jl") +end +@safetestset "Singular Value Decomposition" begin + include("bigfloats/svd.jl") +end +@safetestset "Hermitian Eigenvalue Decomposition" begin + include("bigfloats/eigh.jl") +end +@safetestset "General Eigenvalue Decomposition" begin + include("bigfloats/eig.jl") +end From 5e08f7aa1cd5b58f33f6d6c46ddeeb94bd893a53 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Wed, 29 Oct 2025 14:24:47 +0100 Subject: [PATCH 02/25] include weak dependencies --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 8b50996f..7a6bb6b4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -108,6 +108,7 @@ if AMDGPU.functional() end end +using GenericLinearAlgebra, GenericSchur @safetestset "QR / LQ Decomposition" begin include("bigfloats/qr.jl") include("bigfloats/lq.jl") From 92534bcd190a004339dbc80b406074722a769939 Mon Sep 17 00:00:00 2001 From: Sander De Meyer <74001142+sanderdemeyer@users.noreply.github.com> Date: Thu, 30 Oct 2025 11:24:05 +0100 Subject: [PATCH 03/25] Update ext/MatrixAlgebraKitGenericExt.jl Co-authored-by: Lukas Devos --- ext/MatrixAlgebraKitGenericExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MatrixAlgebraKitGenericExt.jl b/ext/MatrixAlgebraKitGenericExt.jl index b95741ea..07d40380 100644 --- a/ext/MatrixAlgebraKitGenericExt.jl +++ b/ext/MatrixAlgebraKitGenericExt.jl @@ -75,7 +75,7 @@ end function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eigh_Francis)::Tuple{Diagonal{BigFloat}, Matrix{T}} where {T <: Union{BigFloat, Complex{BigFloat}}} check_input(eigh_full!, A, DV, alg) - eigval, eigvec = GenericLinearAlgebra.eigen(A; sortby = λ -> real(λ)) + eigval, eigvec = GenericLinearAlgebra.eigen(A; sortby = real) return Diagonal(eigval), eigvec end From ddaf055eb16d9c0188cc2d93a6e93241596ade85 Mon Sep 17 00:00:00 2001 From: Sander De Meyer <74001142+sanderdemeyer@users.noreply.github.com> Date: Thu, 30 Oct 2025 11:24:12 +0100 Subject: [PATCH 04/25] Update ext/MatrixAlgebraKitGenericExt.jl Co-authored-by: Lukas Devos --- ext/MatrixAlgebraKitGenericExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MatrixAlgebraKitGenericExt.jl b/ext/MatrixAlgebraKitGenericExt.jl index 07d40380..88bfd1a0 100644 --- a/ext/MatrixAlgebraKitGenericExt.jl +++ b/ext/MatrixAlgebraKitGenericExt.jl @@ -81,7 +81,7 @@ end function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::BigFloat_eigh_Francis)::Vector{BigFloat} where {T <: Union{BigFloat, Complex{BigFloat}}} check_input(eigh_vals!, A, D, alg) - D = GenericLinearAlgebra.eigvals(A; sortby = λ -> real(λ)) + D = GenericLinearAlgebra.eigvals(A; sortby = real) return real.(D) end From 88b9902a8fa1f102453167e6809d3258529fb9a5 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Thu, 30 Oct 2025 16:16:42 +0100 Subject: [PATCH 05/25] comments from Lukas --- ext/MatrixAlgebraKitGenericExt.jl | 96 +++++++------------------------ src/interface/decompositions.jl | 1 - 2 files changed, 21 insertions(+), 76 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericExt.jl b/ext/MatrixAlgebraKitGenericExt.jl index 88bfd1a0..1c27f409 100644 --- a/ext/MatrixAlgebraKitGenericExt.jl +++ b/ext/MatrixAlgebraKitGenericExt.jl @@ -19,7 +19,7 @@ function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix{T}, USVᴴ, alg::BigFlo return U, Diagonal(S), V' # conjugation to account for difference in convention end -function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::BigFloat_svd_QRIteration)::Tuple{Matrix{T}, Matrix{BigFloat}, Matrix{T}} where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::BigFloat_svd_QRIteration) where {T <: Union{BigFloat, Complex{BigFloat}}} check_input(svd_full!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ m, n = size(A) @@ -30,15 +30,13 @@ function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::BigFloat_ MatrixAlgebraKit.one!(Vᴴ) return USVᴴ end - S̃ = zeros(eltype(S), size(S)) + S̃ = fill!(S, zero(T)) U_compact, S_compact, V_compact = GenericLinearAlgebra.svd(A) S̃[1:minmn, 1:minmn] .= Diagonal(S_compact) - Ũ = _gram_schmidt(U_compact) - Ṽ = _gram_schmidt(V_compact) - - copyto!(U, Ũ) copyto!(S, S̃) - copyto!(Vᴴ, Ṽ') + + U = _gram_schmidt!(U, U_compact) + Vᴴ = _gram_schmidt!(Vᴴ, V_compact; adjoint = true) return MatrixAlgebraKit.gaugefix!(svd_full!, U, S, Vᴴ, m, n) end @@ -54,7 +52,7 @@ function MatrixAlgebraKit.default_eig_algorithm(::Type{T}; kwargs...) where {T < return BigFloat_eig_Francis(; kwargs...) end -function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eig_Francis)::Tuple{Diagonal{Complex{BigFloat}}, Matrix{Complex{BigFloat}}} where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eig_Francis) where {T <: Union{BigFloat, Complex{BigFloat}}} D, V = DV D̃, Ṽ = GenericSchur.eigen!(A) copyto!(D, Diagonal(D̃)) @@ -62,13 +60,12 @@ function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eig_ return D, V end -function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::BigFloat_eig_Francis)::Vector{Complex{BigFloat}} where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::BigFloat_eig_Francis) where {T <: Union{BigFloat, Complex{BigFloat}}} check_input(eig_vals!, A, D, alg) eigval = GenericSchur.eigvals!(A) return eigval end - function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} return BigFloat_eigh_Francis(; kwargs...) end @@ -79,7 +76,7 @@ function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eig return Diagonal(eigval), eigvec end -function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::BigFloat_eigh_Francis)::Vector{BigFloat} where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::BigFloat_eigh_Francis) where {T <: Union{BigFloat, Complex{BigFloat}}} check_input(eigh_vals!, A, D, alg) D = GenericLinearAlgebra.eigvals(A; sortby = real) return real.(D) @@ -98,33 +95,14 @@ function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::BigFloat_QR_House Q_zero = zeros(eltype(Q), (m, minmn)) R_zero = zeros(eltype(R), (minmn, n)) Q_compact, R_compact = _bigfloat_householder_qr!(A, Q_zero, R_zero; alg.kwargs...) - copyto!(Q, _gram_schmidt(Q_compact[:, 1:min(m, n)])) + Q = _gram_schmidt!(Q, Q_compact[:, 1:min(m, n)]) if computeR - R̃ = zeros(eltype(R), m, n) - R̃[1:minmn, 1:n] .= R_compact - copyto!(R, R̃) + R = fill!(R, zero(eltype(R))) + R[1:minmn, 1:n] .= R_compact end return Q, R end -function MatrixAlgebraKit.lq_full!(A::AbstractMatrix, LQ, alg::BigFloat_LQ_Householder) - L, Q = LQ - m, n = size(A) - minmn = min(m, n) - computeL = length(L) > 0 - - L_zero = zeros(eltype(L), (m, minmn)) - Q_zero = zeros(eltype(Q), (minmn, n)) - L_compact, Q_compact = _bigfloat_householder_lq!(A, L_zero, Q_zero; alg.kwargs...) - copyto!(Q, _gram_schmidt(Q_compact'[:, 1:min(m, n)])') - if computeL - L̃ = zeros(eltype(L), m, n) - L̃[1:m, 1:minmn] .= L_compact - copyto!(L, L̃) - end - return L, Q -end - function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::BigFloat_QR_Householder) check_input(qr_compact!, A, QR, alg) Q, R = QR @@ -181,50 +159,18 @@ function _gram_schmidt(Q_compact) return Q end -function MatrixAlgebraKit.default_lq_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return BigFloat_LQ_Householder(; kwargs...) -end - -function MatrixAlgebraKit.lq_compact!(A::AbstractMatrix, LQ, alg::BigFloat_LQ_Householder) - check_input(lq_compact!, A, LQ, alg) - L, Q = LQ - L, Q = _bigfloat_householder_lq!(A, L, Q; alg.kwargs...) - return L, Q -end - - -function _bigfloat_householder_lq!(A::AbstractMatrix{T}, L, Q; positive = false, blocksize = 1, pivoted = false) where {T <: Union{BigFloat, Complex{BigFloat}}} - pivoted && throw(ArgumentError("Only pivoted = false implemented for BigFloats.")) - (blocksize == 1) || throw(ArgumentError("Only blocksize = 1 implemented for BigFloats.")) - - m, n = size(A) - k = min(m, n) - computeL = length(L) > 0 - - Q̃, R̃ = GenericLinearAlgebra.qr(A') - Q̃ = convert(Array, Q̃) - - if positive - @inbounds for j in 1:k - s = sign_safe(R̃[j, j]) - @simd for i in 1:n - Q̃[i, j] *= s - end - end +function _gram_schmidt!(Q, Q_compact; adjoint = false) + Q̃ = _gram_schmidt(Q_compact) + if adjoint + copyto!(Q, Q̃') + else + copyto!(Q, Q̃) end - copyto!(Q, Q̃') - if computeL - if positive - @inbounds for j in m:-1:1 - for i in 1:min(k, j) - R̃[i, j] = R̃[i, j] * conj(sign_safe(R̃[i, i])) - end - end - end - copyto!(L, R̃') + return Q +end - end - return L, Q +function MatrixAlgebraKit.default_lq_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} + return MatrixAlgebraKit.LQViaTransposedQR(BigFloat_QR_Householder(; kwargs...)) end end diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index 7a1b020c..a37dcb05 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -34,7 +34,6 @@ elements of `L` are non-negative. # TODO: @algdef BigFloat_QR_Householder -@algdef BigFloat_LQ_Householder @algdef LAPACK_HouseholderQL @algdef LAPACK_HouseholderRQ From 3eee758d3de9b2014bfe835af4504f4815dbc142 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Thu, 30 Oct 2025 16:36:01 +0100 Subject: [PATCH 06/25] remove `BigFloat_LQ_Householder` --- src/MatrixAlgebraKit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index a5480a65..943959db 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -33,7 +33,7 @@ export left_orth!, right_orth!, left_null!, right_null! export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert, LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations, LAPACK_DivideAndConquer, LAPACK_Jacobi -export BigFloat_QR_Householder, BigFloat_LQ_Householder, BigFloat_eig_Francis, BigFloat_eigh_Francis, BigFloat_svd_QRIteration +export BigFloat_QR_Householder, BigFloat_eig_Francis, BigFloat_eigh_Francis, BigFloat_svd_QRIteration export LQViaTransposedQR export PolarViaSVD, PolarNewton export DiagonalAlgorithm From f8d331ea8ad78e7f79c9c2f519a6f88e5cf03daf Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Thu, 30 Oct 2025 16:52:54 +0100 Subject: [PATCH 07/25] Change struct names and relax type restrictions --- ext/MatrixAlgebraKitGenericExt.jl | 39 ++++++++++++++++--------------- src/MatrixAlgebraKit.jl | 2 +- src/interface/decompositions.jl | 8 +++---- test/bigfloats/eig.jl | 8 +++---- test/bigfloats/eigh.jl | 10 ++++---- test/bigfloats/lq.jl | 6 ++--- test/bigfloats/svd.jl | 10 ++++---- 7 files changed, 42 insertions(+), 41 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericExt.jl b/ext/MatrixAlgebraKitGenericExt.jl index 1c27f409..56c8af3a 100644 --- a/ext/MatrixAlgebraKitGenericExt.jl +++ b/ext/MatrixAlgebraKitGenericExt.jl @@ -10,16 +10,16 @@ using GenericSchur using LinearAlgebra: I, Diagonal function MatrixAlgebraKit.default_svd_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return BigFloat_svd_QRIteration() + return GLA_svd_QRIteration() end -function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix{T}, USVᴴ, alg::BigFloat_svd_QRIteration) where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_QRIteration) where {T} check_input(svd_compact!, A, USVᴴ, alg) U, S, V = GenericLinearAlgebra.svd(A) return U, Diagonal(S), V' # conjugation to account for difference in convention end -function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::BigFloat_svd_QRIteration) where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_QRIteration) where {T} check_input(svd_full!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ m, n = size(A) @@ -41,7 +41,7 @@ function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::BigFloat_ return MatrixAlgebraKit.gaugefix!(svd_full!, U, S, Vᴴ, m, n) end -function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::BigFloat_svd_QRIteration) where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::GLA_svd_QRIteration) where {T} check_input(svd_vals!, A, S, alg) S̃ = GenericLinearAlgebra.svdvals!(A) copyto!(S, S̃) @@ -49,10 +49,11 @@ function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::BigFloat_svd_Q end function MatrixAlgebraKit.default_eig_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return BigFloat_eig_Francis(; kwargs...) + return GLA_eig_Francis(; kwargs...) end -function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eig_Francis) where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GLA_eig_Francis) where {T} + check_input(eig_full!, A, DV, alg) D, V = DV D̃, Ṽ = GenericSchur.eigen!(A) copyto!(D, Diagonal(D̃)) @@ -60,33 +61,33 @@ function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eig_ return D, V end -function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::BigFloat_eig_Francis) where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::GLA_eig_Francis) where {T} check_input(eig_vals!, A, D, alg) eigval = GenericSchur.eigvals!(A) return eigval end function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return BigFloat_eigh_Francis(; kwargs...) + return GLA_eigh_Francis(; kwargs...) end -function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::BigFloat_eigh_Francis)::Tuple{Diagonal{BigFloat}, Matrix{T}} where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_eigh_Francis)::Tuple{Diagonal{real(T)}, Matrix{T}} where {T} check_input(eigh_full!, A, DV, alg) eigval, eigvec = GenericLinearAlgebra.eigen(A; sortby = real) return Diagonal(eigval), eigvec end -function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::BigFloat_eigh_Francis) where {T <: Union{BigFloat, Complex{BigFloat}}} +function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::GLA_eigh_Francis) where {T} check_input(eigh_vals!, A, D, alg) D = GenericLinearAlgebra.eigvals(A; sortby = real) return real.(D) end function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return BigFloat_QR_Householder(; kwargs...) + return GLA_QR_Householder(; kwargs...) end -function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::BigFloat_QR_Householder) +function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::GLA_QR_Householder) Q, R = QR m, n = size(A) minmn = min(m, n) @@ -94,7 +95,7 @@ function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::BigFloat_QR_House Q_zero = zeros(eltype(Q), (m, minmn)) R_zero = zeros(eltype(R), (minmn, n)) - Q_compact, R_compact = _bigfloat_householder_qr!(A, Q_zero, R_zero; alg.kwargs...) + Q_compact, R_compact = _gla_householder_qr!(A, Q_zero, R_zero; alg.kwargs...) Q = _gram_schmidt!(Q, Q_compact[:, 1:min(m, n)]) if computeR R = fill!(R, zero(eltype(R))) @@ -103,16 +104,16 @@ function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::BigFloat_QR_House return Q, R end -function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::BigFloat_QR_Householder) +function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::GLA_QR_Householder) check_input(qr_compact!, A, QR, alg) Q, R = QR - Q, R = _bigfloat_householder_qr!(A, Q, R; alg.kwargs...) + Q, R = _gla_householder_qr!(A, Q, R; alg.kwargs...) return Q, R end -function _bigfloat_householder_qr!(A::AbstractMatrix{T}, Q, R; positive = false, blocksize = 1, pivoted = false) where {T <: Union{BigFloat, Complex{BigFloat}}} - pivoted && throw(ArgumentError("Only pivoted = false implemented for BigFloats.")) - (blocksize == 1) || throw(ArgumentError("Only blocksize = 1 implemented for BigFloats.")) +function _gla_householder_qr!(A::AbstractMatrix{T}, Q, R; positive = false, blocksize = 1, pivoted = false) where {T} + pivoted && throw(ArgumentError("Only pivoted = false implemented for GLA_QR_Householder.")) + (blocksize == 1) || throw(ArgumentError("Only blocksize = 1 implemented for GLA_QR_Householder.")) m, n = size(A) k = min(m, n) @@ -170,7 +171,7 @@ function _gram_schmidt!(Q, Q_compact; adjoint = false) end function MatrixAlgebraKit.default_lq_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return MatrixAlgebraKit.LQViaTransposedQR(BigFloat_QR_Householder(; kwargs...)) + return MatrixAlgebraKit.LQViaTransposedQR(GLA_QR_Householder(; kwargs...)) end end diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index 943959db..c73f2269 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -33,7 +33,7 @@ export left_orth!, right_orth!, left_null!, right_null! export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert, LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations, LAPACK_DivideAndConquer, LAPACK_Jacobi -export BigFloat_QR_Householder, BigFloat_eig_Francis, BigFloat_eigh_Francis, BigFloat_svd_QRIteration +export GLA_QR_Householder, GLA_eig_Francis, GLA_eigh_Francis, GLA_svd_QRIteration export LQViaTransposedQR export PolarViaSVD, PolarNewton export DiagonalAlgorithm diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index a37dcb05..d138e701 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -33,7 +33,7 @@ elements of `L` are non-negative. @algdef LAPACK_HouseholderLQ # TODO: -@algdef BigFloat_QR_Householder +@algdef GLA_QR_Householder @algdef LAPACK_HouseholderQL @algdef LAPACK_HouseholderRQ @@ -58,7 +58,7 @@ eigenvalue decomposition of a matrix. const LAPACK_EigAlgorithm = Union{LAPACK_Simple, LAPACK_Expert} # TODO: -@algdef BigFloat_eig_Francis +@algdef GLA_eig_Francis # Hermitian Eigenvalue Decomposition # ---------------------------------- @@ -105,7 +105,7 @@ const LAPACK_EighAlgorithm = Union{ } # TODO: -@algdef BigFloat_eigh_Francis +@algdef GLA_eigh_Francis # Singular Value Decomposition # ---------------------------- @@ -125,7 +125,7 @@ const LAPACK_SVDAlgorithm = Union{ } # TODO: -@algdef BigFloat_svd_QRIteration +@algdef GLA_svd_QRIteration # ========================= # Polar decompositions diff --git a/test/bigfloats/eig.jl b/test/bigfloats/eig.jl index ac244b4a..5f347c78 100644 --- a/test/bigfloats/eig.jl +++ b/test/bigfloats/eig.jl @@ -12,7 +12,7 @@ const eltypes = (BigFloat, Complex{BigFloat}) @testset "eig_full! for T = $T" for T in eltypes rng = StableRNG(123) m = 24 - alg = BigFloat_eig_Francis() + alg = GLA_eig_Francis() A = randn(rng, T, m, m) Tc = complex(T) @@ -36,7 +36,7 @@ end @testset "eig_trunc! for T = $T" for T in eltypes rng = StableRNG(123) m = 6 - alg = BigFloat_eig_Francis() + alg = GLA_eig_Francis() A = randn(rng, T, m, m) A *= A' # TODO: deal with eigenvalue ordering etc # eigenvalues are sorted by ascending real component... @@ -81,13 +81,13 @@ end V = randn(rng, T, m, m) D = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) A = V * D * inv(V) - alg = TruncatedAlgorithm(BigFloat_eig_Francis(), truncrank(2)) + alg = TruncatedAlgorithm(GLA_eig_Francis(), truncrank(2)) D2, V2, ϵ2 = @constinferred eig_trunc(A; alg) @test diagview(D2) ≈ diagview(D)[1:2] @test ϵ2 ≈ norm(diagview(D)[3:4]) atol = atol @test_throws ArgumentError eig_trunc(A; alg, trunc = (; maxrank = 2)) - alg = TruncatedAlgorithm(BigFloat_eig_Francis(), truncerror(; atol = 0.2, p = 1)) + alg = TruncatedAlgorithm(GLA_eig_Francis(), truncerror(; atol = 0.2, p = 1)) D3, V3, ϵ3 = @constinferred eig_trunc(A; alg) @test diagview(D3) ≈ diagview(D)[1:2] @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol diff --git a/test/bigfloats/eigh.jl b/test/bigfloats/eigh.jl index 8042df20..cca4677a 100644 --- a/test/bigfloats/eigh.jl +++ b/test/bigfloats/eigh.jl @@ -12,7 +12,7 @@ const eltypes = (BigFloat, Complex{BigFloat}) @testset "eigh_full! for T = $T" for T in eltypes rng = StableRNG(123) m = 54 - alg = BigFloat_eigh_Francis() + alg = GLA_eigh_Francis() A = randn(rng, T, m, m) A = (A + A') / 2 @@ -33,7 +33,7 @@ end @testset "eigh_trunc! for T = $T" for T in eltypes rng = StableRNG(123) m = 54 - alg = BigFloat_eigh_Francis() + alg = GLA_eigh_Francis() A = randn(rng, T, m, m) A = A * A' A = (A + A') / 2 @@ -81,13 +81,13 @@ end D = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) A = V * D * V' A = (A + A') / 2 - alg = TruncatedAlgorithm(BigFloat_eigh_Francis(), truncrank(2)) + alg = TruncatedAlgorithm(GLA_eigh_Francis(), truncrank(2)) D2, V2, ϵ2 = @constinferred eigh_trunc(A; alg) @test diagview(D2) ≈ diagview(D)[1:2] @test_throws ArgumentError eigh_trunc(A; alg, trunc = (; maxrank = 2)) @test ϵ2 ≈ norm(diagview(D)[3:4]) atol = atol - alg = TruncatedAlgorithm(BigFloat_eigh_Francis(), truncerror(; atol = 0.2)) + alg = TruncatedAlgorithm(GLA_eigh_Francis(), truncerror(; atol = 0.2)) D3, V3, ϵ3 = @constinferred eigh_trunc(A; alg) @test diagview(D3) ≈ diagview(D)[1:2] @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol @@ -111,7 +111,7 @@ end @test diagview(D) ≈ D2 A2 = Diagonal(T[0.9, 0.3, 0.1, 0.01]) - alg = TruncatedAlgorithm(BigFloat_eigh_Francis(), truncrank(2)) + alg = TruncatedAlgorithm(GLA_eigh_Francis(), truncrank(2)) D2, V2, ϵ2 = @constinferred eigh_trunc(A2; alg) @test diagview(D2) ≈ diagview(A2)[1:2] @test ϵ2 ≈ norm(diagview(A2)[3:4]) atol = atol diff --git a/test/bigfloats/lq.jl b/test/bigfloats/lq.jl index c0e4faba..9d904da9 100644 --- a/test/bigfloats/lq.jl +++ b/test/bigfloats/lq.jl @@ -33,7 +33,7 @@ eltypes = (BigFloat, Complex{BigFloat}) @test Q == Q2 # Transposed QR algorithm - qr_alg = BigFloat_QR_Householder() + qr_alg = GLA_QR_Householder() lq_alg = LQViaTransposedQR(qr_alg) L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q), lq_alg) @test L2 === L @@ -81,7 +81,7 @@ end @test Q[1:minmn, n] ≈ Q2[1:minmn, n] # Transposed QR algorithm - qr_alg = BigFloat_QR_Householder() + qr_alg = GLA_QR_Householder() lq_alg = LQViaTransposedQR(qr_alg) L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q), lq_alg) @test L2 === L @@ -105,7 +105,7 @@ end lq_full!(copy!(Ac, A), (noL, Q2); positive = true) @test Q[1:minmn, n] ≈ Q2[1:minmn, n] - qr_alg = BigFloat_QR_Householder(; positive = true) + qr_alg = GLA_QR_Householder(; positive = true) lq_alg = LQViaTransposedQR(qr_alg) lq_full!(copy!(Ac, A), (L, Q), lq_alg) @test L * Q ≈ A diff --git a/test/bigfloats/svd.jl b/test/bigfloats/svd.jl index 5e616d37..f52f7e45 100644 --- a/test/bigfloats/svd.jl +++ b/test/bigfloats/svd.jl @@ -14,7 +14,7 @@ eltypes = (BigFloat, Complex{BigFloat}) m = 54 @testset "size ($m, $n)" for n in (37, m, 63, 0) k = min(m, n) - alg = BigFloat_svd_QRIteration() + alg = GLA_svd_QRIteration() minmn = min(m, n) A = randn(rng, T, m, n) @@ -53,7 +53,7 @@ end rng = StableRNG(123) m = 54 @testset "size ($m, $n)" for n in (37, m, 63, 0) - alg = BigFloat_svd_QRIteration() + alg = GLA_svd_QRIteration() A = randn(rng, T, m, n) U, S, Vᴴ = svd_full(A; alg) @test U isa Matrix{T} && size(U) == (m, m) @@ -99,7 +99,7 @@ end rng = StableRNG(123) m = 54 atol = sqrt(eps(real(T))) - alg = BigFloat_svd_QRIteration() + alg = GLA_svd_QRIteration() @testset "size ($m, $n)" for n in (37, m, 63) n > m && alg isa LAPACK_Jacobi && continue # not supported @@ -137,7 +137,7 @@ end @testset "svd_trunc! mix maxrank and tol for T = $T" for T in eltypes rng = StableRNG(123) - alg = BigFloat_svd_QRIteration() + alg = GLA_svd_QRIteration() m = 4 U = qr_compact(randn(rng, T, m, m))[1] S = Diagonal(T[0.9, 0.3, 0.1, 0.01]) @@ -166,7 +166,7 @@ end S = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) Vᴴ = qr_compact(randn(rng, T, m, m))[1] A = U * S * Vᴴ - alg = TruncatedAlgorithm(BigFloat_svd_QRIteration(), trunctol(; atol = 0.2)) + alg = TruncatedAlgorithm(GLA_svd_QRIteration(), trunctol(; atol = 0.2)) U2, S2, V2ᴴ, ϵ2 = @constinferred svd_trunc(A; alg) @test diagview(S2) ≈ diagview(S)[1:2] @test ϵ2 ≈ norm(diagview(S)[3:4]) atol = atol From 03caa562f2480b5971e5d081734d978c2219aee7 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Thu, 30 Oct 2025 17:59:53 +0100 Subject: [PATCH 08/25] Split extensions --- Project.toml | 3 ++- ...atrixAlgebraKitGenericLinearAlgebraExt.jl} | 23 ++----------------- ext/MatrixAlgebraKitGenericSchurExt.jl | 23 +++++++++++++++++++ test/bigfloats/eig.jl | 1 - test/bigfloats/eigh.jl | 1 - test/bigfloats/lq.jl | 1 - test/bigfloats/qr.jl | 1 - test/bigfloats/svd.jl | 1 - 8 files changed, 27 insertions(+), 27 deletions(-) rename ext/{MatrixAlgebraKitGenericExt.jl => MatrixAlgebraKitGenericLinearAlgebraExt.jl} (87%) create mode 100644 ext/MatrixAlgebraKitGenericSchurExt.jl diff --git a/Project.toml b/Project.toml index aa1f3cf1..934c0ceb 100644 --- a/Project.toml +++ b/Project.toml @@ -17,7 +17,8 @@ GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e" MatrixAlgebraKitChainRulesCoreExt = "ChainRulesCore" MatrixAlgebraKitAMDGPUExt = "AMDGPU" MatrixAlgebraKitCUDAExt = "CUDA" -MatrixAlgebraKitGenericExt = ["GenericLinearAlgebra", "GenericSchur"] +MatrixAlgebraKitGenericLinearAlgebraExt = "GenericLinearAlgebra" +MatrixAlgebraKitGenericSchurExt = "GenericSchur" [compat] AMDGPU = "2" diff --git a/ext/MatrixAlgebraKitGenericExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl similarity index 87% rename from ext/MatrixAlgebraKitGenericExt.jl rename to ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 56c8af3a..9c451150 100644 --- a/ext/MatrixAlgebraKitGenericExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -1,12 +1,8 @@ -module MatrixAlgebraKitGenericExt +module MatrixAlgebraKitGenericLinearAlgebraExt using MatrixAlgebraKit -using MatrixAlgebraKit: LAPACK_SVDAlgorithm, LAPACK_EigAlgorithm, LAPACK_EighAlgorithm, LAPACK_QRIteration -using MatrixAlgebraKit: uppertriangular! -using MatrixAlgebraKit: @algdef, Algorithm, check_input -using MatrixAlgebraKit: sign_safe +using MatrixAlgebraKit: sign_safe, check_input using GenericLinearAlgebra -using GenericSchur using LinearAlgebra: I, Diagonal function MatrixAlgebraKit.default_svd_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} @@ -52,21 +48,6 @@ function MatrixAlgebraKit.default_eig_algorithm(::Type{T}; kwargs...) where {T < return GLA_eig_Francis(; kwargs...) end -function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GLA_eig_Francis) where {T} - check_input(eig_full!, A, DV, alg) - D, V = DV - D̃, Ṽ = GenericSchur.eigen!(A) - copyto!(D, Diagonal(D̃)) - copyto!(V, Ṽ) - return D, V -end - -function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::GLA_eig_Francis) where {T} - check_input(eig_vals!, A, D, alg) - eigval = GenericSchur.eigvals!(A) - return eigval -end - function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} return GLA_eigh_Francis(; kwargs...) end diff --git a/ext/MatrixAlgebraKitGenericSchurExt.jl b/ext/MatrixAlgebraKitGenericSchurExt.jl new file mode 100644 index 00000000..36645ec5 --- /dev/null +++ b/ext/MatrixAlgebraKitGenericSchurExt.jl @@ -0,0 +1,23 @@ +module MatrixAlgebraKitGenericSchurExt + +using MatrixAlgebraKit +using MatrixAlgebraKit: check_input +using LinearAlgebra: Diagonal +using GenericSchur + +function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GLA_eig_Francis) where {T} + check_input(eig_full!, A, DV, alg) + D, V = DV + D̃, Ṽ = GenericSchur.eigen!(A) + copyto!(D, Diagonal(D̃)) + copyto!(V, Ṽ) + return D, V +end + +function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::GLA_eig_Francis) where {T} + check_input(eig_vals!, A, D, alg) + eigval = GenericSchur.eigvals!(A) + return eigval +end + +end diff --git a/test/bigfloats/eig.jl b/test/bigfloats/eig.jl index 5f347c78..37150c8f 100644 --- a/test/bigfloats/eig.jl +++ b/test/bigfloats/eig.jl @@ -4,7 +4,6 @@ using TestExtras using StableRNGs using LinearAlgebra: Diagonal using MatrixAlgebraKit: TruncatedAlgorithm, diagview, norm -using GenericLinearAlgebra using GenericSchur const eltypes = (BigFloat, Complex{BigFloat}) diff --git a/test/bigfloats/eigh.jl b/test/bigfloats/eigh.jl index cca4677a..b9bc9e16 100644 --- a/test/bigfloats/eigh.jl +++ b/test/bigfloats/eigh.jl @@ -5,7 +5,6 @@ using StableRNGs using LinearAlgebra: LinearAlgebra, Diagonal, I using MatrixAlgebraKit: TruncatedAlgorithm, diagview, norm using GenericLinearAlgebra -using GenericSchur const eltypes = (BigFloat, Complex{BigFloat}) diff --git a/test/bigfloats/lq.jl b/test/bigfloats/lq.jl index 9d904da9..d071afee 100644 --- a/test/bigfloats/lq.jl +++ b/test/bigfloats/lq.jl @@ -4,7 +4,6 @@ using TestExtras using StableRNGs using LinearAlgebra: diag, I, Diagonal using GenericLinearAlgebra -using GenericSchur eltypes = (BigFloat, Complex{BigFloat}) diff --git a/test/bigfloats/qr.jl b/test/bigfloats/qr.jl index 32c8e375..9e1ba515 100644 --- a/test/bigfloats/qr.jl +++ b/test/bigfloats/qr.jl @@ -4,7 +4,6 @@ using TestExtras using StableRNGs using LinearAlgebra: diag, I, Diagonal using GenericLinearAlgebra -using GenericSchur eltypes = (BigFloat, Complex{BigFloat}) diff --git a/test/bigfloats/svd.jl b/test/bigfloats/svd.jl index f52f7e45..abe369e0 100644 --- a/test/bigfloats/svd.jl +++ b/test/bigfloats/svd.jl @@ -5,7 +5,6 @@ using StableRNGs using LinearAlgebra: LinearAlgebra, Diagonal, I, isposdef, norm using MatrixAlgebraKit: TruncatedAlgorithm, diagview, isisometric using GenericLinearAlgebra -using GenericSchur eltypes = (BigFloat, Complex{BigFloat}) From faa0921465f8d9290d61023cd6b8076c42845a76 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 31 Oct 2025 10:05:47 +0100 Subject: [PATCH 09/25] fix copies --- ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl | 13 ++++++++++--- ext/MatrixAlgebraKitGenericSchurExt.jl | 3 ++- test/bigfloats/eig.jl | 4 ++-- test/bigfloats/eigh.jl | 4 ++-- test/bigfloats/svd.jl | 6 +++--- 5 files changed, 19 insertions(+), 11 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 9c451150..8457c51a 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -11,8 +11,12 @@ end function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_QRIteration) where {T} check_input(svd_compact!, A, USVᴴ, alg) - U, S, V = GenericLinearAlgebra.svd(A) - return U, Diagonal(S), V' # conjugation to account for difference in convention + U, S, Vᴴ = USVᴴ + Ũ, S̃, Ṽ = GenericLinearAlgebra.svd(A) + copyto!(U, Ũ) + copyto!(S, Diagonal(S̃)) + copyto!(Vᴴ, Ṽ') # conjugation to account for difference in convention + return U, S, Vᴴ end function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_QRIteration) where {T} @@ -54,8 +58,11 @@ end function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_eigh_Francis)::Tuple{Diagonal{real(T)}, Matrix{T}} where {T} check_input(eigh_full!, A, DV, alg) + D, V = DV eigval, eigvec = GenericLinearAlgebra.eigen(A; sortby = real) - return Diagonal(eigval), eigvec + copyto!(D, Diagonal(eigval)) + copyto!(V, eigvec) + return D, V end function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::GLA_eigh_Francis) where {T} diff --git a/ext/MatrixAlgebraKitGenericSchurExt.jl b/ext/MatrixAlgebraKitGenericSchurExt.jl index 36645ec5..b576b6c2 100644 --- a/ext/MatrixAlgebraKitGenericSchurExt.jl +++ b/ext/MatrixAlgebraKitGenericSchurExt.jl @@ -17,7 +17,8 @@ end function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::GLA_eig_Francis) where {T} check_input(eig_vals!, A, D, alg) eigval = GenericSchur.eigvals!(A) - return eigval + copyto!(D, eigval) + return D end end diff --git a/test/bigfloats/eig.jl b/test/bigfloats/eig.jl index 37150c8f..2e166231 100644 --- a/test/bigfloats/eig.jl +++ b/test/bigfloats/eig.jl @@ -23,8 +23,8 @@ const eltypes = (BigFloat, Complex{BigFloat}) Ac = similar(A) D2, V2 = @constinferred eig_full!(copy!(Ac, A), (D, V), alg′) - @test D2 ≈ D - @test V2 ≈ V + @test D2 === D + @test V2 === V @test A * V ≈ V * D Dc = @constinferred eig_vals(A, alg′) diff --git a/test/bigfloats/eigh.jl b/test/bigfloats/eigh.jl index b9bc9e16..9fc99408 100644 --- a/test/bigfloats/eigh.jl +++ b/test/bigfloats/eigh.jl @@ -22,8 +22,8 @@ const eltypes = (BigFloat, Complex{BigFloat}) @test all(isreal, D) D2, V2 = eigh_full!(copy(A), (D, V), alg) - @test D2 ≈ D - @test V2 ≈ V + @test D2 === D + @test V2 === V D3 = @constinferred eigh_vals(A, alg) @test D ≈ Diagonal(D3) diff --git a/test/bigfloats/svd.jl b/test/bigfloats/svd.jl index abe369e0..5f7884bd 100644 --- a/test/bigfloats/svd.jl +++ b/test/bigfloats/svd.jl @@ -35,9 +35,9 @@ eltypes = (BigFloat, Complex{BigFloat}) Sc = similar(A, real(T), min(m, n)) alg′ = @constinferred MatrixAlgebraKit.select_algorithm(svd_compact!, A, $alg) U2, S2, V2ᴴ = @constinferred svd_compact!(copy!(Ac, A), (U, S, Vᴴ), alg′) - @test U2 ≈ U - @test S2 ≈ S - @test V2ᴴ ≈ Vᴴ + @test U2 === U + @test S2 === S + @test V2ᴴ === Vᴴ @test U * S * Vᴴ ≈ A @test isisometric(U) @test isisometric(Vᴴ; side = :right) From 9ed02bd9ed25b94a682071428a115ccf140c4e00 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 31 Oct 2025 10:29:09 +0100 Subject: [PATCH 10/25] fix type instability use Hermitian() small cleanup --- ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 8457c51a..25e849a3 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -2,7 +2,7 @@ module MatrixAlgebraKitGenericLinearAlgebraExt using MatrixAlgebraKit using MatrixAlgebraKit: sign_safe, check_input -using GenericLinearAlgebra +using GenericLinearAlgebra: svd, svdvals!, eigen, eigvals, Hermitian, qr using LinearAlgebra: I, Diagonal function MatrixAlgebraKit.default_svd_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} @@ -12,7 +12,7 @@ end function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_QRIteration) where {T} check_input(svd_compact!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ - Ũ, S̃, Ṽ = GenericLinearAlgebra.svd(A) + Ũ, S̃, Ṽ = svd(A) copyto!(U, Ũ) copyto!(S, Diagonal(S̃)) copyto!(Vᴴ, Ṽ') # conjugation to account for difference in convention @@ -31,7 +31,7 @@ function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_Q return USVᴴ end S̃ = fill!(S, zero(T)) - U_compact, S_compact, V_compact = GenericLinearAlgebra.svd(A) + U_compact, S_compact, V_compact = svd(A) S̃[1:minmn, 1:minmn] .= Diagonal(S_compact) copyto!(S, S̃) @@ -43,7 +43,7 @@ end function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::GLA_svd_QRIteration) where {T} check_input(svd_vals!, A, S, alg) - S̃ = GenericLinearAlgebra.svdvals!(A) + S̃ = svdvals!(A) copyto!(S, S̃) return S end @@ -56,10 +56,10 @@ function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T return GLA_eigh_Francis(; kwargs...) end -function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_eigh_Francis)::Tuple{Diagonal{real(T)}, Matrix{T}} where {T} +function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_eigh_Francis) where {T} check_input(eigh_full!, A, DV, alg) D, V = DV - eigval, eigvec = GenericLinearAlgebra.eigen(A; sortby = real) + eigval, eigvec = eigen(Hermitian(A); sortby = real) copyto!(D, Diagonal(eigval)) copyto!(V, eigvec) return D, V @@ -67,7 +67,7 @@ end function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::GLA_eigh_Francis) where {T} check_input(eigh_vals!, A, D, alg) - D = GenericLinearAlgebra.eigvals(A; sortby = real) + D = eigvals(A; sortby = real) return real.(D) end @@ -106,7 +106,7 @@ function _gla_householder_qr!(A::AbstractMatrix{T}, Q, R; positive = false, bloc m, n = size(A) k = min(m, n) computeR = length(R) > 0 - Q̃, R̃ = GenericLinearAlgebra.qr(A) + Q̃, R̃ = qr(A) Q̃ = convert(Array, Q̃) if positive @inbounds for j in 1:k From c3e069ec1c4cde6d2c74ea92b7eb54fe2e7cc71b Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 31 Oct 2025 10:44:29 +0100 Subject: [PATCH 11/25] Name change From GLA to GS, since eig using GenericSchur and not GenericLinearAlgebra --- ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl | 4 ---- ext/MatrixAlgebraKitGenericSchurExt.jl | 8 ++++++-- src/MatrixAlgebraKit.jl | 2 +- src/interface/decompositions.jl | 2 +- test/bigfloats/eig.jl | 8 ++++---- 5 files changed, 12 insertions(+), 12 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 25e849a3..12fcf2ed 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -48,10 +48,6 @@ function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::GLA_svd_QRIter return S end -function MatrixAlgebraKit.default_eig_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return GLA_eig_Francis(; kwargs...) -end - function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} return GLA_eigh_Francis(; kwargs...) end diff --git a/ext/MatrixAlgebraKitGenericSchurExt.jl b/ext/MatrixAlgebraKitGenericSchurExt.jl index b576b6c2..fd06fc1b 100644 --- a/ext/MatrixAlgebraKitGenericSchurExt.jl +++ b/ext/MatrixAlgebraKitGenericSchurExt.jl @@ -5,7 +5,11 @@ using MatrixAlgebraKit: check_input using LinearAlgebra: Diagonal using GenericSchur -function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GLA_eig_Francis) where {T} +function MatrixAlgebraKit.default_eig_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} + return GS_eig_Francis(; kwargs...) +end + +function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GS_eig_Francis) where {T} check_input(eig_full!, A, DV, alg) D, V = DV D̃, Ṽ = GenericSchur.eigen!(A) @@ -14,7 +18,7 @@ function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GLA_eig_Franc return D, V end -function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::GLA_eig_Francis) where {T} +function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::GS_eig_Francis) where {T} check_input(eig_vals!, A, D, alg) eigval = GenericSchur.eigvals!(A) copyto!(D, eigval) diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index c73f2269..71403167 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -33,7 +33,7 @@ export left_orth!, right_orth!, left_null!, right_null! export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert, LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations, LAPACK_DivideAndConquer, LAPACK_Jacobi -export GLA_QR_Householder, GLA_eig_Francis, GLA_eigh_Francis, GLA_svd_QRIteration +export GLA_QR_Householder, GS_eig_Francis, GLA_eigh_Francis, GLA_svd_QRIteration export LQViaTransposedQR export PolarViaSVD, PolarNewton export DiagonalAlgorithm diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index d138e701..944a210d 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -58,7 +58,7 @@ eigenvalue decomposition of a matrix. const LAPACK_EigAlgorithm = Union{LAPACK_Simple, LAPACK_Expert} # TODO: -@algdef GLA_eig_Francis +@algdef GS_eig_Francis # Hermitian Eigenvalue Decomposition # ---------------------------------- diff --git a/test/bigfloats/eig.jl b/test/bigfloats/eig.jl index 2e166231..69a04761 100644 --- a/test/bigfloats/eig.jl +++ b/test/bigfloats/eig.jl @@ -11,7 +11,7 @@ const eltypes = (BigFloat, Complex{BigFloat}) @testset "eig_full! for T = $T" for T in eltypes rng = StableRNG(123) m = 24 - alg = GLA_eig_Francis() + alg = GS_eig_Francis() A = randn(rng, T, m, m) Tc = complex(T) @@ -35,7 +35,7 @@ end @testset "eig_trunc! for T = $T" for T in eltypes rng = StableRNG(123) m = 6 - alg = GLA_eig_Francis() + alg = GS_eig_Francis() A = randn(rng, T, m, m) A *= A' # TODO: deal with eigenvalue ordering etc # eigenvalues are sorted by ascending real component... @@ -80,13 +80,13 @@ end V = randn(rng, T, m, m) D = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) A = V * D * inv(V) - alg = TruncatedAlgorithm(GLA_eig_Francis(), truncrank(2)) + alg = TruncatedAlgorithm(GS_eig_Francis(), truncrank(2)) D2, V2, ϵ2 = @constinferred eig_trunc(A; alg) @test diagview(D2) ≈ diagview(D)[1:2] @test ϵ2 ≈ norm(diagview(D)[3:4]) atol = atol @test_throws ArgumentError eig_trunc(A; alg, trunc = (; maxrank = 2)) - alg = TruncatedAlgorithm(GLA_eig_Francis(), truncerror(; atol = 0.2, p = 1)) + alg = TruncatedAlgorithm(GS_eig_Francis(), truncerror(; atol = 0.2, p = 1)) D3, V3, ϵ3 = @constinferred eig_trunc(A; alg) @test diagview(D3) ≈ diagview(D)[1:2] @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol From 7af4d52c93f247b1a3e5861d04f581f244c55b22 Mon Sep 17 00:00:00 2001 From: Sander De Meyer <74001142+sanderdemeyer@users.noreply.github.com> Date: Tue, 4 Nov 2025 11:25:27 +0100 Subject: [PATCH 12/25] Update ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl Co-authored-by: Jutho --- ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 12fcf2ed..07b35bac 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -102,7 +102,7 @@ function _gla_householder_qr!(A::AbstractMatrix{T}, Q, R; positive = false, bloc m, n = size(A) k = min(m, n) computeR = length(R) > 0 - Q̃, R̃ = qr(A) + Q̃, R̃ = qr!(A) Q̃ = convert(Array, Q̃) if positive @inbounds for j in 1:k From c6ba4e43a71178a9cbd4a0554dcbe9be7cb8e423 Mon Sep 17 00:00:00 2001 From: Sander De Meyer <74001142+sanderdemeyer@users.noreply.github.com> Date: Tue, 4 Nov 2025 11:25:37 +0100 Subject: [PATCH 13/25] Update ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl Co-authored-by: Jutho --- ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 07b35bac..0634ed4d 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -12,7 +12,7 @@ end function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_QRIteration) where {T} check_input(svd_compact!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ - Ũ, S̃, Ṽ = svd(A) + Ũ, S̃, Ṽ = svd!(A) copyto!(U, Ũ) copyto!(S, Diagonal(S̃)) copyto!(Vᴴ, Ṽ') # conjugation to account for difference in convention From fdb1222b8f16323d064b885ec357a835ebfbc48d Mon Sep 17 00:00:00 2001 From: Sander De Meyer <74001142+sanderdemeyer@users.noreply.github.com> Date: Tue, 4 Nov 2025 11:25:46 +0100 Subject: [PATCH 14/25] Update ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl Co-authored-by: Jutho --- ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 0634ed4d..aa4a8af0 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -55,7 +55,7 @@ end function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_eigh_Francis) where {T} check_input(eigh_full!, A, DV, alg) D, V = DV - eigval, eigvec = eigen(Hermitian(A); sortby = real) + eigval, eigvec = eigen!(Hermitian(A); sortby = real) copyto!(D, Diagonal(eigval)) copyto!(V, eigvec) return D, V From 837a7a41b5d56a511d6b52b428b53ae88decba30 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Tue, 4 Nov 2025 15:23:41 +0100 Subject: [PATCH 15/25] Resolve comments Name change some copying issues resolved --- ...MatrixAlgebraKitGenericLinearAlgebraExt.jl | 73 ++++++++----------- ext/MatrixAlgebraKitGenericSchurExt.jl | 6 +- src/MatrixAlgebraKit.jl | 2 +- src/interface/decompositions.jl | 7 +- test/bigfloats/eig.jl | 8 +- test/bigfloats/eigh.jl | 10 +-- test/bigfloats/svd.jl | 10 +-- 7 files changed, 52 insertions(+), 64 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index aa4a8af0..198e5ec7 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -2,24 +2,24 @@ module MatrixAlgebraKitGenericLinearAlgebraExt using MatrixAlgebraKit using MatrixAlgebraKit: sign_safe, check_input -using GenericLinearAlgebra: svd, svdvals!, eigen, eigvals, Hermitian, qr -using LinearAlgebra: I, Diagonal +using GenericLinearAlgebra: svd!, svdvals!, eigen!, eigvals!, Hermitian, qr! +using LinearAlgebra: I, Diagonal, rmul!, lmul!, transpose!, dot function MatrixAlgebraKit.default_svd_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return GLA_svd_QRIteration() + return GLA_QRIteration() end -function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_QRIteration) where {T} +function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix, USVᴴ, alg::GLA_QRIteration) check_input(svd_compact!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ - Ũ, S̃, Ṽ = svd!(A) - copyto!(U, Ũ) - copyto!(S, Diagonal(S̃)) - copyto!(Vᴴ, Ṽ') # conjugation to account for difference in convention + F = svd!(A) + copyto!(U, F.U) + copyto!(S, Diagonal(F.S)) + copyto!(Vᴴ, F.Vt) # conjugation to account for difference in convention return U, S, Vᴴ end -function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_QRIteration) where {T} +function MatrixAlgebraKit.svd_full!(A::AbstractMatrix, USVᴴ, alg::GLA_QRIteration) check_input(svd_full!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ m, n = size(A) @@ -30,18 +30,17 @@ function MatrixAlgebraKit.svd_full!(A::AbstractMatrix{T}, USVᴴ, alg::GLA_svd_Q MatrixAlgebraKit.one!(Vᴴ) return USVᴴ end - S̃ = fill!(S, zero(T)) - U_compact, S_compact, V_compact = svd(A) - S̃[1:minmn, 1:minmn] .= Diagonal(S_compact) - copyto!(S, S̃) + S = MatrixAlgebraKit.zero!(S) + U_compact, S_compact, Vᴴ_compact = svd!(A) + S[1:minmn, 1:minmn] .= Diagonal(S_compact) U = _gram_schmidt!(U, U_compact) - Vᴴ = _gram_schmidt!(Vᴴ, V_compact; adjoint = true) + Vᴴ = _gram_schmidt!(Vᴴ, Vᴴ_compact; adjoint = true) return MatrixAlgebraKit.gaugefix!(svd_full!, U, S, Vᴴ, m, n) end -function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::GLA_svd_QRIteration) where {T} +function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::GLA_QRIteration) where {T} check_input(svd_vals!, A, S, alg) S̃ = svdvals!(A) copyto!(S, S̃) @@ -49,10 +48,10 @@ function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::GLA_svd_QRIter end function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return GLA_eigh_Francis(; kwargs...) + return GLA_QRIteration(; kwargs...) end -function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_eigh_Francis) where {T} +function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_QRIteration) where {T} check_input(eigh_full!, A, DV, alg) D, V = DV eigval, eigvec = eigen!(Hermitian(A); sortby = real) @@ -61,10 +60,9 @@ function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_eigh_Fra return D, V end -function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::GLA_eigh_Francis) where {T} +function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::GLA_QRIteration) where {T} check_input(eigh_vals!, A, D, alg) - D = eigvals(A; sortby = real) - return real.(D) + return eigvals!(Hermitian(A); sortby = real) end function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} @@ -80,10 +78,10 @@ function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::GLA_QR_Householde Q_zero = zeros(eltype(Q), (m, minmn)) R_zero = zeros(eltype(R), (minmn, n)) Q_compact, R_compact = _gla_householder_qr!(A, Q_zero, R_zero; alg.kwargs...) - Q = _gram_schmidt!(Q, Q_compact[:, 1:min(m, n)]) + Q = _gram_schmidt!(Q, view(Q_compact, :, 1:minmn)) if computeR - R = fill!(R, zero(eltype(R))) - R[1:minmn, 1:n] .= R_compact + R[1:minmn, :] .= R_compact + R[(minmn + 1):end, :] .= zero(eltype(R)) end return Q, R end @@ -103,16 +101,20 @@ function _gla_householder_qr!(A::AbstractMatrix{T}, Q, R; positive = false, bloc k = min(m, n) computeR = length(R) > 0 Q̃, R̃ = qr!(A) - Q̃ = convert(Array, Q̃) + + if k < m + copyto!(Q, Q̃) + else + rmul!(MatrixAlgebraKit.one!(Q), Q̃) + end if positive @inbounds for j in 1:k s = sign_safe(R̃[j, j]) @simd for i in 1:m - Q̃[i, j] *= s + Q[i, j] *= s end end end - copyto!(Q, Q̃) if computeR if positive @inbounds for j in n:-1:1 @@ -126,30 +128,19 @@ function _gla_householder_qr!(A::AbstractMatrix{T}, Q, R; positive = false, bloc return Q, R end -function _gram_schmidt(Q_compact) +function _gram_schmidt!(Q, Q_compact; adjoint = false) m, minmn = size(Q_compact) - if minmn >= m - return Q_compact - end - Q = zeros(eltype(Q_compact), (m, m)) Q[:, 1:minmn] .= Q_compact for j in (minmn + 1):m v = rand(eltype(Q), m) for i in 1:(j - 1) - r = sum([v[k] * conj(Q[k, i])] for k in 1:size(v)[1])[1] - v .= v .- r * Q[:, i] + r = dot(view(Q, :, i), v) + v .-= r * view(Q, :, i) end Q[:, j] = v ./ MatrixAlgebraKit.norm(v) end - return Q -end - -function _gram_schmidt!(Q, Q_compact; adjoint = false) - Q̃ = _gram_schmidt(Q_compact) if adjoint - copyto!(Q, Q̃') - else - copyto!(Q, Q̃) + copyto!(Q, Q') end return Q end diff --git a/ext/MatrixAlgebraKitGenericSchurExt.jl b/ext/MatrixAlgebraKitGenericSchurExt.jl index fd06fc1b..04d45769 100644 --- a/ext/MatrixAlgebraKitGenericSchurExt.jl +++ b/ext/MatrixAlgebraKitGenericSchurExt.jl @@ -6,10 +6,10 @@ using LinearAlgebra: Diagonal using GenericSchur function MatrixAlgebraKit.default_eig_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return GS_eig_Francis(; kwargs...) + return GS_QRIteration(; kwargs...) end -function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GS_eig_Francis) where {T} +function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GS_QRIteration) where {T} check_input(eig_full!, A, DV, alg) D, V = DV D̃, Ṽ = GenericSchur.eigen!(A) @@ -18,7 +18,7 @@ function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GS_eig_Franci return D, V end -function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::GS_eig_Francis) where {T} +function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::GS_QRIteration) where {T} check_input(eig_vals!, A, D, alg) eigval = GenericSchur.eigvals!(A) copyto!(D, eigval) diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index 71403167..1fd4fabe 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -33,7 +33,7 @@ export left_orth!, right_orth!, left_null!, right_null! export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert, LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations, LAPACK_DivideAndConquer, LAPACK_Jacobi -export GLA_QR_Householder, GS_eig_Francis, GLA_eigh_Francis, GLA_svd_QRIteration +export GLA_QR_Householder, GLA_QRIteration, GS_QRIteration export LQViaTransposedQR export PolarViaSVD, PolarNewton export DiagonalAlgorithm diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index ddaa13b4..b23fbeab 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -58,7 +58,7 @@ eigenvalue decomposition of a matrix. const LAPACK_EigAlgorithm = Union{LAPACK_Simple, LAPACK_Expert} # TODO: -@algdef GS_eig_Francis +@algdef GS_QRIteration # Hermitian Eigenvalue Decomposition # ---------------------------------- @@ -105,7 +105,7 @@ const LAPACK_EighAlgorithm = Union{ } # TODO: -@algdef GLA_eigh_Francis +@algdef GLA_QRIteration # Singular Value Decomposition # ---------------------------- @@ -124,9 +124,6 @@ const LAPACK_SVDAlgorithm = Union{ LAPACK_Jacobi, } -# TODO: -@algdef GLA_svd_QRIteration - # ========================= # Polar decompositions # ========================= diff --git a/test/bigfloats/eig.jl b/test/bigfloats/eig.jl index 69a04761..88de1b60 100644 --- a/test/bigfloats/eig.jl +++ b/test/bigfloats/eig.jl @@ -11,7 +11,7 @@ const eltypes = (BigFloat, Complex{BigFloat}) @testset "eig_full! for T = $T" for T in eltypes rng = StableRNG(123) m = 24 - alg = GS_eig_Francis() + alg = GS_QRIteration() A = randn(rng, T, m, m) Tc = complex(T) @@ -35,7 +35,7 @@ end @testset "eig_trunc! for T = $T" for T in eltypes rng = StableRNG(123) m = 6 - alg = GS_eig_Francis() + alg = GS_QRIteration() A = randn(rng, T, m, m) A *= A' # TODO: deal with eigenvalue ordering etc # eigenvalues are sorted by ascending real component... @@ -80,13 +80,13 @@ end V = randn(rng, T, m, m) D = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) A = V * D * inv(V) - alg = TruncatedAlgorithm(GS_eig_Francis(), truncrank(2)) + alg = TruncatedAlgorithm(GS_QRIteration(), truncrank(2)) D2, V2, ϵ2 = @constinferred eig_trunc(A; alg) @test diagview(D2) ≈ diagview(D)[1:2] @test ϵ2 ≈ norm(diagview(D)[3:4]) atol = atol @test_throws ArgumentError eig_trunc(A; alg, trunc = (; maxrank = 2)) - alg = TruncatedAlgorithm(GS_eig_Francis(), truncerror(; atol = 0.2, p = 1)) + alg = TruncatedAlgorithm(GS_QRIteration(), truncerror(; atol = 0.2, p = 1)) D3, V3, ϵ3 = @constinferred eig_trunc(A; alg) @test diagview(D3) ≈ diagview(D)[1:2] @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol diff --git a/test/bigfloats/eigh.jl b/test/bigfloats/eigh.jl index 9fc99408..a176981d 100644 --- a/test/bigfloats/eigh.jl +++ b/test/bigfloats/eigh.jl @@ -11,7 +11,7 @@ const eltypes = (BigFloat, Complex{BigFloat}) @testset "eigh_full! for T = $T" for T in eltypes rng = StableRNG(123) m = 54 - alg = GLA_eigh_Francis() + alg = GLA_QRIteration() A = randn(rng, T, m, m) A = (A + A') / 2 @@ -32,7 +32,7 @@ end @testset "eigh_trunc! for T = $T" for T in eltypes rng = StableRNG(123) m = 54 - alg = GLA_eigh_Francis() + alg = GLA_QRIteration() A = randn(rng, T, m, m) A = A * A' A = (A + A') / 2 @@ -80,13 +80,13 @@ end D = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) A = V * D * V' A = (A + A') / 2 - alg = TruncatedAlgorithm(GLA_eigh_Francis(), truncrank(2)) + alg = TruncatedAlgorithm(GLA_QRIteration(), truncrank(2)) D2, V2, ϵ2 = @constinferred eigh_trunc(A; alg) @test diagview(D2) ≈ diagview(D)[1:2] @test_throws ArgumentError eigh_trunc(A; alg, trunc = (; maxrank = 2)) @test ϵ2 ≈ norm(diagview(D)[3:4]) atol = atol - alg = TruncatedAlgorithm(GLA_eigh_Francis(), truncerror(; atol = 0.2)) + alg = TruncatedAlgorithm(GLA_QRIteration(), truncerror(; atol = 0.2)) D3, V3, ϵ3 = @constinferred eigh_trunc(A; alg) @test diagview(D3) ≈ diagview(D)[1:2] @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol @@ -110,7 +110,7 @@ end @test diagview(D) ≈ D2 A2 = Diagonal(T[0.9, 0.3, 0.1, 0.01]) - alg = TruncatedAlgorithm(GLA_eigh_Francis(), truncrank(2)) + alg = TruncatedAlgorithm(DiagonalAlgorithm(), truncrank(2)) D2, V2, ϵ2 = @constinferred eigh_trunc(A2; alg) @test diagview(D2) ≈ diagview(A2)[1:2] @test ϵ2 ≈ norm(diagview(A2)[3:4]) atol = atol diff --git a/test/bigfloats/svd.jl b/test/bigfloats/svd.jl index 5f7884bd..a127a5f6 100644 --- a/test/bigfloats/svd.jl +++ b/test/bigfloats/svd.jl @@ -13,7 +13,7 @@ eltypes = (BigFloat, Complex{BigFloat}) m = 54 @testset "size ($m, $n)" for n in (37, m, 63, 0) k = min(m, n) - alg = GLA_svd_QRIteration() + alg = GLA_QRIteration() minmn = min(m, n) A = randn(rng, T, m, n) @@ -52,7 +52,7 @@ end rng = StableRNG(123) m = 54 @testset "size ($m, $n)" for n in (37, m, 63, 0) - alg = GLA_svd_QRIteration() + alg = GLA_QRIteration() A = randn(rng, T, m, n) U, S, Vᴴ = svd_full(A; alg) @test U isa Matrix{T} && size(U) == (m, m) @@ -98,7 +98,7 @@ end rng = StableRNG(123) m = 54 atol = sqrt(eps(real(T))) - alg = GLA_svd_QRIteration() + alg = GLA_QRIteration() @testset "size ($m, $n)" for n in (37, m, 63) n > m && alg isa LAPACK_Jacobi && continue # not supported @@ -136,7 +136,7 @@ end @testset "svd_trunc! mix maxrank and tol for T = $T" for T in eltypes rng = StableRNG(123) - alg = GLA_svd_QRIteration() + alg = GLA_QRIteration() m = 4 U = qr_compact(randn(rng, T, m, m))[1] S = Diagonal(T[0.9, 0.3, 0.1, 0.01]) @@ -165,7 +165,7 @@ end S = Diagonal(real(T)[0.9, 0.3, 0.1, 0.01]) Vᴴ = qr_compact(randn(rng, T, m, m))[1] A = U * S * Vᴴ - alg = TruncatedAlgorithm(GLA_svd_QRIteration(), trunctol(; atol = 0.2)) + alg = TruncatedAlgorithm(GLA_QRIteration(), trunctol(; atol = 0.2)) U2, S2, V2ᴴ, ϵ2 = @constinferred svd_trunc(A; alg) @test diagview(S2) ≈ diagview(S)[1:2] @test ϵ2 ≈ norm(diagview(S)[3:4]) atol = atol From 78f92f7472d4571eaedd7b46895a5799db9e5ccf Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Tue, 4 Nov 2025 15:27:36 +0100 Subject: [PATCH 16/25] change folder names --- test/{bigfloats => genericlinearalgebra}/eigh.jl | 0 test/{bigfloats => genericlinearalgebra}/lq.jl | 0 test/{bigfloats => genericlinearalgebra}/qr.jl | 0 test/{bigfloats => genericlinearalgebra}/svd.jl | 0 test/{bigfloats => genericschur}/eig.jl | 0 test/runtests.jl | 14 ++++++++------ 6 files changed, 8 insertions(+), 6 deletions(-) rename test/{bigfloats => genericlinearalgebra}/eigh.jl (100%) rename test/{bigfloats => genericlinearalgebra}/lq.jl (100%) rename test/{bigfloats => genericlinearalgebra}/qr.jl (100%) rename test/{bigfloats => genericlinearalgebra}/svd.jl (100%) rename test/{bigfloats => genericschur}/eig.jl (100%) diff --git a/test/bigfloats/eigh.jl b/test/genericlinearalgebra/eigh.jl similarity index 100% rename from test/bigfloats/eigh.jl rename to test/genericlinearalgebra/eigh.jl diff --git a/test/bigfloats/lq.jl b/test/genericlinearalgebra/lq.jl similarity index 100% rename from test/bigfloats/lq.jl rename to test/genericlinearalgebra/lq.jl diff --git a/test/bigfloats/qr.jl b/test/genericlinearalgebra/qr.jl similarity index 100% rename from test/bigfloats/qr.jl rename to test/genericlinearalgebra/qr.jl diff --git a/test/bigfloats/svd.jl b/test/genericlinearalgebra/svd.jl similarity index 100% rename from test/bigfloats/svd.jl rename to test/genericlinearalgebra/svd.jl diff --git a/test/bigfloats/eig.jl b/test/genericschur/eig.jl similarity index 100% rename from test/bigfloats/eig.jl rename to test/genericschur/eig.jl diff --git a/test/runtests.jl b/test/runtests.jl index 7a6bb6b4..ec255538 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -108,17 +108,19 @@ if AMDGPU.functional() end end -using GenericLinearAlgebra, GenericSchur +using GenericLinearAlgebra @safetestset "QR / LQ Decomposition" begin - include("bigfloats/qr.jl") - include("bigfloats/lq.jl") + include("genericlinearalgebra/qr.jl") + include("genericlinearalgebra/lq.jl") end @safetestset "Singular Value Decomposition" begin - include("bigfloats/svd.jl") + include("genericlinearalgebra/svd.jl") end @safetestset "Hermitian Eigenvalue Decomposition" begin - include("bigfloats/eigh.jl") + include("genericlinearalgebra/eigh.jl") end + +using GenericSchur @safetestset "General Eigenvalue Decomposition" begin - include("bigfloats/eig.jl") + include("genericschur/eig.jl") end From fb1994f8b76624d11643120dbcd4633d2cab9181 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Tue, 4 Nov 2025 15:40:04 +0100 Subject: [PATCH 17/25] add docs and namechange GS and GLA algorithms are now defined in the docs `GLA_QR_Householder` is now `GLA_HouseholderQR` to be consistent with the LAPACK algorithms --- ...MatrixAlgebraKitGenericLinearAlgebraExt.jl | 12 ++++----- src/MatrixAlgebraKit.jl | 2 +- src/interface/decompositions.jl | 26 ++++++++++++++++--- test/genericlinearalgebra/lq.jl | 6 ++--- 4 files changed, 33 insertions(+), 13 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 198e5ec7..0bc246b5 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -66,10 +66,10 @@ function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::GLA_QRIterati end function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return GLA_QR_Householder(; kwargs...) + return GLA_HouseholderQR(; kwargs...) end -function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::GLA_QR_Householder) +function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::GLA_HouseholderQR) Q, R = QR m, n = size(A) minmn = min(m, n) @@ -86,7 +86,7 @@ function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::GLA_QR_Householde return Q, R end -function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::GLA_QR_Householder) +function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::GLA_HouseholderQR) check_input(qr_compact!, A, QR, alg) Q, R = QR Q, R = _gla_householder_qr!(A, Q, R; alg.kwargs...) @@ -94,8 +94,8 @@ function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::GLA_QR_Househo end function _gla_householder_qr!(A::AbstractMatrix{T}, Q, R; positive = false, blocksize = 1, pivoted = false) where {T} - pivoted && throw(ArgumentError("Only pivoted = false implemented for GLA_QR_Householder.")) - (blocksize == 1) || throw(ArgumentError("Only blocksize = 1 implemented for GLA_QR_Householder.")) + pivoted && throw(ArgumentError("Only pivoted = false implemented for GLA_HouseholderQR.")) + (blocksize == 1) || throw(ArgumentError("Only blocksize = 1 implemented for GLA_HouseholderQR.")) m, n = size(A) k = min(m, n) @@ -146,7 +146,7 @@ function _gram_schmidt!(Q, Q_compact; adjoint = false) end function MatrixAlgebraKit.default_lq_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} - return MatrixAlgebraKit.LQViaTransposedQR(GLA_QR_Householder(; kwargs...)) + return MatrixAlgebraKit.LQViaTransposedQR(GLA_HouseholderQR(; kwargs...)) end end diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index 1fd4fabe..4a846f85 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -33,7 +33,7 @@ export left_orth!, right_orth!, left_null!, right_null! export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert, LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations, LAPACK_DivideAndConquer, LAPACK_Jacobi -export GLA_QR_Householder, GLA_QRIteration, GS_QRIteration +export GLA_HouseholderQR, GLA_QRIteration, GS_QRIteration export LQViaTransposedQR export PolarViaSVD, PolarNewton export DiagonalAlgorithm diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index b23fbeab..68a46ddc 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -32,8 +32,17 @@ elements of `L` are non-negative. """ @algdef LAPACK_HouseholderLQ +""" + GLA_HouseholderQR(; positive = false) + +Algorithm type to denote the GenericLinearAlgebra.jl implementation for computing the QR decomposition +of a matrix using Householder reflectors. Currenlty, only `blocksize = 1` and `pivoted == false` +are supported. The keyword `positive=true` can be used to ensure that the diagonal elements +of `R` are non-negative. +""" +@algdef GLA_HouseholderQR + # TODO: -@algdef GLA_QR_Householder @algdef LAPACK_HouseholderQL @algdef LAPACK_HouseholderRQ @@ -57,7 +66,12 @@ eigenvalue decomposition of a matrix. const LAPACK_EigAlgorithm = Union{LAPACK_Simple, LAPACK_Expert} -# TODO: +""" + GS_QRIteration() + +Algorithm type to denote the GenericSchur.jl implementation for computing the +eigenvalue decomposition of a non-Hermitian matrix. +""" @algdef GS_QRIteration # Hermitian Eigenvalue Decomposition @@ -104,7 +118,13 @@ const LAPACK_EighAlgorithm = Union{ LAPACK_MultipleRelativelyRobustRepresentations, } -# TODO: +""" + GLA_QRIteration() + +Algorithm type to denote the GenericLinearAlgebra.jl implementation for computing the +eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of +a general matrix. +""" @algdef GLA_QRIteration # Singular Value Decomposition diff --git a/test/genericlinearalgebra/lq.jl b/test/genericlinearalgebra/lq.jl index d071afee..b3ce2379 100644 --- a/test/genericlinearalgebra/lq.jl +++ b/test/genericlinearalgebra/lq.jl @@ -32,7 +32,7 @@ eltypes = (BigFloat, Complex{BigFloat}) @test Q == Q2 # Transposed QR algorithm - qr_alg = GLA_QR_Householder() + qr_alg = GLA_HouseholderQR() lq_alg = LQViaTransposedQR(qr_alg) L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q), lq_alg) @test L2 === L @@ -80,7 +80,7 @@ end @test Q[1:minmn, n] ≈ Q2[1:minmn, n] # Transposed QR algorithm - qr_alg = GLA_QR_Householder() + qr_alg = GLA_HouseholderQR() lq_alg = LQViaTransposedQR(qr_alg) L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q), lq_alg) @test L2 === L @@ -104,7 +104,7 @@ end lq_full!(copy!(Ac, A), (noL, Q2); positive = true) @test Q[1:minmn, n] ≈ Q2[1:minmn, n] - qr_alg = GLA_QR_Householder(; positive = true) + qr_alg = GLA_HouseholderQR(; positive = true) lq_alg = LQViaTransposedQR(qr_alg) lq_full!(copy!(Ac, A), (L, Q), lq_alg) @test L * Q ≈ A From 3a364460dbe942fcb20e0e3f8e3884d44a1ce23e Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 4 Nov 2025 11:36:31 -0500 Subject: [PATCH 18/25] Remove unnecessary type parameters --- ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 0bc246b5..8d7042e2 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -40,8 +40,8 @@ function MatrixAlgebraKit.svd_full!(A::AbstractMatrix, USVᴴ, alg::GLA_QRIterat return MatrixAlgebraKit.gaugefix!(svd_full!, U, S, Vᴴ, m, n) end -function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix{T}, S, alg::GLA_QRIteration) where {T} check_input(svd_vals!, A, S, alg) +function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix, S, ::GLA_QRIteration) S̃ = svdvals!(A) copyto!(S, S̃) return S @@ -51,7 +51,7 @@ function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T return GLA_QRIteration(; kwargs...) end -function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_QRIteration) where {T} +function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix, DV, alg::GLA_QRIteration) check_input(eigh_full!, A, DV, alg) D, V = DV eigval, eigvec = eigen!(Hermitian(A); sortby = real) @@ -60,7 +60,7 @@ function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix{T}, DV, alg::GLA_QRIterat return D, V end -function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix{T}, D, alg::GLA_QRIteration) where {T} +function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix, D, alg::GLA_QRIteration) check_input(eigh_vals!, A, D, alg) return eigvals!(Hermitian(A); sortby = real) end @@ -93,7 +93,7 @@ function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::GLA_Householde return Q, R end -function _gla_householder_qr!(A::AbstractMatrix{T}, Q, R; positive = false, blocksize = 1, pivoted = false) where {T} +function _gla_householder_qr!(A::AbstractMatrix, Q, R; positive = false, blocksize = 1, pivoted = false) pivoted && throw(ArgumentError("Only pivoted = false implemented for GLA_HouseholderQR.")) (blocksize == 1) || throw(ArgumentError("Only blocksize = 1 implemented for GLA_HouseholderQR.")) From e6f64b2e956366620de8ff29e90ee8e95c538802 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 4 Nov 2025 11:44:51 -0500 Subject: [PATCH 19/25] simplify SVD and remove allocations --- ...MatrixAlgebraKitGenericLinearAlgebraExt.jl | 46 ++++++------------- test/genericlinearalgebra/svd.jl | 16 +++---- 2 files changed, 22 insertions(+), 40 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 8d7042e2..e499a80d 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -1,7 +1,7 @@ module MatrixAlgebraKitGenericLinearAlgebraExt using MatrixAlgebraKit -using MatrixAlgebraKit: sign_safe, check_input +using MatrixAlgebraKit: sign_safe, check_input, diagview using GenericLinearAlgebra: svd!, svdvals!, eigen!, eigvals!, Hermitian, qr! using LinearAlgebra: I, Diagonal, rmul!, lmul!, transpose!, dot @@ -9,42 +9,26 @@ function MatrixAlgebraKit.default_svd_algorithm(::Type{T}; kwargs...) where {T < return GLA_QRIteration() end -function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix, USVᴴ, alg::GLA_QRIteration) - check_input(svd_compact!, A, USVᴴ, alg) - U, S, Vᴴ = USVᴴ - F = svd!(A) - copyto!(U, F.U) - copyto!(S, Diagonal(F.S)) - copyto!(Vᴴ, F.Vt) # conjugation to account for difference in convention - return U, S, Vᴴ +for f! in (:svd_compact!, :svd_full!, :svd_vals!) + @eval MatrixAlgebraKit.initialize_output(::typeof($f!), A::AbstractMatrix, ::GLA_QRIteration) = nothing end -function MatrixAlgebraKit.svd_full!(A::AbstractMatrix, USVᴴ, alg::GLA_QRIteration) - check_input(svd_full!, A, USVᴴ, alg) - U, S, Vᴴ = USVᴴ - m, n = size(A) - minmn = min(m, n) - if minmn == 0 - MatrixAlgebraKit.one!(U) - MatrixAlgebraKit.zero!(S) - MatrixAlgebraKit.one!(Vᴴ) - return USVᴴ - end - S = MatrixAlgebraKit.zero!(S) - U_compact, S_compact, Vᴴ_compact = svd!(A) - S[1:minmn, 1:minmn] .= Diagonal(S_compact) - - U = _gram_schmidt!(U, U_compact) - Vᴴ = _gram_schmidt!(Vᴴ, Vᴴ_compact; adjoint = true) +function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix, USVᴴ, ::GLA_QRIteration) + F = svd!(A) + U, S, Vᴴ = F.U, Diagonal(F.S), F.Vt + return MatrixAlgebraKit.gaugefix!(svd_compact!, U, S, Vᴴ, size(A)...) +end - return MatrixAlgebraKit.gaugefix!(svd_full!, U, S, Vᴴ, m, n) +function MatrixAlgebraKit.svd_full!(A::AbstractMatrix, USVᴴ, ::GLA_QRIteration) + F = svd!(A; full = true) + U, Vᴴ = F.U, F.Vt + S = MatrixAlgebraKit.zero!(similar(F.S, (size(U, 2), size(Vᴴ, 1)))) + diagview(S) .= F.S + return MatrixAlgebraKit.gaugefix!(svd_full!, U, S, Vᴴ, size(A)...) end - check_input(svd_vals!, A, S, alg) function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix, S, ::GLA_QRIteration) - S̃ = svdvals!(A) - copyto!(S, S̃) - return S + return svdvals!(A) end function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} diff --git a/test/genericlinearalgebra/svd.jl b/test/genericlinearalgebra/svd.jl index a127a5f6..4e9c6d87 100644 --- a/test/genericlinearalgebra/svd.jl +++ b/test/genericlinearalgebra/svd.jl @@ -35,9 +35,9 @@ eltypes = (BigFloat, Complex{BigFloat}) Sc = similar(A, real(T), min(m, n)) alg′ = @constinferred MatrixAlgebraKit.select_algorithm(svd_compact!, A, $alg) U2, S2, V2ᴴ = @constinferred svd_compact!(copy!(Ac, A), (U, S, Vᴴ), alg′) - @test U2 === U - @test S2 === S - @test V2ᴴ === Vᴴ + @test U2 ≈ U + @test S2 ≈ S + @test V2ᴴ ≈ Vᴴ @test U * S * Vᴴ ≈ A @test isisometric(U) @test isisometric(Vᴴ; side = :right) @@ -65,17 +65,15 @@ end Ac = similar(A) U2, S2, V2ᴴ = @constinferred svd_full!(copy!(Ac, A), (U, S, Vᴴ), alg) - @test U2 === U - @test S2 === S - @test V2ᴴ === Vᴴ + @test U2 ≈ U + @test S2 ≈ S + @test V2ᴴ ≈ Vᴴ @test U * S * Vᴴ ≈ A @test isunitary(U) @test isunitary(Vᴴ) @test all(isposdef, diagview(S)) - Sc = similar(A, real(T), min(m, n)) - Sc2 = svd_vals!(copy!(Ac, A), Sc, alg) - @test Sc === Sc2 + Sc = svd_vals!(copy!(Ac, A), alg) @test diagview(S) ≈ Sc end @testset "size (0, 0)" begin From 4323b17d24396d7dcd7527ce8a0f820faeef2c76 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 4 Nov 2025 12:03:19 -0500 Subject: [PATCH 20/25] simplify Eigh and remove allocations --- ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl | 15 +++++++-------- test/genericlinearalgebra/eigh.jl | 4 ++-- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index e499a80d..96bcc419 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -35,17 +35,16 @@ function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T return GLA_QRIteration(; kwargs...) end -function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix, DV, alg::GLA_QRIteration) - check_input(eigh_full!, A, DV, alg) - D, V = DV +for f! in (:eigh_full!, :eigh_vals!) + @eval MatrixAlgebraKit.initialize_output(::typeof($f!), A::AbstractMatrix, ::GLA_QRIteration) = nothing +end + +function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix, DV, ::GLA_QRIteration) eigval, eigvec = eigen!(Hermitian(A); sortby = real) - copyto!(D, Diagonal(eigval)) - copyto!(V, eigvec) - return D, V + return Diagonal(eigval::AbstractVector{real(eltype(A))}), eigvec::AbstractMatrix{eltype(A)} end -function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix, D, alg::GLA_QRIteration) - check_input(eigh_vals!, A, D, alg) +function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix, D, ::GLA_QRIteration) return eigvals!(Hermitian(A); sortby = real) end diff --git a/test/genericlinearalgebra/eigh.jl b/test/genericlinearalgebra/eigh.jl index a176981d..9d8400ff 100644 --- a/test/genericlinearalgebra/eigh.jl +++ b/test/genericlinearalgebra/eigh.jl @@ -22,8 +22,8 @@ const eltypes = (BigFloat, Complex{BigFloat}) @test all(isreal, D) D2, V2 = eigh_full!(copy(A), (D, V), alg) - @test D2 === D - @test V2 === V + @test D2 ≈ D + @test V2 ≈ V D3 = @constinferred eigh_vals(A, alg) @test D ≈ Diagonal(D3) From f99112b9cce357043691fb3e712473dea35472ef Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 4 Nov 2025 12:17:22 -0500 Subject: [PATCH 21/25] simplify QR --- ...MatrixAlgebraKitGenericLinearAlgebraExt.jl | 51 +++++-------------- 1 file changed, 12 insertions(+), 39 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 96bcc419..9150876d 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -53,27 +53,15 @@ function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T <: end function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::GLA_HouseholderQR) + check_input(qr_full!, A, QR, alg) Q, R = QR - m, n = size(A) - minmn = min(m, n) - computeR = length(R) > 0 - - Q_zero = zeros(eltype(Q), (m, minmn)) - R_zero = zeros(eltype(R), (minmn, n)) - Q_compact, R_compact = _gla_householder_qr!(A, Q_zero, R_zero; alg.kwargs...) - Q = _gram_schmidt!(Q, view(Q_compact, :, 1:minmn)) - if computeR - R[1:minmn, :] .= R_compact - R[(minmn + 1):end, :] .= zero(eltype(R)) - end - return Q, R + return _gla_householder_qr!(A, Q, R; alg.kwargs...) end function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::GLA_HouseholderQR) check_input(qr_compact!, A, QR, alg) Q, R = QR - Q, R = _gla_householder_qr!(A, Q, R; alg.kwargs...) - return Q, R + return _gla_householder_qr!(A, Q, R; alg.kwargs...) end function _gla_householder_qr!(A::AbstractMatrix, Q, R; positive = false, blocksize = 1, pivoted = false) @@ -83,13 +71,10 @@ function _gla_householder_qr!(A::AbstractMatrix, Q, R; positive = false, blocksi m, n = size(A) k = min(m, n) computeR = length(R) > 0 + compact = k < m Q̃, R̃ = qr!(A) - if k < m - copyto!(Q, Q̃) - else - rmul!(MatrixAlgebraKit.one!(Q), Q̃) - end + Q = compact ? copyto!(Q, Q̃) : rmul!(MatrixAlgebraKit.one!(Q), Q̃) if positive @inbounds for j in 1:k s = sign_safe(R̃[j, j]) @@ -102,32 +87,20 @@ function _gla_householder_qr!(A::AbstractMatrix, Q, R; positive = false, blocksi if positive @inbounds for j in n:-1:1 @simd for i in 1:min(k, j) - R̃[i, j] = R̃[i, j] * conj(sign_safe(R̃[i, i])) + R[i, j] = R̃[i, j] * conj(sign_safe(R̃[i, i])) + end + @simd for i in (min(k, j) + 1):size(R, 1) + R[i, j] = zero(eltype(R)) end end + else + R[1:k, :] .= R̃ + MatrixAlgebraKit.zero!(@view(R[(k + 1):end, :])) end - copyto!(R, R̃) end return Q, R end -function _gram_schmidt!(Q, Q_compact; adjoint = false) - m, minmn = size(Q_compact) - Q[:, 1:minmn] .= Q_compact - for j in (minmn + 1):m - v = rand(eltype(Q), m) - for i in 1:(j - 1) - r = dot(view(Q, :, i), v) - v .-= r * view(Q, :, i) - end - Q[:, j] = v ./ MatrixAlgebraKit.norm(v) - end - if adjoint - copyto!(Q, Q') - end - return Q -end - function MatrixAlgebraKit.default_lq_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} return MatrixAlgebraKit.LQViaTransposedQR(GLA_HouseholderQR(; kwargs...)) end From 89fd42b9c8755155e5fde3d759414d7f125723cc Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 4 Nov 2025 13:03:53 -0500 Subject: [PATCH 22/25] simplify Eig and remove allocations --- ext/MatrixAlgebraKitGenericSchurExt.jl | 21 +++++++++------------ test/genericschur/eig.jl | 4 ++-- 2 files changed, 11 insertions(+), 14 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericSchurExt.jl b/ext/MatrixAlgebraKitGenericSchurExt.jl index 04d45769..d278b5c5 100644 --- a/ext/MatrixAlgebraKitGenericSchurExt.jl +++ b/ext/MatrixAlgebraKitGenericSchurExt.jl @@ -9,20 +9,17 @@ function MatrixAlgebraKit.default_eig_algorithm(::Type{T}; kwargs...) where {T < return GS_QRIteration(; kwargs...) end -function MatrixAlgebraKit.eig_full!(A::AbstractMatrix{T}, DV, alg::GS_QRIteration) where {T} - check_input(eig_full!, A, DV, alg) - D, V = DV - D̃, Ṽ = GenericSchur.eigen!(A) - copyto!(D, Diagonal(D̃)) - copyto!(V, Ṽ) - return D, V +for f! in (:eig_full!, :eig_vals!) + @eval MatrixAlgebraKit.initialize_output(::typeof($f!), A::AbstractMatrix, ::GS_QRIteration) = nothing end -function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix{T}, D, alg::GS_QRIteration) where {T} - check_input(eig_vals!, A, D, alg) - eigval = GenericSchur.eigvals!(A) - copyto!(D, eigval) - return D +function MatrixAlgebraKit.eig_full!(A::AbstractMatrix, DV, ::GS_QRIteration) + D, V = GenericSchur.eigen!(A) + return Diagonal(D), V +end + +function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix, D, ::GS_QRIteration) + return GenericSchur.eigvals!(A) end end diff --git a/test/genericschur/eig.jl b/test/genericschur/eig.jl index 88de1b60..ce1e8f1b 100644 --- a/test/genericschur/eig.jl +++ b/test/genericschur/eig.jl @@ -23,8 +23,8 @@ const eltypes = (BigFloat, Complex{BigFloat}) Ac = similar(A) D2, V2 = @constinferred eig_full!(copy!(Ac, A), (D, V), alg′) - @test D2 === D - @test V2 === V + @test D2 ≈ D + @test V2 ≈ V @test A * V ≈ V * D Dc = @constinferred eig_vals(A, alg′) From 35af44d93eea99a17190c2e7ca8bbed5027b6573 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 4 Nov 2025 19:20:06 -0500 Subject: [PATCH 23/25] switch to `lmul!` --- ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 9150876d..7b078381 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -3,7 +3,7 @@ module MatrixAlgebraKitGenericLinearAlgebraExt using MatrixAlgebraKit using MatrixAlgebraKit: sign_safe, check_input, diagview using GenericLinearAlgebra: svd!, svdvals!, eigen!, eigvals!, Hermitian, qr! -using LinearAlgebra: I, Diagonal, rmul!, lmul!, transpose!, dot +using LinearAlgebra: I, Diagonal, lmul! function MatrixAlgebraKit.default_svd_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} return GLA_QRIteration() @@ -70,11 +70,9 @@ function _gla_householder_qr!(A::AbstractMatrix, Q, R; positive = false, blocksi m, n = size(A) k = min(m, n) - computeR = length(R) > 0 - compact = k < m Q̃, R̃ = qr!(A) + lmul!(Q̃, MatrixAlgebraKit.one!(Q)) - Q = compact ? copyto!(Q, Q̃) : rmul!(MatrixAlgebraKit.one!(Q), Q̃) if positive @inbounds for j in 1:k s = sign_safe(R̃[j, j]) @@ -83,6 +81,8 @@ function _gla_householder_qr!(A::AbstractMatrix, Q, R; positive = false, blocksi end end end + + computeR = length(R) > 0 if computeR if positive @inbounds for j in n:-1:1 From 4fe0fe176516f03b779873c9c59a3f3c9907fce7 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 4 Nov 2025 19:22:53 -0500 Subject: [PATCH 24/25] docstring improvements --- src/interface/decompositions.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index 68a46ddc..1bdf1534 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -16,7 +16,7 @@ Algorithm type to denote the standard LAPACK algorithm for computing the QR deco a matrix using Householder reflectors. The specific LAPACK function can be controlled using the keyword arugments, i.e. `?geqrt` will be chosen if `blocksize > 1`. With `blocksize == 1`, `?geqrf` will be chosen if `pivoted == false` and `?geqp3` will be chosen -if `pivoted == true`. The keyword `positive=true` can be used to ensure that the diagonal +if `pivoted == true`. The keyword `positive = true` can be used to ensure that the diagonal elements of `R` are non-negative. """ @algdef LAPACK_HouseholderQR @@ -27,7 +27,7 @@ elements of `R` are non-negative. Algorithm type to denote the standard LAPACK algorithm for computing the LQ decomposition of a matrix using Householder reflectors. The specific LAPACK function can be controlled using the keyword arugments, i.e. `?gelqt` will be chosen if `blocksize > 1` or `?gelqf` will be -chosen if `blocksize == 1`. The keyword `positive=true` can be used to ensure that the diagonal +chosen if `blocksize == 1`. The keyword `positive = true` can be used to ensure that the diagonal elements of `L` are non-negative. """ @algdef LAPACK_HouseholderLQ @@ -35,9 +35,9 @@ elements of `L` are non-negative. """ GLA_HouseholderQR(; positive = false) -Algorithm type to denote the GenericLinearAlgebra.jl implementation for computing the QR decomposition -of a matrix using Householder reflectors. Currenlty, only `blocksize = 1` and `pivoted == false` -are supported. The keyword `positive=true` can be used to ensure that the diagonal elements +Algorithm type to denote the GenericLinearAlgebra.jl implementation for computing the QR decomposition +of a matrix using Householder reflectors. Currently, only `blocksize = 1` and `pivoted == false` +are supported. The keyword `positive = true` can be used to ensure that the diagonal elements of `R` are non-negative. """ @algdef GLA_HouseholderQR @@ -69,7 +69,7 @@ const LAPACK_EigAlgorithm = Union{LAPACK_Simple, LAPACK_Expert} """ GS_QRIteration() -Algorithm type to denote the GenericSchur.jl implementation for computing the +Algorithm type to denote the GenericSchur.jl implementation for computing the eigenvalue decomposition of a non-Hermitian matrix. """ @algdef GS_QRIteration @@ -121,8 +121,8 @@ const LAPACK_EighAlgorithm = Union{ """ GLA_QRIteration() -Algorithm type to denote the GenericLinearAlgebra.jl implementation for computing the -eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of +Algorithm type to denote the GenericLinearAlgebra.jl implementation for computing the +eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix. """ @algdef GLA_QRIteration From c08dad4917ff9cdd57789f00d1aa974935dded28 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 4 Nov 2025 19:39:39 -0500 Subject: [PATCH 25/25] move Diagonal tests --- test/eig.jl | 5 +++-- test/eigh.jl | 5 +++-- test/genericlinearalgebra/eigh.jl | 25 --------------------- test/genericlinearalgebra/lq.jl | 25 --------------------- test/genericlinearalgebra/qr.jl | 25 --------------------- test/genericlinearalgebra/svd.jl | 36 ------------------------------- test/lq.jl | 9 ++++---- test/qr.jl | 9 ++++---- test/svd.jl | 5 +++-- 9 files changed, 19 insertions(+), 125 deletions(-) diff --git a/test/eig.jl b/test/eig.jl index 0ece1b35..6da6d72c 100644 --- a/test/eig.jl +++ b/test/eig.jl @@ -5,7 +5,8 @@ using StableRNGs using LinearAlgebra: Diagonal using MatrixAlgebraKit: TruncatedAlgorithm, diagview, norm -const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, BigFloat, Complex{BigFloat}) @testset "eig_full! for T = $T" for T in BLASFloats rng = StableRNG(123) @@ -91,7 +92,7 @@ end @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol end -@testset "eig for Diagonal{$T}" for T in BLASFloats +@testset "eig for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) rng = StableRNG(123) m = 54 Ad = randn(rng, T, m) diff --git a/test/eigh.jl b/test/eigh.jl index bc04b057..92b0f3a0 100644 --- a/test/eigh.jl +++ b/test/eigh.jl @@ -5,7 +5,8 @@ using StableRNGs using LinearAlgebra: LinearAlgebra, Diagonal, I using MatrixAlgebraKit: TruncatedAlgorithm, diagview, norm -const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, BigFloat, Complex{BigFloat}) @testset "eigh_full! for T = $T" for T in BLASFloats rng = StableRNG(123) @@ -100,7 +101,7 @@ end @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol end -@testset "eigh for Diagonal{$T}" for T in BLASFloats +@testset "eigh for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) rng = StableRNG(123) m = 54 Ad = randn(rng, T, m) diff --git a/test/genericlinearalgebra/eigh.jl b/test/genericlinearalgebra/eigh.jl index 9d8400ff..7e602026 100644 --- a/test/genericlinearalgebra/eigh.jl +++ b/test/genericlinearalgebra/eigh.jl @@ -91,28 +91,3 @@ end @test diagview(D3) ≈ diagview(D)[1:2] @test ϵ3 ≈ norm(diagview(D)[3:4]) atol = atol end - -@testset "eigh for Diagonal{$T}" for T in eltypes - rng = StableRNG(123) - m = 54 - Ad = randn(rng, T, m) - Ad .+= conj.(Ad) - A = Diagonal(Ad) - atol = sqrt(eps(real(T))) - - D, V = @constinferred eigh_full(A) - @test D isa Diagonal{real(T)} && size(D) == size(A) - @test V isa Diagonal{T} && size(V) == size(A) - @test A * V ≈ V * D - - D2 = @constinferred eigh_vals(A) - @test D2 isa AbstractVector{real(T)} && length(D2) == m - @test diagview(D) ≈ D2 - - A2 = Diagonal(T[0.9, 0.3, 0.1, 0.01]) - alg = TruncatedAlgorithm(DiagonalAlgorithm(), truncrank(2)) - D2, V2, ϵ2 = @constinferred eigh_trunc(A2; alg) - @test diagview(D2) ≈ diagview(A2)[1:2] - @test ϵ2 ≈ norm(diagview(A2)[3:4]) atol = atol - -end diff --git a/test/genericlinearalgebra/lq.jl b/test/genericlinearalgebra/lq.jl index b3ce2379..dc186dcb 100644 --- a/test/genericlinearalgebra/lq.jl +++ b/test/genericlinearalgebra/lq.jl @@ -122,28 +122,3 @@ end @test Q[1:minmn, n] ≈ Q2[1:minmn, n] end end - -@testset "lq_compact for Diagonal{$T}" for T in eltypes - rng = StableRNG(123) - atol = eps(real(T))^(3 / 4) - for m in (54, 0) - Ad = randn(rng, T, m) - A = Diagonal(Ad) - - # compact - L, Q = @constinferred lq_compact(A) - @test Q isa Diagonal{T} && size(Q) == (m, m) - @test L isa Diagonal{T} && size(L) == (m, m) - @test L * Q ≈ A - @test isunitary(Q) - - # compact and positive - Lp, Qp = @constinferred lq_compact(A; positive = true) - @test Qp isa Diagonal{T} && size(Qp) == (m, m) - @test Lp isa Diagonal{T} && size(Lp) == (m, m) - @test Lp * Qp ≈ A - @test isunitary(Qp) - @test all(≥(zero(real(T))), real(diag(Lp))) && - all(≈(zero(real(T)); atol), imag(diag(Lp))) - end -end diff --git a/test/genericlinearalgebra/qr.jl b/test/genericlinearalgebra/qr.jl index 9e1ba515..3ce530bb 100644 --- a/test/genericlinearalgebra/qr.jl +++ b/test/genericlinearalgebra/qr.jl @@ -107,28 +107,3 @@ end end end end - -@testset "qr_compact for Diagonal{$T}" for T in eltypes - rng = StableRNG(123) - atol = eps(real(T))^(3 / 4) - for m in (54, 0) - Ad = randn(rng, T, m) - A = Diagonal(Ad) - - # compact - Q, R = @constinferred qr_compact(A) - @test Q isa Diagonal{T} && size(Q) == (m, m) - @test R isa Diagonal{T} && size(R) == (m, m) - @test Q * R ≈ A - @test isunitary(Q) - - # compact and positive - Qp, Rp = @constinferred qr_compact(A; positive = true) - @test Qp isa Diagonal{T} && size(Qp) == (m, m) - @test Rp isa Diagonal{T} && size(Rp) == (m, m) - @test Qp * Rp ≈ A - @test isunitary(Qp) - @test all(≥(zero(real(T))), real(diag(Rp))) && - all(≈(zero(real(T)); atol), imag(diag(Rp))) - end -end diff --git a/test/genericlinearalgebra/svd.jl b/test/genericlinearalgebra/svd.jl index 4e9c6d87..f7177e79 100644 --- a/test/genericlinearalgebra/svd.jl +++ b/test/genericlinearalgebra/svd.jl @@ -169,39 +169,3 @@ end @test ϵ2 ≈ norm(diagview(S)[3:4]) atol = atol @test_throws ArgumentError svd_trunc(A; alg, trunc = (; maxrank = 2)) end - -@testset "svd for Diagonal{$T}" for T in eltypes - rng = StableRNG(123) - atol = sqrt(eps(real(T))) - for m in (54, 0) - Ad = randn(T, m) - A = Diagonal(Ad) - - U, S, Vᴴ = @constinferred svd_compact(A) - @test U isa AbstractMatrix{T} && size(U) == size(A) - @test Vᴴ isa AbstractMatrix{T} && size(Vᴴ) == size(A) - @test S isa Diagonal{real(T)} && size(S) == size(A) - @test isunitary(U) - @test isunitary(Vᴴ) - @test all(≥(0), diagview(S)) - @test A ≈ U * S * Vᴴ - - U, S, Vᴴ = @constinferred svd_full(A) - @test U isa AbstractMatrix{T} && size(U) == size(A) - @test Vᴴ isa AbstractMatrix{T} && size(Vᴴ) == size(A) - @test S isa Diagonal{real(T)} && size(S) == size(A) - @test isunitary(U) - @test isunitary(Vᴴ) - @test all(≥(0), diagview(S)) - @test A ≈ U * S * Vᴴ - - S2 = @constinferred svd_vals(A) - @test S2 isa AbstractVector{real(T)} && length(S2) == m - @test S2 ≈ diagview(S) - - alg = TruncatedAlgorithm(DiagonalAlgorithm(), truncrank(2)) - U3, S3, Vᴴ3, ϵ3 = @constinferred svd_trunc(A; alg) - @test diagview(S3) ≈ S2[1:min(m, 2)] - @test ϵ3 ≈ norm(S2[(min(m, 2) + 1):m]) atol = atol - end -end diff --git a/test/lq.jl b/test/lq.jl index 2c8dfefe..8de5a582 100644 --- a/test/lq.jl +++ b/test/lq.jl @@ -5,9 +5,10 @@ using StableRNGs using LinearAlgebra: diag, I, Diagonal using MatrixAlgebraKit: LQViaTransposedQR, LAPACK_HouseholderQR -eltypes = (Float32, Float64, ComplexF32, ComplexF64) +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, BigFloat, Complex{BigFloat}) -@testset "lq_compact! for T = $T" for T in eltypes +@testset "lq_compact! for T = $T" for T in BLASFloats rng = StableRNG(123) m = 54 for n in (37, m, 63) @@ -114,7 +115,7 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64) end end -@testset "lq_full! for T = $T" for T in eltypes +@testset "lq_full! for T = $T" for T in BLASFloats rng = StableRNG(123) m = 54 for n in (37, m, 63) @@ -208,7 +209,7 @@ end end end -@testset "lq_compact, lq_full and lq_null for Diagonal{$T}" for T in eltypes +@testset "lq_compact, lq_full and lq_null for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) rng = StableRNG(123) atol = eps(real(T))^(3 / 4) for m in (54, 0) diff --git a/test/qr.jl b/test/qr.jl index 826c320b..c4f0c9d6 100644 --- a/test/qr.jl +++ b/test/qr.jl @@ -4,9 +4,10 @@ using TestExtras using StableRNGs using LinearAlgebra: diag, I, Diagonal -eltypes = (Float32, Float64, ComplexF32, ComplexF64) +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, BigFloat, Complex{BigFloat}) -@testset "qr_compact! and qr_null! for T = $T" for T in eltypes +@testset "qr_compact! and qr_null! for T = $T" for T in BLASFloats rng = StableRNG(123) m = 54 for n in (37, m, 63) @@ -99,7 +100,7 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64) end end -@testset "qr_full! for T = $T" for T in eltypes +@testset "qr_full! for T = $T" for T in BLASFloats rng = StableRNG(123) m = 54 for n in (37, m, 63) @@ -176,7 +177,7 @@ end end end -@testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in eltypes +@testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) rng = StableRNG(123) atol = eps(real(T))^(3 / 4) for m in (54, 0) diff --git a/test/svd.jl b/test/svd.jl index acb27946..d055f866 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -5,7 +5,8 @@ using StableRNGs using LinearAlgebra: LinearAlgebra, Diagonal, I, isposdef, norm using MatrixAlgebraKit: TruncatedAlgorithm, diagview, isisometric -const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, BigFloat, Complex{BigFloat}) @testset "svd_compact! for T = $T" for T in BLASFloats rng = StableRNG(123) @@ -202,7 +203,7 @@ end @test_throws ArgumentError svd_trunc(A; alg, trunc = (; maxrank = 2)) end -@testset "svd for Diagonal{$T}" for T in BLASFloats +@testset "svd for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) rng = StableRNG(123) atol = sqrt(eps(real(T))) for m in (54, 0)