diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index 2178d41a..127c6ae2 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -33,7 +33,7 @@ export left_orth!, right_orth!, left_null!, right_null! export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert, LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations, - LAPACK_DivideAndConquer, LAPACK_Jacobi + LAPACK_DivideAndConquer, LAPACK_Jacobi, SafeSVD export GLA_HouseholderQR, GLA_QRIteration, GS_QRIteration export LQViaTransposedQR export PolarViaSVD, PolarNewton diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index 4263c592..1b5c06f6 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -131,6 +131,16 @@ function svd_full!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm) isempty(alg_kwargs) || throw(ArgumentError("invalid keyword arguments for LAPACK_DivideAndConquer")) YALAPACK.gesdd!(A, view(S, 1:minmn, 1), U, Vᴴ) + elseif alg isa SafeSVD + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for SafeSVD")) + # extra copy to avoid modifying input if sdd goes wrong + A′ = copy(A) + try + YALAPACK.gesdd!(A′, view(S, 1:minmn, 1), U, Vᴴ) + catch + YALAPACK.gesvd!(A, view(S, 1:minmn, 1), U, Vᴴ) + end elseif alg isa LAPACK_Bisection throw(ArgumentError("LAPACK_Bisection is not supported for full SVD")) elseif alg isa LAPACK_Jacobi @@ -164,6 +174,16 @@ function svd_compact!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm) isempty(alg_kwargs) || throw(ArgumentError("invalid keyword arguments for LAPACK_DivideAndConquer")) YALAPACK.gesdd!(A, diagview(S), U, Vᴴ) + elseif alg isa SafeSVD + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for SafeSVD")) + # extra copy to avoid modifying input if sdd goes wrong + A′ = copy(A) + try + YALAPACK.gesdd!(A′, diagview(S), U, Vᴴ) + catch + YALAPACK.gesvd!(A, diagview(S), U, Vᴴ) + end elseif alg isa LAPACK_Bisection YALAPACK.gesvdx!(A, diagview(S), U, Vᴴ; alg_kwargs...) elseif alg isa LAPACK_Jacobi diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index 0fc7403d..5dbffbb0 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -153,11 +153,24 @@ see also [`gaugefix!`](@ref). """ @algdef LAPACK_Jacobi +""" + SafeSVD(; fixgauge::Bool = true) + +Algorithm type to denote a combination of LAPACK_DivideAndConquer and LAPACK_QRIteration +for computing the singular value decomposition of a general matrix. +This algorithm first tries to perform the SVD using the faster LAPACK_DivideAndConquer method, +and falls back to the more stable LAPACK_QRIteration method if numerical issues are detected. +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the singular vectors, +see also [`gaugefix!`](@ref), [`LAPACK_DivideAndConquer`](@ref) and [`LAPACK_QRIteration`](@ref). +""" +@algdef SafeSVD + const LAPACK_SVDAlgorithm = Union{ LAPACK_QRIteration, LAPACK_Bisection, LAPACK_DivideAndConquer, LAPACK_Jacobi, + SafeSVD, } # =========================