rework tests and more fixes and cleanup

Author: lkdvos

src/TensorKit.jl

@@ -70,28 +70,30 @@ export inner, dot, norm, normalize, normalize!, tr
 
 # factorizations
 export mul!, lmul!, rmul!, adjoint!, pinv, axpy!, axpby!
-export leftorth, rightorth, leftnull, rightnull,
-       leftorth!, rightorth!, leftnull!, rightnull!,
+export left_orth, right_orth, left_null, right_null,
+       left_orth!, right_orth!, left_null!, right_null!,
        left_polar, left_polar!, right_polar, right_polar!,
        qr_full, qr_compact, qr_null, lq_full, lq_compact, lq_null,
        qr_full!, qr_compact!, qr_null!, lq_full!, lq_compact!, lq_null!,
-       tsvd!, tsvd, eigen, eigen!, eig, eig!, eigh, eigh!, exp, exp!,
-       eigh_full!, eigh_full, eig_full!, eig_full, eigh_vals!, eigh_vals,
-       eig_vals!, eig_vals,
+       svd_compact!, svd_full!, svd_trunc!, svd_compact, svd_full, svd_trunc,
+       exp, exp!,
+       eigh_full!, eigh_full, eigh_trunc!, eigh_trunc, eig_full!, eig_full, eig_trunc!,
+       eig_trunc,
+       eigh_vals!, eigh_vals, eig_vals!, eig_vals,
        isposdef, isposdef!, ishermitian, isisometry, isunitary, sylvester, rank, cond
+# deprecate:
+export eig, eig!, eigh, eigh!, eigen, eigen!, tsvd, tsvd!, leftorth, leftorth!, rightorth,
+       rightorth!, leftnull, leftnull!, rightnull, rightnull!
 export braid, braid!, permute, permute!, transpose, transpose!, twist, twist!, repartition,
        repartition!
 export catdomain, catcodomain, absorb, absorb!
 
-export OrthogonalFactorizationAlgorithm, QR, QRpos, QL, QLpos, LQ, LQpos, RQ, RQpos,
-       SVD, SDD, Polar
-
 # tensor operations
 export @tensor, @tensoropt, @ncon, ncon, @planar, @plansor
 export scalar, add!, contract!
 
 # truncation schemes
-export notrunc, truncerr, truncdim, truncspace, truncbelow
+export notrunc, truncerr, truncrank, truncspace, trunctol
 
 # cache management
 export empty_globalcaches!

src/auxiliary/deprecate.jl

@@ -1,30 +1,5 @@
 import Base: transpose
 
-#! format: off
-# Base.@deprecate(permute(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; copy::Bool=false),
-#                 permute(t, (p1, p2); copy=copy))
-# Base.@deprecate(transpose(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; copy::Bool=false),
-#                 transpose(t, (p1, p2); copy=copy))
-# Base.@deprecate(braid(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple, levels; copy::Bool=false),
-#                 braid(t, (p1, p2), levels; copy=copy))
-
-# Base.@deprecate(tsvd(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...),
-#                 tsvd(t, (p₁, p₂); kwargs...))
-# Base.@deprecate(leftorth(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...),
-#                 leftorth(t, (p₁, p₂); kwargs...))
-# Base.@deprecate(rightorth(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...),
-#                 rightorth(t, (p₁, p₂); kwargs...))
-# Base.@deprecate(leftnull(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...),
-#                 leftnull(t, (p₁, p₂); kwargs...))
-# Base.@deprecate(rightnull(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...),
-#                 rightnull(t, (p₁, p₂); kwargs...))
-# Base.@deprecate(LinearAlgebra.eigen(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...),
-#                 LinearAlgebra.eigen(t, (p₁, p₂); kwargs...), false)
-# Base.@deprecate(eig(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...),
-#                 eig(t, (p₁, p₂); kwargs...))
-# Base.@deprecate(eigh(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...),
-#                 eigh(t, (p₁, p₂); kwargs...))
-
 for f in (:rand, :randn, :zeros, :ones)
     @eval begin
         Base.@deprecate TensorMap(::typeof($f), T::Type, P::HomSpace) $f(T, P)

src/tensors/factorizations/deprecations.jl

@@ -23,58 +23,97 @@ const TruncationScheme = TruncationStrategy
 
 # factorizations
 # --------------
+_kindof(::LAPACK_HouseholderQR) = :qr
+_kindof(::LAPACK_HouseholderLQ) = :lq
+_kindof(::LAPACK_SVDAlgorithm) = :svd
+_kindof(::PolarViaSVD) = :polar
+
+_drop_alg(; alg=nothing, kwargs...) = kwargs
+_drop_p(; p=nothing, kwargs...) = kwargs
+
 # orthogonalization
-@deprecate(leftorth(t::AbstractTensorMap, p::Index2Tuple; kwargs...),
-           leftorth!(permutedcopy_oftype(t, factorization_scalartype(leftorth, t), p);
-                     kwargs...))
-@deprecate(rightorth(t::AbstractTensorMap, p::Index2Tuple; kwargs...),
-           rightorth!(permutedcopy_oftype(t, factorisation_scalartype(rightorth, t), p);
-                      kwargs...))
-function leftorth(t::AbstractTensorMap; kwargs...)
+export leftorth, leftorth!, rightorth, rightorth!
+function leftorth(t::AbstractTensorMap, p::Index2Tuple; kwargs...)
     Base.depwarn("`leftorth` is no longer supported, use `left_orth` instead", :leftorth)
-    return left_orth(t; kwargs...)
+    return leftorth!(permutedcopy_oftype(t, factorisation_scalartype(leftorth, t), p);
+                     kwargs...)
 end
-function leftorth!(t::AbstractTensorMap; kwargs...)
-    Base.depwarn("`leftorth!` is no longer supported, use `left_orth!` instead", :leftorth!)
-    return left_orth!(t; kwargs...)
+function rightorth(t::AbstractTensorMap, p::Index2Tuple; kwargs...)
+    Base.depwarn("`rightorth` is no longer supported, use `right_orth` instead", :rightorth)
+    return rightorth!(permutedcopy_oftype(t, factorisation_scalartype(rightorth, t), p);
+                      kwargs...)
+end
+function leftorth(t::AbstractTensorMap; kwargs...)
+    Base.depwarn("`leftorth` is no longer supported, use `left_orth` instead", :leftorth)
+    return leftorth!(copy_oftype(t, factorisation_scalartype(leftorth, t)); kwargs...)
 end
 function rightorth(t::AbstractTensorMap; kwargs...)
     Base.depwarn("`rightorth` is no longer supported, use `right_orth` instead", :rightorth)
-    return right_orth(t; kwargs...)
+    return rightorth!(copy_oftype(t, factorisation_scalartype(rightorth, t)); kwargs...)
+end
+function leftorth!(t::AbstractTensorMap; kwargs...)
+    Base.depwarn("`leftorth!` is no longer supported, use `left_orth!` instead", :leftorth!)
+    haskey(kwargs, :alg) || return left_orth!(t; kwargs...)
+    alg = kwargs[:alg]
+    kind = _kindof(alg)
+    kind === :svd && return left_orth!(t; kind, alg_svd=alg, _drop_alg(; kwargs...)...)
+    kind === :qr && return left_orth!(t; kind, alg_qr=alg, _drop_alg(; kwargs...)...)
+    kind === :polar && return left_orth!(t; kind, alg_polar=alg, _drop_alg(; kwargs...)...)
+    throw(ArgumentError("invalid leftorth kind"))
 end
 function rightorth!(t::AbstractTensorMap; kwargs...)
     Base.depwarn("`rightorth!` is no longer supported, use `right_orth!` instead",
                  :rightorth!)
-    return right_orth!(t; kwargs...)
+    haskey(kwargs, :alg) || return right_orth!(t; kwargs...)
+    alg = kwargs[:alg]
+    kind = _kindof(alg)
+    kind === :svd && return right_orth!(t; kind, alg_svd=alg, _drop_alg(; kwargs...)...)
+    kind === :lq && return right_orth!(t; kind, alg_lq=alg, _drop_alg(; kwargs...)...)
+    kind === :polar && return right_orth!(t; kind, alg_polar=alg, _drop_alg(; kwargs...)...)
+    throw(ArgumentError("invalid rightorth kind"))
 end
 
 # nullspaces
-@deprecate(leftnull(t::AbstractTensorMap, p::Index2Tuple; kwargs...),
-           leftnull!(permutedcopy_oftype(t, factorization_scalartype(leftnull, t), p);
-                     kwargs...))
-@deprecate(rightnull(t::AbstractTensorMap, p::Index2Tuple; kwargs...),
-           rightnull!(permutedcopy_oftype(t, factorisation_scalartype(rightnull, t), p);
-                      kwargs...))
+export leftnull, leftnull!, rightnull, rightnull!
 function leftnull(t::AbstractTensorMap; kwargs...)
     Base.depwarn("`leftnull` is no longer supported, use `left_null` instead", :leftnull)
-    return left_null(t; kwargs...)
+    return leftnull!(copy_oftype(t, factorisation_scalartype(leftnull, t)); kwargs...)
 end
-function leftnull!(t::AbstractTensorMap; kwargs...)
-    Base.depwarn("`left_null!` is no longer supported, use `left_null!` instead",
-                 :leftnull!)
-    return left_null!(t; kwargs...)
+function leftnull(t::AbstractTensorMap, p::Index2Tuple; kwargs...)
+    Base.depwarn("`leftnull` is no longer supported, use `left_null` instead", :leftnull)
+    return leftnull!(permutedcopy_oftype(t, factorisation_scalartype(leftnull, t), p); kwargs...)
 end
 function rightnull(t::AbstractTensorMap; kwargs...)
     Base.depwarn("`rightnull` is no longer supported, use `right_null` instead", :rightnull)
-    return right_null(t; kwargs...)
+    return rightnull!(copy_oftype(t, factorisation_scalartype(rightnull, t)); kwargs...)
+end
+function rightnull(t::AbstractTensorMap, p::Index2Tuple; kwargs...)
+    Base.depwarn("`rightnull` is no longer supported, use `right_null` instead", :rightnull)
+    return rightnull!(permutedcopy_oftype(t, factorisation_scalartype(rightnull, t), p); kwargs...)
+end
+function leftnull!(t::AbstractTensorMap; kwargs...)
+    Base.depwarn("`left_null!` is no longer supported, use `left_null!` instead",
+                 :leftnull!)
+    haskey(kwargs, :alg) || return left_null!(t; kwargs...)
+    alg = kwargs[:alg]
+    kind = _kindof(alg)
+    kind === :svd && return left_null!(t; kind, alg_svd=alg, _drop_alg(; kwargs...)...)
+    kind === :qr && return left_null!(t; kind, alg_qr=alg, _drop_alg(; kwargs...)...)
+    throw(ArgumentError("invalid leftnull kind"))
 end
 function rightnull!(t::AbstractTensorMap; kwargs...)
     Base.depwarn("`rightnull!` is no longer supported, use `right_null!` instead",
                  :rightnull!)
-    return right_null!(t; kwargs...)
+    haskey(kwargs, :alg) || return right_null!(t; kwargs...)
+    alg = kwargs[:alg]
+    kind = _kindof(alg)
+    kind === :svd && return right_null!(t; kind, alg_svd=alg, _drop_alg(; kwargs...)...)
+    kind === :lq && return right_null!(t; kind, alg_lq=alg, _drop_alg(; kwargs...)...)
+    throw(ArgumentError("invalid rightnull kind"))
 end
 
 # eigen values
+export eig, eig!, eigh, eigh!, eigen, eigen!
 @deprecate(eig(t::AbstractTensorMap, p::Index2Tuple; kwargs...),
            eig!(permutedcopy_oftype(t, factorisation_scalartype(eig, t), p); kwargs...))
 @deprecate(eigh(t::AbstractTensorMap, p::Index2Tuple; kwargs...),
@@ -103,7 +142,7 @@ function eigh!(t::AbstractTensorMap; kwargs...)
 end
 
 # singular values
-_drop_p(; p=nothing, kwargs...) = kwargs
+export tsvd, tsvd!
 @deprecate(tsvd(t::AbstractTensorMap, p::Index2Tuple; kwargs...),
            tsvd!(permutedcopy_oftype(t, factorisation_scalartype(tsvd, t), p); kwargs...))
 function tsvd(t::AbstractTensorMap; kwargs...)

src/tensors/factorizations/factorizations.jl

@@ -65,6 +65,13 @@ const RealOrComplexFloat = Union{AbstractFloat,Complex{<:AbstractFloat}}
 # LinearAlgebra overloads
 #------------------------------#
 
+function eigen(t::AbstractTensorMap; kwargs...)
+    return ishermitian(t) ? eigh_full(t; kwargs...) : eig_full(t; kwargs...)
+end
+function eigen!(t::AbstractTensorMap; kwargs...)
+    return ishermitian(t) ? eigh_full!(t; kwargs...) : eig_full!(t; kwargs...)
+end
+
 function LinearAlgebra.eigvals(t::AbstractTensorMap; kwargs...)
     tcopy = copy_oftype(t, factorisation_scalartype(eigen, t))
     return LinearAlgebra.eigvals!(tcopy; kwargs...)
@@ -89,6 +96,9 @@ function LinearAlgebra.ishermitian(t::AbstractTensorMap)
     return true
 end
 
+function LinearAlgebra.isposdef(t::AbstractTensorMap)
+    return isposdef!(copy_oftype(t, factorisation_scalartype(isposdef, t)))
+end
 function LinearAlgebra.isposdef!(t::AbstractTensorMap)
     domain(t) == codomain(t) ||
         throw(SpaceMismatch("`isposdef` requires domain and codomain to be the same"))

src/tensors/factorizations/matrixalgebrakit.jl

@@ -1,13 +1,17 @@
 # Algorithm selection
 # -------------------
-for f! in
-    [:svd_compact!, :svd_full!, :svd_trunc!, :svd_vals!, :qr_compact!, :qr_full!, :qr_null!,
-     :lq_compact!, :lq_full!, :lq_null!, :eig_full!, :eig_trunc!, :eig_vals!, :eigh_full!,
-     :eigh_trunc!, :eigh_vals!, :left_polar!, :right_polar!]
+for f in
+    [:svd_compact, :svd_full, :svd_trunc, :svd_vals, :qr_compact, :qr_full, :qr_null,
+     :lq_compact, :lq_full, :lq_null, :eig_full, :eig_trunc, :eig_vals, :eigh_full,
+     :eigh_trunc, :eigh_vals, :left_polar, :right_polar]
+     f! = Symbol(f, :!)
     @eval function default_algorithm(::typeof($f!), ::Type{T};
                                      kwargs...) where {T<:AbstractTensorMap}
         return default_algorithm($f!, blocktype(T); kwargs...)
     end
+    @eval function copy_input(::typeof($f), t::AbstractTensorMap)
+        return copy_oftype(t, factorisation_scalartype($f, t))
+    end
 end
 
 function _select_truncation(f, ::AbstractTensorMap,
@@ -388,7 +392,7 @@ function check_input(::typeof(lq_full!), t::AbstractTensorMap, LQ, ::AbstractAlg
 
     # type checks
     @assert L isa AbstractTensorMap
-    @assert R isa AbstractTensorMap
+    @assert Q isa AbstractTensorMap
 
     # scalartype checks
     @check_scalar L t
@@ -407,7 +411,7 @@ function check_input(::typeof(lq_compact!), t::AbstractTensorMap, LQ, ::Abstract
 
     # type checks
     @assert L isa AbstractTensorMap
-    @assert R isa AbstractTensorMap
+    @assert Q isa AbstractTensorMap
 
     # scalartype checks
     @check_scalar L t
@@ -546,7 +550,7 @@ function initialize_output(::typeof(right_polar!), t::AbstractTensorMap,
                            ::AbstractAlgorithm)
     P = similar(t, codomain(t) ← codomain(t))
     Wᴴ = similar(t, space(t))
-    return Wᴴ, P
+    return P, Wᴴ
 end
 
 # Needed to get algorithm selection to behave

src/tensors/factorizations/truncation.jl

@@ -64,48 +64,29 @@ function truncate!(::typeof(left_null!),
     return Ũ
 end
 
-function truncate!(::typeof(eigh_trunc!),
-                   (D, V)::Tuple{AbstractTensorMap,AbstractTensorMap},
-                   strategy::TruncationStrategy)
-    ind = findtruncated(diagview(D), strategy)
-    V_truncated = spacetype(D)(c => length(I) for (c, I) in ind)
-
-    D̃ = DiagonalTensorMap{scalartype(D)}(undef, V_truncated)
-    for (c, b) in blocks(D̃)
-        I = get(ind, c, nothing)
-        @assert !isnothing(I)
-        copy!(b.diag, @view(block(D, c).diag[I]))
-    end
-
-    Ṽ = similar(V, V_truncated ← domain(V))
-    for (c, b) in blocks(Ṽ)
-        I = get(ind, c, nothing)
-        @assert !isnothing(I)
-        copy!(b, @view(block(V, c)[I, :]))
-    end
-
-    return D̃, Ṽ
-end
-function truncate!(::typeof(eig_trunc!), (D, V)::Tuple{AbstractTensorMap,AbstractTensorMap},
-                   strategy::TruncationStrategy)
-    ind = findtruncated(diagview(D), strategy)
-    V_truncated = spacetype(D)(c => length(I) for (c, I) in ind)
+for f! in (:eig_trunc!, :eigh_trunc!)
+    @eval function truncate!(::typeof($f!),
+                             (D, V)::Tuple{AbstractTensorMap,AbstractTensorMap},
+                             strategy::TruncationStrategy)
+        ind = findtruncated(diagview(D), strategy)
+        V_truncated = spacetype(D)(c => length(I) for (c, I) in ind)
+
+        D̃ = DiagonalTensorMap{scalartype(D)}(undef, V_truncated)
+        for (c, b) in blocks(D̃)
+            I = get(ind, c, nothing)
+            @assert !isnothing(I)
+            copy!(b.diag, @view(block(D, c).diag[I]))
+        end
 
-    D̃ = DiagonalTensorMap{scalartype(D)}(undef, V_truncated)
-    for (c, b) in blocks(D̃)
-        I = get(ind, c, nothing)
-        @assert !isnothing(I)
-        copy!(b.diag, @view(block(D, c).diag[I]))
-    end
+        Ṽ = similar(V, codomain(V) ← V_truncated)
+        for (c, b) in blocks(Ṽ)
+            I = get(ind, c, nothing)
+            @assert !isnothing(I)
+            copy!(b, @view(block(V, c)[:, I]))
+        end
 
-    Ṽ = similar(V, V_truncated ← domain(V))
-    for (c, b) in blocks(Ṽ)
-        I = get(ind, c, nothing)
-        @assert !isnothing(I)
-        copy!(b, @view(block(V, c)[I, :]))
+        return D̃, Ṽ
     end
-
-    return D̃, Ṽ
 end
 
 # Find truncation
@@ -201,7 +182,7 @@ function findtruncated(Sd::SectorDict, strategy::TruncationKeepSorted)
         truncdim[cmin] -= 1
         totaldim -= dim(cmin)
         if totaldim < strategy.howmany
-            truncdim[cmin] += 1
+            # truncdim[cmin] += 1
             break
         end
         if truncdim[cmin] == 0

src/tensors/factorizations/utility.jl

@@ -16,7 +16,7 @@ function permutedcopy_oftype(t::AbstractTensorMap, T::Type{<:Number}, p::Index2T
     return permute!(similar(t, T, permute(space(t), p)), t, p)
 end
 function copy_oftype(t::AbstractTensorMap, T::Type{<:Number})
-    return copy!(similar(t, T), t)
+    return copy!(similar(t, T, space(t)), t)
 end
 
 function _reverse!(t::AbstractTensorMap; dims=:)

src/tensors/linalg.jl

@@ -287,12 +287,13 @@ end
 
 _default_rtol(t) = eps(real(float(scalartype(t)))) * min(dim(domain(t)), dim(codomain(t)))
 
-function LinearAlgebra.rank(t::AbstractTensorMap; atol::Real=0,
-                            rtol::Real=atol > 0 ? 0 : _default_rtol(t))
-    dim(t) == 0 && return 0
+function LinearAlgebra.rank(t::AbstractTensorMap;
+                            atol::Real=0, rtol::Real=atol > 0 ? 0 : _default_rtol(t))
+    init = dim(one(sectortype(t))) * 0
+    dim(t) == 0 && return init
     S = LinearAlgebra.svdvals(t)
     tol = max(atol, rtol * maximum(first, values(S)))
-    return sum(cs -> dim(cs[1]) * count(>(tol), cs[2]), S)
+    return sum(((c, b),) -> dim(c) * count(>(tol), b), S; init)
 end
 
 function LinearAlgebra.cond(t::AbstractTensorMap, p::Real=2)