diff --git a/docs/src/user_interface/decompositions.md b/docs/src/user_interface/decompositions.md index b6d0e94e..1ae1fefe 100644 --- a/docs/src/user_interface/decompositions.md +++ b/docs/src/user_interface/decompositions.md @@ -169,23 +169,205 @@ PolarNewton ## Orthogonal Subspaces -Often it is useful to compute orthogonal bases for a particular subspace defined by a matrix. -Given a matrix `A` we can compute an orthonormal basis for its image or coimage, and factorize the matrix accordingly. +Often it is useful to compute orthogonal bases for particular subspaces defined by a matrix. +Given a matrix `A`, we can compute an orthonormal basis for its image or coimage, and factorize the matrix accordingly. These bases are accessible through [`left_orth`](@ref) and [`right_orth`](@ref) respectively. -This is implemented through a combination of the decompositions mentioned above, and serves as a convenient interface to these operations. + +### Overview + +The [`left_orth`](@ref) function computes an orthonormal basis `V` for the image (column space) of `A`, along with a corestriction matrix `C` such that `A = V * C`. +The resulting `V` has orthonormal columns (`V' * V ≈ I` or `isisometric(V)`). + +Similarly, [`right_orth`](@ref) computes an orthonormal basis for the coimage (row space) of `A`, i.e., the image of `A'`. +It returns matrices `C` and `Vᴴ` such that `A = C * Vᴴ`, where `V = (Vᴴ)'` has orthonormal columns (`isisometric(Vᴴ; side = :right)`). + +These functions serve as high-level interfaces that automatically select the most appropriate decomposition based on the specified options, making them convenient for users who want orthonormalization without worrying about the underlying implementation details. ```@docs; canonical=false left_orth right_orth ``` +### Algorithm Selection + +Both functions support multiple decomposition drivers, which can be selected through the `alg` keyword argument: + +**For `left_orth`:** +- `alg = :qr` (default without truncation): Uses QR decomposition via [`qr_compact`](@ref) +- `alg = :polar`: Uses polar decomposition via [`left_polar`](@ref) +- `alg = :svd` (default with truncation): Uses SVD via [`svd_compact`](@ref) or [`svd_trunc`](@ref) + +**For `right_orth`:** +- `alg = :lq` (default without truncation): Uses LQ decomposition via [`lq_compact`](@ref) +- `alg = :polar`: Uses polar decomposition via [`right_polar`](@ref) +- `alg = :svd` (default with truncation): Uses SVD via [`svd_compact`](@ref) or [`svd_trunc`](@ref) + +When `alg` is not specified, the function automatically selects `:qr`/`:lq` for exact orthogonalization, or `:svd` when a truncation strategy is provided. + +### Extending with Custom Algorithms + +To register a custom algorithm type for use with these functions, you need to define the appropriate conversion function, for example: + +```julia +# For left_orth +MatrixAlgebraKit.left_orth_alg(alg::MyCustomAlgorithm) = LeftOrthAlgorithm{:qr}(alg) + +# For right_orth +MatrixAlgebraKit.right_orth_alg(alg::MyCustomAlgorithm) = RightOrthAlgorithm{:lq}(alg) +``` + +The type parameter (`:qr`, `:lq`, `:polar`, or `:svd`) indicates which factorization backend will be used. +The wrapper algorithm types handle the dispatch to the appropriate implementation: + +```@docs; canonical=false +left_orth_alg +right_orth_alg +LeftOrthAlgorithm +RightOrthAlgorithm +``` + +### Examples + +Basic orthogonalization: + +```jldoctest orthnull; output=false +using MatrixAlgebraKit +using LinearAlgebra + +A = [1.0 2.0; 3.0 4.0; 5.0 6.0] +V, C = left_orth(A) +(V' * V) ≈ I && A ≈ V * C + +# output +true +``` + +Using different algorithms: + +```jldoctest orthnull; output=false +A = randn(4, 3) +V1, C1 = left_orth(A; alg = :qr) +V2, C2 = left_orth(A; alg = :polar) +V3, C3 = left_orth(A; alg = :svd) +A ≈ V1 * C1 ≈ V2 * C2 ≈ V3 * C3 + +# output +true +``` + +With truncation: + +```jldoctest orthnull; output=false +A = [1.0 0.0; 0.0 1e-10; 0.0 0.0] +V, C = left_orth(A; trunc = (atol = 1e-8,)) +size(V, 2) == 1 # Only one column retained + +# output +true +``` + + ## Null Spaces Similarly, it can be convenient to obtain an orthogonal basis for the kernel or cokernel of a matrix. -These are the compliments of the coimage and image, respectively, and can be computed using the [`left_null`](@ref) and [`right_null`](@ref) functions. -Again, this is typically implemented through a combination of the decompositions mentioned above, and serves as a convenient interface to these operations. +These are the orthogonal complements of the coimage and image, respectively, and can be computed using the [`left_null`](@ref) and [`right_null`](@ref) functions. + +### Overview + +The [`left_null`](@ref) function computes an orthonormal basis `N` for the cokernel (left nullspace) of `A`, which is the nullspace of `A'`. +This means `A' * N ≈ 0` and `N' * N ≈ I`. + +Similarly, [`right_null`](@ref) computes an orthonormal basis for the kernel (right nullspace) of `A`. +It returns `Nᴴ` such that `A * Nᴴ' ≈ 0` and `Nᴴ * Nᴴ' ≈ I`, where `N = (Nᴴ)'` has orthonormal columns. + +These functions automatically handle rank determination and provide convenient access to nullspace computation without requiring detailed knowledge of the underlying decomposition methods. ```@docs; canonical=false left_null right_null ``` + +### Algorithm Selection + +Both functions support multiple decomposition drivers, which can be selected through the `alg` keyword argument: + +**For `left_null`:** +- `alg = :qr` (default without truncation): Uses QR-based nullspace computation via [`qr_null`](@ref) +- `alg = :svd` (default with truncation): Uses SVD via [`svd_full`](@ref) with appropriate truncation + +**For `right_null`:** +- `alg = :lq` (default without truncation): Uses LQ-based nullspace computation via [`lq_null`](@ref) +- `alg = :svd` (default with truncation): Uses SVD via [`svd_full`](@ref) with appropriate truncation + +When `alg` is not specified, the function automatically selects `:qr`/`:lq` for exact nullspace computation, or `:svd` when a truncation strategy is provided to handle numerical rank determination. + +!!! note + For nullspace functions, [`notrunc`](@ref) has special meaning when used with the default QR/LQ algorithms. + It indicates that the nullspace should be computed from the exact zeros determined by the additional rows/columns of the extended matrix, without any tolerance-based truncation. + +### Extending with Custom Algorithms + +To register a custom algorithm type for use with these functions, you need to define the appropriate conversion function: + +```julia +# For left_null +MatrixAlgebraKit.left_null_alg(alg::MyCustomAlgorithm) = LeftNullAlgorithm{:qr}(alg) + +# For right_null +MatrixAlgebraKit.right_null_alg(alg::MyCustomAlgorithm) = RightNullAlgorithm{:lq}(alg) +``` + +The type parameter (`:qr`, `:lq`, or `:svd`) indicates which factorization backend will be used. +The wrapper algorithm types handle the dispatch to the appropriate implementation: + +```@docs; canonical=false +LeftNullAlgorithm +RightNullAlgorithm +left_null_alg +right_null_alg +``` + +### Examples + +Basic nullspace computation: + +```jldoctest orthnull; output=false +A = [1.0 2.0 3.0; 4.0 5.0 6.0] # Rank 2 matrix +N = left_null(A) +size(N) == (2, 0) + +# output +true +``` + +```jldoctest orthnull; output=false +Nᴴ = right_null(A) +size(Nᴴ) == (1, 3) && norm(A * Nᴴ') < 1e-14 && isisometric(Nᴴ; side = :right) + +# output +true +``` + +Computing nullspace with rank detection: + +```jldoctest orthnull; output=false +A = [1.0 2.0; 2.0 4.0; 3.0 6.0] # Rank 1 matrix (second column = 2*first) +N = left_null(A; alg = :svd, trunc = (atol = 1e-10,)) +size(N) == (3, 2) && norm(A' * N) < 1e-12 && isisometric(N) + +# output +true +``` + +Using different algorithms: + +```jldoctest orthnull; output=false +A = [1.0 0.0 0.0; 0.0 1.0 0.0] +N1 = right_null(A; alg = :lq) +N2 = right_null(A; alg = :svd) +norm(A * N1') < 1e-14 && norm(A * N2') < 1e-14 && + isisometric(N1; side = :right) && isisometric(N2; side = :right) + +# output +true +``` diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index 48bf56fd..04bee40e 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -4,7 +4,7 @@ using MatrixAlgebraKit using MatrixAlgebraKit: @algdef, Algorithm, check_input using MatrixAlgebraKit: one!, zero!, uppertriangular!, lowertriangular! using MatrixAlgebraKit: diagview, sign_safe -using MatrixAlgebraKit: LQViaTransposedQR, TruncationByValue, AbstractAlgorithm +using MatrixAlgebraKit: LQViaTransposedQR, TruncationStrategy, NoTruncation, TruncationByValue, AbstractAlgorithm using MatrixAlgebraKit: default_qr_algorithm, default_lq_algorithm, default_svd_algorithm, default_eigh_algorithm import MatrixAlgebraKit: _gpu_geqrf!, _gpu_ungqr!, _gpu_unmqr!, _gpu_gesvd!, _gpu_Xgesvdp!, _gpu_gesvdj! import MatrixAlgebraKit: _gpu_heevj!, _gpu_heevd!, _gpu_heev!, _gpu_heevx! @@ -161,7 +161,9 @@ 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} +function MatrixAlgebraKit.truncate( + ::typeof(left_null!), US::Tuple{TU, TS}, strategy::TruncationStrategy + ) where {TU <: ROCMatrix, TS} # TODO: avoid allocation? U, S = US extended_S = vcat(diagview(S), zeros(eltype(S), max(0, size(S, 1) - size(S, 2)))) @@ -170,5 +172,32 @@ function MatrixAlgebraKit.truncate(::typeof(MatrixAlgebraKit.left_null!), US::Tu Utrunc = U[:, trunc_cols] return Utrunc, ind end +function MatrixAlgebraKit.truncate( + ::typeof(right_null!), SVᴴ::Tuple{TS, TVᴴ}, strategy::TruncationStrategy + ) where {TS, TVᴴ <: ROCMatrix} + # TODO: avoid allocation? + S, Vᴴ = SVᴴ + extended_S = vcat(diagview(S), zeros(eltype(S), max(0, size(S, 2) - size(S, 1)))) + ind = MatrixAlgebraKit.findtruncated(extended_S, strategy) + trunc_rows = collect(1:size(Vᴴ, 1))[ind] + Vᴴtrunc = Vᴴ[trunc_rows, :] + return Vᴴtrunc, ind +end + +# disambiguate: +function MatrixAlgebraKit.truncate( + ::typeof(left_null!), (U, S)::Tuple{TU, TS}, ::NoTruncation + ) where {TU <: ROCMatrix, TS} + m, n = size(S) + ind = (n + 1):m + return U[:, ind], ind +end +function MatrixAlgebraKit.truncate( + ::typeof(right_null!), (S, Vᴴ)::Tuple{TS, TVᴴ}, ::NoTruncation + ) where {TS, TVᴴ <: ROCMatrix} + m, n = size(S) + ind = (m + 1):n + return Vᴴ[ind, :], ind +end end diff --git a/src/algorithms.jl b/src/algorithms.jl index a47be2c8..e9e7b8e8 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -53,6 +53,18 @@ function _show_alg(io::IO, alg::Algorithm) return print(io, ")") end +# Algorithm traits +# ---------------- +""" + does_truncate(alg::AbstractAlgorithm) -> Bool + +Indicate whether or not an algorithm will compute a truncated decomposition +(such that composing the factors only approximates the input up to some tolerance). +""" +does_truncate(alg::AbstractAlgorithm) = false + +# Algorithm selection +# ------------------- @doc """ MatrixAlgebraKit.select_algorithm(f, A, alg::AbstractAlgorithm) MatrixAlgebraKit.select_algorithm(f, A, alg::Symbol; kwargs...) @@ -83,7 +95,7 @@ function select_algorithm(f::F, A, alg::Alg = nothing; kwargs...) where {F, Alg} return Algorithm{alg}(; kwargs...) elseif alg isa Type return alg(; kwargs...) - elseif alg isa NamedTuple + elseif alg isa NamedTuple || alg isa Base.Pairs isempty(kwargs) || throw(ArgumentError("Additional keyword arguments are not allowed when algorithm parameters are specified.")) return default_algorithm(f, A; alg...) @@ -160,6 +172,24 @@ function select_truncation(trunc) end end +@doc """ + MatrixAlgebraKit.select_null_truncation(trunc) + +Construct a [`TruncationStrategy`](@ref) from the given `NamedTuple` of keywords or input strategy, to implement a nullspace selection. +""" select_null_truncation + +function select_null_truncation(trunc) + if isnothing(trunc) + return NoTruncation() + elseif trunc isa NamedTuple + return null_truncation_strategy(; trunc...) + elseif trunc isa TruncationStrategy + return trunc + else + return throw(ArgumentError("Unknown truncation strategy: $trunc")) + end +end + @doc """ MatrixAlgebraKit.findtruncated(values::AbstractVector, strategy::TruncationStrategy) @@ -200,6 +230,8 @@ struct TruncatedAlgorithm{A, T} <: AbstractAlgorithm trunc::T end +does_truncate(::TruncatedAlgorithm) = true + # Utility macros # -------------- @@ -230,14 +262,14 @@ end function _arg_expr(::Val{1}, f, f!) return quote # out of place to inplace - $f(A; kwargs...) = $f!(copy_input($f, A); kwargs...) + @inline $f(A; alg = nothing, kwargs...) = $f(A, select_algorithm($f, A, alg; kwargs...)) $f(A, alg::AbstractAlgorithm) = $f!(copy_input($f, A), alg) # fill in arguments - function $f!(A; alg = nothing, kwargs...) + @inline function $f!(A; alg = nothing, kwargs...) return $f!(A, select_algorithm($f!, A, alg; kwargs...)) end - function $f!(A, out; alg = nothing, kwargs...) + @inline function $f!(A, out; alg = nothing, kwargs...) return $f!(A, out, select_algorithm($f!, A, alg; kwargs...)) end function $f!(A, alg::AbstractAlgorithm) diff --git a/src/implementations/orthnull.jl b/src/implementations/orthnull.jl index 76e32ce1..2a15259b 100644 --- a/src/implementations/orthnull.jl +++ b/src/implementations/orthnull.jl @@ -5,287 +5,108 @@ copy_input(::typeof(right_orth), A) = copy_input(lq_compact, A) # do we ever nee copy_input(::typeof(left_null), A) = copy_input(qr_null, A) # do we ever need anything else copy_input(::typeof(right_null), A) = copy_input(lq_null, A) # do we ever need anything else -function check_input(::typeof(left_orth!), A::AbstractMatrix, VC, ::AbstractAlgorithm) - m, n = size(A) - minmn = min(m, n) - V, C = VC - @assert V isa AbstractMatrix && C isa AbstractMatrix - @check_size(V, (m, minmn)) - @check_scalar(V, A) - if !isempty(C) - @check_size(C, (minmn, n)) - @check_scalar(C, A) - end - return nothing -end -function check_input(::typeof(right_orth!), A::AbstractMatrix, CVᴴ, ::AbstractAlgorithm) - m, n = size(A) - minmn = min(m, n) - C, Vᴴ = CVᴴ - @assert C isa AbstractMatrix && Vᴴ isa AbstractMatrix - if !isempty(C) - @check_size(C, (m, minmn)) - @check_scalar(C, A) - end - @check_size(Vᴴ, (minmn, n)) - @check_scalar(Vᴴ, A) - return nothing -end +check_input(::typeof(left_orth!), A, VC, alg::AbstractAlgorithm) = + check_input(left_orth!, A, VC, left_orth_alg(alg)) +check_input(::typeof(left_orth!), A, VC, alg::LeftOrthViaQR) = + check_input(qr_compact!, A, VC, alg.alg) +check_input(::typeof(left_orth!), A, VC, alg::LeftOrthViaPolar) = + check_input(left_polar!, A, VC, alg.alg) +check_input(::typeof(left_orth!), A, VC, alg::LeftOrthViaSVD) = nothing -function check_input(::typeof(left_null!), A::AbstractMatrix, N, ::AbstractAlgorithm) - m, n = size(A) - minmn = min(m, n) - @assert N isa AbstractMatrix - @check_size(N, (m, m - minmn)) - @check_scalar(N, A) - return nothing -end -function check_input(::typeof(right_null!), A::AbstractMatrix, Nᴴ, ::AbstractAlgorithm) - m, n = size(A) - minmn = min(m, n) - @assert Nᴴ isa AbstractMatrix - @check_size(Nᴴ, (n - minmn, n)) - @check_scalar(Nᴴ, A) - return nothing -end +check_input(::typeof(right_orth!), A, CVᴴ, alg::AbstractAlgorithm) = + check_input(right_orth!, A, CVᴴ, right_orth_alg(alg)) +check_input(::typeof(right_orth!), A, VC, alg::RightOrthViaLQ) = + check_input(lq_compact!, A, VC, alg.alg) +check_input(::typeof(right_orth!), A, VC, alg::RightOrthViaPolar) = + check_input(right_polar!, A, VC, alg.alg) +check_input(::typeof(right_orth!), A, VC, alg::RightOrthViaSVD) = nothing + +check_input(::typeof(left_null!), A, N, alg::AbstractAlgorithm) = + check_input(left_null!, A, N, left_null_alg(alg)) +check_input(::typeof(left_null!), A, N, alg::LeftNullViaQR) = + check_input(qr_null!, A, N, alg.alg) +check_input(::typeof(left_null!), A, N, alg::LeftNullViaSVD) = nothing + +check_input(::typeof(right_null!), A, Nᴴ, alg::AbstractAlgorithm) = + check_input(right_null!, A, Nᴴ, right_null_alg(alg)) +check_input(::typeof(right_null!), A, Nᴴ, alg::RightNullViaLQ) = + check_input(lq_null!, A, Nᴴ, alg.alg) +check_input(::typeof(right_null!), A, Nᴴ, alg::RightNullViaSVD) = nothing # Outputs # ------- -function initialize_output(::typeof(left_orth!), A::AbstractMatrix) - m, n = size(A) - minmn = min(m, n) - V = similar(A, (m, minmn)) - C = similar(A, (minmn, n)) - return (V, C) -end -function initialize_output(::typeof(right_orth!), A::AbstractMatrix) - m, n = size(A) - minmn = min(m, n) - C = similar(A, (m, minmn)) - Vᴴ = similar(A, (minmn, n)) - return (C, Vᴴ) -end +initialize_output(::typeof(left_orth!), A, alg::AbstractAlgorithm) = + initialize_output(left_orth!, A, left_orth_alg(alg)) +initialize_output(::typeof(left_orth!), A, alg::LeftOrthViaQR) = + initialize_output(qr_compact!, A, alg.alg) +initialize_output(::typeof(left_orth!), A, alg::LeftOrthViaPolar) = + initialize_output(left_polar!, A, alg.alg) +initialize_output(::typeof(left_orth!), A, alg::LeftOrthViaSVD) = nothing -function initialize_output(::typeof(left_null!), A::AbstractMatrix) - m, n = size(A) - minmn = min(m, n) - N = similar(A, (m, m - minmn)) - return N -end -function initialize_output(::typeof(right_null!), A::AbstractMatrix) - m, n = size(A) - minmn = min(m, n) - Nᴴ = similar(A, (n - minmn, n)) - return Nᴴ -end +initialize_output(::typeof(right_orth!), A, alg::AbstractAlgorithm) = + initialize_output(right_orth!, A, right_orth_alg(alg)) +initialize_output(::typeof(right_orth!), A, alg::RightOrthViaLQ) = + initialize_output(lq_compact!, A, alg.alg) +initialize_output(::typeof(right_orth!), A, alg::RightOrthViaPolar) = + initialize_output(right_polar!, A, alg.alg) +initialize_output(::typeof(right_orth!), A, alg::RightOrthViaSVD) = nothing + +initialize_output(::typeof(left_null!), A, alg::AbstractAlgorithm) = + initialize_output(left_null!, A, left_null_alg(alg)) +initialize_output(::typeof(left_null!), A, alg::LeftNullViaQR) = + initialize_output(qr_null!, A, alg.alg) +initialize_output(::typeof(left_null!), A, alg::LeftNullViaSVD) = nothing + +initialize_output(::typeof(right_null!), A, alg::AbstractAlgorithm) = + initialize_output(right_null!, A, right_null_alg(alg)) +initialize_output(::typeof(right_null!), A, alg::RightNullViaLQ) = + initialize_output(lq_null!, A, alg.alg) +initialize_output(::typeof(right_null!), A, alg::RightNullViaSVD) = nothing # Implementation of orth functions # -------------------------------- -function left_orth!( - A, VC; - trunc = nothing, kind = isnothing(trunc) ? :qr : :svd, - alg_qr = (; positive = true), alg_polar = (;), alg_svd = (;) - ) - if !isnothing(trunc) && kind != :svd - throw(ArgumentError("truncation not supported for left_orth with kind=$kind")) - end - if kind == :qr - return left_orth_qr!(A, VC, alg_qr) - elseif kind == :polar - return left_orth_polar!(A, VC, alg_polar) - elseif kind == :svd - return left_orth_svd!(A, VC, alg_svd, trunc) - else - throw(ArgumentError("`left_orth!` received unknown value `kind = $kind`")) - end -end -function left_orth_qr!(A, VC, alg) - alg′ = select_algorithm(qr_compact!, A, alg) - check_input(left_orth!, A, VC, alg′) - return qr_compact!(A, VC, alg′) -end -function left_orth_polar!(A, VC, alg) - alg′ = select_algorithm(left_polar!, A, alg) - check_input(left_orth!, A, VC, alg′) - return left_polar!(A, VC, alg′) -end -function left_orth_svd!(A, VC, alg, trunc::Nothing = nothing) - alg′ = select_algorithm(svd_compact!, A, alg) - check_input(left_orth!, A, VC, alg′) - U, S, Vᴴ = svd_compact!(A, alg′) - V, C = VC - return copy!(V, U), mul!(C, S, Vᴴ) -end -function left_orth_svd!(A::AbstractMatrix, VC, alg, trunc::Nothing = nothing) - alg′ = select_algorithm(svd_compact!, A, alg) - check_input(left_orth!, A, VC, alg′) - V, C = VC - S = Diagonal(initialize_output(svd_vals!, A, alg′)) - U, S, Vᴴ = svd_compact!(A, (V, S, C), alg′) - return U, lmul!(S, Vᴴ) -end -function left_orth_svd!(A, VC, alg, trunc) - alg′ = select_algorithm(svd_compact!, A, alg) - check_input(left_orth!, A, VC, alg′) - alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc) - U, S, Vᴴ = svd_trunc!(A, alg_trunc) - V, C = VC - return copy!(V, U), mul!(C, S, Vᴴ) -end -function left_orth_svd!(A::AbstractMatrix, VC, alg, trunc) - alg′ = select_algorithm(svd_compact!, A, alg) - check_input(left_orth!, A, VC, alg′) - alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc) - V, C = VC - S = Diagonal(initialize_output(svd_vals!, A, alg_trunc.alg)) - U, S, Vᴴ = svd_trunc!(A, (V, S, C), alg_trunc) - return U, lmul!(S, Vᴴ) +left_orth!(A, VC, alg::AbstractAlgorithm) = left_orth!(A, VC, left_orth_alg(alg)) +left_orth!(A, VC, alg::LeftOrthViaQR) = qr_compact!(A, VC, alg.alg) +left_orth!(A, VC, alg::LeftOrthViaPolar) = left_polar!(A, VC, alg.alg) +function left_orth!(A, VC, alg::LeftOrthViaSVD) + check_input(left_orth!, A, VC, alg) + V, S, C = does_truncate(alg.alg) ? svd_trunc!(A, alg.alg) : svd_compact!(A, alg.alg) + lmul!(S, C) + return V, C end -function right_orth!( - A, CVᴴ; - trunc = nothing, kind = isnothing(trunc) ? :lq : :svd, - alg_lq = (; positive = true), alg_polar = (;), alg_svd = (;) - ) - if !isnothing(trunc) && kind != :svd - throw(ArgumentError("truncation not supported for right_orth with kind=$kind")) - end - if kind == :lq - return right_orth_lq!(A, CVᴴ, alg_lq) - elseif kind == :polar - return right_orth_polar!(A, CVᴴ, alg_polar) - elseif kind == :svd - return right_orth_svd!(A, CVᴴ, alg_svd, trunc) - else - throw(ArgumentError("`right_orth!` received unknown value `kind = $kind`")) - end -end -function right_orth_lq!(A, CVᴴ, alg) - alg′ = select_algorithm(lq_compact!, A, alg) - check_input(right_orth!, A, CVᴴ, alg′) - return lq_compact!(A, CVᴴ, alg′) -end -function right_orth_polar!(A, CVᴴ, alg) - alg′ = select_algorithm(right_polar!, A, alg) - check_input(right_orth!, A, CVᴴ, alg′) - return right_polar!(A, CVᴴ, alg′) -end -function right_orth_svd!(A, CVᴴ, alg, trunc::Nothing = nothing) - alg′ = select_algorithm(svd_compact!, A, alg) - check_input(right_orth!, A, CVᴴ, alg′) - U, S, Vᴴ′ = svd_compact!(A, alg′) - C, Vᴴ = CVᴴ - return mul!(C, U, S), copy!(Vᴴ, Vᴴ′) -end -function right_orth_svd!(A::AbstractMatrix, CVᴴ, alg, trunc::Nothing = nothing) - alg′ = select_algorithm(svd_compact!, A, alg) - check_input(right_orth!, A, CVᴴ, alg′) - C, Vᴴ = CVᴴ - S = Diagonal(initialize_output(svd_vals!, A, alg′)) - U, S, Vᴴ = svd_compact!(A, (C, S, Vᴴ), alg′) - return rmul!(U, S), Vᴴ -end -function right_orth_svd!(A, CVᴴ, alg, trunc) - alg′ = select_algorithm(svd_compact!, A, alg) - check_input(right_orth!, A, CVᴴ, alg′) - alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc) - U, S, Vᴴ′ = svd_trunc!(A, alg_trunc) - C, Vᴴ = CVᴴ - return mul!(C, U, S), copy!(Vᴴ, Vᴴ′) -end -function right_orth_svd!(A::AbstractMatrix, CVᴴ, alg, trunc) - alg′ = select_algorithm(svd_compact!, A, alg) - check_input(right_orth!, A, CVᴴ, alg′) - alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc) - C, Vᴴ = CVᴴ - S = Diagonal(initialize_output(svd_vals!, A, alg_trunc.alg)) - U, S, Vᴴ = svd_trunc!(A, (C, S, Vᴴ), alg_trunc) - return rmul!(U, S), Vᴴ +right_orth!(A, CVᴴ, alg::AbstractAlgorithm) = right_orth!(A, CVᴴ, right_orth_alg(alg)) +right_orth!(A, CVᴴ, alg::RightOrthViaLQ) = lq_compact!(A, CVᴴ, alg.alg) +right_orth!(A, CVᴴ, alg::RightOrthViaPolar) = right_polar!(A, CVᴴ, alg.alg) +function right_orth!(A, CVᴴ, alg::RightOrthViaSVD) + check_input(right_orth!, A, CVᴴ, alg) + C, S, Vᴴ = does_truncate(alg.alg) ? svd_trunc!(A, alg.alg) : svd_compact!(A, alg.alg) + rmul!(C, S) + return C, Vᴴ end # Implementation of null functions # -------------------------------- -function null_truncation_strategy(; atol = nothing, rtol = nothing, maxnullity = nothing) - if isnothing(maxnullity) && isnothing(atol) && isnothing(rtol) - return notrunc() - end - atol = @something atol 0 - rtol = @something rtol 0 - trunc = trunctol(; atol, rtol, keep_below = true) - return !isnothing(maxnullity) ? trunc & truncrank(maxnullity; rev = false) : trunc +left_null!(A, N, alg::AbstractAlgorithm) = left_null!(A, N, left_null_alg(alg)) +left_null!(A, N, alg::LeftNullViaQR) = qr_null!(A, N, alg.alg) +function left_null!(A, N, alg::LeftNullViaSVD{<:TruncatedAlgorithm}) + check_input(left_null!, A, N, alg) + U, S, _ = svd_full!(A, alg.alg.alg) + N, _ = truncate(left_null!, (U, S), alg.alg.trunc) + return N end -function left_null!( - A, N; - trunc = nothing, kind = isnothing(trunc) ? :qr : :svd, - alg_qr = (; positive = true), alg_svd = (;) - ) - if !isnothing(trunc) && kind != :svd - throw(ArgumentError("truncation not supported for left_null with kind=$kind")) - end - return if kind == :qr - left_null_qr!(A, N, alg_qr) - elseif kind == :svd - left_null_svd!(A, N, alg_svd, trunc) - else - throw(ArgumentError("`left_null!` received unknown value `kind = $kind`")) - end -end -function left_null_qr!(A, N, alg) - alg′ = select_algorithm(qr_null!, A, alg) - check_input(left_null!, A, N, alg′) - return qr_null!(A, N, alg′) -end -function left_null_svd!(A, N, alg, trunc::Nothing = nothing) - alg′ = select_algorithm(svd_full!, A, alg) - check_input(left_null!, A, N, alg′) - U, _, _ = svd_full!(A, alg′) - (m, n) = size(A) - return copy!(N, view(U, 1:m, (n + 1):m)) -end -function left_null_svd!(A, N, alg, trunc) - alg′ = select_algorithm(svd_full!, A, alg) - U, S, _ = svd_full!(A, alg′) - trunc′ = trunc isa TruncationStrategy ? trunc : - trunc isa NamedTuple ? null_truncation_strategy(; trunc...) : - throw(ArgumentError("Unknown truncation strategy: $trunc")) - return first(truncate(left_null!, (U, S), trunc′)) +right_null!(A, Nᴴ, alg::AbstractAlgorithm) = right_null!(A, Nᴴ, right_null_alg(alg)) +right_null!(A, Nᴴ, alg::RightNullViaLQ) = lq_null!(A, Nᴴ, alg.alg) +function right_null!(A, Nᴴ, alg::RightNullViaSVD{<:TruncatedAlgorithm}) + check_input(right_null!, A, Nᴴ, alg) + _, S, Vᴴ = svd_full!(A, alg.alg.alg) + Nᴴ, _ = truncate(right_null!, (S, Vᴴ), alg.alg.trunc) + return Nᴴ end -function right_null!( - A, Nᴴ; - trunc = nothing, kind = isnothing(trunc) ? :lq : :svd, - alg_lq = (; positive = true), alg_svd = (;) - ) - if !isnothing(trunc) && kind != :svd - throw(ArgumentError("truncation not supported for right_null with kind=$kind")) - end - if kind == :lq - return right_null_lq!(A, Nᴴ, alg_lq) - elseif kind == :svd - return right_null_svd!(A, Nᴴ, alg_svd, trunc) - else - throw(ArgumentError("`right_null!` received unknown value `kind = $kind`")) - end -end -function right_null_lq!(A, Nᴴ, alg) - alg′ = select_algorithm(lq_null!, A, alg) - check_input(right_null!, A, Nᴴ, alg′) - return lq_null!(A, Nᴴ, alg′) -end -function right_null_svd!(A, Nᴴ, alg, trunc::Nothing = nothing) - alg′ = select_algorithm(svd_full!, A, alg) - check_input(right_null!, A, Nᴴ, alg′) - _, _, Vᴴ = svd_full!(A, alg′) - (m, n) = size(A) - return copy!(Nᴴ, view(Vᴴ, (m + 1):n, 1:n)) -end -function right_null_svd!(A, Nᴴ, alg, trunc) - alg′ = select_algorithm(svd_full!, A, alg) - check_input(right_null!, A, Nᴴ, alg′) - _, S, Vᴴ = svd_full!(A, alg′) - trunc′ = trunc isa TruncationStrategy ? trunc : - trunc isa NamedTuple ? null_truncation_strategy(; trunc...) : - throw(ArgumentError("Unknown truncation strategy: $trunc")) - return first(truncate(right_null!, (S, Vᴴ), trunc′)) -end +# randomized algorithms don't currently work for smallest values: +left_null!(A, N, alg::LeftNullViaSVD{<:TruncatedAlgorithm{<:GPU_Randomized}}) = + throw(ArgumentError("Randomized SVD ($alg) cannot be used for null spaces yet")) +right_null!(A, Nᴴ, alg::RightNullViaSVD{<:TruncatedAlgorithm{<:GPU_Randomized}}) = + throw(ArgumentError("Randomized SVD ($alg) cannot be used for null spaces yet")) diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index 9ebc7319..fed36cd1 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -285,20 +285,6 @@ end # placed here to avoid code duplication since much of the logic is replicable across # CUDA and AMDGPU ### -const CUSOLVER_SVDAlgorithm = Union{ - CUSOLVER_QRIteration, - CUSOLVER_SVDPolar, - CUSOLVER_Jacobi, - CUSOLVER_Randomized, -} -const ROCSOLVER_SVDAlgorithm = Union{ - ROCSOLVER_QRIteration, - ROCSOLVER_Jacobi, -} -const GPU_SVDAlgorithm = Union{CUSOLVER_SVDAlgorithm, ROCSOLVER_SVDAlgorithm} - -const GPU_SVDPolar = Union{CUSOLVER_SVDPolar} -const GPU_Randomized = Union{CUSOLVER_Randomized} function check_input( ::typeof(svd_trunc!), A::AbstractMatrix, USVᴴ, alg::CUSOLVER_Randomized diff --git a/src/implementations/truncation.jl b/src/implementations/truncation.jl index f6201b07..883c7759 100644 --- a/src/implementations/truncation.jl +++ b/src/implementations/truncation.jl @@ -26,6 +26,19 @@ function truncate(::typeof(right_null!), (S, Vᴴ), strategy::TruncationStrategy return Vᴴ[ind, :], ind end +# special case `NoTruncation` for null: should keep exact zeros due to rectangularity +function truncate(::typeof(left_null!), (U, S), strategy::NoTruncation) + m, n = size(S) + ind = (n + 1):m + return U[:, ind], ind +end +function truncate(::typeof(right_null!), (S, Vᴴ), strategy::NoTruncation) + m, n = size(S) + ind = (m + 1):n + return Vᴴ[ind, :], ind +end + + # findtruncated # ------------- # Generic fallback @@ -95,16 +108,20 @@ function _truncerr_impl(values::AbstractVector, I; atol::Real = 0, rtol::Real = end function findtruncated(values::AbstractVector, strategy::TruncationIntersection) - return mapreduce( - Base.Fix1(findtruncated, values), _ind_intersect, strategy.components; - init = trues(length(values)) - ) + length(strategy.components) == 0 && return eachindex(values) + length(strategy.components) == 1 && return findtruncated(values, only(strategy.components)) + + ind1 = findtruncated(values, strategy.components[1]) + ind2 = findtruncated(values, TruncationIntersection(Base.tail(strategy.components))) + return _ind_intersect(ind1, ind2) end function findtruncated_svd(values::AbstractVector, strategy::TruncationIntersection) - return mapreduce( - Base.Fix1(findtruncated_svd, values), _ind_intersect, - strategy.components; init = trues(length(values)) - ) + length(strategy.components) == 0 && return eachindex(values) + length(strategy.components) == 1 && return findtruncated_svd(values, only(strategy.components)) + + ind1 = findtruncated_svd(values, strategy.components[1]) + ind2 = findtruncated_svd(values, TruncationIntersection(Base.tail(strategy.components))) + return _ind_intersect(ind1, ind2) end # when one of the ind selections is a bitvector, have to handle differently diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index bdda6612..c36a6b02 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -140,7 +140,7 @@ until convergence up to tolerance `tol`. @algdef PolarNewton # ========================= -# DIAGONAL ALGORITHMS +# Varia # ========================= """ DiagonalAlgorithm(; kwargs...) @@ -150,6 +150,21 @@ the diagonal structure of the input and outputs. """ @algdef DiagonalAlgorithm +""" + LQViaTransposedQR(qr_alg) + +Algorithm type to denote finding the LQ decomposition of `A` by computing the QR decomposition of `Aᵀ`. +The `qr_alg` specifies which QR-decomposition implementation to use. +""" +struct LQViaTransposedQR{A <: AbstractAlgorithm} <: AbstractAlgorithm + qr_alg::A +end +function Base.show(io::IO, alg::LQViaTransposedQR) + print(io, "LQViaTransposedQR(") + _show_alg(io, alg.qr_alg) + return print(io, ")") +end + # ========================= # CUSOLVER ALGORITHMS # ========================= @@ -203,6 +218,8 @@ for more information. """ @algdef CUSOLVER_Randomized +does_truncate(::TruncatedAlgorithm{<:CUSOLVER_Randomized}) = true + """ CUSOLVER_Simple() @@ -222,6 +239,10 @@ Divide and Conquer algorithm. """ @algdef CUSOLVER_DivideAndConquer +const CUSOLVER_SVDAlgorithm = Union{ + CUSOLVER_QRIteration, CUSOLVER_SVDPolar, CUSOLVER_Jacobi, CUSOLVER_Randomized, +} + # ========================= # ROCSOLVER ALGORITHMS # ========================= @@ -269,6 +290,10 @@ Divide and Conquer algorithm. """ @algdef ROCSOLVER_DivideAndConquer +const ROCSOLVER_SVDAlgorithm = Union{ROCSOLVER_QRIteration, ROCSOLVER_Jacobi} + +# Various consts and unions +# ------------------------- const GPU_Simple = Union{CUSOLVER_Simple} const GPU_EigAlgorithm = Union{GPU_Simple} @@ -277,8 +302,132 @@ const GPU_Jacobi = Union{CUSOLVER_Jacobi, ROCSOLVER_Jacobi} const GPU_DivideAndConquer = Union{CUSOLVER_DivideAndConquer, ROCSOLVER_DivideAndConquer} const GPU_Bisection = Union{ROCSOLVER_Bisection} const GPU_EighAlgorithm = Union{ - GPU_QRIteration, - GPU_Jacobi, - GPU_DivideAndConquer, - GPU_Bisection, + GPU_QRIteration, GPU_Jacobi, GPU_DivideAndConquer, GPU_Bisection, } +const GPU_SVDAlgorithm = Union{CUSOLVER_SVDAlgorithm, ROCSOLVER_SVDAlgorithm} + +const GPU_SVDPolar = Union{CUSOLVER_SVDPolar} +const GPU_Randomized = Union{CUSOLVER_Randomized} + +const QRAlgorithms = Union{LAPACK_HouseholderQR, CUSOLVER_HouseholderQR, ROCSOLVER_HouseholderQR} +const LQAlgorithms = Union{LAPACK_HouseholderLQ, LQViaTransposedQR} +const SVDAlgorithms = Union{LAPACK_SVDAlgorithm, GPU_SVDAlgorithm} +const PolarAlgorithms = Union{PolarViaSVD, PolarNewton} + +# ================================ +# ORTHOGONALIZATION ALGORITHMS +# ================================ + +""" + LeftOrthAlgorithm{Kind, Alg <: AbstractAlgorithm}(alg) + +Wrapper type to denote the `Kind` of factorization that is used as a backend for [`left_orth`](@ref). +By default `Kind` is a symbol, which can be either `:qr`, `:polar` or `:svd`. +""" +struct LeftOrthAlgorithm{Kind, Alg <: AbstractAlgorithm} <: AbstractAlgorithm + alg::Alg +end +LeftOrthAlgorithm{Kind}(alg::Alg) where {Kind, Alg <: AbstractAlgorithm} = LeftOrthAlgorithm{Kind, Alg}(alg) + +# Note: specific algorithm selection is handled by `left_orth_alg` in orthnull.jl +LeftOrthAlgorithm(alg::AbstractAlgorithm) = error( + """ + Unkown or invalid `left_orth` algorithm type `$(typeof(alg))`. + To register the algorithm type for `left_orth`, define + + MatrixAlgebraKit.left_orth_alg(alg::CustomAlgorithm) = LeftOrthAlgorithm{kind}(alg) + + where `kind` selects the factorization type that will be used. + By default, this is either `:qr`, `:polar` or `:svd`, to select [`qr_compact!`](@ref), + [`left_polar!`](@ref), [`svd_compact!`](@ref) or [`svd_trunc!`](@ref) respectively. + """ +) + +const LeftOrthViaQR = LeftOrthAlgorithm{:qr} +const LeftOrthViaPolar = LeftOrthAlgorithm{:polar} +const LeftOrthViaSVD = LeftOrthAlgorithm{:svd} + +""" + RightOrthAlgorithm{Kind, Alg <: AbstractAlgorithm}(alg) + +Wrapper type to denote the `Kind` of factorization that is used as a backend for [`right_orth`](@ref). +By default `Kind` is a symbol, which can be either `:lq`, `:polar` or `:svd`. +""" +struct RightOrthAlgorithm{Kind, Alg <: AbstractAlgorithm} <: AbstractAlgorithm + alg::Alg +end +RightOrthAlgorithm{Kind}(alg::Alg) where {Kind, Alg <: AbstractAlgorithm} = RightOrthAlgorithm{Kind, Alg}(alg) + +# Note: specific algorithm selection is handled by `right_orth_alg` in orthnull.jl +RightOrthAlgorithm(alg::AbstractAlgorithm) = error( + """ + Unkown or invalid `right_orth` algorithm type `$(typeof(alg))`. + To register the algorithm type for `right_orth`, define + + MatrixAlgebraKit.right_orth_alg(alg::CustomAlgorithm) = RightOrthAlgorithm{kind}(alg) + + where `kind` selects the factorization type that will be used. + By default, this is either `:lq`, `:polar` or `:svd`, to select [`lq_compact!`](@ref), + [`right_polar!`](@ref), [`svd_compact!`](@ref) or [`svd_trunc!`](@ref) respectively. + """ +) + +const RightOrthViaLQ = RightOrthAlgorithm{:lq} +const RightOrthViaPolar = RightOrthAlgorithm{:polar} +const RightOrthViaSVD = RightOrthAlgorithm{:svd} + +""" + LeftNullAlgorithm{Kind, Alg <: AbstractAlgorithm}(alg) + +Wrapper type to denote the `Kind` of factorization that is used as a backend for [`left_null`](@ref). +By default `Kind` is a symbol, which can be either `:qr` or `:svd`. +""" +struct LeftNullAlgorithm{Kind, Alg <: AbstractAlgorithm} <: AbstractAlgorithm + alg::Alg +end +LeftNullAlgorithm{Kind}(alg::Alg) where {Kind, Alg <: AbstractAlgorithm} = LeftNullAlgorithm{Kind, Alg}(alg) + +# Note: specific algorithm selection is handled by `left_null_alg` in orthnull.jl +LeftNullAlgorithm(alg::AbstractAlgorithm) = error( + """ + Unkown or invalid `left_null` algorithm type `$(typeof(alg))`. + To register the algorithm type for `left_null`, define + + MatrixAlgebraKit.left_null_alg(alg::CustomAlgorithm) = LeftNullAlgorithm{kind}(alg) + + where `kind` selects the factorization type that will be used. + By default, this is either `:qr` or `:svd`, to select [`qr_null!`](@ref), + [`svd_compact!`](@ref) or [`svd_trunc!`](@ref) respectively. + """ +) + +const LeftNullViaQR = LeftNullAlgorithm{:qr} +const LeftNullViaSVD = LeftNullAlgorithm{:svd} + +""" + RightNullAlgorithm{Kind, Alg <: AbstractAlgorithm}(alg) + +Wrapper type to denote the `Kind` of factorization that is used as a backend for [`right_null`](@ref). +By default `Kind` is a symbol, which can be either `:lq` or `:svd`. +""" +struct RightNullAlgorithm{Kind, Alg <: AbstractAlgorithm} <: AbstractAlgorithm + alg::Alg +end +RightNullAlgorithm{Kind}(alg::Alg) where {Kind, Alg <: AbstractAlgorithm} = RightNullAlgorithm{Kind, Alg}(alg) + +# Note: specific algorithm selection is handled by `right_null_alg` in orthnull.jl +RightNullAlgorithm(alg::AbstractAlgorithm) = error( + """ + Unkown or invalid `right_null` algorithm type `$(typeof(alg))`. + To register the algorithm type for `right_null`, define + + MatrixAlgebraKit.right_null_alg(alg::CustomAlgorithm) = RightNullAlgorithm{kind}(alg) + + where `kind` selects the factorization type that will be used. + By default, this is either `:lq` or `:svd`, to select [`lq_null!`](@ref), + [`svd_compact!`](@ref) or [`svd_trunc!`](@ref) respectively. + """ +) + +const RightNullViaLQ = RightNullAlgorithm{:lq} +const RightNullViaSVD = RightNullAlgorithm{:svd} diff --git a/src/interface/eig.jl b/src/interface/eig.jl index 867d69eb..28b5c69c 100644 --- a/src/interface/eig.jl +++ b/src/interface/eig.jl @@ -45,22 +45,28 @@ selected according to a truncation strategy. The function also returns `ϵ`, the truncation error defined as the 2-norm of the discarded eigenvalues. -## Keyword arguments -The behavior of this function is controlled by the following keyword arguments: - -- `trunc`: Specifies the truncation strategy. This can be: - - A `NamedTuple` with fields `atol`, `rtol`, and/or `maxrank`, which will be converted to - a [`TruncationStrategy`](@ref). For details on available truncation strategies, see - [Truncations](@ref). - - A `TruncationStrategy` object directly (e.g., `truncrank(10)`, `trunctol(atol=1e-6)`, or - combinations using `&`). - - `nothing` (default), which keeps all eigenvalues. - -- Other keyword arguments are passed to the algorithm selection procedure. If no explicit - `alg` is provided, these keywords are used to select and configure the algorithm through - [`MatrixAlgebraKit.select_algorithm`](@ref). The remaining keywords after algorithm - selection are passed to the algorithm constructor. See [`MatrixAlgebraKit.default_algorithm`](@ref) - for the default algorithm selection behavior. +## Truncation +The truncation strategy can be controlled via the `trunc` keyword argument. This can be +either a `NamedTuple` or a [`TruncationStrategy`](@ref). If `trunc` is not provided or +nothing, all values will be kept. + +### `trunc::NamedTuple` +The supported truncation keyword arguments are: + +$docs_truncation_kwargs + +### `trunc::TruncationStrategy` +For more control, a truncation strategy can be supplied directly. +By default, MatrixAlgebraKit supplies the following: + +$docs_truncation_strategies + +## Keyword Arguments +Other keyword arguments are passed to the algorithm selection procedure. If no explicit +`alg` is provided, these keywords are used to select and configure the algorithm through +[`MatrixAlgebraKit.select_algorithm`](@ref). The remaining keywords after algorithm +selection are passed to the algorithm constructor. See [`MatrixAlgebraKit.default_algorithm`](@ref) +for the default algorithm selection behavior. When `alg` is a [`TruncatedAlgorithm`](@ref), the `trunc` keyword cannot be specified as the truncation strategy is already embedded in the algorithm. @@ -71,7 +77,7 @@ truncation strategy is already embedded in the algorithm. as it may not always be possible to use the provided `DV` as output. !!! note -$(docs_eig_note) +$docs_eig_note See also [`eig_full(!)`](@ref eig_full), [`eig_vals(!)`](@ref eig_vals), and [Truncations](@ref) for more information on truncation strategies. @@ -102,7 +108,7 @@ See also [`eig_full(!)`](@ref eig_full) and [`eig_trunc(!)`](@ref eig_trunc). # ------------------- default_eig_algorithm(A; kwargs...) = default_eig_algorithm(typeof(A); kwargs...) default_eig_algorithm(T::Type; kwargs...) = throw(MethodError(default_eig_algorithm, (T,))) -function default_eig_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.BlasMat} +function default_eig_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.MaybeBlasMat} return LAPACK_Expert(; kwargs...) end function default_eig_algorithm(::Type{T}; kwargs...) where {T <: Diagonal} diff --git a/src/interface/eigh.jl b/src/interface/eigh.jl index 314cb934..97f8f95c 100644 --- a/src/interface/eigh.jl +++ b/src/interface/eigh.jl @@ -50,22 +50,28 @@ selected according to a truncation strategy. The function also returns `ϵ`, the truncation error defined as the 2-norm of the discarded eigenvalues. +## Truncation +The truncation strategy can be controlled via the `trunc` keyword argument. This can be +either a `NamedTuple` or a [`TruncationStrategy`](@ref). If `trunc` is not provided or +nothing, all values will be kept. + +### `trunc::NamedTuple` +The supported truncation keyword arguments are: + +$docs_truncation_kwargs + +### `trunc::TruncationStrategy` +For more control, a truncation strategy can be supplied directly. +By default, MatrixAlgebraKit supplies the following: + +$docs_truncation_strategies + ## Keyword arguments -The behavior of this function is controlled by the following keyword arguments: - -- `trunc`: Specifies the truncation strategy. This can be: - - A `NamedTuple` with fields `atol`, `rtol`, and/or `maxrank`, which will be converted to - a [`TruncationStrategy`](@ref). For details on available truncation strategies, see - [Truncations](@ref). - - A `TruncationStrategy` object directly (e.g., `truncrank(10)`, `trunctol(atol=1e-6)`, or - combinations using `&`). - - `nothing` (default), which keeps all eigenvalues. - -- Other keyword arguments are passed to the algorithm selection procedure. If no explicit - `alg` is provided, these keywords are used to select and configure the algorithm through - [`MatrixAlgebraKit.select_algorithm`](@ref). The remaining keywords after algorithm - selection are passed to the algorithm constructor. See [`MatrixAlgebraKit.default_algorithm`](@ref) - for the default algorithm selection behavior. +Other keyword arguments are passed to the algorithm selection procedure. If no explicit +`alg` is provided, these keywords are used to select and configure the algorithm through +[`MatrixAlgebraKit.select_algorithm`](@ref). The remaining keywords after algorithm +selection are passed to the algorithm constructor. See [`MatrixAlgebraKit.default_algorithm`](@ref) +for the default algorithm selection behavior. When `alg` is a [`TruncatedAlgorithm`](@ref), the `trunc` keyword cannot be specified as the truncation strategy is already embedded in the algorithm. @@ -109,7 +115,7 @@ default_eigh_algorithm(A; kwargs...) = default_eigh_algorithm(typeof(A); kwargs. function default_eigh_algorithm(T::Type; kwargs...) throw(MethodError(default_eigh_algorithm, (T,))) end -function default_eigh_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.BlasMat} +function default_eigh_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.MaybeBlasMat} return LAPACK_MultipleRelativelyRobustRepresentations(; kwargs...) end function default_eigh_algorithm(::Type{T}; kwargs...) where {T <: Diagonal} diff --git a/src/interface/gen_eig.jl b/src/interface/gen_eig.jl index 92cd3477..1903837f 100644 --- a/src/interface/gen_eig.jl +++ b/src/interface/gen_eig.jl @@ -60,7 +60,7 @@ See also [`gen_eig_full(!)`](@ref gen_eig_full). default_gen_eig_algorithm(A, B; kwargs...) = default_gen_eig_algorithm(typeof(A), typeof(B); kwargs...) default_gen_eig_algorithm(::Type{TA}, ::Type{TB}; kwargs...) where {TA, TB} = throw(MethodError(default_gen_eig_algorithm, (TA, TB))) -default_gen_eig_algorithm(::Type{TA}, ::Type{TB}; kwargs...) where {TA <: YALAPACK.BlasMat, TB <: YALAPACK.BlasMat} = +default_gen_eig_algorithm(::Type{TA}, ::Type{TB}; kwargs...) where {TA <: YALAPACK.MaybeBlasMat, TB <: YALAPACK.MaybeBlasMat} = LAPACK_Simple(; kwargs...) for f in (:gen_eig_full!, :gen_eig_vals!) diff --git a/src/interface/lq.jl b/src/interface/lq.jl index b0d138d4..a1efa39b 100644 --- a/src/interface/lq.jl +++ b/src/interface/lq.jl @@ -72,7 +72,7 @@ default_lq_algorithm(A; kwargs...) = default_lq_algorithm(typeof(A); kwargs...) function default_lq_algorithm(T::Type; kwargs...) throw(MethodError(default_lq_algorithm, (T,))) end -function default_lq_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.BlasMat} +function default_lq_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.MaybeBlasMat} return LAPACK_HouseholderLQ(; kwargs...) end function default_lq_algorithm(::Type{T}; kwargs...) where {T <: Diagonal} @@ -84,13 +84,3 @@ for f in (:lq_full!, :lq_compact!, :lq_null!) return default_lq_algorithm(A; kwargs...) end end - -# Alternative algorithm (necessary for CUDA) -struct LQViaTransposedQR{A <: AbstractAlgorithm} <: AbstractAlgorithm - qr_alg::A -end -function Base.show(io::IO, alg::LQViaTransposedQR) - print(io, "LQViaTransposedQR(") - _show_alg(io, alg.qr_alg) - return print(io, ")") -end diff --git a/src/interface/orthnull.jl b/src/interface/orthnull.jl index c27179f8..64acb509 100644 --- a/src/interface/orthnull.jl +++ b/src/interface/orthnull.jl @@ -1,220 +1,549 @@ # Orth functions # -------------- """ - left_orth(A; [kind::Symbol, trunc, alg_qr, alg_polar, alg_svd]) -> V, C - left_orth!(A, [VC]; [kind::Symbol, trunc, alg_qr, alg_polar, alg_svd]) -> V, C - -Compute an orthonormal basis `V` for the image of the matrix `A` of size `(m, n)`, -as well as a matrix `C` (the corestriction) such that `A` factors as `A = V * C`. -The keyword argument `kind` can be used to specify the specific orthogonal decomposition -that should be used to factor `A`, whereas `trunc` can be used to control the -precision in determining the rank of `A` via its singular values. - -`trunc` can either be a truncation strategy object or a NamedTuple with fields -`atol`, `rtol`, and `maxrank`. - -This is a high-level wrapper and will use one of the decompositions -[`qr_compact!`](@ref), [`svd_compact!`](@ref)/[`svd_trunc!`](@ref), and[`left_polar!`](@ref) -to compute the orthogonal basis `V`, as controlled by the keyword arguments. - -When `kind` is provided, its possible values are - -* `kind == :qr`: `V` and `C` are computed using the QR decomposition. - This requires `isnothing(trunc)` and `left_orth!(A, [VC])` is equivalent to - `qr_compact!(A, [VC], alg)` with a default value `alg = select_algorithm(qr_compact!, A; positive=true)` - -* `kind == :polar`: `V` and `C` are computed using the polar decomposition, - This requires `isnothing(trunc)` and `left_orth!(A, [VC])` is equivalent to - `left_polar!(A, [VC], alg)` with a default value `alg = select_algorithm(left_polar!, A)` - -* `kind == :svd`: `V` and `C` are computed using the singular value decomposition `svd_trunc!` when a - truncation strategy is specified using the `trunc` keyword argument, and using `svd_compact!` otherwise. - `V` will contain the left singular vectors and `C` is computed as the product of the singular - values and the right singular vectors, i.e. with `U, S, Vᴴ = svd(A)`, we have - `V = U` and `C = S * Vᴴ`. - -When `kind` is not provided, the default value is `:qr` when `isnothing(trunc)` -and `:svd` otherwise. Finally, finer control is obtained by providing an explicit algorithm -for backend factorizations through the `alg_qr`, `alg_polar`, and `alg_svd` keyword arguments, -which will only be used if the corresponding factorization is called based on the other inputs. -If NamedTuples are passed as `alg_qr`, `alg_polar`, or `alg_svd`, a default algorithm is chosen -with `select_algorithm` and the NamedTuple is passed as keyword arguments to that algorithm. -`alg_qr` defaults to `(; positive=true)` so that by default a positive QR decomposition will -be used. + left_orth(A; [alg], [trunc], kwargs...) -> V, C + left_orth!(A, [VC], [alg]; [trunc], kwargs...) -> V, C + +Compute an orthonormal basis `V` for the image of the matrix `A`, as well as a matrix `C` +(the corestriction) such that `A` factors as `A = V * C`. + +This is a high-level wrapper where the keyword arguments can be used to specify and control +the specific orthogonal decomposition that should be used to factor `A`, whereas `trunc` +can optionally be used to control the precision in determining the rank of `A`, typically +via its singular values. + +## Truncation +The optional truncation strategy can be controlled via the `trunc` keyword argument, and +any non-trivial strategy typically requires an SVD-based decompositions. This keyword can +be either a `NamedTuple` or a [`TruncationStrategy`](@ref). + +### `trunc::NamedTuple` +The supported truncation keyword arguments are: + +$(docs_truncation_kwargs) + +### `trunc::TruncationStrategy` +For more control, a truncation strategy can be supplied directly. +By default, MatrixAlgebraKit supplies the following: + +$(docs_truncation_strategies) + +## Keyword arguments +There are 3 major modes of operation, based on the `alg` keyword, with slightly different +application purposes. + +### `alg::Nothing` +This default mode uses the presence of a truncation strategy `trunc` to determine an optimal +decomposition type, which will typically be QR-based for no truncation, or SVD-based for +truncation. The remaining keyword arguments are passed on directly to the algorithm selection +procedure of the chosen decomposition type. + +### `alg::Symbol` +Here, the driving selector is `alg`, which is used to select the kind of decomposition. The +remaining keyword arguments are passed on directly to the algorithm selection procedure of +the chosen decomposition type. By default, the supported kinds are: + +* `:qr` : Factorize via QR decomposition, with further customizations through the other + keywords. This mode requires `isnothing(trunc)`, and is equivalent to +```julia + V, C = qr_compact(A; kwargs...) +``` + +* `:polar` : Factorize via polar decomposition, with further customizations through the + other keywords. This mode requires `isnothing(trunc)`, and is equivalent to +```julia + V, C = left_polar(A; kwargs...) +``` + +* `:svd` : Factorize via SVD, with further customizations through the other keywords. + This mode further allows truncation, which can be selected through the `trunc` argument, + and is roughly equivalent to: +```julia + V, S, C = svd_trunc(A; trunc, kwargs...) + C = lmul!(S, C) +``` + +### `alg::AbstractAlgorithm` +In this expert mode the algorithm is supplied directly, and the kind of decomposition is +deduced from that. This is achieved either directly by providing a +[`LeftOrthAlgorithm{kind}`](@ref LeftOrthAlgorithm), or automatically by attempting to +deduce the decomposition kind with [`left_orth_alg(alg)`](@ref left_orth_alg). + +--- !!! note - The bang method `left_orth!` optionally accepts the output structure and possibly destroys - the input matrix `A`. Always use the return value of the function as it may not always be - possible to use the provided `CV` as output. + The bang method `left_orth!` optionally accepts the output structure and possibly + destroys the input matrix `A`. Always use the return value of the function as it may + not always be possible to use the provided `CV` as output. -See also [`right_orth(!)`](@ref right_orth), [`left_null(!)`](@ref left_null), [`right_null(!)`](@ref right_null) +See also [`right_orth(!)`](@ref right_orth), [`left_null(!)`](@ref left_null) and +[`right_null(!)`](@ref right_null). """ -function left_orth end -function left_orth! end -function left_orth!(A; kwargs...) - return left_orth!(A, initialize_output(left_orth!, A); kwargs...) -end -function left_orth(A; kwargs...) - return left_orth!(copy_input(left_orth, A); kwargs...) -end +@functiondef left_orth + +""" + right_orth(A; [alg], [trunc], kwargs...) -> C, Vᴴ + right_orth!(A, [CVᴴ], [alg]; [trunc], kwargs...) -> C, Vᴴ + +Compute an orthonormal basis `V = adjoint(Vᴴ)` for the coimage of the matrix `A`, i.e. for +the image of `adjoint(A)`, as well as a matrix `C` such that `A` factors as `A = C * Vᴴ`. + +This is a high-level wrapper where the keyword arguments can be used to specify and control +the specific orthogonal decomposition that should be used to factor `A`, whereas `trunc` +can optionally be used to control the precision in determining the rank of `A`, typically +via its singular values. + +## Truncation +The optional truncation strategy can be controlled via the `trunc` keyword argument, and +any non-trivial strategy typically requires an SVD-based decompositions. This keyword can +be either a `NamedTuple` or a [`TruncationStrategy`](@ref). + +### `trunc::NamedTuple` +The supported truncation keyword arguments are: + +$(docs_truncation_kwargs) + +### `trunc::TruncationStrategy` +For more control, a truncation strategy can be supplied directly. +By default, MatrixAlgebraKit supplies the following: + +$(docs_truncation_strategies) + +## Keyword arguments +There are 3 major modes of operation, based on the `alg` keyword, with slightly different +application purposes. + +### `alg::Nothing` +This default mode uses the presence of a truncation strategy `trunc` to determine an optimal +decomposition type, which will typicall be LQ-based for no truncation, or SVD-based for +truncation. The remaining keyword arguments are passed on directly to the algorithm selection +procedure of the chosen decomposition type. + +### `alg::Symbol` +Here, the driving selector is `alg`, which is used to select the kind of decomposition. The +remaining keyword arguments are passed on directly to the algorithm selection procedure of +the chosen decomposition type. By default, the supported kinds are: + +* `:lq` : Factorize via LQ decomposition, with further customizations through the other + keywords. This mode requires `isnothing(trunc)`, and is equivalent to +```julia + C, Vᴴ = lq_compact(A; kwargs...) +``` + +* `:polar` : Factorize via polar decomposition, with further customizations through the + other keywords. This mode requires `isnothing(trunc)`, and is equivalent to +```julia + C, Vᴴ = right_polar(A; kwargs...) +``` + +* `:svd` : Factorize via SVD, with further customizations through the other keywords. + This mode further allows truncation, which can be selected through the `trunc` argument, + and is roughly equivalent to: +```julia + C, S, Vᴴ = svd_trunc(A; trunc, kwargs...) + C = rmul!(C, S) +``` + +### `alg::AbstractAlgorithm` +In this expert mode the algorithm is supplied directly, and the kind of decomposition is +deduced from that. This is achieved either directly by providing a +[`RightOrthAlgorithm{kind}`](@ref RightOrthAlgorithm), or automatically by attempting to +deduce the decomposition kind with [`right_orth_alg`](@ref). + +--- + +!!! note + The bang method `right_orth!` optionally accepts the output structure and possibly + destroys the input matrix `A`. Always use the return value of the function as it may not + always be possible to use the provided `CVᴴ` as output. + +See also [`left_orth(!)`](@ref left_orth), [`left_null(!)`](@ref left_null) and +[`right_null(!)`](@ref right_null). +""" +@functiondef right_orth +# Null functions +# -------------- """ - right_orth(A; [kind::Symbol, trunc, alg_lq, alg_polar, alg_svd]) -> C, Vᴴ - right_orth!(A, [CVᴴ]; [kind::Symbol, trunc, alg_lq, alg_polar, alg_svd]) -> C, Vᴴ + left_null(A; [alg], [trunc], kwargs...) -> N + left_null!(A, [N], [alg]; [trunc], kwargs...) -> N -Compute an orthonormal basis `V = adjoint(Vᴴ)` for the coimage of the matrix `A`, i.e. -for the image of `adjoint(A)`, as well as a matrix `C` such that `A = C * Vᴴ`. -The keyword argument `kind` can be used to specify the specific orthogonal decomposition -that should be used to factor `A`, whereas `trunc` can be used to control the -precision in determining the rank of `A` via its singular values. +Compute an orthonormal basis `N` for the cokernel of the matrix `A`, i.e. the nullspace of +`adjoint(A)`, such that `adjoint(A) * N ≈ 0` and `N' * N ≈ I`. -`trunc` can either be a truncation strategy object or a NamedTuple with fields -`atol`, `rtol`, and `maxrank`. +This is a high-level wrapper where the keyword arguments can be used to specify and control +the underlying orthogonal decomposition that should be used to find the null space of `A'`, +whereas `trunc` can optionally be used to control the precision in determining the rank of +`A`, typically via its singular values. -This is a high-level wrapper and will use one of the decompositions -[`lq_compact!`](@ref), [`svd_compact!`](@ref)/[`svd_trunc!`](@ref), and -[`right_polar!`](@ref) to compute the orthogonal basis `V`, as controlled by the -keyword arguments. +## Truncation +The optional truncation strategy can be controlled via the `trunc` keyword argument, and any +non-trivial strategy typically requires an SVD-based decomposition. This keyword can be +either a `NamedTuple` or a [`TruncationStrategy`](@ref). -When `kind` is provided, its possible values are +### `trunc::NamedTuple` +The supported truncation keyword arguments are: -* `kind == :lq`: `C` and `Vᴴ` are computed using the QR decomposition, - This requires `isnothing(trunc)` and `right_orth!(A, [CVᴴ])` is equivalent to - `lq_compact!(A, [CVᴴ], alg)` with a default value `alg = select_algorithm(lq_compact!, A; positive=true)` +$(docs_null_truncation_kwargs) -* `kind == :polar`: `C` and `Vᴴ` are computed using the polar decomposition, - This requires `isnothing(trunc)` and `right_orth!(A, [CVᴴ])` is equivalent to - `right_polar!(A, [CVᴴ], alg)` with a default value `alg = select_algorithm(right_polar!, A)` +### `trunc::TruncationStrategy` +For more control, a truncation strategy can be supplied directly. By default, +MatrixAlgebraKit supplies the following: -* `kind == :svd`: `C` and `Vᴴ` are computed using the singular value decomposition `svd_trunc!` when - a truncation strategy is specified using the `trunc` keyword argument, and using `svd_compact!` otherwise. - `V = adjoint(Vᴴ)` will contain the right singular vectors corresponding to the singular - values and `C` is computed as the product of the singular values and the right singular vectors, - i.e. with `U, S, Vᴴ = svd(A)`, we have `C = rmul!(U, S)` and `Vᴴ = Vᴴ`. +$(docs_truncation_strategies) -When `kind` is not provided, the default value is `:lq` when `isnothing(trunc)` -and `:svd` otherwise. Finally, finer control is obtained by providing an explicit algorithm -for backend factorizations through the `alg_lq`, `alg_polar`, and `alg_svd` keyword arguments, -which will only be used if the corresponding factorization is called based on the other inputs. -If `alg_lq`, `alg_polar`, or `alg_svd` are NamedTuples, a default algorithm is chosen -with `select_algorithm` and the NamedTuple is passed as keyword arguments to that algorithm. -`alg_lq` defaults to `(; positive=true)` so that by default a positive LQ decomposition will -be used. +!!! note + Here [`notrunc`](@ref) has special meaning, and signifies keeping the values that + correspond to the exact zeros determined from the additional columns of `A`. + +## Keyword arguments +There are 3 major modes of operation, based on the `alg` keyword, with slightly different +application purposes. + +### `alg::Nothing` +This default mode uses the presence of a truncation strategy `trunc` to determine an optimal +decomposition type, which will be QR-based for no truncation, or SVD-based for truncation. +The remaining keyword arguments are passed on directly to the algorithm selection procedure +of the chosen decomposition type. + +### `alg::Symbol` +Here, the driving selector is `alg`, which is used to select the kind of decomposition. The +remaining keyword arguments are passed on directly to the algorithm selection procedure of +the chosen decomposition type. By default, the supported kinds are: + +* `:qr` : Factorize via QR nullspace, with further customizations through the other + keywords. This mode requires `isnothing(trunc)`, and is equivalent to +```julia + N = qr_null(A; kwargs...) +``` + +* `:svd` : Factorize via SVD, with further customizations through the other keywords. + This mode further allows truncation, which can be selected through the `trunc` argument. + It is roughly equivalent to: +```julia + U, S, _ = svd_full(A; kwargs...) + N = truncate(left_null, (U, S), trunc) +``` + +### `alg::AbstractAlgorithm` +In this expert mode the algorithm is supplied directly, and the kind of decomposition is +deduced from that. This is achieved either directly by providing a +[`LeftNullAlgorithm{kind}`](@ref LeftNullAlgorithm), or automatically by attempting to +deduce the decomposition kind with [`left_null_alg(alg)`](@ref left_null_alg). + +--- !!! note - The bang method `right_orth!` optionally accepts the output structure and possibly destroys - the input matrix `A`. Always use the return value of the function as it may not always be - possible to use the provided `CVᴴ` as output. + The bang method `left_null!` optionally accepts the output structure and possibly + destroys the input matrix `A`. Always use the return value of the function as it may not + always be possible to use the provided `N` as output. -See also [`left_orth(!)`](@ref left_orth), [`left_null(!)`](@ref left_null), [`right_null(!)`](@ref right_null) +See also [`right_null(!)`](@ref right_null), [`left_orth(!)`](@ref left_orth) and +[`right_orth(!)`](@ref right_orth). """ -function right_orth end -function right_orth! end -function right_orth!(A; kwargs...) - return right_orth!(A, initialize_output(right_orth!, A); kwargs...) -end -function right_orth(A; kwargs...) - return right_orth!(copy_input(right_orth, A); kwargs...) -end +@functiondef left_null -# Null functions -# -------------- """ - left_null(A; [kind::Symbol, trunc, alg_qr, alg_svd]) -> N - left_null!(A, [N]; [kind::Symbol, alg_qr, alg_svd]) -> N + right_null(A; [alg], [trunc], kwargs...) -> Nᴴ + right_null!(A, [Nᴴ], [alg]; [trunc], kwargs...) -> Nᴴ + +Compute an orthonormal basis `N = adjoint(Nᴴ)` for the kernel of the matrix `A`, i.e. the +nullspace of `A`, such that `A * Nᴴ' ≈ 0` and `Nᴴ * Nᴴ' ≈ I`. -Compute an orthonormal basis `N` for the cokernel of the matrix `A` of size `(m, n)`, i.e. -the nullspace of `adjoint(A)`, such that `adjoint(A)*N ≈ 0` and `N'*N ≈ I`. -The keyword argument `kind` can be used to specify the specific orthogonal decomposition -that should be used to factor `A`, whereas `trunc` can be used to control the -the rank of `A` via its singular values. +This is a high-level wrapper where the keyword arguments can be used to specify and control +the underlying orthogonal decomposition that should be used to find the null space of `A`, +whereas `trunc` can optionally be used to control the precision in determining the rank of +`A`, typically via its singular values. -`trunc` can either be a truncation strategy object or a NamedTuple with fields -`atol`, `rtol`, and `maxnullity`. +## Truncation +The optional truncation strategy can be controlled via the `trunc` keyword argument, and any +non-trivial strategy typically requires an SVD-based decomposition. This keyword can be +either a `NamedTuple` or a [`TruncationStrategy`](@ref). -This is a high-level wrapper and will use one of the decompositions `qr!` or `svd!` -to compute the orthogonal basis `N`, as controlled by the keyword arguments. +### `trunc::NamedTuple` +The supported truncation keyword arguments are: -When `kind` is provided, its possible values are +$(docs_null_truncation_kwargs) -* `kind == :qr`: `N` is computed using the QR decomposition. - This requires `isnothing(trunc)` and `left_null!(A, [N], kind=:qr)` is equivalent to - `qr_null!(A, [N], alg)` with a default value `alg = select_algorithm(qr_compact!, A; positive=true)` +### `trunc::TruncationStrategy` +For more control, a truncation strategy can be supplied directly. By default, +MatrixAlgebraKit supplies the following: -* `kind == :svd`: `N` is computed using the singular value decomposition and will contain - the left singular vectors corresponding to either the zero singular values if `trunc` - isn't specified or the singular values specified by `trunc`. +$(docs_truncation_strategies) -When `kind` is not provided, the default value is `:qr` when `isnothing(trunc)` -and `:svd` otherwise. Finally, finer control is obtained by providing an explicit algorithm -using the `alg_qr` and `alg_svd` keyword arguments, which will only be used by the corresponding -factorization backend. If `alg_qr` or `alg_svd` are NamedTuples, a default algorithm is chosen -with `select_algorithm` and the NamedTuple is passed as keyword arguments to that algorithm. -`alg_qr` defaults to `(; positive=true)` so that by default a positive QR decomposition will -be used. +!!! note + Here [`notrunc`](@ref) has special meaning, and signifies keeping the values that + correspond to the exact zeros determined from the additional rows of `A`. + +## Keyword arguments +There are 3 major modes of operation, based on the `alg` keyword, with slightly different +application purposes. + +### `alg::Nothing` +This default mode uses the presence of a truncation strategy `trunc` to determine an optimal +decomposition type, which will be LQ-based for no truncation, or SVD-based for truncation. +The remaining keyword arguments are passed on directly to the algorithm selection procedure +of the chosen decomposition type. + +### `alg::Symbol` +Here, the driving selector is `alg`, which is used to select the kind of decomposition. The +remaining keyword arguments are passed on directly to the algorithm selection procedure of +the chosen decomposition type. By default, the supported kinds are: + +* `:lq` : Factorize via LQ nullspace, with further customizations through the other + keywords. This mode requires `isnothing(trunc)`, and is equivalent to +```julia + Nᴴ = lq_null(A; kwargs...) +``` + +* `:svd` : Factorize via SVD, with further customizations through the other keywords. + This mode further allows truncation, which can be selected through the `trunc` argument. + It is roughly equivalent to: +```julia + _, S, Vᴴ = svd_full(A; kwargs...) + Nᴴ = truncate(right_null, (S, Vᴴ), trunc) +``` + +### `alg::AbstractAlgorithm` +In this expert mode the algorithm is supplied directly, and the kind of decomposition is +deduced from that. This is achieved either directly by providing a +[`RightNullAlgorithm{kind}`](@ref RightNullAlgorithm), or automatically by attempting to +deduce the decomposition kind with [`right_null_alg(alg)`](@ref right_null_alg). + +--- !!! note - The bang method `left_null!` optionally accepts the output structure and possibly destroys - the input matrix `A`. Always use the return value of the function as it may not always be - possible to use the provided `N` as output. + The bang method `right_null!` optionally accepts the output structure and possibly + destroys the input matrix `A`. Always use the return value of the function as it may not + always be possible to use the provided `Nᴴ` as output. -See also [`right_null(!)`](@ref right_null), [`left_orth(!)`](@ref left_orth), [`right_orth(!)`](@ref right_orth) +See also [`left_null(!)`](@ref left_null), [`left_orth(!)`](@ref left_orth) and +[`right_orth(!)`](@ref right_orth). """ -function left_null end -function left_null! end -function left_null!(A; kwargs...) - return left_null!(A, initialize_output(left_null!, A); kwargs...) +@functiondef right_null + +# Algorithm selection +# ------------------- +# specific override for `alg::Symbol` case, to allow for choosing the kind of factorization. +@inline select_algorithm(::typeof(left_orth!), A, alg::Symbol; kwargs...) = + select_algorithm(left_orth!, A, Val(alg); kwargs...) +@inline select_algorithm(::typeof(right_orth!), A, alg::Symbol; kwargs...) = + select_algorithm(right_orth!, A, Val(alg); kwargs...) +@inline select_algorithm(::typeof(left_null!), A, alg::Symbol; kwargs...) = + select_algorithm(left_null!, A, Val(alg); kwargs...) +@inline select_algorithm(::typeof(right_null!), A, alg::Symbol; kwargs...) = + select_algorithm(right_null!, A, Val(alg); kwargs...) + +function select_algorithm(::typeof(left_orth!), A, ::Val{:qr}; trunc = nothing, kwargs...) + isnothing(trunc) || + throw(ArgumentError("QR-based `left_orth` is incompatible with specifying `trunc`")) + alg′ = select_algorithm(qr_compact!, A; kwargs...) + return LeftOrthViaQR(alg′) +end +function select_algorithm(::typeof(left_orth!), A, ::Val{:polar}; trunc = nothing, kwargs...) + isnothing(trunc) || + throw(ArgumentError("Polar-based `left_orth` is incompatible with specifying `trunc`")) + alg′ = select_algorithm(left_polar!, A; kwargs...) + return LeftOrthViaPolar(alg′) +end +function select_algorithm(::typeof(left_orth!), A, ::Val{:svd}; trunc = nothing, kwargs...) + alg′ = isnothing(trunc) ? select_algorithm(svd_compact!, A; kwargs...) : + select_algorithm(svd_trunc!, A; trunc, kwargs...) + return LeftOrthViaSVD(alg′) +end + +function select_algorithm(::typeof(right_orth!), A, ::Val{:lq}; trunc = nothing, kwargs...) + isnothing(trunc) || + throw(ArgumentError("LQ-based `right_orth` is incompatible with specifying `trunc`")) + alg = select_algorithm(lq_compact!, A; kwargs...) + return RightOrthViaLQ(alg) +end +function select_algorithm(::typeof(right_orth!), A, ::Val{:polar}; trunc = nothing, kwargs...) + isnothing(trunc) || + throw(ArgumentError("Polar-based `right_orth` is incompatible with specifying `trunc`")) + alg = select_algorithm(right_polar!, A; kwargs...) + return RightOrthViaPolar(alg) +end +function select_algorithm(::typeof(right_orth!), A, ::Val{:svd}; trunc = nothing, kwargs...) + alg′ = isnothing(trunc) ? select_algorithm(svd_compact!, A; kwargs...) : + select_algorithm(svd_trunc!, A; trunc, kwargs...) + return RightOrthViaSVD(alg′) +end + +function select_algorithm(::typeof(left_null!), A, ::Val{:qr}; trunc = nothing, kwargs...) + isnothing(trunc) || + throw(ArgumentError("QR-based `left_null` is incompatible with specifying `trunc`")) + alg = select_algorithm(qr_null!, A; kwargs...) + return LeftNullViaQR(alg) +end +function select_algorithm(::typeof(left_null!), A, ::Val{:svd}; trunc = nothing, kwargs...) + alg_svd = select_algorithm(svd_full!, A, get(kwargs, :svd, nothing)) + alg = TruncatedAlgorithm(alg_svd, select_null_truncation(trunc)) + return LeftNullViaSVD(alg) +end + +function select_algorithm(::typeof(right_null!), A, ::Val{:lq}; trunc = nothing, kwargs...) + isnothing(trunc) || + throw(ArgumentError("LQ-based `right_null` is incompatible with specifying `trunc`")) + alg = select_algorithm(lq_null!, A; kwargs...) + return RightNullViaLQ(alg) end -function left_null(A; kwargs...) - return left_null!(copy_input(left_null, A); kwargs...) +function select_algorithm(::typeof(right_null!), A, ::Val{:svd}; trunc = nothing, kwargs...) + alg_svd = select_algorithm(svd_full!, A; kwargs...) + alg = TruncatedAlgorithm(alg_svd, select_null_truncation(trunc)) + return RightNullViaSVD(alg) end +default_algorithm(::typeof(left_orth!), ::Type{A}; trunc = nothing, kwargs...) where {A} = + isnothing(trunc) ? select_algorithm(left_orth!, A, Val(:qr); kwargs...) : + select_algorithm(left_orth!, A, Val(:svd); trunc, kwargs...) + +default_algorithm(::typeof(right_orth!), ::Type{A}; trunc = nothing, kwargs...) where {A} = + isnothing(trunc) ? select_algorithm(right_orth!, A, Val(:lq); kwargs...) : + select_algorithm(right_orth!, A, Val(:svd); trunc, kwargs...) + +default_algorithm(::typeof(left_null!), ::Type{A}; trunc = nothing, kwargs...) where {A} = + isnothing(trunc) ? select_algorithm(left_null!, A, Val(:qr); kwargs...) : + select_algorithm(left_null!, A, Val(:svd); trunc, kwargs...) + +default_algorithm(::typeof(right_null!), ::Type{A}; trunc = nothing, kwargs...) where {A} = + isnothing(trunc) ? select_algorithm(right_null!, A, Val(:lq); kwargs...) : + select_algorithm(right_null!, A, Val(:svd); trunc, kwargs...) + """ - right_null(A; [kind::Symbol, alg_lq, alg_svd]) -> Nᴴ - right_null!(A, [Nᴴ]; [kind::Symbol, alg_lq, alg_svd]) -> Nᴴ + left_orth_alg(alg::AbstractAlgorithm) -> LeftOrthAlgorithm -Compute an orthonormal basis `N = adjoint(Nᴴ)` for the kernel or nullspace of the matrix `A` -of size `(m, n)`, such that `A*adjoint(Nᴴ) ≈ 0` and `Nᴴ*adjoint(Nᴴ) ≈ I`. -The keyword argument `kind` can be used to specify the specific orthogonal decomposition -that should be used to factor `A`, whereas `trunc` can be used to control the -the rank of `A` via its singular values. +Convert an algorithm to a [`LeftOrthAlgorithm`](@ref) wrapper for use with [`left_orth`](@ref). -`trunc` can either be a truncation strategy object or a NamedTuple with fields -`atol`, `rtol`, and `maxnullity`. +This function attempts to deduce the appropriate factorization kind (`:qr`, `:polar`, or `:svd`) +from the algorithm type and wraps it in a `LeftOrthAlgorithm`. Custom algorithm types can be +registered by defining: -This is a high-level wrapper and will use one of the decompositions `lq!` or `svd!` -to compute the orthogonal basis `Nᴴ`, as controlled by the keyword arguments. +```julia +MatrixAlgebraKit.left_orth_alg(alg::CustomAlgorithm) = LeftOrthAlgorithm{kind}(alg) +``` -When `kind` is provided, its possible values are +where `kind` specifies the factorization backend to use. -* `kind == :lq`: `Nᴴ` is computed using the (nonpositive) LQ decomposition. - This requires `isnothing(trunc)` and `right_null!(A, [Nᴴ], kind=:lq)` is equivalent to - `lq_null!(A, [Nᴴ], alg)` with a default value `alg = select_algorithm(lq_compact!, A; positive=true)` +See also [`LeftOrthAlgorithm`](@ref), [`left_orth`](@ref). +""" +left_orth_alg(alg::AbstractAlgorithm) = error( + """ + Unkown or invalid `left_orth` algorithm type `$(typeof(alg))`. + To register the algorithm type for `left_orth`, define + + MatrixAlgebraKit.left_orth_alg(alg::CustomAlgorithm) = LeftOrthAlgorithm{kind}(alg) + + where `kind` selects the factorization type that will be used. + By default, this is either `:qr`, `:polar` or `:svd`, to select [`qr_compact!`](@ref), + [`left_polar!`](@ref), [`svd_compact!`](@ref) or [`svd_trunc!`](@ref) respectively. + """ +) +left_orth_alg(alg::LeftOrthAlgorithm) = alg +left_orth_alg(alg::QRAlgorithms) = LeftOrthViaQR(alg) +left_orth_alg(alg::PolarAlgorithms) = LeftOrthViaPolar(alg) +left_orth_alg(alg::SVDAlgorithms) = LeftOrthViaSVD(alg) +left_orth_alg(alg::TruncatedAlgorithm{<:SVDAlgorithms}) = LeftOrthViaSVD(alg) -* `kind == :svd`: `N` is computed using the singular value decomposition and will contain - the left singular vectors corresponding to the singular values that - are smaller than `max(atol, rtol * σ₁)`, where `σ₁` is the largest singular value of `A`. +""" + right_orth_alg(alg::AbstractAlgorithm) -> RightOrthAlgorithm -When `kind` is not provided, the default value is `:lq` when `isnothing(trunc)` -and `:svd` otherwise. Finally, finer control is obtained by providing an explicit algorithm -using the `alg_lq` and `alg_svd` keyword arguments, which will only be used by the corresponding -factorization backend. If `alg_lq` or `alg_svd` are NamedTuples, a default algorithm is chosen -with `select_algorithm` and the NamedTuple is passed as keyword arguments to that algorithm. -`alg_lq` defaults to `(; positive=true)` so that by default a positive LQ decomposition will -be used. +Convert an algorithm to a [`RightOrthAlgorithm`](@ref) wrapper for use with [`right_orth`](@ref). -!!! note - The bang method `right_null!` optionally accepts the output structure and possibly destroys - the input matrix `A`. Always use the return value of the function as it may not always be - possible to use the provided `Nᴴ` as output. +This function attempts to deduce the appropriate factorization kind (`:lq`, `:polar`, or `:svd`) +from the algorithm type and wraps it in a `RightOrthAlgorithm`. Custom algorithm types can be +registered by defining: -See also [`left_null(!)`](@ref left_null), [`left_orth(!)`](@ref left_orth), [`right_orth(!)`](@ref right_orth) +```julia +MatrixAlgebraKit.right_orth_alg(alg::CustomAlgorithm) = RightOrthAlgorithm{kind}(alg) +``` + +where `kind` specifies the factorization backend to use. + +See also [`RightOrthAlgorithm`](@ref), [`right_orth`](@ref). """ -function right_null end -function right_null! end -function right_null!(A; kwargs...) - return right_null!(A, initialize_output(right_null!, A); kwargs...) -end -function right_null(A; kwargs...) - return right_null!(copy_input(right_null, A); kwargs...) -end +right_orth_alg(alg::AbstractAlgorithm) = error( + """ + Unkown or invalid `right_orth` algorithm type `$(typeof(alg))`. + To register the algorithm type for `right_orth`, define + + MatrixAlgebraKit.right_orth_alg(alg::CustomAlgorithm) = RightOrthAlgorithm{kind}(alg) + + where `kind` selects the factorization type that will be used. + By default, this is either `:lq`, `:polar` or `:svd`, to select [`lq_compact!`](@ref), + [`right_polar!`](@ref), [`svd_compact!`](@ref) or [`svd_trunc!`](@ref) respectively. + """ +) +right_orth_alg(alg::RightOrthAlgorithm) = alg +right_orth_alg(alg::LQAlgorithms) = RightOrthViaLQ(alg) +right_orth_alg(alg::PolarAlgorithms) = RightOrthViaPolar(alg) +right_orth_alg(alg::SVDAlgorithms) = RightOrthViaSVD(alg) +right_orth_alg(alg::TruncatedAlgorithm{<:SVDAlgorithms}) = RightOrthViaSVD(alg) + +""" + left_null_alg(alg::AbstractAlgorithm) -> LeftNullAlgorithm + +Convert an algorithm to a [`LeftNullAlgorithm`](@ref) wrapper for use with [`left_null`](@ref). + +This function attempts to deduce the appropriate factorization kind (`:qr` or `:svd`) from +the algorithm type and wraps it in a `LeftNullAlgorithm`. Custom algorithm types can be +registered by defining: + +```julia +MatrixAlgebraKit.left_null_alg(alg::CustomAlgorithm) = LeftNullAlgorithm{kind}(alg) +``` + +where `kind` specifies the factorization backend to use. + +See also [`LeftNullAlgorithm`](@ref), [`left_null`](@ref). +""" +left_null_alg(alg::AbstractAlgorithm) = error( + """ + Unkown or invalid `left_null` algorithm type `$(typeof(alg))`. + To register the algorithm type for `left_null`, define + + MatrixAlgebraKit.left_null_alg(alg::CustomAlgorithm) = LeftNullAlgorithm{kind}(alg) + + where `kind` selects the factorization type that will be used. + By default, this is either `:qr` or `:svd`, to select [`qr_null!`](@ref), + [`svd_compact!`](@ref) or [`svd_trunc!`](@ref) respectively. + """ +) +left_null_alg(alg::LeftNullAlgorithm) = alg +left_null_alg(alg::QRAlgorithms) = LeftNullViaQR(alg) +left_null_alg(alg::SVDAlgorithms) = LeftNullViaSVD(alg) +left_null_alg(alg::TruncatedAlgorithm{<:SVDAlgorithms}) = LeftNullViaSVD(alg) + +""" + right_null_alg(alg::AbstractAlgorithm) -> RightNullAlgorithm + +Convert an algorithm to a [`RightNullAlgorithm`](@ref) wrapper for use with [`right_null`](@ref). + +This function attempts to deduce the appropriate factorization kind (`:lq` or `:svd`) from +the algorithm type and wraps it in a `RightNullAlgorithm`. Custom algorithm types can be +registered by defining: + +```julia +MatrixAlgebraKit.right_null_alg(alg::CustomAlgorithm) = RightNullAlgorithm{kind}(alg) +``` + +where `kind` specifies the factorization backend to use. + +See also [`RightNullAlgorithm`](@ref), [`right_null`](@ref). +""" +right_null_alg(alg::AbstractAlgorithm) = error( + """ + Unkown or invalid `right_null` algorithm type `$(typeof(alg))`. + To register the algorithm type for `right_null`, define + + MatrixAlgebraKit.right_null_alg(alg::CustomAlgorithm) = RightNullAlgorithm{kind}(alg) + + where `kind` selects the factorization type that will be used. + By default, this is either `:lq` or `:svd`, to select [`lq_null!`](@ref), + [`svd_compact!`](@ref) or [`svd_trunc!`](@ref) respectively. + """ +) +right_null_alg(alg::RightNullAlgorithm) = alg +right_null_alg(alg::LQAlgorithms) = RightNullViaLQ(alg) +right_null_alg(alg::SVDAlgorithms) = RightNullViaSVD(alg) +right_null_alg(alg::TruncatedAlgorithm{<:SVDAlgorithms}) = RightNullViaSVD(alg) diff --git a/src/interface/qr.jl b/src/interface/qr.jl index 7d5f8dfb..56624282 100644 --- a/src/interface/qr.jl +++ b/src/interface/qr.jl @@ -72,7 +72,7 @@ default_qr_algorithm(A; kwargs...) = default_qr_algorithm(typeof(A); kwargs...) function default_qr_algorithm(T::Type; kwargs...) throw(MethodError(default_qr_algorithm, (T,))) end -function default_qr_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.BlasMat} +function default_qr_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.MaybeBlasMat} return LAPACK_HouseholderQR(; kwargs...) end function default_qr_algorithm(::Type{T}; kwargs...) where {T <: Diagonal} diff --git a/src/interface/svd.jl b/src/interface/svd.jl index 606a1c4e..2ea26204 100644 --- a/src/interface/svd.jl +++ b/src/interface/svd.jl @@ -55,22 +55,28 @@ square diagonal matrix of size `(k, k)`, with `k` is set by the truncation strat The function also returns `ϵ`, the truncation error defined as the 2-norm of the discarded singular values. +## Truncation +The truncation strategy can be controlled via the `trunc` keyword argument. This can be +either a `NamedTuple` or a [`TruncationStrategy`](@ref). If `trunc` is not provided or +nothing, all values will be kept. + +### `trunc::NamedTuple` +The supported truncation keyword arguments are: + +$docs_truncation_kwargs + +### `trunc::TruncationStrategy` +For more control, a truncation strategy can be supplied directly. +By default, MatrixAlgebraKit supplies the following: + +$docs_truncation_strategies + ## Keyword arguments -The behavior of this function is controlled by the following keyword arguments: - -- `trunc`: Specifies the truncation strategy. This can be: - - A `NamedTuple` with fields `atol`, `rtol`, and/or `maxrank`, which will be converted to - a [`TruncationStrategy`](@ref). For details on available truncation strategies, see - [Truncations](@ref). - - A `TruncationStrategy` object directly (e.g., `truncrank(10)`, `trunctol(atol=1e-6)`, or - combinations using `&`). - - `nothing` (default), which keeps all singular values. - -- Other keyword arguments are passed to the algorithm selection procedure. If no explicit - `alg` is provided, these keywords are used to select and configure the algorithm through - [`MatrixAlgebraKit.select_algorithm`](@ref). The remaining keywords after algorithm - selection are passed to the algorithm constructor. See [`MatrixAlgebraKit.default_algorithm`](@ref) - for the default algorithm selection behavior. +Other keyword arguments are passed to the algorithm selection procedure. If no explicit +`alg` is provided, these keywords are used to select and configure the algorithm through +[`MatrixAlgebraKit.select_algorithm`](@ref). The remaining keywords after algorithm +selection are passed to the algorithm constructor. See [`MatrixAlgebraKit.default_algorithm`](@ref) +for the default algorithm selection behavior. When `alg` is a [`TruncatedAlgorithm`](@ref), the `trunc` keyword cannot be specified as the truncation strategy is already embedded in the algorithm. @@ -106,7 +112,7 @@ default_svd_algorithm(A; kwargs...) = default_svd_algorithm(typeof(A); kwargs... function default_svd_algorithm(T::Type; kwargs...) throw(MethodError(default_svd_algorithm, (T,))) end -function default_svd_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.BlasMat} +function default_svd_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.MaybeBlasMat} return LAPACK_DivideAndConquer(; kwargs...) end function default_svd_algorithm(::Type{T}; kwargs...) where {T <: Diagonal} diff --git a/src/interface/truncation.jl b/src/interface/truncation.jl index db417edb..a1e993de 100644 --- a/src/interface/truncation.jl +++ b/src/interface/truncation.jl @@ -1,3 +1,25 @@ +const docs_truncation_kwargs = """ +* `atol::Real` : Absolute tolerance for the truncation +* `rtol::Real` : Relative tolerance for the truncation +* `maxrank::Real` : Maximal rank for the truncation +* `maxerror::Real` : Maximal truncation error. +* `filter` : Custom filter to select truncated values. +""" + +const docs_truncation_strategies = """ +- [`notrunc`](@ref) +- [`truncrank`](@ref) +- [`trunctol`](@ref) +- [`truncerror`](@ref) +- [`truncfilter`](@ref) +""" + +const docs_null_truncation_kwargs = """ +* `atol::Real` : Absolute tolerance for the truncation +* `rtol::Real` : Relative tolerance for the truncation +* `maxnullity::Real` : Maximal rank for the truncation +""" + """ TruncationStrategy(; kwargs...) @@ -8,11 +30,7 @@ The following keyword arguments are all optional, and their default value (`noth will be ignored. It is also allowed to combine multiple of these, in which case the kept values will consist of the intersection of the different truncated strategies. -- `atol::Real` : Absolute tolerance for the truncation -- `rtol::Real` : Relative tolerance for the truncation -- `maxrank::Real` : Maximal rank for the truncation -- `maxerror::Real` : Maximal truncation error. -- `filter` : Custom filter to select truncated values. +$docs_truncation_kwargs """ function TruncationStrategy(; atol::Union{Real, Nothing} = nothing, @@ -36,6 +54,28 @@ function TruncationStrategy(; return strategy end +""" + null_truncation_strategy(; kwargs...) + +Select a nullspace truncation strategy based on the provided keyword arguments. + +## Keyword arguments +The following keyword arguments are all optional, and their default value (`nothing`) +will be ignored. It is also allowed to combine multiple of these, in which case the +discarded values will consist of the intersection of the different truncated strategies. + +$docs_null_truncation_kwargs +""" +function null_truncation_strategy(; atol = nothing, rtol = nothing, maxnullity = nothing) + if isnothing(maxnullity) && isnothing(atol) && isnothing(rtol) + return notrunc() + end + atol = @something atol 0 + rtol = @something rtol 0 + trunc = trunctol(; atol, rtol, keep_below = true) + return !isnothing(maxnullity) ? trunc & truncrank(maxnullity; rev = false) : trunc +end + """ NoTruncation() diff --git a/src/yalapack.jl b/src/yalapack.jl index 5bc87073..b2d4c27e 100644 --- a/src/yalapack.jl +++ b/src/yalapack.jl @@ -17,6 +17,8 @@ using LinearAlgebra.LAPACK: chkfinite, chktrans, chkside, chkuplofinite, chklapa # type alias for matrices that are definitely supported by YALAPACK const BlasMat{T <: BlasFloat} = StridedMatrix{T} +# type alias for matrices that are possibly supported by YALAPACK, after conversion +const MaybeBlasMat = Union{BlasMat, AbstractMatrix{<:Integer}} # LU factorisation for (getrf, getrs, elty) in ( diff --git a/test/amd/orthnull.jl b/test/amd/orthnull.jl index 97f5d915..6ed44228 100644 --- a/test/amd/orthnull.jl +++ b/test/amd/orthnull.jl @@ -10,7 +10,9 @@ 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) +eltypes = (Float32, Float64, ComplexF32, ComplexF64) + +@testset "left_orth and left_null for T = $T" for T in eltypes rng = StableRNG(123) m = 54 @testset for n in (37, m, 63) @@ -30,7 +32,7 @@ include(joinpath("..", "linearmap.jl")) @test hV * hV' + hN * hN' ≈ I M = LinearMap(A) - VM, CM = @constinferred left_orth(M; kind = :svd) + VM, CM = @constinferred left_orth(M; alg = :svd) @test parent(VM) * parent(CM) ≈ A if m > n @@ -48,20 +50,33 @@ include(joinpath("..", "linearmap.jl")) @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 + # passing a kind and some kwargs + V, C = @constinferred left_orth(A; alg = :qr, positive = true) + N = @constinferred left_null(A; alg = :qr, positive = true) + @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 LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + hV = collect(V) + hN = collect(N) + @test hV * hV' + hN * hN' ≈ I + + # passing an algorithm + V, C = @constinferred left_orth(A; alg = CUSOLVER_HouseholderQR()) + N = @constinferred left_null(A; alg = :qr, positive = true) + @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 LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + hV = collect(V) + hN = collect(N) + @test hV * hV' + hN * hN' ≈ I Ac = similar(A) V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C)) @@ -82,9 +97,6 @@ include(joinpath("..", "linearmap.jl")) 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) @@ -102,9 +114,6 @@ include(joinpath("..", "linearmap.jl")) 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) @@ -114,16 +123,13 @@ include(joinpath("..", "linearmap.jl")) @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 + @testset for alg in (:qr, :polar, :svd) # explicit alg kwarg + m < n && alg == :polar && continue + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); alg = $(QuoteNode(alg))) @test V2 * C2 ≈ A @test isisometric(V2) - if kind != :polar - N2 = @constinferred left_null!(copy!(Ac, A), N; kind = kind) - @test N2 === N + if alg != :polar + N2 = @constinferred left_null!(copy!(Ac, A), N; alg) @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @test isisometric(N2) hV2 = collect(V2) @@ -131,42 +137,24 @@ include(joinpath("..", "linearmap.jl")) @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) - ) + # with alg and tol kwargs + if alg == :svd + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); alg = $(QuoteNode(alg)), trunc = (; atol)) AMDGPU.@allowscalar begin - N2 = @constinferred left_null!( - copy!(Ac, A), N; kind = kind, - trunc = (; atol = atol) - ) + N2 = @constinferred left_null!(copy!(Ac, A), N; alg, trunc = (; atol)) end - @test V2 !== V - @test C2 !== C - @test N2 !== C @test V2 * C2 ≈ A - @test V2' * V2 ≈ I + @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 - V2, C2 = @constinferred left_orth!( - copy!(Ac, A), (V, C); kind = kind, - trunc = (; rtol = rtol) - ) + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); alg = $(QuoteNode(alg)), trunc = (; rtol)) AMDGPU.@allowscalar begin - N2 = @constinferred left_null!( - copy!(Ac, A), N; kind = kind, - trunc = (; rtol = rtol) - ) + N2 = @constinferred left_null!(copy!(Ac, A), N; alg, trunc = (; 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) @@ -175,31 +163,17 @@ include(joinpath("..", "linearmap.jl")) 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) - ) + @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; atol)) + @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; rtol)) + alg == :polar && continue + @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; atol)) + @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; rtol)) end end end end -@testset "right_orth and right_null for T = $T" for T in ( - Float32, Float64, ComplexF32, - ComplexF64, - ) +@testset "right_orth and right_null for T = $T" for T in eltypes rng = StableRNG(123) m = 54 @testset for n in (37, m, 63) @@ -219,15 +193,12 @@ end @test hVᴴ' * hVᴴ + hNᴴ' * hNᴴ ≈ I M = LinearMap(A) - CM, VMᴴ = @constinferred right_orth(M; kind = :svd) + CM, VMᴴ = @constinferred right_orth(M; alg = :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) @@ -239,10 +210,9 @@ end 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ᴴ + AMDGPU.@allowscalar begin + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; atol = atol)) + end @test C2 * Vᴴ2 ≈ A @test isisometric(Vᴴ2; side = :right) @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @@ -252,10 +222,9 @@ end @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ᴴ + AMDGPU.@allowscalar begin + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; rtol = rtol)) + end @test C2 * Vᴴ2 ≈ A @test isisometric(Vᴴ2; side = :right) @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @@ -264,16 +233,14 @@ end 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ᴴ + + @testset "alg = $alg" for alg in (:lq, :polar, :svd) + n < m && alg == :polar && continue + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); alg = $(QuoteNode(alg))) @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ᴴ + if alg != :polar + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; alg = $(QuoteNode(alg))) @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @test isisometric(Nᴴ2; side = :right) hVᴴ2 = collect(Vᴴ2) @@ -281,18 +248,11 @@ end @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ᴴ + if alg == :svd + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); alg = $(QuoteNode(alg)), trunc = (; atol)) + AMDGPU.@allowscalar begin + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; alg = $(QuoteNode(alg)), trunc = (; atol)) + end @test C2 * Vᴴ2 ≈ A @test isisometric(Vᴴ2; side = :right) @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @@ -301,42 +261,25 @@ end 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ᴴ + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); alg = $(QuoteNode(alg)), trunc = (; rtol)) + AMDGPU.@allowscalar begin + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; alg = $(QuoteNode(alg)), trunc = (; rtol)) + end @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) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I 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) - ) + @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; atol)) + @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; rtol)) + alg == :polar && continue + @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; atol)) + @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; rtol)) end end + end end diff --git a/test/chainrules.jl b/test/chainrules.jl index 441a17fb..ba3f0681 100644 --- a/test/chainrules.jl +++ b/test/chainrules.jl @@ -541,18 +541,18 @@ end ) test_rrule( config, left_orth, A; - fkwargs = (; kind = :qr), atol = atol, rtol = rtol, rrule_f = rrule_via_ad, check_inferred = false + fkwargs = (; alg = :qr), atol = atol, rtol = rtol, rrule_f = rrule_via_ad, check_inferred = false ) m >= n && test_rrule( config, left_orth, A; - fkwargs = (; kind = :polar), atol = atol, rtol = rtol, rrule_f = rrule_via_ad, check_inferred = false + fkwargs = (; alg = :polar), atol = atol, rtol = rtol, rrule_f = rrule_via_ad, check_inferred = false ) - ΔN = left_orth(A; kind = :qr)[1] * randn(rng, T, min(m, n), m - min(m, n)) + ΔN = left_orth(A; alg = :qr)[1] * randn(rng, T, min(m, n), m - min(m, n)) test_rrule( config, left_null, A; - fkwargs = (; kind = :qr), output_tangent = ΔN, atol = atol, rtol = rtol, + fkwargs = (; alg = :qr), output_tangent = ΔN, atol = atol, rtol = rtol, rrule_f = rrule_via_ad, check_inferred = false ) @@ -561,19 +561,19 @@ end atol = atol, rtol = rtol, rrule_f = rrule_via_ad, check_inferred = false ) test_rrule( - config, right_orth, A; fkwargs = (; kind = :lq), + config, right_orth, A; fkwargs = (; alg = :lq), atol = atol, rtol = rtol, rrule_f = rrule_via_ad, check_inferred = false ) m <= n && test_rrule( - config, right_orth, A; fkwargs = (; kind = :polar), + config, right_orth, A; fkwargs = (; alg = :polar), atol = atol, rtol = rtol, rrule_f = rrule_via_ad, check_inferred = false ) - ΔNᴴ = randn(rng, T, n - min(m, n), min(m, n)) * right_orth(A; kind = :lq)[2] + ΔNᴴ = randn(rng, T, n - min(m, n), min(m, n)) * right_orth(A; alg = :lq)[2] test_rrule( config, right_null, A; - fkwargs = (; kind = :lq), output_tangent = ΔNᴴ, + fkwargs = (; alg = :lq), output_tangent = ΔNᴴ, atol = atol, rtol = rtol, rrule_f = rrule_via_ad, check_inferred = false ) end diff --git a/test/cuda/orthnull.jl b/test/cuda/orthnull.jl index cfe7d9e3..2a2a26f6 100644 --- a/test/cuda/orthnull.jl +++ b/test/cuda/orthnull.jl @@ -10,7 +10,9 @@ 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) +eltypes = (Float32, Float64, ComplexF32, ComplexF64) + +@testset "left_orth and left_null for T = $T" for T in eltypes rng = StableRNG(123) m = 54 @testset for n in (37, m, 63) @@ -30,7 +32,7 @@ include(joinpath("..", "linearmap.jl")) @test hV * hV' + hN * hN' ≈ I M = LinearMap(A) - VM, CM = @constinferred left_orth(M; kind = :svd) + VM, CM = @constinferred left_orth(M; alg = :svd) @test parent(VM) * parent(CM) ≈ A if m > n @@ -48,27 +50,37 @@ include(joinpath("..", "linearmap.jl")) @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 + # passing a kind and some kwargs + V, C = @constinferred left_orth(A; alg = :qr, positive = true) + N = @constinferred left_null(A; alg = :qr, positive = true) + @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 LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + hV = collect(V) + hN = collect(N) + @test hV * hV' + hN * hN' ≈ I + + # passing an algorithm + V, C = @constinferred left_orth(A; alg = CUSOLVER_HouseholderQR()) + N = @constinferred left_null(A; alg = :qr, positive = true) + @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 LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + hV = collect(V) + hN = collect(N) + @test hV * hV' + hN * hN' ≈ I 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) @@ -80,9 +92,6 @@ include(joinpath("..", "linearmap.jl")) 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) @@ -98,9 +107,6 @@ include(joinpath("..", "linearmap.jl")) ) 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) @@ -110,16 +116,13 @@ include(joinpath("..", "linearmap.jl")) @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 + @testset for alg in (:qr, :polar, :svd) # explicit alg kwarg + m < n && alg == :polar && continue + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); alg = $(QuoteNode(alg))) @test V2 * C2 ≈ A @test isisometric(V2) - if kind != :polar - N2 = @constinferred left_null!(copy!(Ac, A), N; kind = kind) - @test N2 === N + if alg != :polar + N2 = @constinferred left_null!(copy!(Ac, A), N; alg) @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @test isisometric(N2) hV2 = collect(V2) @@ -127,38 +130,20 @@ include(joinpath("..", "linearmap.jl")) @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 + # with alg and tol kwargs + if alg == :svd + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); alg = $(QuoteNode(alg)), trunc = (; atol)) + N2 = @constinferred left_null!(copy!(Ac, A), N; alg, trunc = (; atol)) @test V2 * C2 ≈ A - @test V2' * V2 ≈ I + @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 - 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 + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); alg = $(QuoteNode(alg)), trunc = (; rtol)) + N2 = @constinferred left_null!(copy!(Ac, A), N; alg, trunc = (; rtol)) @test V2 * C2 ≈ A @test isisometric(V2) @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @@ -167,31 +152,17 @@ include(joinpath("..", "linearmap.jl")) 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) - ) + @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; atol)) + @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; rtol)) + alg == :polar && continue + @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; atol)) + @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; rtol)) end end end end -@testset "right_orth and right_null for T = $T" for T in ( - Float32, Float64, ComplexF32, - ComplexF64, - ) +@testset "right_orth and right_null for T = $T" for T in eltypes rng = StableRNG(123) m = 54 @testset for n in (37, m, 63) @@ -211,15 +182,12 @@ end @test hVᴴ' * hVᴴ + hNᴴ' * hNᴴ ≈ I M = LinearMap(A) - CM, VMᴴ = @constinferred right_orth(M; kind = :svd) + CM, VMᴴ = @constinferred right_orth(M; alg = :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) @@ -232,9 +200,6 @@ end 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) @@ -245,9 +210,6 @@ end 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) @@ -256,16 +218,13 @@ end 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ᴴ + @testset "alg = $alg" for alg in (:lq, :polar, :svd) + n < m && alg == :polar && continue + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); alg = $(QuoteNode(alg))) @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ᴴ + if alg != :polar + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; alg = $(QuoteNode(alg))) @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @test isisometric(Nᴴ2; side = :right) hVᴴ2 = collect(Vᴴ2) @@ -273,18 +232,9 @@ end @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ᴴ + if alg == :svd + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); alg = $(QuoteNode(alg)), trunc = (; atol)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; alg = $(QuoteNode(alg)), trunc = (; atol)) @test C2 * Vᴴ2 ≈ A @test isisometric(Vᴴ2; side = :right) @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @@ -293,41 +243,21 @@ end 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ᴴ + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); alg = $(QuoteNode(alg)), trunc = (; rtol)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; alg = $(QuoteNode(alg)), trunc = (; rtol)) @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) + @test hVᴴ2' * hVᴴ2 + hNᴴ2' * hNᴴ2 ≈ I 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) - ) + @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; atol)) + @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; rtol)) + alg == :polar && continue + @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; atol)) + @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; rtol)) end end end diff --git a/test/linearmap.jl b/test/linearmap.jl index 61753cde..bbf84e5e 100644 --- a/test/linearmap.jl +++ b/test/linearmap.jl @@ -34,13 +34,6 @@ module LinearMaps 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...) diff --git a/test/orthnull.jl b/test/orthnull.jl index cee54252..ce742e8f 100644 --- a/test/orthnull.jl +++ b/test/orthnull.jl @@ -2,15 +2,12 @@ using MatrixAlgebraKit using Test using TestExtras using StableRNGs -using LinearAlgebra: LinearAlgebra, I, mul! -using MatrixAlgebraKit: LAPACK_SVDAlgorithm, check_input, copy_input, default_svd_algorithm, - initialize_output, AbstractAlgorithm +using LinearAlgebra: LinearAlgebra, I # testing non-AbstractArray codepaths: include("linearmap.jl") eltypes = (Float32, Float64, ComplexF32, ComplexF64) - @testset "left_orth and left_null for T = $T" for T in eltypes rng = StableRNG(123) m = 54 @@ -29,7 +26,7 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64) @test V * V' + N * N' ≈ I M = LinearMap(A) - VM, CM = @constinferred left_orth(M; kind = :svd) + VM, CM = @constinferred left_orth(M; alg = :svd) @test parent(VM) * parent(CM) ≈ A if m > n @@ -45,25 +42,33 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64) @test isisometric(N) end - for alg_qr in ((; positive = true), (; positive = false), LAPACK_HouseholderQR()) - V, C = @constinferred left_orth(A; alg_qr) - N = @constinferred left_null(A; alg_qr) - @test V isa Matrix{T} && size(V) == (m, minmn) - @test C isa Matrix{T} && size(C) == (minmn, n) - @test N isa Matrix{T} && size(N) == (m, m - minmn) - @test V * C ≈ A - @test isisometric(V) - @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) - @test isisometric(N) - @test V * V' + N * N' ≈ I - end + # passing a kind and some kwargs + V, C = @constinferred left_orth(A; alg = :qr, positive = true) + N = @constinferred left_null(A; alg = :qr, positive = true) + @test V isa Matrix{T} && size(V) == (m, minmn) + @test C isa Matrix{T} && size(C) == (minmn, n) + @test N isa Matrix{T} && size(N) == (m, m - minmn) + @test V * C ≈ A + @test isisometric(V) + @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + @test V * V' + N * N' ≈ I + + # passing an algorithm + V, C = @constinferred left_orth(A; alg = LAPACK_HouseholderQR()) + N = @constinferred left_null(A; alg = :qr, positive = true) + @test V isa Matrix{T} && size(V) == (m, minmn) + @test C isa Matrix{T} && size(C) == (minmn, n) + @test N isa Matrix{T} && size(N) == (m, m - minmn) + @test V * C ≈ A + @test isisometric(V) + @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + @test V * V' + N * N' ≈ I 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) @@ -73,9 +78,6 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64) 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) @@ -89,9 +91,6 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64) ) 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) @@ -99,49 +98,41 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64) @test V2 * V2' + N2 * N2' ≈ I end - 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 + for alg in (:qr, :polar, :svd) # explicit kind kwarg + m < n && alg === :polar && continue + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); alg = $(QuoteNode(alg))) @test V2 * C2 ≈ A @test isisometric(V2) - if kind != :polar - N2 = @constinferred left_null!(copy!(Ac, A), N; kind = kind) - @test N2 === N + if alg != :polar + N2 = @constinferred left_null!(copy!(Ac, A), N; alg) @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @test isisometric(N2) @test V2 * V2' + N2 * N2' ≈ I end # with kind and tol kwargs - if kind == :svd - V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); kind, trunc = (; atol)) - N2 = @constinferred left_null!(copy!(Ac, A), N; kind, trunc = (; atol)) - @test V2 !== V - @test C2 !== C - @test N2 !== C + if alg == :svd + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); alg = $(QuoteNode(alg)), trunc = (; atol)) + N2 = @constinferred left_null!(copy!(Ac, A), N; alg, trunc = (; atol)) @test V2 * C2 ≈ A - @test V2' * V2 ≈ I + @test isisometric(V2) @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) - @test N2' * N2 ≈ I + @test isisometric(N2) @test V2 * V2' + N2 * N2' ≈ I - V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); kind, trunc = (; rtol)) - N2 = @constinferred left_null!(copy!(Ac, A), N; kind, trunc = (; rtol)) - @test V2 !== V - @test C2 !== C - @test N2 !== C + V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); alg = $(QuoteNode(alg)), trunc = (; rtol)) + N2 = @constinferred left_null!(copy!(Ac, A), N; alg, trunc = (; rtol)) @test V2 * C2 ≈ A @test isisometric(V2) @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @test isisometric(N2) @test V2 * V2' + N2 * N2' ≈ I else - @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); kind, trunc = (; atol)) - @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); kind, trunc = (; rtol)) - @test_throws ArgumentError left_null!(copy!(Ac, A), N; kind, trunc = (; atol)) - @test_throws ArgumentError left_null!(copy!(Ac, A), N; kind, trunc = (; rtol)) + @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; atol)) + @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; rtol)) + alg == :polar && continue + @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; atol)) + @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; rtol)) end end end @@ -165,15 +156,12 @@ end @test Vᴴ' * Vᴴ + Nᴴ' * Nᴴ ≈ I M = LinearMap(A) - CM, VMᴴ = @constinferred right_orth(M; kind = :svd) + CM, VMᴴ = @constinferred right_orth(M; alg = :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) @@ -183,9 +171,6 @@ end atol = eps(real(T)) C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; atol)) Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; 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) @@ -195,57 +180,46 @@ end rtol = eps(real(T)) C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; rtol)) Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; 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) @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I - for kind in (:lq, :polar, :svd) - n < m && kind == :polar && continue - C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); kind) - @test C2 === C - @test Vᴴ2 === Vᴴ + for alg in (:lq, :polar, :svd) + n < m && alg == :polar && continue + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); alg = $(QuoteNode(alg))) @test C2 * Vᴴ2 ≈ A @test isisometric(Vᴴ2; side = :right) - if kind != :polar - Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind) - @test Nᴴ2 === Nᴴ + if alg != :polar + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; alg = $(QuoteNode(alg))) @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) @test isisometric(Nᴴ2; side = :right) @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I end - if kind == :svd - C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); kind, trunc = (; atol)) - Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind, trunc = (; atol)) - @test C2 !== C - @test Vᴴ2 !== Vᴴ - @test Nᴴ2 !== Nᴴ + if alg == :svd + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); alg = $(QuoteNode(alg)), trunc = (; atol)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; alg = $(QuoteNode(alg)), trunc = (; atol)) @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) @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I - C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); kind, trunc = (; rtol)) - Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind, trunc = (; rtol)) - @test C2 !== C - @test Vᴴ2 !== Vᴴ - @test Nᴴ2 !== Nᴴ + C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); alg = $(QuoteNode(alg)), trunc = (; rtol)) + Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; alg = $(QuoteNode(alg)), trunc = (; rtol)) @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) @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I else - @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); kind, trunc = (; atol)) - @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); kind, trunc = (; rtol)) - @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; kind, trunc = (; atol)) - @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; kind, trunc = (; rtol)) + @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; atol)) + @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; rtol)) + alg == :polar && continue + @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; atol)) + @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; rtol)) end end end