diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index 8f495f46..48bf56fd 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -161,4 +161,14 @@ 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 = U[:, trunc_cols] + return Utrunc, ind +end + end diff --git a/test/amd/orthnull.jl b/test/amd/orthnull.jl new file mode 100644 index 00000000..97f5d915 --- /dev/null +++ b/test/amd/orthnull.jl @@ -0,0 +1,342 @@ +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 + +# 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) + 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..cfe7d9e3 --- /dev/null +++ b/test/cuda/orthnull.jl @@ -0,0 +1,334 @@ +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 + +# 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) + 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/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) 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