From 6e83d2d1764242e4cc80ab388e5819e4f379b466 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Wed, 15 Oct 2025 18:19:19 -0400 Subject: [PATCH 1/3] Add tests for image and null space for GPU --- src/implementations/truncation.jl | 5 +- test/amd/orthnull.jl | 384 ++++++++++++++++++++++++++++++ test/cuda/orthnull.jl | 376 +++++++++++++++++++++++++++++ test/runtests.jl | 6 + 4 files changed, 770 insertions(+), 1 deletion(-) create mode 100644 test/amd/orthnull.jl create mode 100644 test/cuda/orthnull.jl diff --git a/src/implementations/truncation.jl b/src/implementations/truncation.jl index f6201b07..ebe52825 100644 --- a/src/implementations/truncation.jl +++ b/src/implementations/truncation.jl @@ -17,7 +17,10 @@ function truncate(::typeof(left_null!), (U, S), strategy::TruncationStrategy) # TODO: avoid allocation? extended_S = vcat(diagview(S), zeros(eltype(S), max(0, size(S, 1) - size(S, 2)))) ind = findtruncated(extended_S, strategy) - return U[:, ind], ind + trunc_cols = collect(1:size(U, 2))[ind] + Utrunc = similar(U, (size(U, 1), length(trunc_cols))) + Utrunc .= U[:, trunc_cols] + return Utrunc, ind end function truncate(::typeof(right_null!), (S, Vᴴ), strategy::TruncationStrategy) # TODO: avoid allocation? diff --git a/test/amd/orthnull.jl b/test/amd/orthnull.jl new file mode 100644 index 00000000..5f81ffe0 --- /dev/null +++ b/test/amd/orthnull.jl @@ -0,0 +1,384 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using LinearAlgebra: LinearAlgebra, I, mul!, diagm, norm +using MatrixAlgebraKit: GPU_SVDAlgorithm, check_input, copy_input, default_svd_algorithm, + initialize_output, AbstractAlgorithm +using AMDGPU + +# Used to test non-AbstractMatrix codepaths. +struct LinearMap{P <: AbstractMatrix} + parent::P +end +Base.parent(A::LinearMap) = getfield(A, :parent) +function Base.copy!(dest::LinearMap, src::LinearMap) + copy!(parent(dest), parent(src)) + return dest +end +function LinearAlgebra.mul!(C::LinearMap, A::LinearMap, B::LinearMap) + mul!(parent(C), parent(A), parent(B)) + return C +end + +function MatrixAlgebraKit.copy_input(::typeof(qr_compact), A::LinearMap) + return LinearMap(copy_input(qr_compact, parent(A))) +end +function MatrixAlgebraKit.copy_input(::typeof(lq_compact), A::LinearMap) + return LinearMap(copy_input(lq_compact, parent(A))) +end +function MatrixAlgebraKit.initialize_output(::typeof(left_orth!), A::LinearMap) + return LinearMap.(initialize_output(left_orth!, parent(A))) +end +function MatrixAlgebraKit.initialize_output(::typeof(right_orth!), A::LinearMap) + return LinearMap.(initialize_output(right_orth!, parent(A))) +end +function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::LinearMap, VC, alg::AbstractAlgorithm) + return check_input(left_orth!, parent(A), parent.(VC), alg) +end +function MatrixAlgebraKit.check_input(::typeof(right_orth!), A::LinearMap, VC, alg::AbstractAlgorithm) + return check_input(right_orth!, parent(A), parent.(VC), alg) +end +function MatrixAlgebraKit.default_svd_algorithm(::Type{LinearMap{A}}; kwargs...) where {A} + return default_svd_algorithm(A; kwargs...) +end +function MatrixAlgebraKit.initialize_output( + ::typeof(svd_compact!), A::LinearMap, + alg::GPU_SVDAlgorithm + ) + return LinearMap.(initialize_output(svd_compact!, parent(A), alg)) +end +function MatrixAlgebraKit.svd_compact!(A::LinearMap, USVᴴ, alg::GPU_SVDAlgorithm) + return LinearMap.(svd_compact!(parent(A), parent.(USVᴴ), alg)) +end + +@testset "left_orth and left_null for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) + rng = StableRNG(123) + m = 54 + @testset for n in (37, m, 63) + minmn = min(m, n) + A = ROCArray(randn(rng, T, m, n)) + V, C = @constinferred left_orth(A) + N = @constinferred left_null(A) + @test V isa ROCMatrix{T} && size(V) == (m, minmn) + @test C isa ROCMatrix{T} && size(C) == (minmn, n) + @test N isa ROCMatrix{T} && size(N) == (m, m - minmn) + @test V * C ≈ A + @test isisometric(V) + @test norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + hV = collect(V) + hN = collect(N) + @test hV * hV' + hN * hN' ≈ I + + M = LinearMap(A) + VM, CM = @constinferred left_orth(M; kind = :svd) + @test parent(VM) * parent(CM) ≈ A + + if m > n + nullity = 5 + V, C = @constinferred left_orth(A) + AMDGPU.@allowscalar begin + N = @constinferred left_null(A; trunc = (; maxnullity = nullity)) + end + @test V isa ROCMatrix{T} && size(V) == (m, minmn) + @test C isa ROCMatrix{T} && size(C) == (minmn, n) + @test N isa ROCMatrix{T} && size(N) == (m, nullity) + @test V * C ≈ A + @test isisometric(V) + @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + end + + for alg_qr in ((; positive = true), (; positive = false), ROCSOLVER_HouseholderQR()) + V, C = @constinferred left_orth(A; alg_qr) + N = @constinferred left_null(A; alg_qr) + @test V isa ROCMatrix{T} && size(V) == (m, minmn) + @test C isa ROCMatrix{T} && size(C) == (minmn, n) + @test N isa ROCMatrix{T} && size(N) == (m, m - minmn) + @test V * C ≈ A + @test isisometric(V) + @test norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + hV = collect(V) + hN = collect(N) + @test hV * hV' + hN * hN' ≈ I + end + + Ac = similar(A) + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C)) + N2 = @constinferred left_null!(copy!(Ac, A), N) + @test V2 === V + @test C2 === C + @test N2 === N + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + + atol = eps(real(T)) + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); trunc = (; atol = atol)) + AMDGPU.@allowscalar begin + N2 = @constinferred left_null!(copy!(Ac, A), N; trunc = (; atol = atol)) + end + @test V2 !== V + @test C2 !== C + @test N2 !== C + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + + rtol = eps(real(T)) + for (trunc_orth, trunc_null) in ( + ((; rtol = rtol), (; rtol = rtol)), + (trunctol(; rtol), trunctol(; rtol, keep_below = true)), + ) + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); trunc = trunc_orth) + AMDGPU.@allowscalar begin + N2 = @constinferred left_null!(copy!(Ac, A), N; trunc = trunc_null) + end + @test V2 !== V + @test C2 !== C + @test N2 !== C + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + end + + @testset for kind in (:qr, :polar, :svd) # explicit kind kwarg + m < n && kind == :polar && continue + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); kind = kind) + @test V2 === V + @test C2 === C + @test V2 * C2 ≈ A + @test isisometric(V2) + if kind != :polar + N2 = @constinferred left_null!(copy!(Ac, A), N; kind = kind) + @test N2 === N + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + end + + # with kind and tol kwargs + if kind == :svd + V2, C2 = @constinferred left_orth!( + copy!(Ac, A), (V, C); kind = kind, + trunc = (; atol = atol) + ) + AMDGPU.@allowscalar begin + N2 = @constinferred left_null!( + copy!(Ac, A), N; kind = kind, + trunc = (; atol = atol) + ) + end + @test V2 !== V + @test C2 !== C + @test N2 !== C + @test V2 * C2 ≈ A + @test V2' * V2 ≈ I + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + + V2, C2 = @constinferred left_orth!( + copy!(Ac, A), (V, C); kind = kind, + trunc = (; rtol = rtol) + ) + AMDGPU.@allowscalar begin + N2 = @constinferred left_null!( + copy!(Ac, A), N; kind = kind, + trunc = (; rtol = rtol) + ) + end + @test V2 !== V + @test C2 !== C + @test N2 !== C + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + else + @test_throws ArgumentError left_orth!( + copy!(Ac, A), (V, C); kind = kind, + trunc = (; atol = atol) + ) + @test_throws ArgumentError left_orth!( + copy!(Ac, A), (V, C); kind = kind, + trunc = (; rtol = rtol) + ) + @test_throws ArgumentError left_null!( + copy!(Ac, A), N; kind = kind, + trunc = (; atol = atol) + ) + @test_throws ArgumentError left_null!( + copy!(Ac, A), N; kind = kind, + trunc = (; rtol = rtol) + ) + end + end + end +end + +@testset "right_orth and right_null for T = $T" for T in ( + Float32, Float64, ComplexF32, + ComplexF64, + ) + rng = StableRNG(123) + m = 54 + @testset for n in (37, m, 63) + minmn = min(m, n) + A = ROCArray(randn(rng, T, m, n)) + C, Vᴴ = @constinferred right_orth(A) + Nᴴ = @constinferred right_null(A) + @test C isa ROCMatrix{T} && size(C) == (m, minmn) + @test Vᴴ isa ROCMatrix{T} && size(Vᴴ) == (minmn, n) + @test Nᴴ isa ROCMatrix{T} && size(Nᴴ) == (n - minmn, n) + @test C * Vᴴ ≈ A + @test isisometric(Vᴴ; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ; side = :right) + hVᴴ = collect(Vᴴ) + hNᴴ = collect(Nᴴ) + @test hVᴴ' * hVᴴ + hNᴴ' * hNᴴ ≈ I + + M = LinearMap(A) + CM, VMᴴ = @constinferred right_orth(M; kind = :svd) + @test parent(CM) * parent(VMᴴ) ≈ A + + Ac = similar(A) + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ) + @test C2 === C + @test Vᴴ2 === Vᴴ + @test Nᴴ2 === Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + + atol = eps(real(T)) + rtol = eps(real(T)) + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; atol = atol)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; atol = atol)) + @test C2 !== C + @test Vᴴ2 !== Vᴴ + @test Nᴴ2 !== Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; rtol = rtol)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; rtol = rtol)) + @test C2 !== C + @test Vᴴ2 !== Vᴴ + @test Nᴴ2 !== Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + + @testset "kind = $kind" for kind in (:lq, :polar, :svd) + n < m && kind == :polar && continue + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); kind = kind) + @test C2 === C + @test Vᴴ2 === Vᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + if kind != :polar + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind = kind) + @test Nᴴ2 === Nᴴ + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + end + + if kind == :svd + C2, Vᴴ2 = @constinferred right_orth!( + copy!(Ac, A), (C, Vᴴ); kind = kind, + trunc = (; atol = atol) + ) + Nᴴ2 = @constinferred right_null!( + copy!(Ac, A), Nᴴ; kind = kind, + trunc = (; atol = atol) + ) + @test C2 !== C + @test Vᴴ2 !== Vᴴ + @test Nᴴ2 !== Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + + C2, Vᴴ2 = @constinferred right_orth!( + copy!(Ac, A), (C, Vᴴ); kind = kind, + trunc = (; rtol = rtol) + ) + Nᴴ2 = @constinferred right_null!( + copy!(Ac, A), Nᴴ; kind = kind, + trunc = (; rtol = rtol) + ) + @test C2 !== C + @test Vᴴ2 !== Vᴴ + @test Nᴴ2 !== Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ diagm(ones(T, size(Vᴴ2, 2))) atol = m * n * MatrixAlgebraKit.defaulttol(T) + else + @test_throws ArgumentError right_orth!( + copy!(Ac, A), (C, Vᴴ); kind = kind, + trunc = (; atol = atol) + ) + @test_throws ArgumentError right_orth!( + copy!(Ac, A), (C, Vᴴ); kind = kind, + trunc = (; rtol = rtol) + ) + @test_throws ArgumentError right_null!( + copy!(Ac, A), Nᴴ; kind = kind, + trunc = (; atol = atol) + ) + @test_throws ArgumentError right_null!( + copy!(Ac, A), Nᴴ; kind = kind, + trunc = (; rtol = rtol) + ) + end + end + end +end diff --git a/test/cuda/orthnull.jl b/test/cuda/orthnull.jl new file mode 100644 index 00000000..44ce6388 --- /dev/null +++ b/test/cuda/orthnull.jl @@ -0,0 +1,376 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using LinearAlgebra: LinearAlgebra, I, mul!, diagm, norm +using MatrixAlgebraKit: GPU_SVDAlgorithm, check_input, copy_input, default_svd_algorithm, + initialize_output, AbstractAlgorithm +using CUDA + +# Used to test non-AbstractMatrix codepaths. +struct LinearMap{P <: AbstractMatrix} + parent::P +end +Base.parent(A::LinearMap) = getfield(A, :parent) +function Base.copy!(dest::LinearMap, src::LinearMap) + copy!(parent(dest), parent(src)) + return dest +end +function LinearAlgebra.mul!(C::LinearMap, A::LinearMap, B::LinearMap) + mul!(parent(C), parent(A), parent(B)) + return C +end + +function MatrixAlgebraKit.copy_input(::typeof(qr_compact), A::LinearMap) + return LinearMap(copy_input(qr_compact, parent(A))) +end +function MatrixAlgebraKit.copy_input(::typeof(lq_compact), A::LinearMap) + return LinearMap(copy_input(lq_compact, parent(A))) +end +function MatrixAlgebraKit.initialize_output(::typeof(left_orth!), A::LinearMap) + return LinearMap.(initialize_output(left_orth!, parent(A))) +end +function MatrixAlgebraKit.initialize_output(::typeof(right_orth!), A::LinearMap) + return LinearMap.(initialize_output(right_orth!, parent(A))) +end +function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::LinearMap, VC, alg::AbstractAlgorithm) + return check_input(left_orth!, parent(A), parent.(VC), alg) +end +function MatrixAlgebraKit.check_input(::typeof(right_orth!), A::LinearMap, VC, alg::AbstractAlgorithm) + return check_input(right_orth!, parent(A), parent.(VC), alg) +end +function MatrixAlgebraKit.default_svd_algorithm(::Type{LinearMap{A}}; kwargs...) where {A} + return default_svd_algorithm(A; kwargs...) +end +function MatrixAlgebraKit.initialize_output( + ::typeof(svd_compact!), A::LinearMap, + alg::GPU_SVDAlgorithm + ) + return LinearMap.(initialize_output(svd_compact!, parent(A), alg)) +end +function MatrixAlgebraKit.svd_compact!(A::LinearMap, USVᴴ, alg::GPU_SVDAlgorithm) + return LinearMap.(svd_compact!(parent(A), parent.(USVᴴ), alg)) +end + +@testset "left_orth and left_null for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) + rng = StableRNG(123) + m = 54 + @testset for n in (37, m, 63) + minmn = min(m, n) + A = CuArray(randn(rng, T, m, n)) + V, C = @constinferred left_orth(A) + N = @constinferred left_null(A) + @test V isa CuMatrix{T} && size(V) == (m, minmn) + @test C isa CuMatrix{T} && size(C) == (minmn, n) + @test N isa CuMatrix{T} && size(N) == (m, m - minmn) + @test V * C ≈ A + @test isisometric(V) + @test norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + hV = collect(V) + hN = collect(N) + @test hV * hV' + hN * hN' ≈ I + + M = LinearMap(A) + VM, CM = @constinferred left_orth(M; kind = :svd) + @test parent(VM) * parent(CM) ≈ A + + if m > n + nullity = 5 + V, C = @constinferred left_orth(A) + CUDA.@allowscalar begin + N = @constinferred left_null(A; trunc = (; maxnullity = nullity)) + end + @test V isa CuMatrix{T} && size(V) == (m, minmn) + @test C isa CuMatrix{T} && size(C) == (minmn, n) + @test N isa CuMatrix{T} && size(N) == (m, nullity) + @test V * C ≈ A + @test isisometric(V) + @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + end + + for alg_qr in ((; positive = true), (; positive = false), CUSOLVER_HouseholderQR()) + V, C = @constinferred left_orth(A; alg_qr) + N = @constinferred left_null(A; alg_qr) + @test V isa CuMatrix{T} && size(V) == (m, minmn) + @test C isa CuMatrix{T} && size(C) == (minmn, n) + @test N isa CuMatrix{T} && size(N) == (m, m - minmn) + @test V * C ≈ A + @test isisometric(V) + @test norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + hV = collect(V) + hN = collect(N) + @test hV * hV' + hN * hN' ≈ I + end + + Ac = similar(A) + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C)) + N2 = @constinferred left_null!(copy!(Ac, A), N) + @test V2 === V + @test C2 === C + @test N2 === N + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + + atol = eps(real(T)) + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); trunc = (; atol = atol)) + N2 = @constinferred left_null!(copy!(Ac, A), N; trunc = (; atol = atol)) + @test V2 !== V + @test C2 !== C + @test N2 !== C + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + + rtol = eps(real(T)) + for (trunc_orth, trunc_null) in ( + ((; rtol = rtol), (; rtol = rtol)), + (trunctol(; rtol), trunctol(; rtol, keep_below = true)), + ) + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); trunc = trunc_orth) + N2 = @constinferred left_null!(copy!(Ac, A), N; trunc = trunc_null) + @test V2 !== V + @test C2 !== C + @test N2 !== C + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + end + + @testset for kind in (:qr, :polar, :svd) # explicit kind kwarg + m < n && kind == :polar && continue + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); kind = kind) + @test V2 === V + @test C2 === C + @test V2 * C2 ≈ A + @test isisometric(V2) + if kind != :polar + N2 = @constinferred left_null!(copy!(Ac, A), N; kind = kind) + @test N2 === N + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + end + + # with kind and tol kwargs + if kind == :svd + V2, C2 = @constinferred left_orth!( + copy!(Ac, A), (V, C); kind = kind, + trunc = (; atol = atol) + ) + N2 = @constinferred left_null!( + copy!(Ac, A), N; kind = kind, + trunc = (; atol = atol) + ) + @test V2 !== V + @test C2 !== C + @test N2 !== C + @test V2 * C2 ≈ A + @test V2' * V2 ≈ I + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + + V2, C2 = @constinferred left_orth!( + copy!(Ac, A), (V, C); kind = kind, + trunc = (; rtol = rtol) + ) + N2 = @constinferred left_null!( + copy!(Ac, A), N; kind = kind, + trunc = (; rtol = rtol) + ) + @test V2 !== V + @test C2 !== C + @test N2 !== C + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + hV2 = collect(V2) + hN2 = collect(N2) + @test hV2 * hV2' + hN2 * hN2' ≈ I + else + @test_throws ArgumentError left_orth!( + copy!(Ac, A), (V, C); kind = kind, + trunc = (; atol = atol) + ) + @test_throws ArgumentError left_orth!( + copy!(Ac, A), (V, C); kind = kind, + trunc = (; rtol = rtol) + ) + @test_throws ArgumentError left_null!( + copy!(Ac, A), N; kind = kind, + trunc = (; atol = atol) + ) + @test_throws ArgumentError left_null!( + copy!(Ac, A), N; kind = kind, + trunc = (; rtol = rtol) + ) + end + end + end +end + +@testset "right_orth and right_null for T = $T" for T in ( + Float32, Float64, ComplexF32, + ComplexF64, + ) + rng = StableRNG(123) + m = 54 + @testset for n in (37, m, 63) + minmn = min(m, n) + A = CuArray(randn(rng, T, m, n)) + C, Vᴴ = @constinferred right_orth(A) + Nᴴ = @constinferred right_null(A) + @test C isa CuMatrix{T} && size(C) == (m, minmn) + @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (minmn, n) + @test Nᴴ isa CuMatrix{T} && size(Nᴴ) == (n - minmn, n) + @test C * Vᴴ ≈ A + @test isisometric(Vᴴ; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ; side = :right) + hVᴴ = collect(Vᴴ) + hNᴴ = collect(Nᴴ) + @test hVᴴ' * hVᴴ + hNᴴ' * hNᴴ ≈ I + + M = LinearMap(A) + CM, VMᴴ = @constinferred right_orth(M; kind = :svd) + @test parent(CM) * parent(VMᴴ) ≈ A + + Ac = similar(A) + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ) + @test C2 === C + @test Vᴴ2 === Vᴴ + @test Nᴴ2 === Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + + atol = eps(real(T)) + rtol = eps(real(T)) + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; atol = atol)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; atol = atol)) + @test C2 !== C + @test Vᴴ2 !== Vᴴ + @test Nᴴ2 !== Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; rtol = rtol)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; rtol = rtol)) + @test C2 !== C + @test Vᴴ2 !== Vᴴ + @test Nᴴ2 !== Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + + @testset "kind = $kind" for kind in (:lq, :polar, :svd) + n < m && kind == :polar && continue + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); kind = kind) + @test C2 === C + @test Vᴴ2 === Vᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + if kind != :polar + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind = kind) + @test Nᴴ2 === Nᴴ + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + end + + if kind == :svd + C2, Vᴴ2 = @constinferred right_orth!( + copy!(Ac, A), (C, Vᴴ); kind = kind, + trunc = (; atol = atol) + ) + Nᴴ2 = @constinferred right_null!( + copy!(Ac, A), Nᴴ; kind = kind, + trunc = (; atol = atol) + ) + @test C2 !== C + @test Vᴴ2 !== Vᴴ + @test Nᴴ2 !== Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I + + C2, Vᴴ2 = @constinferred right_orth!( + copy!(Ac, A), (C, Vᴴ); kind = kind, + trunc = (; rtol = rtol) + ) + Nᴴ2 = @constinferred right_null!( + copy!(Ac, A), Nᴴ; kind = kind, + trunc = (; rtol = rtol) + ) + @test C2 !== C + @test Vᴴ2 !== Vᴴ + @test Nᴴ2 !== Nᴴ + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + hVᴴ2 = collect(Vᴴ2) + hNᴴ2 = collect(Nᴴ2) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ diagm(ones(T, size(Vᴴ2, 2))) atol = m * n * MatrixAlgebraKit.defaulttol(T) + else + @test_throws ArgumentError right_orth!( + copy!(Ac, A), (C, Vᴴ); kind = kind, + trunc = (; atol = atol) + ) + @test_throws ArgumentError right_orth!( + copy!(Ac, A), (C, Vᴴ); kind = kind, + trunc = (; rtol = rtol) + ) + @test_throws ArgumentError right_null!( + copy!(Ac, A), Nᴴ; kind = kind, + trunc = (; atol = atol) + ) + @test_throws ArgumentError right_null!( + copy!(Ac, A), Nᴴ; kind = kind, + trunc = (; rtol = rtol) + ) + end + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index c0d72578..af4996ed 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -78,6 +78,9 @@ if CUDA.functional() @safetestset "CUDA Polar Decomposition" begin include("cuda/polar.jl") end + @safetestset "CUDA Image and Null Space" begin + include("cuda/orthnull.jl") + end end using AMDGPU @@ -100,4 +103,7 @@ if AMDGPU.functional() @safetestset "AMDGPU Polar Decomposition" begin include("amd/polar.jl") end + @safetestset "AMDGPU Image and Null Space" begin + include("amd/orthnull.jl") + end end From 3f08b07cf1d646345d17dcb6d519c0870482cf2c Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Thu, 16 Oct 2025 14:58:02 +0200 Subject: [PATCH 2/3] Only use the scalar method for AMDGPU --- .../MatrixAlgebraKitAMDGPUExt.jl | 11 +++++++++++ src/implementations/truncation.jl | 5 +---- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index 8f495f46..e7f03c29 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -161,4 +161,15 @@ function MatrixAlgebraKit._avgdiff!(A::StridedROCMatrix, B::StridedROCMatrix) return A, B end +function MatrixAlgebraKit.truncate(::typeof(MatrixAlgebraKit.left_null!), US::Tuple{TU, TS}, strategy::MatrixAlgebraKit.TruncationStrategy) where {TU <: ROCArray, TS} + # TODO: avoid allocation? + U, S = US + extended_S = vcat(diagview(S), zeros(eltype(S), max(0, size(S, 1) - size(S, 2)))) + ind = MatrixAlgebraKit.findtruncated(extended_S, strategy) + trunc_cols = collect(1:size(U, 2))[ind] + Utrunc = similar(U, (size(U, 1), length(trunc_cols))) + Utrunc .= U[:, trunc_cols] + return Utrunc, ind +end + end diff --git a/src/implementations/truncation.jl b/src/implementations/truncation.jl index ebe52825..f6201b07 100644 --- a/src/implementations/truncation.jl +++ b/src/implementations/truncation.jl @@ -17,10 +17,7 @@ function truncate(::typeof(left_null!), (U, S), strategy::TruncationStrategy) # TODO: avoid allocation? extended_S = vcat(diagview(S), zeros(eltype(S), max(0, size(S, 1) - size(S, 2)))) ind = findtruncated(extended_S, strategy) - trunc_cols = collect(1:size(U, 2))[ind] - Utrunc = similar(U, (size(U, 1), length(trunc_cols))) - Utrunc .= U[:, trunc_cols] - return Utrunc, ind + return U[:, ind], ind end function truncate(::typeof(right_null!), (S, Vᴴ), strategy::TruncationStrategy) # TODO: avoid allocation? From ff7d5208e6b7dc4afe7238a60f8d8bc73e8943b6 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 17 Oct 2025 11:53:41 +0200 Subject: [PATCH 3/3] Comments --- .../MatrixAlgebraKitAMDGPUExt.jl | 3 +- test/amd/orthnull.jl | 46 +---------------- test/cuda/orthnull.jl | 46 +---------------- test/linearmap.jl | 50 ++++++++++++++++++ test/orthnull.jl | 51 ++----------------- 5 files changed, 58 insertions(+), 138 deletions(-) create mode 100644 test/linearmap.jl diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index e7f03c29..48bf56fd 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -167,8 +167,7 @@ function MatrixAlgebraKit.truncate(::typeof(MatrixAlgebraKit.left_null!), US::Tu extended_S = vcat(diagview(S), zeros(eltype(S), max(0, size(S, 1) - size(S, 2)))) ind = MatrixAlgebraKit.findtruncated(extended_S, strategy) trunc_cols = collect(1:size(U, 2))[ind] - Utrunc = similar(U, (size(U, 1), length(trunc_cols))) - Utrunc .= U[:, trunc_cols] + Utrunc = U[:, trunc_cols] return Utrunc, ind end diff --git a/test/amd/orthnull.jl b/test/amd/orthnull.jl index 5f81ffe0..97f5d915 100644 --- a/test/amd/orthnull.jl +++ b/test/amd/orthnull.jl @@ -7,50 +7,8 @@ using MatrixAlgebraKit: GPU_SVDAlgorithm, check_input, copy_input, default_svd_a initialize_output, AbstractAlgorithm using AMDGPU -# Used to test non-AbstractMatrix codepaths. -struct LinearMap{P <: AbstractMatrix} - parent::P -end -Base.parent(A::LinearMap) = getfield(A, :parent) -function Base.copy!(dest::LinearMap, src::LinearMap) - copy!(parent(dest), parent(src)) - return dest -end -function LinearAlgebra.mul!(C::LinearMap, A::LinearMap, B::LinearMap) - mul!(parent(C), parent(A), parent(B)) - return C -end - -function MatrixAlgebraKit.copy_input(::typeof(qr_compact), A::LinearMap) - return LinearMap(copy_input(qr_compact, parent(A))) -end -function MatrixAlgebraKit.copy_input(::typeof(lq_compact), A::LinearMap) - return LinearMap(copy_input(lq_compact, parent(A))) -end -function MatrixAlgebraKit.initialize_output(::typeof(left_orth!), A::LinearMap) - return LinearMap.(initialize_output(left_orth!, parent(A))) -end -function MatrixAlgebraKit.initialize_output(::typeof(right_orth!), A::LinearMap) - return LinearMap.(initialize_output(right_orth!, parent(A))) -end -function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::LinearMap, VC, alg::AbstractAlgorithm) - return check_input(left_orth!, parent(A), parent.(VC), alg) -end -function MatrixAlgebraKit.check_input(::typeof(right_orth!), A::LinearMap, VC, alg::AbstractAlgorithm) - return check_input(right_orth!, parent(A), parent.(VC), alg) -end -function MatrixAlgebraKit.default_svd_algorithm(::Type{LinearMap{A}}; kwargs...) where {A} - return default_svd_algorithm(A; kwargs...) -end -function MatrixAlgebraKit.initialize_output( - ::typeof(svd_compact!), A::LinearMap, - alg::GPU_SVDAlgorithm - ) - return LinearMap.(initialize_output(svd_compact!, parent(A), alg)) -end -function MatrixAlgebraKit.svd_compact!(A::LinearMap, USVᴴ, alg::GPU_SVDAlgorithm) - return LinearMap.(svd_compact!(parent(A), parent.(USVᴴ), alg)) -end +# testing non-AbstractArray codepaths: +include(joinpath("..", "linearmap.jl")) @testset "left_orth and left_null for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) rng = StableRNG(123) diff --git a/test/cuda/orthnull.jl b/test/cuda/orthnull.jl index 44ce6388..cfe7d9e3 100644 --- a/test/cuda/orthnull.jl +++ b/test/cuda/orthnull.jl @@ -7,50 +7,8 @@ using MatrixAlgebraKit: GPU_SVDAlgorithm, check_input, copy_input, default_svd_a initialize_output, AbstractAlgorithm using CUDA -# Used to test non-AbstractMatrix codepaths. -struct LinearMap{P <: AbstractMatrix} - parent::P -end -Base.parent(A::LinearMap) = getfield(A, :parent) -function Base.copy!(dest::LinearMap, src::LinearMap) - copy!(parent(dest), parent(src)) - return dest -end -function LinearAlgebra.mul!(C::LinearMap, A::LinearMap, B::LinearMap) - mul!(parent(C), parent(A), parent(B)) - return C -end - -function MatrixAlgebraKit.copy_input(::typeof(qr_compact), A::LinearMap) - return LinearMap(copy_input(qr_compact, parent(A))) -end -function MatrixAlgebraKit.copy_input(::typeof(lq_compact), A::LinearMap) - return LinearMap(copy_input(lq_compact, parent(A))) -end -function MatrixAlgebraKit.initialize_output(::typeof(left_orth!), A::LinearMap) - return LinearMap.(initialize_output(left_orth!, parent(A))) -end -function MatrixAlgebraKit.initialize_output(::typeof(right_orth!), A::LinearMap) - return LinearMap.(initialize_output(right_orth!, parent(A))) -end -function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::LinearMap, VC, alg::AbstractAlgorithm) - return check_input(left_orth!, parent(A), parent.(VC), alg) -end -function MatrixAlgebraKit.check_input(::typeof(right_orth!), A::LinearMap, VC, alg::AbstractAlgorithm) - return check_input(right_orth!, parent(A), parent.(VC), alg) -end -function MatrixAlgebraKit.default_svd_algorithm(::Type{LinearMap{A}}; kwargs...) where {A} - return default_svd_algorithm(A; kwargs...) -end -function MatrixAlgebraKit.initialize_output( - ::typeof(svd_compact!), A::LinearMap, - alg::GPU_SVDAlgorithm - ) - return LinearMap.(initialize_output(svd_compact!, parent(A), alg)) -end -function MatrixAlgebraKit.svd_compact!(A::LinearMap, USVᴴ, alg::GPU_SVDAlgorithm) - return LinearMap.(svd_compact!(parent(A), parent.(USVᴴ), alg)) -end +# testing non-AbstractArray codepaths: +include(joinpath("..", "linearmap.jl")) @testset "left_orth and left_null for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) rng = StableRNG(123) diff --git a/test/linearmap.jl b/test/linearmap.jl new file mode 100644 index 00000000..61753cde --- /dev/null +++ b/test/linearmap.jl @@ -0,0 +1,50 @@ +module LinearMaps + + export LinearMap + + using MatrixAlgebraKit + using MatrixAlgebraKit: AbstractAlgorithm + import MatrixAlgebraKit as MAK + + using LinearAlgebra: LinearAlgebra, lmul!, rmul! + + # Used to test non-AbstractMatrix codepaths. + struct LinearMap{P <: AbstractMatrix} + parent::P + end + Base.parent(A::LinearMap) = A.parent + + Base.copy!(dest::LinearMap, src::LinearMap) = (copy!(parent(dest), parent(src)); dest) + + # necessary for orth_svd default implementations + LinearAlgebra.lmul!(D::LinearMap, A::LinearMap) = (lmul!(parent(D), parent(A)); A) + LinearAlgebra.rmul!(A::LinearMap, D::LinearMap) = (rmul!(parent(A), parent(D)); A) + LinearAlgebra.mul!(C::LinearMap, A::LinearMap, B::LinearMap, α, β) = LinearAlgebra.mul!(parent(C), parent(A), parent(B), α, β) + + for f in (:qr_compact, :lq_compact) + @eval MAK.copy_input(::typeof($f), A::LinearMap) = LinearMap(MAK.copy_input($f, parent(A))) + end + + for f! in (:qr_compact!, :lq_compact!, :svd_compact!, :svd_full!, :svd_trunc!) + @eval MAK.check_input(::typeof($f!), A::LinearMap, F, alg::AbstractAlgorithm) = + MAK.check_input($f!, parent(A), parent.(F), alg) + @eval MAK.initialize_output(::typeof($f!), A::LinearMap, alg::AbstractAlgorithm) = + LinearMap.(MAK.initialize_output($f!, parent(A), alg)) + @eval MAK.$f!(A::LinearMap, F, alg::AbstractAlgorithm) = + LinearMap.(MAK.$f!(parent(A), parent.(F), alg)) + end + + for f! in (:left_orth!, :right_orth!) + @eval MAK.check_input(::typeof($f!), A::LinearMap, F, alg) = + MAK.check_input($f!, parent(A), parent.(F), alg) + @eval MAK.initialize_output(::typeof($f!), A::LinearMap) = + LinearMap.(MAK.initialize_output($f!, parent(A))) + end + + for f in (:qr, :lq, :svd) + default_f = Symbol(:default_, f, :_algorithm) + @eval MAK.$default_f(::Type{LinearMap{A}}; kwargs...) where {A} = MAK.$default_f(A; kwargs...) + end +end + +using .LinearMaps diff --git a/test/orthnull.jl b/test/orthnull.jl index fc364916..cee54252 100644 --- a/test/orthnull.jl +++ b/test/orthnull.jl @@ -6,55 +6,10 @@ using LinearAlgebra: LinearAlgebra, I, mul! using MatrixAlgebraKit: LAPACK_SVDAlgorithm, check_input, copy_input, default_svd_algorithm, initialize_output, AbstractAlgorithm -eltypes = (Float32, Float64, ComplexF32, ComplexF64) - -# Used to test non-AbstractMatrix codepaths. -struct LinearMap{P <: AbstractMatrix} - parent::P -end -Base.parent(A::LinearMap) = getfield(A, :parent) -function Base.copy!(dest::LinearMap, src::LinearMap) - copy!(parent(dest), parent(src)) - return dest -end -function LinearAlgebra.mul!(C::LinearMap, A::LinearMap, B::LinearMap) - mul!(parent(C), parent(A), parent(B)) - return C -end +# testing non-AbstractArray codepaths: +include("linearmap.jl") -function MatrixAlgebraKit.copy_input(::typeof(qr_compact), A::LinearMap) - return LinearMap(copy_input(qr_compact, parent(A))) -end -function MatrixAlgebraKit.copy_input(::typeof(lq_compact), A::LinearMap) - return LinearMap(copy_input(lq_compact, parent(A))) -end -function MatrixAlgebraKit.initialize_output(::typeof(left_orth!), A::LinearMap) - return LinearMap.(initialize_output(left_orth!, parent(A))) -end -function MatrixAlgebraKit.initialize_output(::typeof(right_orth!), A::LinearMap) - return LinearMap.(initialize_output(right_orth!, parent(A))) -end -function MatrixAlgebraKit.check_input( - ::typeof(left_orth!), A::LinearMap, VC, alg::AbstractAlgorithm - ) - return check_input(left_orth!, parent(A), parent.(VC), alg) -end -function MatrixAlgebraKit.check_input( - ::typeof(right_orth!), A::LinearMap, VC, alg::AbstractAlgorithm - ) - return check_input(right_orth!, parent(A), parent.(VC), alg) -end -function MatrixAlgebraKit.default_svd_algorithm(::Type{LinearMap{A}}; kwargs...) where {A} - return default_svd_algorithm(A; kwargs...) -end -function MatrixAlgebraKit.initialize_output( - ::typeof(svd_compact!), A::LinearMap, alg::LAPACK_SVDAlgorithm - ) - return LinearMap.(initialize_output(svd_compact!, parent(A), alg)) -end -function MatrixAlgebraKit.svd_compact!(A::LinearMap, USVᴴ, alg::LAPACK_SVDAlgorithm) - return LinearMap.(svd_compact!(parent(A), parent.(USVᴴ), alg)) -end +eltypes = (Float32, Float64, ComplexF32, ComplexF64) @testset "left_orth and left_null for T = $T" for T in eltypes rng = StableRNG(123)