diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index b1483f5..925a508 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -11,7 +11,7 @@ on: workflow_dispatch: jobs: - build: + documentation: runs-on: ${{ matrix.os }} strategy: matrix: diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index 69d8c33..562307b 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -12,4 +12,4 @@ permissions: jobs: runic: - uses: "QuantumKitHub/QuantumKitHubActions/.github/workflows/FormatCheck.yml@main"mKitHub/QuantumKitHubActions/.github/workflows/FormatCheck.yml@main" + uses: "QuantumKitHub/QuantumKitHubActions/.github/workflows/FormatCheck.yml@main" diff --git a/Project.toml b/Project.toml index 05fcaeb..9a7a204 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,13 @@ name = "BlockTensorKit" uuid = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd" authors = ["Lukas Devos and contributors"] -version = "0.2.0" +version = "0.3.0" [deps] BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" @@ -20,10 +21,11 @@ BlockArrays = "1" Combinatorics = "1" Compat = "4.13" LinearAlgebra = "1" +MatrixAlgebraKit = "0.5" Random = "1" SafeTestsets = "0.1" Strided = "2" -TensorKit = "0.14" +TensorKit = "0.15" TensorOperations = "5" Test = "1" TestExtras = "0.2, 0.3" diff --git a/src/BlockTensorKit.jl b/src/BlockTensorKit.jl index bdc4042..4f23d61 100644 --- a/src/BlockTensorKit.jl +++ b/src/BlockTensorKit.jl @@ -15,25 +15,12 @@ export dropzeros!, droptol! export undef_blocks using TensorKit -using TensorKit: - OneOrNoneIterator, - HomSpace, - MatrixAlgebra, - SectorDict, - AdjointTensorMap, - adjointtensorindices, - compose, - sectorscalartype +using TensorKit: OneOrNoneIterator, HomSpace, SectorDict, AdjointTensorMap, + adjointtensorindices, compose, sectorscalartype using VectorInterface using TensorOperations -using TensorOperations: - dimcheck_tensoradd, - dimcheck_tensorcontract, - dimcheck_tensortrace, - argcheck_tensoradd, - argcheck_tensorcontract, - argcheck_tensortrace, - AbstractBackend +using TensorOperations: dimcheck_tensoradd, dimcheck_tensorcontract, dimcheck_tensortrace, + argcheck_tensoradd, argcheck_tensorcontract, argcheck_tensortrace, AbstractBackend using LinearAlgebra using Strided using BlockArrays @@ -48,6 +35,7 @@ import VectorInterface as VI import TensorKit as TK import TensorOperations as TO import TupleTools as TT +import MatrixAlgebraKit as MAK # Spaces include("vectorspaces/sumspace.jl") @@ -64,7 +52,6 @@ include("tensors/vectorinterface.jl") include("tensors/tensoroperations.jl") include("linalg/linalg.jl") -include("linalg/matrixalgebra.jl") include("linalg/factorizations.jl") include("auxiliary/sparsetensorarray.jl") diff --git a/src/linalg/factorizations.jl b/src/linalg/factorizations.jl index 1962219..73b5462 100644 --- a/src/linalg/factorizations.jl +++ b/src/linalg/factorizations.jl @@ -1,267 +1,128 @@ -using TensorKit: QR, QRpos, QL, QLpos, SVD, SDD, Polar, LQ, LQpos, RQ, RQpos +using MatrixAlgebraKit +using MatrixAlgebraKit: AbstractAlgorithm, YALAPACK.BlasMat, Algorithm +import MatrixAlgebraKit as MAK -function TK.leftorth!( - t::AbstractBlockTensorMap; - alg::Union{QR, QRpos, QL, QLpos, SVD, SDD, Polar} = QRpos(), - atol::Real = zero(float(real(scalartype(t)))), - rtol::Real = if (alg ∉ (SVD(), SDD())) - zero(float(real(scalartype(t)))) - else - eps(real(float(one(scalartype(t))))) * iszero(atol) - end, - ) - InnerProductStyle(t) === EuclideanInnerProduct() || throw_invalid_innerproduct(:leftorth!) - if !iszero(rtol) - atol = max(atol, rtol * norm(t)) - end - I = sectortype(t) - dims = TK.SectorDict{I, Int}() - - # compute QR factorization for each block - if !isempty(TK.blocks(t)) - generator = Base.Iterators.map(TK.blocks(t)) do (c, b) - Qc, Rc = TK.MatrixAlgebra.leftorth!(b, alg, atol) - dims[c] = size(Qc, 2) - return c => (Qc, Rc) - end - QRdata = SectorDict(generator) - end +# Type piracy for defining the MAK rules on BlockArrays! +# ----------------------------------------------------- - # construct new space - S = spacetype(t) - V = S(dims) - if alg isa Polar - @assert V ≅ domain(t) - W = domain(t) - else - W = ProductSpace(V) - end +MAK.@algdef BlockAlgorithm - # construct output tensors - T = float(scalartype(t)) - Q = similar(t, T, codomain(t) ← W) - R = similar(t, T, W ← domain(t)) - if !isempty(blocksectors(domain(t))) - for (c, (Qc, Rc)) in QRdata - block(Q, c) .= Qc - block(R, c) .= Rc - end - end - return Q, R -end -function TK.leftorth!(t::SparseBlockTensorMap; kwargs...) - return leftorth!(BlockTensorMap(t); kwargs...) -end +const BlockBlasMat{T <: MAK.BlasFloat} = BlockMatrix{T} -function TK.leftnull!( - t::BlockTensorMap; - alg::Union{QR, QRpos, SVD, SDD} = QRpos(), - atol::Real = zero(float(real(scalartype(t)))), - rtol::Real = if (alg ∉ (SVD(), SDD())) - zero(float(real(scalartype(t)))) - else - eps(real(float(one(scalartype(t))))) * iszero(atol) - end, +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, ) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:leftnull!) - if !iszero(rtol) - atol = max(atol, rtol * norm(t)) - end - I = sectortype(t) - dims = SectorDict{I, Int}() + f! = Symbol(f, :!) - # compute QR factorization for each block - V = codomain(t) - if !isempty(blocksectors(V)) - generator = Base.Iterators.map(blocksectors(V)) do c - Nc = MatrixAlgebra.leftnull!(block(t, c), alg, atol) - dims[c] = size(Nc, 2) - return c => Nc - end - Ndata = SectorDict(generator) - end + @eval MAK.copy_input(::typeof($f), A::BlockBlasMat) = Array(A) + @eval MAK.$f!(t::BlockBlasMat, F, alg::BlockAlgorithm) = $f(t; alg.kwargs...) +end - # construct new space - S = spacetype(t) - W = S(dims) +for f in (:qr, :lq, :eig, :eigh, :gen_eig, :svd, :polar) + default_f_algorithm = Symbol(:default_, f, :_algorithm) + @eval MAK.$default_f_algorithm(::Type{<:BlockBlasMat{T}}; kwargs...) where {T} = + BlockAlgorithm(; kwargs...) +end - # construct output tensor - T = float(scalartype(t)) - N = similar(t, T, V ← W) - if !isempty(blocksectors(V)) - for (c, Nc) in Ndata - copy!(block(N, c), Nc) - end - end +# Make sure sparse blocktensormaps have dense outputs +function MAK.initialize_output(::typeof(qr_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_Q = fuse(codomain(t)) + Q = dense_similar(t, codomain(t) ← V_Q) + R = dense_similar(t, V_Q ← domain(t)) + return Q, R +end +function MAK.initialize_output(::typeof(qr_compact!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + Q = dense_similar(t, codomain(t) ← V_Q) + R = dense_similar(t, V_Q ← domain(t)) + return Q, R +end +function MAK.initialize_output(::typeof(qr_null!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(codomain(t)), V_Q) + N = dense_similar(t, codomain(t) ← V_N) return N end -TK.leftnull!(t::SparseBlockTensorMap; kwargs...) = leftnull!(BlockTensorMap(t); kwargs...) - -function TK.rightorth!( - t::AbstractBlockTensorMap; - alg::Union{LQ, LQpos, RQ, RQpos, SVD, SDD, Polar} = LQpos(), - atol::Real = zero(float(real(scalartype(t)))), - rtol::Real = if (alg ∉ (SVD(), SDD())) - zero(float(real(scalartype(t)))) - else - eps(real(float(one(scalartype(t))))) * iszero(atol) - end, - ) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:rightorth!) - if !iszero(rtol) - atol = max(atol, rtol * norm(t)) - end - I = sectortype(t) - dims = TK.SectorDict{I, Int}() - - # compute LQ factorization for each block - if !isempty(TK.blocks(t)) - generator = Base.Iterators.map(TK.blocks(t)) do (c, b) - Lc, Qc = TK.MatrixAlgebra.rightorth!(b, alg, atol) - dims[c] = size(Qc, 1) - return c => (Lc, Qc) - end - LQdata = SectorDict(generator) - end - # construct new space - S = spacetype(t) - V = S(dims) - if alg isa Polar - @assert V ≅ codomain(t) - W = codomain(t) - else - W = ProductSpace(V) - end - - # construct output tensors - T = float(scalartype(t)) - L = similar(t, T, codomain(t) ← W) - Q = similar(t, T, W ← domain(t)) - if !isempty(blocksectors(codomain(t))) - for (c, (Lc, Qc)) in LQdata - block(L, c) .= Lc - block(Q, c) .= Qc - end - end +function MAK.initialize_output(::typeof(lq_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_Q = fuse(domain(t)) + L = dense_similar(t, codomain(t) ← V_Q) + Q = dense_similar(t, V_Q ← domain(t)) return L, Q end -function TK.rightorth!(t::SparseBlockTensorMap; kwargs...) - return rightorth!(BlockTensorMap(t); kwargs...) +function MAK.initialize_output(::typeof(lq_compact!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + L = dense_similar(t, codomain(t) ← V_Q) + Q = dense_similar(t, V_Q ← domain(t)) + return L, Q end - -function TK.rightnull!( - t::BlockTensorMap; - alg::Union{LQ, LQpos, SVD, SDD} = LQpos(), - atol::Real = zero(float(real(scalartype(t)))), - rtol::Real = if (alg ∉ (SVD(), SDD())) - zero(float(real(scalartype(t)))) - else - eps(real(float(one(scalartype(t))))) * iszero(atol) - end, - ) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:rightnull!) - if !iszero(rtol) - atol = max(atol, rtol * norm(t)) - end - I = sectortype(t) - dims = SectorDict{I, Int}() - - # compute LQ factorization for each block - V = domain(t) - if !isempty(blocksectors(V)) - generator = Base.Iterators.map(blocksectors(V)) do c - Nc = MatrixAlgebra.rightnull!(block(t, c), alg, atol) - dims[c] = size(Nc, 1) - return c => Nc - end - Ndata = SectorDict(generator) - end - - # construct new space - S = spacetype(t) - W = S(dims) - - # construct output tensor - T = float(scalartype(t)) - N = similar(t, T, W ← V) - if !isempty(blocksectors(V)) - for (c, Nc) in Ndata - copy!(block(N, c), Nc) - end - end +function MAK.initialize_output(::typeof(lq_null!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(domain(t)), V_Q) + N = dense_similar(t, V_N ← domain(t)) return N end -TK.rightnull!(t::SparseBlockTensorMap; kwargs...) = rightnull!(BlockTensorMap(t); kwargs...) -function TK.tsvd!(t::AbstractBlockTensorMap; trunc = TK.NoTruncation(), p::Real = 2, alg = SDD()) - return TK._tsvd!(t, alg, trunc, p) +function MAK.initialize_output(::typeof(left_polar!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + W = dense_similar(t, space(t)) + P = dense_similar(t, domain(t) ← domain(t)) + return W, P end -function TK.tsvd!(t::SparseBlockTensorMap; kwargs...) - return tsvd!(BlockTensorMap(t); kwargs...) +function MAK.initialize_output(::typeof(right_polar!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + P = dense_similar(t, codomain(t) ← codomain(t)) + Wᴴ = dense_similar(t, space(t)) + return P, Wᴴ end -function TK._tsvd!( - t::BlockTensorMap, alg::Union{SVD, SDD}, trunc::TruncationScheme, p::Real = 2 - ) - # early return - if isempty(blocksectors(t)) - truncerr = zero(real(scalartype(t))) - return TK._empty_svdtensors(t)..., truncerr - end - - # compute SVD factorization for each block - S = spacetype(t) - SVDdata, dims = TK._compute_svddata!(t, alg) - Σdata = SectorDict(c => Σ for (c, (U, Σ, V)) in SVDdata) - truncdim = TK._compute_truncdim(Σdata, trunc, p) - truncerr = TK._compute_truncerr(Σdata, truncdim, p) - - # construct output tensors - U, Σ, V⁺ = TK._create_svdtensors(t, SVDdata, truncdim) - return U, Σ, V⁺, truncerr +function MAK.initialize_output(::typeof(left_null!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(codomain(t)), V_Q) + return dense_similar(t, codomain(t) ← V_N) +end +function MAK.initialize_output(::typeof(right_null!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(domain(t)), V_Q) + return dense_similar(t, V_N ← domain(t)) end -function TK._compute_svddata!(t::AbstractBlockTensorMap, alg::Union{SVD, SDD}) - InnerProductStyle(t) === EuclideanInnerProduct() || throw_invalid_innerproduct(:tsvd!) - I = sectortype(t) - dims = SectorDict{I, Int}() - generator = Base.Iterators.map(TK.blocks(t)) do (c, b) - U, Σ, V = TK.MatrixAlgebra.svd!(b, alg) - dims[c] = length(Σ) - return c => (U, Σ, V) - end - SVDdata = SectorDict(generator) - return SVDdata, dims +function MAK.initialize_output(::typeof(eigh_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_D = fuse(domain(t)) + T = real(scalartype(t)) + D = DiagonalTensorMap{T}(undef, V_D) + V = dense_similar(t, codomain(t) ← V_D) + return D, V end -function TK._create_svdtensors(t::AbstractBlockTensorMap, SVDdata, dims) - S = spacetype(t) - W = S(dims) - T = float(scalartype(t)) - U = similar(t, T, codomain(t) ← W) - Σ = similar(t, real(T), W ← W) - V⁺ = similar(t, T, W ← domain(t)) - for (c, (Uc, Σc, V⁺c)) in SVDdata - r = Base.OneTo(dims[c]) - block(U, c) .= view(Uc, :, r) - block(Σ, c) .= Diagonal(view(Σc, r)) - block(V⁺, c) .= view(V⁺c, r, :) - end - return U, Σ, V⁺ +function MAK.initialize_output(::typeof(eig_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_D = fuse(domain(t)) + Tc = complex(scalartype(t)) + D = DiagonalTensorMap{Tc}(undef, V_D) + V = dense_similar(t, Tc, codomain(t) ← V_D) + return D, V end -function TK._empty_svdtensors(t::AbstractBlockTensorMap) - T = scalartype(t) - S = spacetype(t) - I = sectortype(t) - dims = SectorDict{I, Int}() - W = S(dims) +function MAK.initialize_output(::typeof(svd_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_cod = fuse(codomain(t)) + V_dom = fuse(domain(t)) + U = dense_similar(t, codomain(t) ← V_cod) + S = similar(t, real(scalartype(t)), V_cod ← V_dom) + Vᴴ = dense_similar(t, V_dom ← domain(t)) + return U, S, Vᴴ +end +function MAK.initialize_output(::typeof(svd_compact!), t::AbstractBlockTensorMap, ::AbstractAlgorithm) + V_cod = V_dom = infimum(fuse(codomain(t)), fuse(domain(t))) + U = dense_similar(t, codomain(t) ← V_cod) + S = DiagonalTensorMap{real(scalartype(t))}(undef, V_cod) + Vᴴ = dense_similar(t, V_dom ← domain(t)) + return U, S, Vᴴ +end - U = similar(t, codomain(t) ← W) - Σ = similar(t, real(T), W ← W) - V⁺ = similar(t, W ← domain(t)) - return U, Σ, V⁺ +# Disambiguate Diagonal implementations +# ------------------------------------- +# these shouldn't ever happen as blocktensors aren't diagonal +for f! in (:eig_full!, :eigh_full!, :lq_compact!, :lq_full!, :qr_compact!, :qr_full!, :svd_full!) + @eval function MAK.initialize_output(::typeof($f!), t::AbstractBlockTensorMap, alg::DiagonalAlgorithm) + error("Blocktensors are incompatible with diagonal algorithm") + end end diff --git a/src/linalg/matrixalgebra.jl b/src/linalg/matrixalgebra.jl deleted file mode 100644 index 10cfe08..0000000 --- a/src/linalg/matrixalgebra.jl +++ /dev/null @@ -1,26 +0,0 @@ -# TODO: is it possible to avoid making this contiguous? -function MatrixAlgebra.leftorth!(A::BlockMatrix, alg, atol::Real) - return MatrixAlgebra.leftorth!(Array(A), alg, atol) -end -function MatrixAlgebra.leftnull!(A::BlockMatrix, alg, atol::Real) - return MatrixAlgebra.leftnull!(Array(A), alg, atol) -end - -function MatrixAlgebra.rightorth!(A::BlockMatrix, alg, atol::Real) - return MatrixAlgebra.rightorth!(Array(A), alg, atol) -end -function MatrixAlgebra.rightnull!(A::BlockMatrix, alg, atol::Real) - return MatrixAlgebra.rightnull!(Array(A), alg, atol) -end - -function MatrixAlgebra.one!(A::BlockMatrix) - @inbounds for i in axes(A, 1), j in axes(A, 2) - A[i, j] = i == j ? one(eltype(A)) : zero(eltype(A)) - end - return A -end - -MatrixAlgebra.svd!(A::BlockMatrix, alg) = MatrixAlgebra.svd!(Array(A), alg) - -# piracy -LinearAlgebra.isposdef!(A::BlockMatrix) = LinearAlgebra.isposdef!(Array(A)) diff --git a/src/tensors/abstractblocktensor/abstractarray.jl b/src/tensors/abstractblocktensor/abstractarray.jl index 8b997b3..9211fa9 100644 --- a/src/tensors/abstractblocktensor/abstractarray.jl +++ b/src/tensors/abstractblocktensor/abstractarray.jl @@ -241,8 +241,18 @@ Base.similar(t::AbstractBlockTensorMap, P::TensorMapSumSpace) = similar(t, eltyp function Base.similar( t::AbstractTensorMap, ::Type{TorA}, P::TensorMapSumSpace{S} ) where {S, TorA} + return issparse(t) ? sparse_similar(t, TorA, P) : dense_similar(t, TorA, P) +end + +dense_similar(t::AbstractTensorMap, P::TensorMapSumSpace) = dense_similar(t, TK.similarstoragetype(t), P) +function dense_similar(t::AbstractTensorMap, ::Type{TorA}, P::TensorMapSumSpace) where {TorA} + TT = similar_tensormaptype(t, TorA, P) + return BlockTensorMap{TT}(undef, P) +end +sparse_similar(t::AbstractTensorMap, P::TensorMapSumSpace) = sparse_similar(t, TK.similarstoragetype(t), P) +function sparse_similar(t::AbstractTensorMap, ::Type{TorA}, P::TensorMapSumSpace) where {TorA} TT = similar_tensormaptype(t, TorA, P) - return issparse(t) ? SparseBlockTensorMap{TT}(undef, P) : BlockTensorMap{TT}(undef, P) + return SparseBlockTensorMap{TT}(undef, P) end function similar_tensormaptype( diff --git a/src/tensors/abstractblocktensor/abstracttensormap.jl b/src/tensors/abstractblocktensor/abstracttensormap.jl index bac400b..6253c14 100644 --- a/src/tensors/abstractblocktensor/abstracttensormap.jl +++ b/src/tensors/abstractblocktensor/abstracttensormap.jl @@ -50,22 +50,59 @@ end return setindex!(t, v′, f₁, f₂) end -function TensorKit.block(t::AbstractBlockTensorMap, c::Sector) +function TensorKit.block(t::AbstractBlockTensorMap, c::Sector)::TK.blocktype(t) sectortype(t) == typeof(c) || throw(SectorMismatch()) rows = prod(TT.getindices(size(t), codomainind(t))) cols = prod(TT.getindices(size(t), domainind(t))) - @assert rows != 0 && cols != 0 "to be added" + + if rows == 0 || cols == 0 + allblocks = Matrix{TK.blocktype(eltype(t))}(undef, rows, cols) + + rowaxes = Int[] + if rows != 0 + W′ = codomain(t) ← zero(spacetype(t)) + for V in eachspace(W′) + push!(rowaxes, blockdim(codomain(V), c)) + end + end + + colaxes = Int[] + if cols != 0 + W′ = zero(spacetype(t)) ← domain(t) + for V in eachspace(W′) + push!(colaxes, blockdim(domain(V), c)) + end + end + + return mortar(allblocks, rowaxes, colaxes) + end allblocks = map(Base.Fix2(block, c), parent(t)) return mortar(reshape(allblocks, rows, cols)) end # TODO: this might get fixed once new tensormap is implemented -TensorKit.blocks(t::AbstractBlockTensorMap) = ((c => block(t, c)) for c in blocksectors(t)) TensorKit.blocksectors(t::AbstractBlockTensorMap) = blocksectors(space(t)) TensorKit.hasblock(t::AbstractBlockTensorMap, c::Sector) = c in blocksectors(t) +TensorKit.blocks(t::AbstractBlockTensorMap) = TK.BlockIterator(t, blocksectors(t)) +Base.@assume_effects :foldable function TensorKit.blocktype(::Type{TT}) where {TT <: AbstractBlockTensorMap} + T = scalartype(TT) + B = TK.blocktype(eltype(TT)) + (B <: AbstractMatrix{T}) || (B = AbstractMatrix{T}) # safeguard against type-instability + BS = NTuple{2, BlockedOneTo{Int, Vector{Int}}} + return BlockMatrix{T, Matrix{B}, BS} +end + +function Base.iterate(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, state...) + next = iterate(iter.structure, state...) + isnothing(next) && return next + c, newstate = next + return c => block(iter.t, c), newstate +end +Base.getindex(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, c::Sector) = block(iter.t, c) + function TensorKit.storagetype(::Type{TT}) where {TT <: AbstractBlockTensorMap} return if isconcretetype(eltype(TT)) storagetype(eltype(TT)) diff --git a/src/vectorspaces/sumspace.jl b/src/vectorspaces/sumspace.jl index e805ede..5c6d37c 100644 --- a/src/vectorspaces/sumspace.jl +++ b/src/vectorspaces/sumspace.jl @@ -114,6 +114,8 @@ TensorKit.compose(V, W) = TensorKit.compose(promote(V, W)...) # bit of a hack to make spacechecks happy? Base.:(==)(V::SumSpace{S}, W::S) where {S} = ==(promote(V, W)...) Base.:(==)(V::S, W::SumSpace{S}) where {S} = ==(promote(V, W)...) +Base.:(==)(V::ProductSumSpace{S}, W::ProductSpace{S}) where {S <: ElementarySpace} = ==(promote(V, W)...) +Base.:(==)(V::ProductSpace{S}, W::ProductSumSpace{S}) where {S <: ElementarySpace} = ==(promote(V, W)...) function Base.:(==)(V::TensorMapSumSpace{S}, W::TensorMapSpace{S}) where {S <: IndexSpace} return ==(promote(V, W)...) end @@ -124,6 +126,11 @@ end function Base.:(==)(V::TensorMapSumSpace{S}, W::TensorMapSumSpace{S}) where {S <: IndexSpace} return @invoke ==(V::HomSpace, W::HomSpace) end + + +TensorKit.infimum(V::S, W::S) where {S <: SumSpace} = infimum(TensorKit.oplus(V), TensorKit.oplus(W)) +TensorKit.supremum(V::S, W::S) where {S <: SumSpace} = supremum(TensorKit.oplus(V), TensorKit.oplus(W)) +TensorKit.ominus(V::S, W::S) where {S <: SumSpace} = ominus(TensorKit.oplus(V), TensorKit.oplus(W)) # this conflicts with the definition in TensorKit, so users always need to specify # ⊕(Vs::IndexSpace...) = SumSpace(Vs...) @@ -148,21 +155,16 @@ function ⊕(V₁::SumSpace{S}, V₂::SumSpace{S}) where {S} end #! format: off -function TensorKit.:⊕(S::SumSpace) - if length(S) == 1 - return only(S.spaces) - else - return TensorKit.oplus(S.spaces...) - end -end -TensorKit.:⊕(V1::SumSpace, V2::SumSpace...) = TensorKit.oplus(⊕(V1, V2...)) +TensorKit.:⊕(V::SumSpace{S}) where {S} = reduce(TK.oplus, V.spaces; init = isdual(V) ? zero(S)' : zero(S)) +TensorKit.:⊕(V1::SumSpace{S}, V2::SumSpace{S}...) where {S} = TensorKit.oplus(⊕(V1, V2...)) #! format: on -function TensorKit.fuse(V1::S, V2::S) where {S <: SumSpace} - return SumSpace(vec([fuse(v1, v2) for (v1, v2) in Base.product(V1.spaces, V2.spaces)])) -end +TensorKit.fuse(V1::SumSpace) = fuse(TK.oplus(V1)) +TK.fuse(V1::S, V2::S) where {S <: SumSpace} = fuse(TK.oplus(V1), TK.oplus(V2)) Base.oneunit(S::Type{<:SumSpace}) = SumSpace(oneunit(eltype(S))) +Base.zero(V::SumSpace{S}) where {S} = SumSpace{S}(; dual = isdual(V)) +Base.zero(::Type{SumSpace{S}}) where {S} = SumSpace{S}() # Promotion and conversion # ------------------------ diff --git a/test/linalg/factorizations.jl b/test/linalg/factorizations.jl index 3f7ee05..2321bc1 100644 --- a/test/linalg/factorizations.jl +++ b/test/linalg/factorizations.jl @@ -4,6 +4,7 @@ using TensorKit using BlockTensorKit using Random using Combinatorics +using LinearAlgebra Vtr = ( SumSpace(ℂ^3), @@ -20,122 +21,455 @@ for V in (Vtr,) end spacelist = (Vtr,) -scalartypes = (Float64, ComplexF64) -V = first(spacelist) -# @testset "Tensors with symmetry: $(TensorKit.type_repr(sectortype(first(V))))" verbose = true failfast=true for V in -# spacelist -I = sectortype(first(V)) -V1, V2, V3, V4, V5 = V -W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 -W_empty = V1 ⊗ V2 ← eltype(V1)() - -const leftorth_algs = ( - TensorKit.QR(), - TensorKit.QRpos(), - TensorKit.QL(), - TensorKit.QLpos(), - TensorKit.Polar(), - TensorKit.SVD(), - TensorKit.SDD(), -) +eltypes = (Float64, ComplexF64) + +for V in spacelist + I = sectortype(first(V)) + Istr = TensorKit.type_repr(I) + println("---------------------------------------") + println("Factorizations with symmetry: $Istr") + println("---------------------------------------") + @timedtestset "Factorizations with symmetry: $Istr" verbose = true begin + V1, V2, V3, V4, V5 = V + W = V1 ⊗ V2 + + @testset "QR decomposition" begin + for T in eltypes, + t in ( + rand(T, W, W), rand(T, W, W)', rand(T, W, V1), rand(T, V1, W)', + ) + + Q, R = @constinferred qr_full(t) + @test Q * R ≈ t + @test isunitary(Q) + + Q, R = @constinferred qr_compact(t) + @test Q * R ≈ t + @test isisometry(Q) + + Q, R = @constinferred left_orth(t; kind = :qr) + @test Q * R ≈ t + @test isisometry(Q) + + N = @constinferred qr_null(t) + @test isisometry(N) + @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) + + N = @constinferred left_null(t; kind = :qr) + @test isisometry(N) + @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) + end -@testset "leftorth with $alg" for alg in leftorth_algs - for T in (Float32, ComplexF64), isadjoint in (false, true) - t = isadjoint ? rand(T, W)' : rand(T, W) - Q, R = @inferred leftorth(t, ((3, 4, 2), (1, 5)); alg) - QdQ = Q' * Q - @test QdQ ≈ one(QdQ) - @test Q * R ≈ permute(t, ((3, 4, 2), (1, 5))) - if alg isa Polar - @test isposdef(R) - @test domain(R) == codomain(R) == space(t, 1)' ⊗ space(t, 5)' + # empty tensor + for T in eltypes + t = rand(T, V1 ⊗ V2, zero(V1)) + + Q, R = @constinferred qr_full(t) + @test Q * R ≈ t + @test isunitary(Q) + @test dim(R) == dim(t) == 0 + + Q, R = @constinferred qr_compact(t) + @test Q * R ≈ t + @test isisometry(Q) + @test dim(Q) == dim(R) == dim(t) + + Q, R = @constinferred left_orth(t; kind = :qr) + @test Q * R ≈ t + @test isisometry(Q) + @test dim(Q) == dim(R) == dim(t) + + N = @constinferred qr_null(t) + @test isunitary(N) + @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) + end end - t_empty = isadjoint ? rand(T, W_empty')' : rand(T, W_empty) - Q, R = @constinferred leftorth(t_empty; alg) - @test Q == t_empty - @test dim(Q) == dim(R) == 0 - end -end + @testset "LQ decomposition" begin + for T in eltypes, + t in ( + rand(T, W, W), rand(T, W, W)', rand(T, W, V1), rand(T, V1, W)', + ) -const leftnull_algs = (TensorKit.QR(), TensorKit.SVD(), TensorKit.SDD()) -@testset "leftnull with $alg" for alg in leftnull_algs - for T in (Float32, ComplexF64), isadjoint in (false, true) - t = isadjoint ? rand(T, W)' : rand(T, W) - N = @constinferred leftnull(t, ((3, 4, 2), (1, 5)); alg = alg) - NdN = N' * N - @test NdN ≈ one(NdN) - @test norm(N' * permute(t, ((3, 4, 2), (1, 5)))) < 100 * eps(norm(t)) - - t_empty = isadjoint ? rand(T, W_empty')' : rand(T, W_empty) - N = @constinferred leftnull(t_empty; alg = alg) - @test N' * N ≈ id(domain(N)) - @test N * N' ≈ id(codomain(N)) - end -end + L, Q = @constinferred lq_full(t) + @test L * Q ≈ t + @test isunitary(Q) -const rightorth_algs = ( - TensorKit.RQ(), - TensorKit.RQpos(), - TensorKit.LQ(), - TensorKit.LQpos(), - TensorKit.Polar(), - TensorKit.SVD(), - TensorKit.SDD(), -) -@testset "rightorth with $alg" for alg in rightorth_algs - for T in (Float32, ComplexF64), isadjoint in (false, true) - t = isadjoint ? rand(T, W)' : rand(T, W) - L, Q = @constinferred rightorth(t, ((3, 4), (2, 1, 5)); alg = alg) - QQd = Q * Q' - @test QQd ≈ one(QQd) - @test L * Q ≈ permute(t, ((3, 4), (2, 1, 5))) - if alg isa Polar - @test isposdef(L) - @test domain(L) == codomain(L) == space(t, 3) ⊗ space(t, 4) + L, Q = @constinferred lq_compact(t) + @test L * Q ≈ t + @test isisometry(Q; side = :right) + + L, Q = @constinferred right_orth(t; kind = :lq) + @test L * Q ≈ t + @test isisometry(Q; side = :right) + + Nᴴ = @constinferred lq_null(t) + @test isisometry(Nᴴ; side = :right) + @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) + end + + for T in eltypes + # empty tensor + t = rand(T, zero(V1), V1 ⊗ V2) + + L, Q = @constinferred lq_full(t) + @test L * Q ≈ t + @test isunitary(Q) + @test dim(L) == dim(t) == 0 + + L, Q = @constinferred lq_compact(t) + @test L * Q ≈ t + @test isisometry(Q; side = :right) + @test dim(Q) == dim(L) == dim(t) + + L, Q = @constinferred right_orth(t; kind = :lq) + @test L * Q ≈ t + @test isisometry(Q; side = :right) + @test dim(Q) == dim(L) == dim(t) + + Nᴴ = @constinferred lq_null(t) + @test isunitary(Nᴴ) + @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) + end end - t_empty = isadjoint ? rand(T, W_empty)' : rand(T, W_empty') - L, Q = @constinferred rightorth(t_empty; alg = alg) - @test Q == t_empty - @test dim(Q) == dim(L) == 0 - end -end + @testset "Polar decomposition" begin + for T in eltypes, + t in ( + rand(T, W, W), rand(T, W, W)', rand(T, W, V1), rand(T, V1, W)', + ) -const rightnull_algs = (TensorKit.LQ(), TensorKit.SVD(), TensorKit.SDD()) -@testset "rightnull with $alg" for alg in rightnull_algs - for T in (Float32, ComplexF64), isadjoint in (false, true) - t = isadjoint ? rand(T, W)' : rand(T, W) - M = @constinferred rightnull(t, ((3, 4), (2, 1, 5)); alg = alg) - MMd = M * M' - @test MMd ≈ one(MMd) - @test norm(permute(t, ((3, 4), (2, 1, 5))) * M') < 100 * eps(norm(t)) - - t_empty = isadjoint ? rand(T, W_empty)' : rand(T, W_empty') - M = @constinferred rightnull(t_empty; alg) - @test M * M' ≈ id(codomain(M)) - @test M' * M ≈ id(domain(M)) - end -end + @assert domain(t) ≾ codomain(t) + w, p = @constinferred left_polar(t) + @test w * p ≈ t + @test isisometry(w) + # @test isposdef(p) + + w, p = @constinferred left_orth(t; kind = :polar) + @test w * p ≈ t + @test isisometry(w) + end + + for T in eltypes, + t in (rand(T, W, W), rand(T, W, W)', rand(T, V1, W), rand(T, W, V1)') + + @assert codomain(t) ≾ domain(t) + p, wᴴ = @constinferred right_polar(t) + @test p * wᴴ ≈ t + @test isisometry(wᴴ; side = :right) + # @test isposdef(p) + + p, wᴴ = @constinferred right_orth(t; kind = :polar) + @test p * wᴴ ≈ t + @test isisometry(wᴴ; side = :right) + end + end + + @testset "SVD" begin + for T in eltypes, + t in ( + rand(T, W, W), rand(T, W, W)', + rand(T, W, V1), rand(T, V1, W), + rand(T, W, V1)', rand(T, V1, W)', + ) + + u, s, vᴴ = @constinferred svd_full(t) + @test u * s * vᴴ ≈ t + @test isunitary(u) + @test isunitary(vᴴ) + + u, s, vᴴ = @constinferred svd_compact(t) + @test u * s * vᴴ ≈ t + @test isisometry(u) + # @test isposdef(s) + @test isisometry(vᴴ; side = :right) + + s′ = LinearAlgebra.diag(s) + for (c, b) in LinearAlgebra.svdvals(t) + @test b ≈ s′[c] + end -const svd_algs = (TensorKit.SVD(), TensorKit.SDD()) -@testset "tsvd with $alg" for alg in svd_algs - for T in (Float32, ComplexF64), isadjoint in (false, true) - t = isadjoint ? rand(T, W)' : rand(T, W) - U, S, V = @constinferred tsvd(t, ((3, 4, 2), (1, 5)); alg) - UdU = U' * U - @test UdU ≈ one(UdU) - VVd = V * V' - @test VVd ≈ one(VVd) - @test U * S * V ≈ permute(t, ((3, 4, 2), (1, 5))) - - t_empty = isadjoint ? rand(T, W_empty')' : rand(T, W_empty) - U, S, V = @inferred tsvd(t_empty; alg = alg) - @test U == t_empty - @test dim(U) == dim(S) == dim(V) + v, c = @constinferred left_orth(t; kind = :svd) + @test v * c ≈ t + @test isisometry(v) + + N = @constinferred left_null(t; kind = :svd) + @test isisometry(N) + @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) + + Nᴴ = @constinferred right_null(t; kind = :svd) + @test isisometry(Nᴴ; side = :right) + @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) + end + + # empty tensor + for T in eltypes, t in (rand(T, W, zero(V1)), rand(T, zero(V1), W)) + U, S, Vᴴ = @constinferred svd_full(t) + @test U * S * Vᴴ ≈ t + @test isunitary(U) + @test isunitary(Vᴴ) + + U, S, Vᴴ = @constinferred svd_compact(t) + @test U * S * Vᴴ ≈ t + @test dim(U) == dim(S) == dim(Vᴴ) == dim(t) == 0 + end + end + + @testset "truncated SVD" begin + for T in eltypes, + t in ( + randn(T, W, W), randn(T, W, W)', + randn(T, W, V1), randn(T, V1, W), + randn(T, W, V1)', randn(T, V1, W)', + ) + + @constinferred normalize!(t) + + U, S, Vᴴ = @constinferred svd_trunc(t; trunc = notrunc()) + @test U * S * Vᴴ ≈ t + @test isisometry(U) + @test isisometry(Vᴴ; side = :right) + + trunc = truncrank(dim(domain(S)) ÷ 2) + U1, S1, Vᴴ1 = @constinferred svd_trunc(t; trunc) + @test t * Vᴴ1' ≈ U1 * S1 + @test isisometry(U1) + @test isisometry(Vᴴ1; side = :right) + @test dim(domain(S1)) <= trunc.howmany + + λ = minimum(minimum, values(LinearAlgebra.diag(S1))) + trunc = trunctol(; atol = λ - 10eps(λ)) + U2, S2, Vᴴ2 = @constinferred svd_trunc(t; trunc) + @test t * Vᴴ2' ≈ U2 * S2 + @test isisometry(U2) + @test isisometry(Vᴴ2; side = :right) + @test minimum(minimum, values(LinearAlgebra.diag(S1))) >= λ + @test U2 ≈ U1 + @test S2 ≈ S1 + @test Vᴴ2 ≈ Vᴴ1 + + trunc = truncspace(space(S2, 1)) + U3, S3, Vᴴ3 = @constinferred svd_trunc(t; trunc) + @test t * Vᴴ3' ≈ U3 * S3 + @test isisometry(U3) + @test isisometry(Vᴴ3; side = :right) + @test space(S3, 1) ≾ space(S2, 1) + + trunc = truncerror(; atol = 0.5) + U4, S4, Vᴴ4 = @constinferred svd_trunc(t; trunc) + @test t * Vᴴ4' ≈ U4 * S4 + @test isisometry(U4) + @test isisometry(Vᴴ4; side = :right) + @test norm(t - U4 * S4 * Vᴴ4) <= 0.5 + end + end + + @testset "Eigenvalue decomposition" begin + for T in eltypes, + t in ( + rand(T, V1, V1), rand(T, W, W), rand(T, W, W)', + ) + + d, v = @constinferred eig_full(t) + @test t * v ≈ v * d + + d′ = LinearAlgebra.diag(d) + for (c, b) in LinearAlgebra.eigvals(t) + @test sort(b; by = abs) ≈ sort(d′[c]; by = abs) + end + + # vdv = v' * v + # vdv = (vdv + vdv') / 2 + # @test @constinferred isposdef(vdv) + # t isa DiagonalTensorMap || @test !isposdef(t) # unlikely for non-hermitian map + + d, v = @constinferred eig_trunc(t; trunc = truncrank(dim(domain(t)) ÷ 2)) + @test t * v ≈ v * d + @test dim(domain(d)) ≤ dim(domain(t)) ÷ 2 + + t2 = (t + t') + D, V = eigen(t2) + @test isisometry(V) + D̃, Ṽ = @constinferred eigh_full(t2) + @test D ≈ D̃ + @test V ≈ Ṽ + # λ = minimum( + # minimum(real(LinearAlgebra.diag(b))) + # for (c, b) in blocks(D) + # ) + # @test cond(Ṽ) ≈ one(real(T)) + # @test isposdef(t2) == isposdef(λ) + # @test isposdef(t2 - λ * one(t2) + 0.1 * one(t2)) + # @test !isposdef(t2 - λ * one(t2) - 0.1 * one(t2)) + + d, v = @constinferred eigh_full(t2) + @test t2 * v ≈ v * d + @test isunitary(v) + + # λ = minimum(minimum(real(LinearAlgebra.diag(b))) for (c, b) in blocks(d)) + # @test cond(v) ≈ one(real(T)) + # @test isposdef(t2) == isposdef(λ) + # @test isposdef(t2 - λ * one(t2) + 0.1 * one(t2)) + # @test !isposdef(t2 - λ * one(t2) - 0.1 * one(t2)) + + d, v = @constinferred eigh_trunc(t2; trunc = truncrank(dim(domain(t2)) ÷ 2)) + @test t2 * v ≈ v * d + @test dim(domain(d)) ≤ dim(domain(t2)) ÷ 2 + end + end + + # TODO: currently not supported + # @testset "Condition number and rank" begin + # for T in eltypes, + # t in ( + # rand(T, W, W), rand(T, W, W)', + # rand(T, W, V1), rand(T, V1, W), + # rand(T, W, V1)', rand(T, V1, W)', + # ) + # + # d1, d2 = dim(codomain(t)), dim(domain(t)) + # @test rank(t) == min(d1, d2) + # M = left_null(t) + # @test @constinferred(rank(M)) + rank(t) == d1 + # Mᴴ = right_null(t) + # @test rank(Mᴴ) + rank(t) == d2 + # end + # for T in eltypes + # u = unitary(T, V1 ⊗ V2, V1 ⊗ V2) + # @test @constinferred(cond(u)) ≈ one(real(T)) + # @test @constinferred(rank(u)) == dim(V1 ⊗ V2) + # + # t = rand(T, zero(V1), W) + # @test rank(t) == 0 + # t2 = rand(T, zero(V1) * zero(V2), zero(V1) * zero(V2)) + # @test rank(t2) == 0 + # @test cond(t2) == 0.0 + # end + # for T in eltypes, t in (rand(T, W, W), rand(T, W, W)') + # t += t' + # vals = @constinferred LinearAlgebra.eigvals(t) + # λmax = maximum(s -> maximum(abs, s), values(vals)) + # λmin = minimum(s -> minimum(abs, s), values(vals)) + # @test cond(t) ≈ λmax / λmin + # end + # end end end +# @testset "Tensors with symmetry: $(TensorKit.type_repr(sectortype(first(V))))" verbose = true failfast=true for V in +# spacelist +# I = sectortype(first(V)) +# V1, V2, V3, V4, V5 = V +# W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 +# W_empty = V1 ⊗ V2 ← eltype(V1)() +# +# const leftorth_algs = ( +# TensorKit.QR(), +# TensorKit.QRpos(), +# TensorKit.QL(), +# TensorKit.QLpos(), +# TensorKit.Polar(), +# TensorKit.SVD(), +# TensorKit.SDD(), +# ) +# +# @testset "leftorth with $alg" for alg in leftorth_algs +# for T in (Float32, ComplexF64), isadjoint in (false, true) +# t = isadjoint ? rand(T, W)' : rand(T, W) +# Q, R = @inferred leftorth(t, ((3, 4, 2), (1, 5)); alg) +# QdQ = Q' * Q +# @test QdQ ≈ one(QdQ) +# @test Q * R ≈ permute(t, ((3, 4, 2), (1, 5))) +# if alg isa Polar +# @test isposdef(R) +# @test domain(R) == codomain(R) == space(t, 1)' ⊗ space(t, 5)' +# end +# +# t_empty = isadjoint ? rand(T, W_empty')' : rand(T, W_empty) +# Q, R = @constinferred leftorth(t_empty; alg) +# @test Q == t_empty +# @test dim(Q) == dim(R) == 0 +# end +# end +# +# const leftnull_algs = (TensorKit.QR(), TensorKit.SVD(), TensorKit.SDD()) +# @testset "leftnull with $alg" for alg in leftnull_algs +# for T in (Float32, ComplexF64), isadjoint in (false, true) +# t = isadjoint ? rand(T, W)' : rand(T, W) +# N = @constinferred leftnull(t, ((3, 4, 2), (1, 5)); alg = alg) +# NdN = N' * N +# @test NdN ≈ one(NdN) +# @test norm(N' * permute(t, ((3, 4, 2), (1, 5)))) < 100 * eps(norm(t)) +# +# t_empty = isadjoint ? rand(T, W_empty')' : rand(T, W_empty) +# N = @constinferred leftnull(t_empty; alg = alg) +# @test N' * N ≈ id(domain(N)) +# @test N * N' ≈ id(codomain(N)) +# end +# end +# +# const rightorth_algs = ( +# TensorKit.RQ(), +# TensorKit.RQpos(), +# TensorKit.LQ(), +# TensorKit.LQpos(), +# TensorKit.Polar(), +# TensorKit.SVD(), +# TensorKit.SDD(), +# ) +# @testset "rightorth with $alg" for alg in rightorth_algs +# for T in (Float32, ComplexF64), isadjoint in (false, true) +# t = isadjoint ? rand(T, W)' : rand(T, W) +# L, Q = @constinferred rightorth(t, ((3, 4), (2, 1, 5)); alg = alg) +# QQd = Q * Q' +# @test QQd ≈ one(QQd) +# @test L * Q ≈ permute(t, ((3, 4), (2, 1, 5))) +# if alg isa Polar +# @test isposdef(L) +# @test domain(L) == codomain(L) == space(t, 3) ⊗ space(t, 4) +# end +# +# t_empty = isadjoint ? rand(T, W_empty)' : rand(T, W_empty') +# L, Q = @constinferred rightorth(t_empty; alg = alg) +# @test Q == t_empty +# @test dim(Q) == dim(L) == 0 +# end +# end +# +# const rightnull_algs = (TensorKit.LQ(), TensorKit.SVD(), TensorKit.SDD()) +# @testset "rightnull with $alg" for alg in rightnull_algs +# for T in (Float32, ComplexF64), isadjoint in (false, true) +# t = isadjoint ? rand(T, W)' : rand(T, W) +# M = @constinferred rightnull(t, ((3, 4), (2, 1, 5)); alg = alg) +# MMd = M * M' +# @test MMd ≈ one(MMd) +# @test norm(permute(t, ((3, 4), (2, 1, 5))) * M') < 100 * eps(norm(t)) +# +# t_empty = isadjoint ? rand(T, W_empty)' : rand(T, W_empty') +# M = @constinferred rightnull(t_empty; alg) +# @test M * M' ≈ id(codomain(M)) +# @test M' * M ≈ id(domain(M)) +# end +# end +# +# const svd_algs = (TensorKit.SVD(), TensorKit.SDD()) +# @testset "tsvd with $alg" for alg in svd_algs +# for T in (Float32, ComplexF64), isadjoint in (false, true) +# t = isadjoint ? rand(T, W)' : rand(T, W) +# U, S, V = @constinferred tsvd(t, ((3, 4, 2), (1, 5)); alg) +# UdU = U' * U +# @test UdU ≈ one(UdU) +# VVd = V * V' +# @test VVd ≈ one(VVd) +# @test U * S * V ≈ permute(t, ((3, 4, 2), (1, 5))) +# +# t_empty = isadjoint ? rand(T, W_empty')' : rand(T, W_empty) +# U, S, V = @inferred tsvd(t_empty; alg = alg) +# @test U == t_empty +# @test dim(U) == dim(S) == dim(V) +# end +# end + # t = Tensor(rand, T, V1 ⊗ V1' ⊗ V2 ⊗ V2') # @testset "eig and isposdef" begin # D, V = eigen(t, (1, 3), (2, 4))