From f50f00eecb62aede264f714d3ddca06c7bf0fea1 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 10 Oct 2025 17:12:49 -0400 Subject: [PATCH 1/2] Vendor TypeParameterAccessors --- NDTensors/Project.toml | 4 +- NDTensors/src/NDTensors.jl | 5 + .../generic_array_constructors.jl | 48 +- NDTensors/src/abstractarray/iscu.jl | 4 +- NDTensors/src/abstractarray/set_types.jl | 4 +- NDTensors/src/abstractarray/similar.jl | 28 +- NDTensors/src/adapt.jl | 8 +- NDTensors/src/blocksparse/blockdims.jl | 94 +-- .../src/blocksparse/blocksparsetensor.jl | 6 +- NDTensors/src/blocksparse/diagblocksparse.jl | 2 +- NDTensors/src/blocksparse/linearalgebra.jl | 712 +++++++++--------- NDTensors/src/blocksparse/similar.jl | 22 +- NDTensors/src/default_kwargs.jl | 2 +- NDTensors/src/dense/dense.jl | 120 +-- .../src/dense/generic_array_constructors.jl | 48 +- NDTensors/src/dense/set_types.jl | 10 +- NDTensors/src/diag/diagtensor.jl | 222 +++--- NDTensors/src/diag/set_types.jl | 14 +- NDTensors/src/diag/similar.jl | 14 +- NDTensors/src/dims.jl | 6 +- NDTensors/src/empty/empty.jl | 4 +- .../AMDGPUExtensions/src/AMDGPUExtensions.jl | 8 + NDTensors/src/lib/AMDGPUExtensions/src/roc.jl | 4 +- .../lib/CUDAExtensions/src/CUDAExtensions.jl | 8 + NDTensors/src/lib/CUDAExtensions/src/cuda.jl | 4 +- NDTensors/src/lib/Expose/src/Expose.jl | 7 + NDTensors/src/lib/Expose/src/exposed.jl | 12 +- .../src/GPUArraysCoreExtensions.jl | 8 + .../src/gpuarrayscore.jl | 8 +- .../MetalExtensions/src/MetalExtensions.jl | 8 + .../src/lib/MetalExtensions/src/metal.jl | 4 +- .../src/RankFactorization.jl | 9 + .../RankFactorization/src/default_kwargs.jl | 2 +- .../src/truncate_spectrum.jl | 2 +- NDTensors/src/linearalgebra/linearalgebra.jl | 662 ++++++++-------- NDTensors/src/linearalgebra/svd.jl | 2 +- NDTensors/src/tensor/set_types.jl | 28 +- NDTensors/src/tensor/similar.jl | 50 +- .../src/tensorstorage/default_storage.jl | 12 +- NDTensors/src/tensorstorage/set_types.jl | 10 +- NDTensors/src/tensorstorage/similar.jl | 68 +- NDTensors/src/truncate.jl | 4 +- .../src/TypeParameterAccessors.jl | 23 + .../src/base/abstractarray.jl | 94 +++ .../src/base/linearalgebra.jl | 25 + .../src/base/similartype.jl | 93 +++ .../TypeParameterAccessors/src/ndims.jl | 6 + .../src/type_parameters.jl | 331 ++++++++ .../TypeParameterAccessors/src/type_utils.jl | 14 + NDTensors/test/Project.toml | 1 - NDTensors/vendor/Project.toml | 5 + NDTensors/vendor/README.md | 13 + NDTensors/vendor/run.jl | 18 + 53 files changed, 1796 insertions(+), 1124 deletions(-) create mode 100755 NDTensors/src/vendored/TypeParameterAccessors/src/TypeParameterAccessors.jl create mode 100755 NDTensors/src/vendored/TypeParameterAccessors/src/base/abstractarray.jl create mode 100755 NDTensors/src/vendored/TypeParameterAccessors/src/base/linearalgebra.jl create mode 100755 NDTensors/src/vendored/TypeParameterAccessors/src/base/similartype.jl create mode 100755 NDTensors/src/vendored/TypeParameterAccessors/src/ndims.jl create mode 100755 NDTensors/src/vendored/TypeParameterAccessors/src/type_parameters.jl create mode 100755 NDTensors/src/vendored/TypeParameterAccessors/src/type_utils.jl create mode 100644 NDTensors/vendor/Project.toml create mode 100644 NDTensors/vendor/README.md create mode 100644 NDTensors/vendor/run.jl diff --git a/NDTensors/Project.toml b/NDTensors/Project.toml index a6f620e419..48409c3a3e 100644 --- a/NDTensors/Project.toml +++ b/NDTensors/Project.toml @@ -1,7 +1,7 @@ name = "NDTensors" uuid = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf" authors = ["Matthew Fishman "] -version = "0.4.12" +version = "0.4.13" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" @@ -28,7 +28,6 @@ Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" -TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [weakdeps] @@ -88,7 +87,6 @@ StridedViews = "0.2.2, 0.3, 0.4" TBLIS = "0.2" TimerOutputs = "0.5.5" TupleTools = "1.2.0" -TypeParameterAccessors = "0.3" VectorInterface = "0.4.2, 0.5" cuTENSOR = "2" julia = "1.10" diff --git a/NDTensors/src/NDTensors.jl b/NDTensors/src/NDTensors.jl index d984c09758..fbb1a320bd 100644 --- a/NDTensors/src/NDTensors.jl +++ b/NDTensors/src/NDTensors.jl @@ -1,4 +1,9 @@ module NDTensors + +module Vendored + include(joinpath("vendored", "TypeParameterAccessors", "src", "TypeParameterAccessors.jl")) +end + ##################################### # Imports and exports # diff --git a/NDTensors/src/abstractarray/generic_array_constructors.jl b/NDTensors/src/abstractarray/generic_array_constructors.jl index 8912c81efa..2bf99eba92 100644 --- a/NDTensors/src/abstractarray/generic_array_constructors.jl +++ b/NDTensors/src/abstractarray/generic_array_constructors.jl @@ -1,8 +1,8 @@ -using TypeParameterAccessors: - unwrap_array_type, - specify_default_type_parameters, - specify_type_parameters, - type_parameters +using .Vendored.TypeParameterAccessors: + unwrap_array_type, + specify_default_type_parameters, + specify_type_parameters, + type_parameters # Convert to Array, avoiding copying if possible array(a::AbstractArray) = a @@ -12,33 +12,33 @@ vector(a::AbstractVector) = a ## Warning to use these functions it is necessary to define `TypeParameterAccessors.position(::Type{<:YourArrayType}, ::typeof(ndims)))` # Implementation, catches if `ndims(arraytype) != length(dims)`. ## TODO convert ndims to `type_parameters(::, typeof(ndims))` -function generic_randn(arraytype::Type{<:AbstractArray}, dims...; rng=Random.default_rng()) - arraytype_specified = specify_type_parameters( - unwrap_array_type(arraytype), ndims, length(dims) - ) - arraytype_specified = specify_default_type_parameters(arraytype_specified) - @assert length(dims) == ndims(arraytype_specified) - data = similar(arraytype_specified, dims...) - return randn!(rng, data) +function generic_randn(arraytype::Type{<:AbstractArray}, dims...; rng = Random.default_rng()) + arraytype_specified = specify_type_parameters( + unwrap_array_type(arraytype), ndims, length(dims) + ) + arraytype_specified = specify_default_type_parameters(arraytype_specified) + @assert length(dims) == ndims(arraytype_specified) + data = similar(arraytype_specified, dims...) + return randn!(rng, data) end function generic_randn( - arraytype::Type{<:AbstractArray}, dims::Tuple; rng=Random.default_rng() -) - return generic_randn(arraytype, dims...; rng) + arraytype::Type{<:AbstractArray}, dims::Tuple; rng = Random.default_rng() + ) + return generic_randn(arraytype, dims...; rng) end # Implementation, catches if `ndims(arraytype) != length(dims)`. function generic_zeros(arraytype::Type{<:AbstractArray}, dims...) - arraytype_specified = specify_type_parameters( - unwrap_array_type(arraytype), ndims, length(dims) - ) - arraytype_specified = specify_default_type_parameters(arraytype_specified) - @assert length(dims) == ndims(arraytype_specified) - ElT = eltype(arraytype_specified) - return fill!(similar(arraytype_specified, dims...), zero(ElT)) + arraytype_specified = specify_type_parameters( + unwrap_array_type(arraytype), ndims, length(dims) + ) + arraytype_specified = specify_default_type_parameters(arraytype_specified) + @assert length(dims) == ndims(arraytype_specified) + ElT = eltype(arraytype_specified) + return fill!(similar(arraytype_specified, dims...), zero(ElT)) end function generic_zeros(arraytype::Type{<:AbstractArray}, dims::Tuple) - return generic_zeros(arraytype, dims...) + return generic_zeros(arraytype, dims...) end diff --git a/NDTensors/src/abstractarray/iscu.jl b/NDTensors/src/abstractarray/iscu.jl index fc7560c1e2..185f7295a9 100644 --- a/NDTensors/src/abstractarray/iscu.jl +++ b/NDTensors/src/abstractarray/iscu.jl @@ -1,7 +1,7 @@ -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type # TODO: Make `isgpu`, `ismtl`, etc. # For `isgpu`, will require a `NDTensorsGPUArrayCoreExt`. iscu(A::AbstractArray) = iscu(typeof(A)) function iscu(A::Type{<:AbstractArray}) - return (unwrap_array_type(A) == A ? false : iscu(unwrap_array_type(A))) + return (unwrap_array_type(A) == A ? false : iscu(unwrap_array_type(A))) end diff --git a/NDTensors/src/abstractarray/set_types.jl b/NDTensors/src/abstractarray/set_types.jl index 85a470ad50..2e393392f6 100644 --- a/NDTensors/src/abstractarray/set_types.jl +++ b/NDTensors/src/abstractarray/set_types.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: TypeParameterAccessors +using .Vendored.TypeParameterAccessors: TypeParameterAccessors """ # Do we still want to define things like this? @@ -13,5 +13,5 @@ TODO: Use `Accessors.jl` notation: # `FillArray` instead. This is a stand-in # to make things work with the current design. function TypeParameterAccessors.set_ndims(numbertype::Type{<:Number}, ndims) - return numbertype + return numbertype end diff --git a/NDTensors/src/abstractarray/similar.jl b/NDTensors/src/abstractarray/similar.jl index dbc8c223ff..9e0cea080e 100644 --- a/NDTensors/src/abstractarray/similar.jl +++ b/NDTensors/src/abstractarray/similar.jl @@ -1,5 +1,5 @@ using Base: DimOrInd, Dims, OneTo -using TypeParameterAccessors: IsWrappedArray, unwrap_array_type, set_eltype, similartype +using .Vendored.TypeParameterAccessors: IsWrappedArray, unwrap_array_type, set_eltype, similartype ## Custom `NDTensors.similar` implementation. ## More extensive than `Base.similar`. @@ -7,8 +7,8 @@ using TypeParameterAccessors: IsWrappedArray, unwrap_array_type, set_eltype, sim # This function actually allocates the data. # NDTensors.similar function similar(arraytype::Type{<:AbstractArray}, dims::Tuple) - shape = NDTensors.to_shape(arraytype, dims) - return similartype(arraytype, shape)(undef, NDTensors.to_shape(arraytype, shape)) + shape = NDTensors.to_shape(arraytype, dims) + return similartype(arraytype, shape)(undef, NDTensors.to_shape(arraytype, shape)) end # This function actually allocates the data. @@ -16,27 +16,27 @@ end # dimensions specified by integers with `Base.to_shape`. # NDTensors.similar function similar(arraytype::Type{<:AbstractArray}, dims::Dims) - return similartype(arraytype, dims)(undef, dims) + return similartype(arraytype, dims)(undef, dims) end # NDTensors.similar function similar(arraytype::Type{<:AbstractArray}, dims::DimOrInd...) - return similar(arraytype, NDTensors.to_shape(dims)) + return similar(arraytype, NDTensors.to_shape(dims)) end # Handles range inputs, `Base.to_shape` converts them to integer dimensions. # See Julia's `base/abstractarray.jl`. # NDTensors.similar function similar( - arraytype::Type{<:AbstractArray}, - shape::Tuple{Union{Integer,OneTo},Vararg{Union{Integer,OneTo}}}, -) - return NDTensors.similar(arraytype, NDTensors.to_shape(shape)) + arraytype::Type{<:AbstractArray}, + shape::Tuple{Union{Integer, OneTo}, Vararg{Union{Integer, OneTo}}}, + ) + return NDTensors.similar(arraytype, NDTensors.to_shape(shape)) end # NDTensors.similar function similar(arraytype::Type{<:AbstractArray}, eltype::Type, dims::Tuple) - return NDTensors.similar(similartype(arraytype, eltype, dims), dims) + return NDTensors.similar(similartype(arraytype, eltype, dims), dims) end # TODO: Add an input `structure` which can store things like the nonzero @@ -70,19 +70,19 @@ end # TODO: Maybe makes an empty array, i.e. `similartype(arraytype, eltype)()`? # NDTensors.similar function similar(arraytype::Type{<:AbstractArray}, eltype::Type) - return error("Must specify dimensions.") + return error("Must specify dimensions.") end ## NDTensors.similar for instances # NDTensors.similar function similar(array::AbstractArray, eltype::Type, dims::Tuple) - return NDTensors.similar(similartype(typeof(array), eltype), dims) + return NDTensors.similar(similartype(typeof(array), eltype), dims) end # NDTensors.similar function similar(array::AbstractArray, eltype::Type, dims::Int) - return NDTensors.similar(similartype(typeof(array), eltype), dims) + return NDTensors.similar(similartype(typeof(array), eltype), dims) end # NDTensors.similar @@ -91,7 +91,7 @@ similar(array::AbstractArray, dims::Tuple) = NDTensors.similar(typeof(array), di # Use the `size` to determine the dimensions # NDTensors.similar function similar(array::AbstractArray, eltype::Type) - return NDTensors.similar(typeof(array), eltype, size(array)) + return NDTensors.similar(typeof(array), eltype, size(array)) end # Use the `size` to determine the dimensions diff --git a/NDTensors/src/adapt.jl b/NDTensors/src/adapt.jl index fb68b40186..ba449bcee4 100644 --- a/NDTensors/src/adapt.jl +++ b/NDTensors/src/adapt.jl @@ -3,7 +3,7 @@ adapt_structure(to, x::TensorStorage) = setdata(x, adapt(to, data(x))) adapt_structure(to, x::Tensor) = setstorage(x, adapt(to, storage(x))) function GPUArraysCoreExtensions.cpu(eltype::Type{<:Number}, x) - return fmap(x -> adapt(Array{eltype}, x), x) + return fmap(x -> adapt(Array{eltype}, x), x) end GPUArraysCoreExtensions.cpu(x) = fmap(x -> adapt(Array, x), x) @@ -27,11 +27,11 @@ double_precision(x) = fmap(x -> adapt(double_precision(eltype(x)), x), x) # Used to adapt `EmptyStorage` types # -using TypeParameterAccessors: specify_type_parameters +using .Vendored.TypeParameterAccessors: specify_type_parameters function adapt_storagetype(to::Type{<:AbstractVector}, x::Type{<:TensorStorage}) - return set_datatype(x, specify_type_parameters(to, eltype, eltype(x))) + return set_datatype(x, specify_type_parameters(to, eltype, eltype(x))) end function adapt_storagetype(to::Type{<:AbstractArray}, x::Type{<:TensorStorage}) - return set_datatype(x, specify_type_parameters(to, (ndims, eltype), (1, eltype(x)))) + return set_datatype(x, specify_type_parameters(to, (ndims, eltype), (1, eltype(x)))) end diff --git a/NDTensors/src/blocksparse/blockdims.jl b/NDTensors/src/blocksparse/blockdims.jl index 8c60ec1b78..f4cfe89de4 100644 --- a/NDTensors/src/blocksparse/blockdims.jl +++ b/NDTensors/src/blocksparse/blockdims.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: TypeParameterAccessors +using .Vendored.TypeParameterAccessors: TypeParameterAccessors """ BlockDim @@ -16,12 +16,12 @@ dim(d::BlockDim) = sum(d) Dimensions used for BlockSparse NDTensors. Each entry lists the block sizes in each dimension. """ -const BlockDims{N} = NTuple{N,BlockDim} +const BlockDims{N} = NTuple{N, BlockDim} Base.ndims(ds::Type{<:BlockDims{N}}) where {N} = N function TypeParameterAccessors.similartype(::Type{<:BlockDims}, ::Type{Val{N}}) where {N} - return BlockDims{N} + return BlockDims{N} end Base.copy(ds::BlockDims) = ds @@ -32,7 +32,7 @@ dim(::BlockDims,::Integer) Return the total extent of the specified dimensions. """ function dim(ds::BlockDims{N}, i::Integer) where {N} - return sum(ds[i]) + return sum(ds[i]) end """ @@ -42,7 +42,7 @@ Return the total extents of the dense space the block dimensions live in. """ function dims(ds::BlockDims{N}) where {N} - return ntuple(i -> dim(ds, i), Val(N)) + return ntuple(i -> dim(ds, i), Val(N)) end """ @@ -52,7 +52,7 @@ Return the total extent of the dense space the block dimensions live in. """ function dim(ds::BlockDims{N}) where {N} - return prod(dims(ds)) + return prod(dims(ds)) end """ @@ -61,7 +61,7 @@ end The number of blocks of the BlockDim. """ function nblocks(ind::BlockDim) - return length(ind) + return length(ind) end """ @@ -70,7 +70,7 @@ end The number of blocks along the diagonal. """ function ndiagblocks(x) - return minimum(nblocks(x)) + return minimum(nblocks(x)) end """ @@ -79,7 +79,7 @@ end The number of blocks in the specified dimension. """ function nblocks(inds::Tuple, i::Integer) - return nblocks(inds[i]) + return nblocks(inds[i]) end """ @@ -87,8 +87,8 @@ end The number of blocks in the specified dimensions. """ -function nblocks(inds::Tuple, is::NTuple{N,Int}) where {N} - return ntuple(i -> nblocks(inds, is[i]), Val(N)) +function nblocks(inds::Tuple, is::NTuple{N, Int}) where {N} + return ntuple(i -> nblocks(inds, is[i]), Val(N)) end """ @@ -97,16 +97,16 @@ end A tuple of the number of blocks in each dimension. """ -function nblocks(inds::NTuple{N,<:Any}) where {N} - return ntuple(i -> nblocks(inds, i), Val(N)) +function nblocks(inds::NTuple{N, <:Any}) where {N} + return ntuple(i -> nblocks(inds, i), Val(N)) end function eachblock(inds::Tuple) - return (Block(b) for b in CartesianIndices(_Tuple(nblocks(inds)))) + return (Block(b) for b in CartesianIndices(_Tuple(nblocks(inds)))) end function eachdiagblock(inds::Tuple) - return (Block(ntuple(_ -> i, length(inds))) for i in 1:ndiagblocks(inds)) + return (Block(ntuple(_ -> i, length(inds))) for i in 1:ndiagblocks(inds)) end """ @@ -116,13 +116,13 @@ The size of the specified block in the specified dimension. """ function blockdim(ind::BlockDim, i::Integer) - return ind[i] + return ind[i] end function blockdim(ind::Integer, i) - return error( - "`blockdim(i::Integer, b)` not currently defined for non-block index $i of type `$(typeof(i))`. In the future this may be defined for `b == Block(1)` or `b == 1` as `dim(i)` and error otherwise.", - ) + return error( + "`blockdim(i::Integer, b)` not currently defined for non-block index $i of type `$(typeof(i))`. In the future this may be defined for `b == Block(1)` or `b == 1` as `dim(i)` and error otherwise.", + ) end """ @@ -132,7 +132,7 @@ The size of the specified block in the specified dimension. """ function blockdim(inds, block, i::Integer) - return blockdim(inds[i], block[i]) + return blockdim(inds[i], block[i]) end """ @@ -141,7 +141,7 @@ end The size of the specified block. """ function blockdims(inds, block) - return ntuple(i -> blockdim(inds, block, i), ValLength(inds)) + return ntuple(i -> blockdim(inds, block, i), ValLength(inds)) end """ @@ -150,7 +150,7 @@ end The total size of the specified block. """ function blockdim(inds, block) - return prod(blockdims(inds, block)) + return prod(blockdims(inds, block)) end """ @@ -159,45 +159,45 @@ end The length of the diagonal of the specified block. """ function blockdiaglength(inds, block) - return minimum(blockdims(inds, block)) + return minimum(blockdims(inds, block)) end function outer(dim1, dim2, dim3, dims...; kwargs...) - return outer(outer(dim1, dim2), dim3, dims...; kwargs...) + return outer(outer(dim1, dim2), dim3, dims...; kwargs...) end function outer(dim1::BlockDim, dim2::BlockDim) - dimR = BlockDim(undef, nblocks(dim1) * nblocks(dim2)) - for (i, t) in enumerate(Iterators.product(dim1, dim2)) - dimR[i] = prod(t) - end - return dimR + dimR = BlockDim(undef, nblocks(dim1) * nblocks(dim2)) + for (i, t) in enumerate(Iterators.product(dim1, dim2)) + dimR[i] = prod(t) + end + return dimR end function permuteblocks(dim::BlockDim, perm) - return dim[perm] + return dim[perm] end # Given a CartesianIndex in the range dims(T), get the block it is in # and the index within that block -function blockindex(T, i::Vararg{Integer,N}) where {N} - # Bounds check. - # Do something more robust like: - # @boundscheck Base.checkbounds_indices(Bool, map(Base.oneto, dims(T)), i) || throw_boundserror(T, i) - @boundscheck any(iszero, i) && Base.throw_boundserror(T, i) - - # Start in the (1,1,...,1) block - current_block_loc = @MVector ones(Int, N) - current_block_dims = blockdims(T, Tuple(current_block_loc)) - block_index = MVector(i) - for dim in 1:N - while block_index[dim] > current_block_dims[dim] - block_index[dim] -= current_block_dims[dim] - current_block_loc[dim] += 1 - current_block_dims = blockdims(T, Tuple(current_block_loc)) +function blockindex(T, i::Vararg{Integer, N}) where {N} + # Bounds check. + # Do something more robust like: + # @boundscheck Base.checkbounds_indices(Bool, map(Base.oneto, dims(T)), i) || throw_boundserror(T, i) + @boundscheck any(iszero, i) && Base.throw_boundserror(T, i) + + # Start in the (1,1,...,1) block + current_block_loc = @MVector ones(Int, N) + current_block_dims = blockdims(T, Tuple(current_block_loc)) + block_index = MVector(i) + for dim in 1:N + while block_index[dim] > current_block_dims[dim] + block_index[dim] -= current_block_dims[dim] + current_block_loc[dim] += 1 + current_block_dims = blockdims(T, Tuple(current_block_loc)) + end end - end - return Tuple(block_index), Block{N}(current_block_loc) + return Tuple(block_index), Block{N}(current_block_loc) end blockindex(T) = (), Block{0}() diff --git a/NDTensors/src/blocksparse/blocksparsetensor.jl b/NDTensors/src/blocksparse/blocksparsetensor.jl index 1bb9fe69fc..89a8011d68 100644 --- a/NDTensors/src/blocksparse/blocksparsetensor.jl +++ b/NDTensors/src/blocksparse/blocksparsetensor.jl @@ -1,5 +1,5 @@ using SparseArrays: nnz -using TypeParameterAccessors: similartype +using .Vendored.TypeParameterAccessors: similartype # # BlockSparseTensor (Tensor using BlockSparse storage) @@ -258,7 +258,7 @@ end # Defaults to adding zeros. # Returns the offset of the new block added. # XXX rename to insertblock!, no need to return offset -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type using .Expose: Exposed, expose, unexpose function insertblock_offset!(T::BlockSparseTensor{ElT, N}, newblock::Block{N}) where {ElT, N} newdim = blockdim(T, newblock) @@ -792,7 +792,7 @@ end # permfactor(perm, block, inds) = 1 -using TypeParameterAccessors: set_type_parameters, parenttype +using .Vendored.TypeParameterAccessors: set_type_parameters, parenttype function permutedims!( R::BlockSparseTensor{<:Number, N}, T::BlockSparseTensor{<:Number, N}, diff --git a/NDTensors/src/blocksparse/diagblocksparse.jl b/NDTensors/src/blocksparse/diagblocksparse.jl index 4fe47ec603..2b52c6434b 100644 --- a/NDTensors/src/blocksparse/diagblocksparse.jl +++ b/NDTensors/src/blocksparse/diagblocksparse.jl @@ -1,5 +1,5 @@ using LinearAlgebra: LinearAlgebra -using TypeParameterAccessors: similartype +using .Vendored.TypeParameterAccessors: similartype export DiagBlockSparse, DiagBlockSparseTensor diff --git a/NDTensors/src/blocksparse/linearalgebra.jl b/NDTensors/src/blocksparse/linearalgebra.jl index dad1312f81..6bcafb8c9b 100644 --- a/NDTensors/src/blocksparse/linearalgebra.jl +++ b/NDTensors/src/blocksparse/linearalgebra.jl @@ -1,29 +1,29 @@ -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type using .Expose: expose -const BlockSparseMatrix{ElT,StoreT,IndsT} = BlockSparseTensor{ElT,2,StoreT,IndsT} -const DiagBlockSparseMatrix{ElT,StoreT,IndsT} = DiagBlockSparseTensor{ElT,2,StoreT,IndsT} -const DiagMatrix{ElT,StoreT,IndsT} = DiagTensor{ElT,2,StoreT,IndsT} +const BlockSparseMatrix{ElT, StoreT, IndsT} = BlockSparseTensor{ElT, 2, StoreT, IndsT} +const DiagBlockSparseMatrix{ElT, StoreT, IndsT} = DiagBlockSparseTensor{ElT, 2, StoreT, IndsT} +const DiagMatrix{ElT, StoreT, IndsT} = DiagTensor{ElT, 2, StoreT, IndsT} function _truncated_blockdim( - S::DiagMatrix, docut::Real; singular_values=false, truncate=true, min_blockdim=nothing -) - min_blockdim = replace_nothing(min_blockdim, 0) - # TODO: Replace `cpu` with `Expose` dispatch. - S = cpu(S) - full_dim = diaglength(S) - !truncate && return full_dim - min_blockdim = min(min_blockdim, full_dim) - newdim = 0 - val = singular_values ? getdiagindex(S, newdim + 1)^2 : abs(getdiagindex(S, newdim + 1)) - while newdim + 1 ≤ full_dim && val > docut - newdim += 1 - if newdim + 1 ≤ full_dim - val = - singular_values ? getdiagindex(S, newdim + 1)^2 : abs(getdiagindex(S, newdim + 1)) + S::DiagMatrix, docut::Real; singular_values = false, truncate = true, min_blockdim = nothing + ) + min_blockdim = replace_nothing(min_blockdim, 0) + # TODO: Replace `cpu` with `Expose` dispatch. + S = cpu(S) + full_dim = diaglength(S) + !truncate && return full_dim + min_blockdim = min(min_blockdim, full_dim) + newdim = 0 + val = singular_values ? getdiagindex(S, newdim + 1)^2 : abs(getdiagindex(S, newdim + 1)) + while newdim + 1 ≤ full_dim && val > docut + newdim += 1 + if newdim + 1 ≤ full_dim + val = + singular_values ? getdiagindex(S, newdim + 1)^2 : abs(getdiagindex(S, newdim + 1)) + end end - end - (newdim >= min_blockdim) || (newdim = min_blockdim) - return newdim + (newdim >= min_blockdim) || (newdim = min_blockdim) + return newdim end """ @@ -37,391 +37,391 @@ This assumption makes it so the result can be computed from the dense svds of seperate blocks. """ function svd( - T::Tensor{ElT,2,<:BlockSparse}; - min_blockdim=nothing, - mindim=nothing, - maxdim=nothing, - cutoff=nothing, - alg=nothing, - use_absolute_cutoff=nothing, - use_relative_cutoff=nothing, -) where {ElT} - Us = Vector{DenseTensor{ElT,2}}(undef, nnzblocks(T)) - Ss = Vector{DiagTensor{real(ElT),2}}(undef, nnzblocks(T)) - Vs = Vector{DenseTensor{ElT,2}}(undef, nnzblocks(T)) - - # Sorted eigenvalues - d = Vector{real(ElT)}() - - for (n, b) in enumerate(eachnzblock(T)) - blockT = blockview(T, b) - USVb = svd(blockT; alg) - if isnothing(USVb) - return nothing + T::Tensor{ElT, 2, <:BlockSparse}; + min_blockdim = nothing, + mindim = nothing, + maxdim = nothing, + cutoff = nothing, + alg = nothing, + use_absolute_cutoff = nothing, + use_relative_cutoff = nothing, + ) where {ElT} + Us = Vector{DenseTensor{ElT, 2}}(undef, nnzblocks(T)) + Ss = Vector{DiagTensor{real(ElT), 2}}(undef, nnzblocks(T)) + Vs = Vector{DenseTensor{ElT, 2}}(undef, nnzblocks(T)) + + # Sorted eigenvalues + d = Vector{real(ElT)}() + + for (n, b) in enumerate(eachnzblock(T)) + blockT = blockview(T, b) + USVb = svd(blockT; alg) + if isnothing(USVb) + return nothing + end + Ub, Sb, Vb = USVb + Us[n] = Ub + Ss[n] = Sb + Vs[n] = Vb + # Previously this was: + # vector(diag(Sb)) + # But it broke, did `diag(::Tensor)` change types? + # TODO: call this a function `diagonal`, i.e.: + # https://github.com/JuliaLang/julia/issues/30250 + # or make `diag(::Tensor)` return a view by default. + append!(expose(d), data(Sb)) end - Ub, Sb, Vb = USVb - Us[n] = Ub - Ss[n] = Sb - Vs[n] = Vb - # Previously this was: - # vector(diag(Sb)) - # But it broke, did `diag(::Tensor)` change types? - # TODO: call this a function `diagonal`, i.e.: - # https://github.com/JuliaLang/julia/issues/30250 - # or make `diag(::Tensor)` return a view by default. - append!(expose(d), data(Sb)) - end - - # Square the singular values to get - # the eigenvalues - d .= d .^ 2 - sort!(d; rev=true) - - # Get the list of blocks of T - # that are not dropped - nzblocksT = nzblocks(T) - - dropblocks = Int[] - if any(!isnothing, (maxdim, cutoff)) - d, truncerr, docut = truncate!!( - d; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff - ) - for n in 1:nnzblocks(T) - blockdim = _truncated_blockdim( - Ss[n], docut; min_blockdim, singular_values=true, truncate=true - ) - if blockdim == 0 - push!(dropblocks, n) - else - # TODO: Replace call to `data` with `diagview`. - Strunc = tensor(Diag(data(Ss[n])[1:blockdim]), (blockdim, blockdim)) - Us[n] = Us[n][1:dim(Us[n], 1), 1:blockdim] - Ss[n] = Strunc - Vs[n] = Vs[n][1:dim(Vs[n], 1), 1:blockdim] - end + + # Square the singular values to get + # the eigenvalues + d .= d .^ 2 + sort!(d; rev = true) + + # Get the list of blocks of T + # that are not dropped + nzblocksT = nzblocks(T) + + dropblocks = Int[] + if any(!isnothing, (maxdim, cutoff)) + d, truncerr, docut = truncate!!( + d; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff + ) + for n in 1:nnzblocks(T) + blockdim = _truncated_blockdim( + Ss[n], docut; min_blockdim, singular_values = true, truncate = true + ) + if blockdim == 0 + push!(dropblocks, n) + else + # TODO: Replace call to `data` with `diagview`. + Strunc = tensor(Diag(data(Ss[n])[1:blockdim]), (blockdim, blockdim)) + Us[n] = Us[n][1:dim(Us[n], 1), 1:blockdim] + Ss[n] = Strunc + Vs[n] = Vs[n][1:dim(Vs[n], 1), 1:blockdim] + end + end + deleteat!(Us, dropblocks) + deleteat!(Ss, dropblocks) + deleteat!(Vs, dropblocks) + deleteat!(nzblocksT, dropblocks) + else + truncerr, docut = 0.0, 0.0 end - deleteat!(Us, dropblocks) - deleteat!(Ss, dropblocks) - deleteat!(Vs, dropblocks) - deleteat!(nzblocksT, dropblocks) - else - truncerr, docut = 0.0, 0.0 - end - - # The number of non-zero blocks of T remaining - nnzblocksT = length(nzblocksT) - - # - # Make indices of U and V - # that connect to S - # - i1 = ind(T, 1) - i2 = ind(T, 2) - uind = dag(sim(i1)) - vind = dag(sim(i2)) - resize!(uind, nnzblocksT) - resize!(vind, nnzblocksT) - for (n, blockT) in enumerate(nzblocksT) - Udim = size(Us[n], 2) - b1 = block(i1, blockT[1]) - setblock!(uind, resize(b1, Udim), n) - Vdim = size(Vs[n], 2) - b2 = block(i2, blockT[2]) - setblock!(vind, resize(b2, Vdim), n) - end - - # - # Put the blocks into U,S,V - # - - nzblocksU = Vector{Block{2}}(undef, nnzblocksT) - nzblocksS = Vector{Block{2}}(undef, nnzblocksT) - nzblocksV = Vector{Block{2}}(undef, nnzblocksT) - - for (n, blockT) in enumerate(nzblocksT) - blockU = (blockT[1], UInt(n)) - nzblocksU[n] = blockU - - blockS = (n, n) - nzblocksS[n] = blockS - - blockV = (blockT[2], UInt(n)) - nzblocksV[n] = blockV - end - - indsU = setindex(inds(T), uind, 2) - - indsV = setindex(inds(T), vind, 1) - indsV = permute(indsV, (2, 1)) - - indsS = setindex(inds(T), dag(uind), 1) - indsS = setindex(indsS, dag(vind), 2) - - U = BlockSparseTensor(unwrap_array_type(T), undef, nzblocksU, indsU) - S = DiagBlockSparseTensor( - set_eltype(unwrap_array_type(T), real(ElT)), undef, nzblocksS, indsS - ) - V = BlockSparseTensor(unwrap_array_type(T), undef, nzblocksV, indsV) - - for n in 1:nnzblocksT - Ub, Sb, Vb = Us[n], Ss[n], Vs[n] - - blockU = nzblocksU[n] - blockS = nzblocksS[n] - blockV = nzblocksV[n] - - if VERSION < v"1.5" - # In v1.3 and v1.4 of Julia, Ub has - # a very complicated view wrapper that - # can't be handled efficiently - Ub = copy(Ub) - Vb = copy(Vb) + + # The number of non-zero blocks of T remaining + nnzblocksT = length(nzblocksT) + + # + # Make indices of U and V + # that connect to S + # + i1 = ind(T, 1) + i2 = ind(T, 2) + uind = dag(sim(i1)) + vind = dag(sim(i2)) + resize!(uind, nnzblocksT) + resize!(vind, nnzblocksT) + for (n, blockT) in enumerate(nzblocksT) + Udim = size(Us[n], 2) + b1 = block(i1, blockT[1]) + setblock!(uind, resize(b1, Udim), n) + Vdim = size(Vs[n], 2) + b2 = block(i2, blockT[2]) + setblock!(vind, resize(b2, Vdim), n) end - # - sU = right_arrow_sign(uind, blockU[2]) + # + # Put the blocks into U,S,V + # - if sU == -1 - Ub *= -1 - end - copyto!(expose(blockview(U, blockU)), expose(Ub)) - - blockviewS = blockview(S, blockS) - # TODO: Replace `data` with `diagview`. - copyto!(expose(data(blockviewS)), expose(data(Sb))) - - # - sV = left_arrow_sign(vind, blockV[2]) - # This sign (sVP) accounts for the fact that - # V is transposed, i.e. the index connecting to S - # is the second index: - sVP = 1 - if using_auto_fermion() - sVP = -block_sign(vind, blockV[2]) + nzblocksU = Vector{Block{2}}(undef, nnzblocksT) + nzblocksS = Vector{Block{2}}(undef, nnzblocksT) + nzblocksV = Vector{Block{2}}(undef, nnzblocksT) + + for (n, blockT) in enumerate(nzblocksT) + blockU = (blockT[1], UInt(n)) + nzblocksU[n] = blockU + + blockS = (n, n) + nzblocksS[n] = blockS + + blockV = (blockT[2], UInt(n)) + nzblocksV[n] = blockV end - if (sV * sVP) == -1 - Vb *= -1 + indsU = setindex(inds(T), uind, 2) + + indsV = setindex(inds(T), vind, 1) + indsV = permute(indsV, (2, 1)) + + indsS = setindex(inds(T), dag(uind), 1) + indsS = setindex(indsS, dag(vind), 2) + + U = BlockSparseTensor(unwrap_array_type(T), undef, nzblocksU, indsU) + S = DiagBlockSparseTensor( + set_eltype(unwrap_array_type(T), real(ElT)), undef, nzblocksS, indsS + ) + V = BlockSparseTensor(unwrap_array_type(T), undef, nzblocksV, indsV) + + for n in 1:nnzblocksT + Ub, Sb, Vb = Us[n], Ss[n], Vs[n] + + blockU = nzblocksU[n] + blockS = nzblocksS[n] + blockV = nzblocksV[n] + + if VERSION < v"1.5" + # In v1.3 and v1.4 of Julia, Ub has + # a very complicated view wrapper that + # can't be handled efficiently + Ub = copy(Ub) + Vb = copy(Vb) + end + + # + sU = right_arrow_sign(uind, blockU[2]) + + if sU == -1 + Ub *= -1 + end + copyto!(expose(blockview(U, blockU)), expose(Ub)) + + blockviewS = blockview(S, blockS) + # TODO: Replace `data` with `diagview`. + copyto!(expose(data(blockviewS)), expose(data(Sb))) + + # + sV = left_arrow_sign(vind, blockV[2]) + # This sign (sVP) accounts for the fact that + # V is transposed, i.e. the index connecting to S + # is the second index: + sVP = 1 + if using_auto_fermion() + sVP = -block_sign(vind, blockV[2]) + end + + if (sV * sVP) == -1 + Vb *= -1 + end + copyto!(blockview(V, blockV), Vb) end - copyto!(blockview(V, blockV), Vb) - end - return U, S, V, Spectrum(d, truncerr) + return U, S, V, Spectrum(d, truncerr) end -_eigen_eltypes(T::Hermitian{ElT,<:BlockSparseMatrix{ElT}}) where {ElT} = real(ElT), ElT +_eigen_eltypes(T::Hermitian{ElT, <:BlockSparseMatrix{ElT}}) where {ElT} = real(ElT), ElT _eigen_eltypes(T::BlockSparseMatrix{ElT}) where {ElT} = complex(ElT), complex(ElT) function LinearAlgebra.eigen( - T::Union{Hermitian{ElT,<:Tensor{ElT,2,<:BlockSparse}},Tensor{ElT,2,<:BlockSparse}}; - min_blockdim=nothing, - mindim=nothing, - maxdim=nothing, - cutoff=nothing, - use_absolute_cutoff=nothing, - use_relative_cutoff=nothing, -) where {ElT<:Union{Real,Complex}} - ElD, ElV = _eigen_eltypes(T) - - # Sorted eigenvalues - d = Vector{real(ElT)}() - - for b in eachnzblock(T) - all(==(b[1]), b) || error("Eigen currently only supports block diagonal matrices.") - end - - b = first(eachnzblock(T)) - blockT = blockview(T, b) - Db, Vb = eigen(expose(blockT)) - Ds = [Db] - Vs = [Vb] - append!(expose(d), abs.(data(Db))) - for (n, b) in enumerate(eachnzblock(T)) - n == 1 && continue + T::Union{Hermitian{ElT, <:Tensor{ElT, 2, <:BlockSparse}}, Tensor{ElT, 2, <:BlockSparse}}; + min_blockdim = nothing, + mindim = nothing, + maxdim = nothing, + cutoff = nothing, + use_absolute_cutoff = nothing, + use_relative_cutoff = nothing, + ) where {ElT <: Union{Real, Complex}} + ElD, ElV = _eigen_eltypes(T) + + # Sorted eigenvalues + d = Vector{real(ElT)}() + + for b in eachnzblock(T) + all(==(b[1]), b) || error("Eigen currently only supports block diagonal matrices.") + end + + b = first(eachnzblock(T)) blockT = blockview(T, b) Db, Vb = eigen(expose(blockT)) - push!(Ds, Db) - push!(Vs, Vb) + Ds = [Db] + Vs = [Vb] append!(expose(d), abs.(data(Db))) - end - - dropblocks = Int[] - sort!(d; rev=true, by=abs) + for (n, b) in enumerate(eachnzblock(T)) + n == 1 && continue + blockT = blockview(T, b) + Db, Vb = eigen(expose(blockT)) + push!(Ds, Db) + push!(Vs, Vb) + append!(expose(d), abs.(data(Db))) + end - if any(!isnothing, (maxdim, cutoff)) - d, truncerr, docut = truncate!!( - d; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff - ) - for n in 1:nnzblocks(T) - blockdim = _truncated_blockdim( - Ds[n], docut; min_blockdim, singular_values=false, truncate=true - ) - if blockdim == 0 - push!(dropblocks, n) - else - # TODO: Replace call to `data` with `diagview`. - Dtrunc = tensor(Diag(data(Ds[n])[1:blockdim]), (blockdim, blockdim)) - Ds[n] = Dtrunc - new_size = (dim(Vs[n], 1), blockdim) - new_data = array(Vs[n])[1:new_size[1], 1:new_size[2]] - Vs[n] = tensor(Dense(new_data), new_size) - end + dropblocks = Int[] + sort!(d; rev = true, by = abs) + + if any(!isnothing, (maxdim, cutoff)) + d, truncerr, docut = truncate!!( + d; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff + ) + for n in 1:nnzblocks(T) + blockdim = _truncated_blockdim( + Ds[n], docut; min_blockdim, singular_values = false, truncate = true + ) + if blockdim == 0 + push!(dropblocks, n) + else + # TODO: Replace call to `data` with `diagview`. + Dtrunc = tensor(Diag(data(Ds[n])[1:blockdim]), (blockdim, blockdim)) + Ds[n] = Dtrunc + new_size = (dim(Vs[n], 1), blockdim) + new_data = array(Vs[n])[1:new_size[1], 1:new_size[2]] + Vs[n] = tensor(Dense(new_data), new_size) + end + end + deleteat!(Ds, dropblocks) + deleteat!(Vs, dropblocks) + else + truncerr = 0.0 end - deleteat!(Ds, dropblocks) - deleteat!(Vs, dropblocks) - else - truncerr = 0.0 - end - # Get the list of blocks of T - # that are not dropped - nzblocksT = nzblocks(T) - deleteat!(nzblocksT, dropblocks) + # Get the list of blocks of T + # that are not dropped + nzblocksT = nzblocks(T) + deleteat!(nzblocksT, dropblocks) - # The number of blocks of T remaining - nnzblocksT = nnzblocks(T) - length(dropblocks) + # The number of blocks of T remaining + nnzblocksT = nnzblocks(T) - length(dropblocks) - # - # Put the blocks into D, V - # + # + # Put the blocks into D, V + # - i1, i2 = inds(T) - l = sim(i1) + i1, i2 = inds(T) + l = sim(i1) - lkeepblocks = Int[bT[1] for bT in nzblocksT] - ldropblocks = setdiff(1:nblocks(l), lkeepblocks) - deleteat!(l, ldropblocks) + lkeepblocks = Int[bT[1] for bT in nzblocksT] + ldropblocks = setdiff(1:nblocks(l), lkeepblocks) + deleteat!(l, ldropblocks) - # l may have too many blocks - (nblocks(l) > nnzblocksT) && error("New index l in eigen has too many blocks") + # l may have too many blocks + (nblocks(l) > nnzblocksT) && error("New index l in eigen has too many blocks") - # Truncation may have changed - # some block sizes - for n in 1:nnzblocksT - setblockdim!(l, minimum(dims(Ds[n])), n) - end + # Truncation may have changed + # some block sizes + for n in 1:nnzblocksT + setblockdim!(l, minimum(dims(Ds[n])), n) + end - r = dag(sim(l)) + r = dag(sim(l)) - indsD = (l, r) - indsV = (dag(i2), r) + indsD = (l, r) + indsV = (dag(i2), r) - nzblocksD = Vector{Block{2}}(undef, nnzblocksT) - nzblocksV = Vector{Block{2}}(undef, nnzblocksT) - for n in 1:nnzblocksT - blockT = nzblocksT[n] + nzblocksD = Vector{Block{2}}(undef, nnzblocksT) + nzblocksV = Vector{Block{2}}(undef, nnzblocksT) + for n in 1:nnzblocksT + blockT = nzblocksT[n] - blockD = (n, n) - nzblocksD[n] = blockD + blockD = (n, n) + nzblocksD[n] = blockD - blockV = (blockT[1], n) - nzblocksV[n] = blockV - end + blockV = (blockT[1], n) + nzblocksV[n] = blockV + end - D = DiagBlockSparseTensor( - set_ndims(set_eltype(unwrap_array_type(T), ElD), 1), undef, nzblocksD, indsD - ) - V = BlockSparseTensor(set_eltype(unwrap_array_type(T), ElV), undef, nzblocksV, indsV) + D = DiagBlockSparseTensor( + set_ndims(set_eltype(unwrap_array_type(T), ElD), 1), undef, nzblocksD, indsD + ) + V = BlockSparseTensor(set_eltype(unwrap_array_type(T), ElV), undef, nzblocksV, indsV) - for n in 1:nnzblocksT - Db, Vb = Ds[n], Vs[n] + for n in 1:nnzblocksT + Db, Vb = Ds[n], Vs[n] - blockD = nzblocksD[n] - blockviewD = blockview(D, blockD) - # TODO: Replace `data` with `diagview`. - copyto!(expose(data(blockviewD)), expose(data(Db))) + blockD = nzblocksD[n] + blockviewD = blockview(D, blockD) + # TODO: Replace `data` with `diagview`. + copyto!(expose(data(blockviewD)), expose(data(Db))) - blockV = nzblocksV[n] - copyto!(blockview(V, blockV), Vb) - end + blockV = nzblocksV[n] + copyto!(blockview(V, blockV), Vb) + end - return D, V, Spectrum(d, truncerr) + return D, V, Spectrum(d, truncerr) end -Expose.ql(T::BlockSparseTensor{<:Any,2}; kwargs...) = qx(ql, T; kwargs...) -qr(T::BlockSparseTensor{<:Any,2}; kwargs...) = qx(qr, T; kwargs...) +Expose.ql(T::BlockSparseTensor{<:Any, 2}; kwargs...) = qx(ql, T; kwargs...) +qr(T::BlockSparseTensor{<:Any, 2}; kwargs...) = qx(qr, T; kwargs...) # # Generic function to implelement blocks sparse qr/ql decomposition. It calls # the dense qr or ql for each block. The X tensor = R or L. # This code thanks to Niklas Tausendpfund # https://github.com/ntausend/variance_iTensor/blob/main/Hubig_variance_test.ipynb # -function qx(qx::Function, T::BlockSparseTensor{<:Any,2}; positive=nothing) - ElT = eltype(T) - # getting total number of blocks - nnzblocksT = nnzblocks(T) - nzblocksT = nzblocks(T) +function qx(qx::Function, T::BlockSparseTensor{<:Any, 2}; positive = nothing) + ElT = eltype(T) + # getting total number of blocks + nnzblocksT = nnzblocks(T) + nzblocksT = nzblocks(T) + + Qs = Vector{DenseTensor{ElT, 2}}(undef, nnzblocksT) + Xs = Vector{DenseTensor{ElT, 2}}(undef, nnzblocksT) + + for (jj, b) in enumerate(eachnzblock(T)) + blockT = blockview(T, b) + QXb = qx(blockT; positive) + + if (isnothing(QXb)) + return nothing + end + + Q, X = QXb + Qs[jj] = Q + Xs[jj] = X + end - Qs = Vector{DenseTensor{ElT,2}}(undef, nnzblocksT) - Xs = Vector{DenseTensor{ElT,2}}(undef, nnzblocksT) + # + # Make the new index connecting Q and R + # + itl = ind(T, 1) #left index of T + iq = dag(sim(itl)) #start with similar to the left index of T + resize!(iq, nnzblocksT) #adjust the size to match the block count + for (n, blockT) in enumerate(nzblocksT) + Qdim = size(Qs[n], 2) #get the block dim on right side of Q. + b1 = block(itl, blockT[1]) + setblock!(iq, resize(b1, Qdim), n) + end - for (jj, b) in enumerate(eachnzblock(T)) - blockT = blockview(T, b) - QXb = qx(blockT; positive) + indsQ = setindex(inds(T), iq, 2) + indsX = setindex(inds(T), dag(iq), 1) + + nzblocksQ = Vector{Block{2}}(undef, nnzblocksT) + nzblocksX = Vector{Block{2}}(undef, nnzblocksT) - if (isnothing(QXb)) - return nothing + for n in 1:nnzblocksT + blockT = nzblocksT[n] + nzblocksQ[n] = (blockT[1], UInt(n)) + nzblocksX[n] = (UInt(n), blockT[2]) end - Q, X = QXb - Qs[jj] = Q - Xs[jj] = X - end - - # - # Make the new index connecting Q and R - # - itl = ind(T, 1) #left index of T - iq = dag(sim(itl)) #start with similar to the left index of T - resize!(iq, nnzblocksT) #adjust the size to match the block count - for (n, blockT) in enumerate(nzblocksT) - Qdim = size(Qs[n], 2) #get the block dim on right side of Q. - b1 = block(itl, blockT[1]) - setblock!(iq, resize(b1, Qdim), n) - end - - indsQ = setindex(inds(T), iq, 2) - indsX = setindex(inds(T), dag(iq), 1) - - nzblocksQ = Vector{Block{2}}(undef, nnzblocksT) - nzblocksX = Vector{Block{2}}(undef, nnzblocksT) - - for n in 1:nnzblocksT - blockT = nzblocksT[n] - nzblocksQ[n] = (blockT[1], UInt(n)) - nzblocksX[n] = (UInt(n), blockT[2]) - end - - Q = BlockSparseTensor(unwrap_array_type(T), undef, nzblocksQ, indsQ) - X = BlockSparseTensor(unwrap_array_type(T), undef, nzblocksX, indsX) - - for n in 1:nnzblocksT - copyto!(blockview(Q, nzblocksQ[n]), Qs[n]) - copyto!(blockview(X, nzblocksX[n]), Xs[n]) - end - - Q = adapt(unwrap_array_type(T), Q) - X = adapt(unwrap_array_type(T), X) - return Q, X + Q = BlockSparseTensor(unwrap_array_type(T), undef, nzblocksQ, indsQ) + X = BlockSparseTensor(unwrap_array_type(T), undef, nzblocksX, indsX) + + for n in 1:nnzblocksT + copyto!(blockview(Q, nzblocksQ[n]), Qs[n]) + copyto!(blockview(X, nzblocksX[n]), Xs[n]) + end + + Q = adapt(unwrap_array_type(T), Q) + X = adapt(unwrap_array_type(T), X) + return Q, X end function exp( - T::Union{BlockSparseMatrix{ElT},Hermitian{ElT,<:BlockSparseMatrix{ElT}}} -) where {ElT<:Union{Real,Complex}} - expT = BlockSparseTensor(ElT, undef, nzblocks(T), inds(T)) - for b in eachnzblock(T) - all(==(b[1]), b) || error("exp currently supports only block-diagonal matrices") - end - for b in eachdiagblock(T) - blockT = blockview(T, b) - if isnothing(blockT) - # Block was not found in the list, treat as 0 - id_block = Matrix{ElT}(I, blockdims(T, b)) - insertblock!(expT, b) - blockview(expT, b) .= id_block - else - blockview(expT, b) .= exp(blockT) + T::Union{BlockSparseMatrix{ElT}, Hermitian{ElT, <:BlockSparseMatrix{ElT}}} + ) where {ElT <: Union{Real, Complex}} + expT = BlockSparseTensor(ElT, undef, nzblocks(T), inds(T)) + for b in eachnzblock(T) + all(==(b[1]), b) || error("exp currently supports only block-diagonal matrices") + end + for b in eachdiagblock(T) + blockT = blockview(T, b) + if isnothing(blockT) + # Block was not found in the list, treat as 0 + id_block = Matrix{ElT}(I, blockdims(T, b)) + insertblock!(expT, b) + blockview(expT, b) .= id_block + else + blockview(expT, b) .= exp(blockT) + end end - end - return expT + return expT end diff --git a/NDTensors/src/blocksparse/similar.jl b/NDTensors/src/blocksparse/similar.jl index d9dc3b5e37..a850710e0d 100644 --- a/NDTensors/src/blocksparse/similar.jl +++ b/NDTensors/src/blocksparse/similar.jl @@ -1,35 +1,35 @@ using SparseArrays: nnz -using TypeParameterAccessors: similartype +using .Vendored.TypeParameterAccessors: similartype # NDTensors.similar function similar(storagetype::Type{<:BlockSparse}, blockoffsets::BlockOffsets, dims::Tuple) - data = similar(datatype(storagetype), nnz(blockoffsets, dims)) - return BlockSparse(data, blockoffsets) + data = similar(datatype(storagetype), nnz(blockoffsets, dims)) + return BlockSparse(data, blockoffsets) end # NDTensors.similar function similar(storagetype::Type{<:BlockSparse}, dims::Tuple) - # Create an empty BlockSparse storage - return similartype(storagetype, dims)() + # Create an empty BlockSparse storage + return similartype(storagetype, dims)() end # NDTensors.similar function similar(storagetype::Type{<:BlockSparse}, dims::Dims) - # Create an empty BlockSparse storage - return similartype(storagetype, dims)() + # Create an empty BlockSparse storage + return similartype(storagetype, dims)() end ## TODO: Is there a way to make this generic? # NDTensors.similar function similar( - tensortype::Type{<:BlockSparseTensor}, blockoffsets::BlockOffsets, dims::Tuple -) - return Tensor(similar(storagetype(tensortype), blockoffsets, dims), dims) + tensortype::Type{<:BlockSparseTensor}, blockoffsets::BlockOffsets, dims::Tuple + ) + return Tensor(similar(storagetype(tensortype), blockoffsets, dims), dims) end # NDTensors.similar function similar(tensor::BlockSparseTensor, blockoffsets::BlockOffsets, dims::Tuple) - return similar(typeof(tensor), blockoffsets, dims) + return similar(typeof(tensor), blockoffsets, dims) end ## ## TODO: Determine if the methods below are needed. diff --git a/NDTensors/src/default_kwargs.jl b/NDTensors/src/default_kwargs.jl index fda4fe9786..9cb170b0e5 100644 --- a/NDTensors/src/default_kwargs.jl +++ b/NDTensors/src/default_kwargs.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type replace_nothing(::Nothing, replacement) = replacement replace_nothing(value, replacement) = value diff --git a/NDTensors/src/dense/dense.jl b/NDTensors/src/dense/dense.jl index 0cd079f2b6..abf18e4a47 100644 --- a/NDTensors/src/dense/dense.jl +++ b/NDTensors/src/dense/dense.jl @@ -2,94 +2,94 @@ # Dense storage # -struct Dense{ElT,DataT<:AbstractArray} <: TensorStorage{ElT} - data::DataT - function Dense{ElT,DataT}(data::DataT) where {ElT,DataT<:AbstractVector} - @assert ElT == eltype(DataT) - return new{ElT,DataT}(data) - end +struct Dense{ElT, DataT <: AbstractArray} <: TensorStorage{ElT} + data::DataT + function Dense{ElT, DataT}(data::DataT) where {ElT, DataT <: AbstractVector} + @assert ElT == eltype(DataT) + return new{ElT, DataT}(data) + end - function Dense{ElT,DataT}(data::DataT) where {ElT,DataT<:AbstractArray} - println("Only Vector-based datatypes are currently supported.") - throw(TypeError) - end + function Dense{ElT, DataT}(data::DataT) where {ElT, DataT <: AbstractArray} + println("Only Vector-based datatypes are currently supported.") + throw(TypeError) + end end #Start with high information constructors and move to low information constructors -function Dense{ElT,DataT}() where {ElT,DataT<:AbstractArray} - return Dense{ElT,DataT}(DataT()) +function Dense{ElT, DataT}() where {ElT, DataT <: AbstractArray} + return Dense{ElT, DataT}(DataT()) end # Construct from a set of indices # This will fail if zero(ElT) is not defined for the ElT -function Dense{ElT,DataT}(inds::Tuple) where {ElT,DataT<:AbstractArray} - return Dense{ElT,DataT}(generic_zeros(DataT, dim(inds))) +function Dense{ElT, DataT}(inds::Tuple) where {ElT, DataT <: AbstractArray} + return Dense{ElT, DataT}(generic_zeros(DataT, dim(inds))) end -function Dense{ElT,DataT}(dim::Integer) where {ElT,DataT<:AbstractArray} - return Dense{ElT,DataT}(generic_zeros(DataT, dim)) +function Dense{ElT, DataT}(dim::Integer) where {ElT, DataT <: AbstractArray} + return Dense{ElT, DataT}(generic_zeros(DataT, dim)) end -function Dense{ElT,DataT}(::UndefInitializer, inds::Tuple) where {ElT,DataT<:AbstractArray} - return Dense{ElT,DataT}(similar(DataT, dim(inds))) +function Dense{ElT, DataT}(::UndefInitializer, inds::Tuple) where {ElT, DataT <: AbstractArray} + return Dense{ElT, DataT}(similar(DataT, dim(inds))) end -function Dense{ElT,DataT}(x, dim::Integer) where {ElT,DataT<:AbstractVector} - return Dense{ElT,DataT}(fill!(similar(DataT, dim), ElT(x))) +function Dense{ElT, DataT}(x, dim::Integer) where {ElT, DataT <: AbstractVector} + return Dense{ElT, DataT}(fill!(similar(DataT, dim), ElT(x))) end -function Dense{ElR,DataT}(data::AbstractArray) where {ElR,DataT<:AbstractArray} - data = convert(DataT, data) - return Dense{ElR,DataT}(data) +function Dense{ElR, DataT}(data::AbstractArray) where {ElR, DataT <: AbstractArray} + data = convert(DataT, data) + return Dense{ElR, DataT}(data) end # This function is ill-defined. It cannot transform a complex type to real... -function Dense{ElR}(data::AbstractArray{ElT}) where {ElR,ElT} - return Dense{ElR}(convert(similartype(typeof(data), ElR), data)) +function Dense{ElR}(data::AbstractArray{ElT}) where {ElR, ElT} + return Dense{ElR}(convert(similartype(typeof(data), ElR), data)) end function Dense{ElT}(data::AbstractArray{ElT}) where {ElT} - return Dense{ElT,typeof(data)}(data) + return Dense{ElT, typeof(data)}(data) end function Dense{ElT}(inds::Tuple) where {ElT} - return Dense{ElT}(dim(inds)) + return Dense{ElT}(dim(inds)) end function Dense{ElT}(dim::Integer) where {ElT} - return Dense{ElT,default_datatype(ElT)}(dim) + return Dense{ElT, default_datatype(ElT)}(dim) end -Dense{ElT}() where {ElT} = Dense{ElT,default_datatype(ElT)}() +Dense{ElT}() where {ElT} = Dense{ElT, default_datatype(ElT)}() function Dense(data::AbstractVector) - return Dense{eltype(data)}(data) + return Dense{eltype(data)}(data) end -function Dense(data::DataT) where {DataT<:AbstractArray{<:Any,N}} where {N} - #println("Warning: Only vector based datatypes are currenlty supported by Dense. The data structure provided will be vectorized.") - return Dense(vec(data)) +function Dense(data::DataT) where {DataT <: AbstractArray{<:Any, N}} where {N} + #println("Warning: Only vector based datatypes are currenlty supported by Dense. The data structure provided will be vectorized.") + return Dense(vec(data)) end function Dense(DataT::Type{<:AbstractArray}, dim::Integer) - ElT = eltype(DataT) - return Dense{ElT,DataT}(dim) + ElT = eltype(DataT) + return Dense{ElT, DataT}(dim) end Dense(ElT::Type{<:Number}, dim::Integer) = Dense{ElT}(dim) function Dense(ElT::Type{<:Number}, ::UndefInitializer, dim::Integer) - return Dense{ElT,default_datatype(ElT)}(undef, (dim,)) + return Dense{ElT, default_datatype(ElT)}(undef, (dim,)) end function Dense(::UndefInitializer, dim::Integer) - datatype = default_datatype() - return Dense{eltype(datatype),datatype}(undef, (dim,)) + datatype = default_datatype() + return Dense{eltype(datatype), datatype}(undef, (dim,)) end function Dense(x::Number, dim::Integer) - ElT = typeof(x) - return Dense{ElT,default_datatype(ElT)}(x, dim) + ElT = typeof(x) + return Dense{ElT, default_datatype(ElT)}(x, dim) end Dense(dim::Integer) = Dense(default_eltype(), dim) @@ -102,47 +102,47 @@ setdata(D::Dense, ndata) = Dense(ndata) setdata(storagetype::Type{<:Dense}, data) = Dense(data) function copy(D::Dense) - return Dense(copy(expose(data(D)))) + return Dense(copy(expose(data(D)))) end function Base.copyto!(R::Dense, T::Dense) - copyto!(expose(data(R)), expose(data(T))) - return R + copyto!(expose(data(R)), expose(data(T))) + return R end function Base.real(T::Type{<:Dense}) - return set_datatype(T, similartype(datatype(T), real(eltype(T)))) + return set_datatype(T, similartype(datatype(T), real(eltype(T)))) end function complex(T::Type{<:Dense}) - return set_datatype(T, similartype(datatype(T), complex(eltype(T)))) + return set_datatype(T, similartype(datatype(T), complex(eltype(T)))) end # TODO: Define a generic `dense` for `Tensor`, `TensorStorage`. dense(storagetype::Type{<:Dense}) = storagetype # TODO: make these more general, move to tensorstorage.jl -datatype(storetype::Type{<:Dense{<:Any,DataT}}) where {DataT} = DataT +datatype(storetype::Type{<:Dense{<:Any, DataT}}) where {DataT} = DataT -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type function promote_rule( - ::Type{<:Dense{ElT1,DataT1}}, ::Type{<:Dense{ElT2,DataT2}} -) where {ElT1,DataT1,ElT2,DataT2} - ElR = promote_type(ElT1, ElT2) - VecR = promote_type(unwrap_array_type(DataT1), unwrap_array_type(DataT2)) - VecR = similartype(VecR, ElR) - return Dense{ElR,VecR} + ::Type{<:Dense{ElT1, DataT1}}, ::Type{<:Dense{ElT2, DataT2}} + ) where {ElT1, DataT1, ElT2, DataT2} + ElR = promote_type(ElT1, ElT2) + VecR = promote_type(unwrap_array_type(DataT1), unwrap_array_type(DataT2)) + VecR = similartype(VecR, ElR) + return Dense{ElR, VecR} end # This is for type promotion for Scalar*Dense function promote_rule( - ::Type{<:Dense{ElT1,DataT}}, ::Type{ElT2} -) where {DataT,ElT1,ElT2<:Number} - ElR = promote_type(ElT1, ElT2) - DataR = set_eltype(DataT, ElR) - return Dense{ElR,DataR} + ::Type{<:Dense{ElT1, DataT}}, ::Type{ElT2} + ) where {DataT, ElT1, ElT2 <: Number} + ElR = promote_type(ElT1, ElT2) + DataR = set_eltype(DataT, ElR) + return Dense{ElR, DataR} end -function convert(::Type{<:Dense{ElR,DataT}}, D::Dense) where {ElR,DataT} - return Dense(convert(DataT, data(D))) +function convert(::Type{<:Dense{ElR, DataT}}, D::Dense) where {ElR, DataT} + return Dense(convert(DataT, data(D))) end diff --git a/NDTensors/src/dense/generic_array_constructors.jl b/NDTensors/src/dense/generic_array_constructors.jl index c9ef550e63..a34a97dcba 100644 --- a/NDTensors/src/dense/generic_array_constructors.jl +++ b/NDTensors/src/dense/generic_array_constructors.jl @@ -1,36 +1,36 @@ -using TypeParameterAccessors: - default_type_parameters, - parenttype, - set_eltype, - specify_default_type_parameters, - specify_type_parameters, - type_parameters +using .Vendored.TypeParameterAccessors: + default_type_parameters, + parenttype, + set_eltype, + specify_default_type_parameters, + specify_type_parameters, + type_parameters ##TODO replace randn in ITensors with generic_randn ## and replace zeros with generic_zeros # This is a file to write generic fills for NDTensors. # This includes random fills, zeros, ... -function generic_randn(StoreT::Type{<:Dense}, dims::Integer; rng=Random.default_rng()) - StoreT = specify_default_type_parameters(StoreT) - DataT = specify_type_parameters( - type_parameters(StoreT, parenttype), eltype, eltype(StoreT) - ) - @assert eltype(StoreT) == eltype(DataT) +function generic_randn(StoreT::Type{<:Dense}, dims::Integer; rng = Random.default_rng()) + StoreT = specify_default_type_parameters(StoreT) + DataT = specify_type_parameters( + type_parameters(StoreT, parenttype), eltype, eltype(StoreT) + ) + @assert eltype(StoreT) == eltype(DataT) - data = generic_randn(DataT, dims; rng=rng) - StoreT = set_datatype(StoreT, typeof(data)) - return StoreT(data) + data = generic_randn(DataT, dims; rng = rng) + StoreT = set_datatype(StoreT, typeof(data)) + return StoreT(data) end function generic_zeros(StoreT::Type{<:Dense}, dims::Integer) - StoreT = specify_default_type_parameters(StoreT) - DataT = specify_type_parameters( - type_parameters(StoreT, parenttype), eltype, eltype(StoreT) - ) - @assert eltype(StoreT) == eltype(DataT) + StoreT = specify_default_type_parameters(StoreT) + DataT = specify_type_parameters( + type_parameters(StoreT, parenttype), eltype, eltype(StoreT) + ) + @assert eltype(StoreT) == eltype(DataT) - data = generic_zeros(DataT, dims) - StoreT = set_datatype(StoreT, typeof(data)) - return StoreT(data) + data = generic_zeros(DataT, dims) + StoreT = set_datatype(StoreT, typeof(data)) + return StoreT(data) end diff --git a/NDTensors/src/dense/set_types.jl b/NDTensors/src/dense/set_types.jl index e4a4fa4ea1..8ac590f003 100644 --- a/NDTensors/src/dense/set_types.jl +++ b/NDTensors/src/dense/set_types.jl @@ -1,13 +1,13 @@ -using TypeParameterAccessors: TypeParameterAccessors, Position, parenttype +using .Vendored.TypeParameterAccessors: TypeParameterAccessors, Position, parenttype function set_datatype(storagetype::Type{<:Dense}, datatype::Type{<:AbstractVector}) - return Dense{eltype(datatype),datatype} + return Dense{eltype(datatype), datatype} end function set_datatype(storagetype::Type{<:Dense}, datatype::Type{<:AbstractArray}) - return error( - "Setting the `datatype` of the storage type `$storagetype` to a $(ndims(datatype))-dimsional array of type `$datatype` is not currently supported, use an `AbstractVector` instead.", - ) + return error( + "Setting the `datatype` of the storage type `$storagetype` to a $(ndims(datatype))-dimsional array of type `$datatype` is not currently supported, use an `AbstractVector` instead.", + ) end TypeParameterAccessors.default_type_parameters(::Type{<:Dense}) = (Float64, Vector) diff --git a/NDTensors/src/diag/diagtensor.jl b/NDTensors/src/diag/diagtensor.jl index ef114eb3c4..de9db5b3e9 100644 --- a/NDTensors/src/diag/diagtensor.jl +++ b/NDTensors/src/diag/diagtensor.jl @@ -1,48 +1,48 @@ -const DiagTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:Diag} -const NonuniformDiagTensor{ElT,N,StoreT,IndsT} = - Tensor{ElT,N,StoreT,IndsT} where {StoreT<:NonuniformDiag} -const UniformDiagTensor{ElT,N,StoreT,IndsT} = - Tensor{ElT,N,StoreT,IndsT} where {StoreT<:UniformDiag} +const DiagTensor{ElT, N, StoreT, IndsT} = Tensor{ElT, N, StoreT, IndsT} where {StoreT <: Diag} +const NonuniformDiagTensor{ElT, N, StoreT, IndsT} = + Tensor{ElT, N, StoreT, IndsT} where {StoreT <: NonuniformDiag} +const UniformDiagTensor{ElT, N, StoreT, IndsT} = + Tensor{ElT, N, StoreT, IndsT} where {StoreT <: UniformDiag} function diag(tensor::DiagTensor) - tensor_diag = NDTensors.similar(dense(typeof(tensor)), (diaglength(tensor),)) - # TODO: Define `eachdiagindex`. - diagview(tensor_diag) .= diagview(tensor) - return tensor_diag + tensor_diag = NDTensors.similar(dense(typeof(tensor)), (diaglength(tensor),)) + # TODO: Define `eachdiagindex`. + diagview(tensor_diag) .= diagview(tensor) + return tensor_diag end IndexStyle(::Type{<:DiagTensor}) = IndexCartesian() # TODO: this needs to be better (promote element type, check order compatibility, # etc. -function convert(::Type{<:DenseTensor{ElT,N}}, T::DiagTensor{ElT,N}) where {ElT<:Number,N} - return dense(T) +function convert(::Type{<:DenseTensor{ElT, N}}, T::DiagTensor{ElT, N}) where {ElT <: Number, N} + return dense(T) end -convert(::Type{Diagonal}, D::DiagTensor{<:Number,2}) = Diagonal(data(D)) +convert(::Type{Diagonal}, D::DiagTensor{<:Number, 2}) = Diagonal(data(D)) -function Array{ElT,N}(T::DiagTensor{ElT,N}) where {ElT,N} - return array(T) +function Array{ElT, N}(T::DiagTensor{ElT, N}) where {ElT, N} + return array(T) end -function Array(T::DiagTensor{ElT,N}) where {ElT,N} - return Array{ElT,N}(T) +function Array(T::DiagTensor{ElT, N}) where {ElT, N} + return Array{ElT, N}(T) end function diagview(T::NonuniformDiagTensor) - return data(T) + return data(T) end function zeros(tensortype::Type{<:DiagTensor}, inds) - return tensor(generic_zeros(storagetype(tensortype), mindim(inds)), inds) + return tensor(generic_zeros(storagetype(tensortype), mindim(inds)), inds) end function zeros(tensortype::Type{<:DiagTensor}, inds::Dims) - return tensor(generic_zeros(storagetype(tensortype), mindim(inds)), inds) + return tensor(generic_zeros(storagetype(tensortype), mindim(inds)), inds) end function zeros(tensortype::Type{<:DiagTensor}, inds::Tuple{}) - return tensor(generic_zeros(storagetype(tensortype), mindim(inds)), inds) + return tensor(generic_zeros(storagetype(tensortype), mindim(inds)), inds) end # Compute the norm of Uniform diagonal tensor @@ -71,151 +71,151 @@ Set the entire diagonal of a uniform DiagTensor. setdiag(T::UniformDiagTensor, val) = tensor(Diag(val), inds(T)) function Base.copyto!(R::DenseTensor, T::DiagTensor) - diagview(R) .= diagview(T) - return R + diagview(R) .= diagview(T) + return R end @propagate_inbounds function getindex( - T::DiagTensor{ElT,N}, inds::Vararg{Int,N} -) where {ElT,N} - if all(==(inds[1]), inds) - return getdiagindex(T, inds[1]) - else - return zero(eltype(ElT)) - end -end -@propagate_inbounds getindex(T::DiagTensor{<:Number,1}, ind::Int) = storage(T)[ind] + T::DiagTensor{ElT, N}, inds::Vararg{Int, N} + ) where {ElT, N} + if all(==(inds[1]), inds) + return getdiagindex(T, inds[1]) + else + return zero(eltype(ElT)) + end +end +@propagate_inbounds getindex(T::DiagTensor{<:Number, 1}, ind::Int) = storage(T)[ind] using .Expose: expose -@propagate_inbounds getindex(T::DiagTensor{<:Number,0}) = getindex(expose(storage(T))) +@propagate_inbounds getindex(T::DiagTensor{<:Number, 0}) = getindex(expose(storage(T))) # Set diagonal elements # Throw error for off-diagonal @propagate_inbounds function setindex!( - T::DiagTensor{<:Number,N}, val, inds::Vararg{Int,N} -) where {N} - all(==(inds[1]), inds) || error("Cannot set off-diagonal element of Diag storage") - setdiagindex!(T, val, inds[1]) - return T + T::DiagTensor{<:Number, N}, val, inds::Vararg{Int, N} + ) where {N} + all(==(inds[1]), inds) || error("Cannot set off-diagonal element of Diag storage") + setdiagindex!(T, val, inds[1]) + return T end -@propagate_inbounds function setindex!(T::DiagTensor{<:Number,1}, val, ind::Int) - return (storage(T)[ind] = val) +@propagate_inbounds function setindex!(T::DiagTensor{<:Number, 1}, val, ind::Int) + return (storage(T)[ind] = val) end -@propagate_inbounds setindex!(T::DiagTensor{<:Number,0}, val) = (storage(T)[1] = val) +@propagate_inbounds setindex!(T::DiagTensor{<:Number, 0}, val) = (storage(T)[1] = val) -function setindex!(T::UniformDiagTensor{<:Number,N}, val, inds::Vararg{Int,N}) where {N} - return error("Cannot set elements of a uniform Diag storage") +function setindex!(T::UniformDiagTensor{<:Number, N}, val, inds::Vararg{Int, N}) where {N} + return error("Cannot set elements of a uniform Diag storage") end # TODO: make a fill!! that works for uniform and non-uniform #fill!(T::DiagTensor,v) = fill!(storage(T),v) -function dense(::Type{<:Tensor{ElT,N,StoreT,IndsT}}) where {ElT,N,StoreT<:Diag,IndsT} - return Tensor{ElT,N,dense(StoreT),IndsT} +function dense(::Type{<:Tensor{ElT, N, StoreT, IndsT}}) where {ElT, N, StoreT <: Diag, IndsT} + return Tensor{ElT, N, dense(StoreT), IndsT} end -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type # convert to Dense function dense(T::DiagTensor) - R = zeros(dense(typeof(T)), inds(T)) - diagview(R) .= diagview(T) - return R + R = zeros(dense(typeof(T)), inds(T)) + diagview(R) .= diagview(T) + return R end denseblocks(T::DiagTensor) = dense(T) function permutedims!( - R::DiagTensor{<:Number,N}, - T::DiagTensor{<:Number,N}, - perm::NTuple{N,Int}, - f::Function=(r, t) -> t, -) where {N} - # TODO: check that inds(R)==permute(inds(T),perm)? - diagview(R) .= f.(diagview(R), diagview(T)) - return R + R::DiagTensor{<:Number, N}, + T::DiagTensor{<:Number, N}, + perm::NTuple{N, Int}, + f::Function = (r, t) -> t, + ) where {N} + # TODO: check that inds(R)==permute(inds(T),perm)? + diagview(R) .= f.(diagview(R), diagview(T)) + return R end function permutedims( - T::DiagTensor{<:Number,N}, perm::NTuple{N,Int}, f::Function=identity -) where {N} - R = NDTensors.similar(T) - g(r, t) = f(t) - permutedims!(R, T, perm, g) - return R + T::DiagTensor{<:Number, N}, perm::NTuple{N, Int}, f::Function = identity + ) where {N} + R = NDTensors.similar(T) + g(r, t) = f(t) + permutedims!(R, T, perm, g) + return R end function permutedims( - T::UniformDiagTensor{<:Number,N}, perm::NTuple{N,Int}, f::Function=identity -) where {N} - R = tensor(Diag(f(getdiagindex(T, 1))), permute(inds(T), perm)) - return R + T::UniformDiagTensor{<:Number, N}, perm::NTuple{N, Int}, f::Function = identity + ) where {N} + R = tensor(Diag(f(getdiagindex(T, 1))), permute(inds(T), perm)) + return R end # Version that may overwrite in-place or may return the result function permutedims!!( - R::NonuniformDiagTensor{<:Number,N}, - T::NonuniformDiagTensor{<:Number,N}, - perm::NTuple{N,Int}, - f::Function=(r, t) -> t, -) where {N} - R = convert(promote_type(typeof(R), typeof(T)), R) - permutedims!(R, T, perm, f) - return R + R::NonuniformDiagTensor{<:Number, N}, + T::NonuniformDiagTensor{<:Number, N}, + perm::NTuple{N, Int}, + f::Function = (r, t) -> t, + ) where {N} + R = convert(promote_type(typeof(R), typeof(T)), R) + permutedims!(R, T, perm, f) + return R end function permutedims!!( - R::UniformDiagTensor{ElR,N}, - T::UniformDiagTensor{ElT,N}, - perm::NTuple{N,Int}, - f::Function=(r, t) -> t, -) where {ElR,ElT,N} - R = convert(promote_type(typeof(R), typeof(T)), R) - R = tensor(Diag(f(getdiagindex(R, 1), getdiagindex(T, 1))), inds(R)) - return R + R::UniformDiagTensor{ElR, N}, + T::UniformDiagTensor{ElT, N}, + perm::NTuple{N, Int}, + f::Function = (r, t) -> t, + ) where {ElR, ElT, N} + R = convert(promote_type(typeof(R), typeof(T)), R) + R = tensor(Diag(f(getdiagindex(R, 1), getdiagindex(T, 1))), inds(R)) + return R end function permutedims!( - R::DenseTensor{ElR,N}, T::DiagTensor{ElT,N}, perm::NTuple{N,Int}, f::Function=(r, t) -> t -) where {ElR,ElT,N} - diagview(R) .= f.(diagview(R), diagview(T)) - return R + R::DenseTensor{ElR, N}, T::DiagTensor{ElT, N}, perm::NTuple{N, Int}, f::Function = (r, t) -> t + ) where {ElR, ElT, N} + diagview(R) .= f.(diagview(R), diagview(T)) + return R end function permutedims!!( - R::DenseTensor{ElR,N}, T::DiagTensor{ElT,N}, perm::NTuple{N,Int}, f::Function=(r, t) -> t -) where {ElR,ElT,N} - RR = convert(promote_type(typeof(R), typeof(T)), R) - permutedims!(RR, T, perm, f) - return RR + R::DenseTensor{ElR, N}, T::DiagTensor{ElT, N}, perm::NTuple{N, Int}, f::Function = (r, t) -> t + ) where {ElR, ElT, N} + RR = convert(promote_type(typeof(R), typeof(T)), R) + permutedims!(RR, T, perm, f) + return RR end # TODO: make a single implementation since this is # the same as the version with the input types # swapped. function permutedims!!( - R::DiagTensor{ElR,N}, T::DenseTensor{ElT,N}, perm::NTuple{N,Int}, f::Function=(r, t) -> t -) where {ElR,ElT,N} - RR = convert(promote_type(typeof(R), typeof(T)), R) - permutedims!(RR, T, perm, f) - return RR + R::DiagTensor{ElR, N}, T::DenseTensor{ElT, N}, perm::NTuple{N, Int}, f::Function = (r, t) -> t + ) where {ElR, ElT, N} + RR = convert(promote_type(typeof(R), typeof(T)), R) + permutedims!(RR, T, perm, f) + return RR end function Base.mapreduce(f, op, t1::DiagTensor, t_tail::DiagTensor...; kwargs...) - elt = mapreduce(eltype, promote_type, (t1, t_tail...)) - if !iszero(f(zero(elt))) - return mapreduce(f, op, array(t1), array.(t_tail)...; kwargs...) - end - if length(t1) > diaglength(t1) - # Some elements are zero, account for that - # with the initial value. - init_kwargs = (; init=zero(elt)) - else - init_kwargs = (;) - end - return mapreduce(f, op, diagview(t1), diagview.(t_tail)...; kwargs..., init_kwargs...) + elt = mapreduce(eltype, promote_type, (t1, t_tail...)) + if !iszero(f(zero(elt))) + return mapreduce(f, op, array(t1), array.(t_tail)...; kwargs...) + end + if length(t1) > diaglength(t1) + # Some elements are zero, account for that + # with the initial value. + init_kwargs = (; init = zero(elt)) + else + init_kwargs = (;) + end + return mapreduce(f, op, diagview(t1), diagview.(t_tail)...; kwargs..., init_kwargs...) end function Base.show(io::IO, mime::MIME"text/plain", T::DiagTensor) - summary(io, T) - print_tensor(io, T) - return nothing + summary(io, T) + print_tensor(io, T) + return nothing end diff --git a/NDTensors/src/diag/set_types.jl b/NDTensors/src/diag/set_types.jl index bbb798a5a4..4527136471 100644 --- a/NDTensors/src/diag/set_types.jl +++ b/NDTensors/src/diag/set_types.jl @@ -1,20 +1,20 @@ -using TypeParameterAccessors: TypeParameterAccessors +using .Vendored.TypeParameterAccessors: TypeParameterAccessors function TypeParameterAccessors.set_eltype(storagetype::Type{<:UniformDiag}, eltype::Type) - return Diag{eltype,eltype} + return Diag{eltype, eltype} end function TypeParameterAccessors.set_eltype( - storagetype::Type{<:NonuniformDiag}, eltype::Type{<:AbstractArray} -) - return Diag{eltype,similartype(storagetype, eltype)} + storagetype::Type{<:NonuniformDiag}, eltype::Type{<:AbstractArray} + ) + return Diag{eltype, similartype(storagetype, eltype)} end # TODO: Remove this once uniform diagonal tensors use FillArrays for the data. function set_datatype(storagetype::Type{<:UniformDiag}, datatype::Type) - return Diag{datatype,datatype} + return Diag{datatype, datatype} end function set_datatype(storagetype::Type{<:NonuniformDiag}, datatype::Type{<:AbstractArray}) - return Diag{eltype(datatype),datatype} + return Diag{eltype(datatype), datatype} end diff --git a/NDTensors/src/diag/similar.jl b/NDTensors/src/diag/similar.jl index 7231caffe2..9dcec8dd3c 100644 --- a/NDTensors/src/diag/similar.jl +++ b/NDTensors/src/diag/similar.jl @@ -1,22 +1,22 @@ -using TypeParameterAccessors: TypeParameterAccessors +using .Vendored.TypeParameterAccessors: TypeParameterAccessors # NDTensors.similar function similar(storagetype::Type{<:Diag}, dims::Dims) - return setdata(storagetype, similar(datatype(storagetype), mindim(dims))) + return setdata(storagetype, similar(datatype(storagetype), mindim(dims))) end # TODO: Redesign UniformDiag to make it handled better # by generic code. function TypeParameterAccessors.similartype(storagetype::Type{<:UniformDiag}, eltype::Type) - # This will also set the `datatype`. - return set_eltype(storagetype, eltype) + # This will also set the `datatype`. + return set_eltype(storagetype, eltype) end # Needed to get slice of DiagTensor like T[1:3,1:3] function similar( - T::DiagTensor{<:Number,N}, ::Type{ElR}, inds::Dims{N} -) where {ElR<:Number,N} - return tensor(similar(storage(T), ElR, minimum(inds)), inds) + T::DiagTensor{<:Number, N}, ::Type{ElR}, inds::Dims{N} + ) where {ElR <: Number, N} + return tensor(similar(storage(T), ElR, minimum(inds)), inds) end similar(storage::NonuniformDiag) = setdata(storage, similar(data(storage))) diff --git a/NDTensors/src/dims.jl b/NDTensors/src/dims.jl index 38db99f6d9..14948841a2 100644 --- a/NDTensors/src/dims.jl +++ b/NDTensors/src/dims.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: TypeParameterAccessors +using .Vendored.TypeParameterAccessors: TypeParameterAccessors export dense, dims, dim, mindim, diaglength @@ -18,7 +18,7 @@ dim(::Tuple{}) = 1 dense(ds::Dims) = ds -dense(::Type{DimsT}) where {DimsT<:Dims} = DimsT +dense(::Type{DimsT}) where {DimsT <: Dims} = DimsT dim(ds::Dims) = prod(ds) @@ -72,7 +72,7 @@ sim(i::Int) = i # More complicated definition makes Order(Ref(2)[]) faster @eval struct Order{N} - (OrderT::Type{<:Order})() = $(Expr(:new, :OrderT)) + (OrderT::Type{<:Order})() = $(Expr(:new, :OrderT)) end @doc """ diff --git a/NDTensors/src/empty/empty.jl b/NDTensors/src/empty/empty.jl index 78aaea0e91..37ca1dcda2 100644 --- a/NDTensors/src/empty/empty.jl +++ b/NDTensors/src/empty/empty.jl @@ -1,5 +1,5 @@ using SparseArrays: SparseArrays -using TypeParameterAccessors: TypeParameterAccessors, set_eltype, similartype +using .Vendored.TypeParameterAccessors: TypeParameterAccessors, set_eltype, similartype # # Represents a tensor order that could be set to any order. @@ -95,6 +95,6 @@ function Base.show(io::IO, mime::MIME"text/plain", S::EmptyStorage) return println(io, typeof(S)) end -using TypeParameterAccessors: TypeParameterAccessors +using .Vendored.TypeParameterAccessors: TypeParameterAccessors TypeParameterAccessors.parenttype(empty::Type{<:EmptyStorage}) = storagetype(empty) zero(empty::EmptyStorage) = empty diff --git a/NDTensors/src/lib/AMDGPUExtensions/src/AMDGPUExtensions.jl b/NDTensors/src/lib/AMDGPUExtensions/src/AMDGPUExtensions.jl index e9e77e6cc5..c7ec035a90 100644 --- a/NDTensors/src/lib/AMDGPUExtensions/src/AMDGPUExtensions.jl +++ b/NDTensors/src/lib/AMDGPUExtensions/src/AMDGPUExtensions.jl @@ -1,4 +1,12 @@ module AMDGPUExtensions + +module Vendored + include(joinpath( + "..", "..", "..", "vendored", "TypeParameterAccessors", "src", + "TypeParameterAccessors.jl", + )) +end + include("roc.jl") end diff --git a/NDTensors/src/lib/AMDGPUExtensions/src/roc.jl b/NDTensors/src/lib/AMDGPUExtensions/src/roc.jl index 03c28c613e..0aa1a58693 100644 --- a/NDTensors/src/lib/AMDGPUExtensions/src/roc.jl +++ b/NDTensors/src/lib/AMDGPUExtensions/src/roc.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: TypeParameterAccessors, Position +using .Vendored.TypeParameterAccessors: TypeParameterAccessors, Position using ..GPUArraysCoreExtensions: storagemode # Implemented in NDTensorsAMDGPUExt function roc end @@ -10,5 +10,5 @@ function roc end struct ROCArrayAdaptor{B} end function TypeParameterAccessors.position(::Type{<:ROCArrayAdaptor}, ::typeof(storagemode)) - return Position(1) + return Position(1) end diff --git a/NDTensors/src/lib/CUDAExtensions/src/CUDAExtensions.jl b/NDTensors/src/lib/CUDAExtensions/src/CUDAExtensions.jl index 63a030f9cb..1e101699ad 100644 --- a/NDTensors/src/lib/CUDAExtensions/src/CUDAExtensions.jl +++ b/NDTensors/src/lib/CUDAExtensions/src/CUDAExtensions.jl @@ -1,4 +1,12 @@ module CUDAExtensions + +module Vendored + include(joinpath( + "..", "..", "..", "vendored", "TypeParameterAccessors", "src", + "TypeParameterAccessors.jl", + )) +end + include("cuda.jl") end diff --git a/NDTensors/src/lib/CUDAExtensions/src/cuda.jl b/NDTensors/src/lib/CUDAExtensions/src/cuda.jl index 948f0b4b17..60dad8d5d3 100644 --- a/NDTensors/src/lib/CUDAExtensions/src/cuda.jl +++ b/NDTensors/src/lib/CUDAExtensions/src/cuda.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: TypeParameterAccessors, Position +using .Vendored.TypeParameterAccessors: TypeParameterAccessors, Position using ..GPUArraysCoreExtensions: storagemode # Implemented in NDTensorsCUDAExt function cu end @@ -10,5 +10,5 @@ function cu end struct CuArrayAdaptor{B} end function TypeParameterAccessors.position(::Type{<:CuArrayAdaptor}, ::typeof(storagemode)) - return Position(1) + return Position(1) end diff --git a/NDTensors/src/lib/Expose/src/Expose.jl b/NDTensors/src/lib/Expose/src/Expose.jl index fb64927d7b..e0c38b56b8 100644 --- a/NDTensors/src/lib/Expose/src/Expose.jl +++ b/NDTensors/src/lib/Expose/src/Expose.jl @@ -5,6 +5,13 @@ using Base: ReshapedArray using StridedViews using Adapt: Adapt, adapt, adapt_structure +module Vendored + include(joinpath( + "..", "..", "..", "vendored", "TypeParameterAccessors", "src", + "TypeParameterAccessors.jl", + )) +end + include("exposed.jl") include("import.jl") diff --git a/NDTensors/src/lib/Expose/src/exposed.jl b/NDTensors/src/lib/Expose/src/exposed.jl index b57f19e34c..3ebc88e416 100644 --- a/NDTensors/src/lib/Expose/src/exposed.jl +++ b/NDTensors/src/lib/Expose/src/exposed.jl @@ -1,17 +1,17 @@ -using TypeParameterAccessors: - TypeParameterAccessors, unwrap_array_type, parenttype, type_parameters -struct Exposed{Unwrapped,Object} - object::Object +using .Vendored.TypeParameterAccessors: + TypeParameterAccessors, unwrap_array_type, parenttype, type_parameters +struct Exposed{Unwrapped, Object} + object::Object end -expose(object) = Exposed{unwrap_array_type(object),typeof(object)}(object) +expose(object) = Exposed{unwrap_array_type(object), typeof(object)}(object) unexpose(E::Exposed) = E.object ## TODO remove TypeParameterAccessors when SetParameters is removed TypeParameterAccessors.parenttype(type::Type{<:Exposed}) = type_parameters(type, parenttype) function TypeParameterAccessors.position(::Type{<:Exposed}, ::typeof(parenttype)) - return TypeParameterAccessors.Position(1) + return TypeParameterAccessors.Position(1) end TypeParameterAccessors.unwrap_array_type(type::Type{<:Exposed}) = parenttype(type) TypeParameterAccessors.unwrap_array_type(E::Exposed) = unwrap_array_type(typeof(E)) diff --git a/NDTensors/src/lib/GPUArraysCoreExtensions/src/GPUArraysCoreExtensions.jl b/NDTensors/src/lib/GPUArraysCoreExtensions/src/GPUArraysCoreExtensions.jl index 6845816fca..059718e5bb 100644 --- a/NDTensors/src/lib/GPUArraysCoreExtensions/src/GPUArraysCoreExtensions.jl +++ b/NDTensors/src/lib/GPUArraysCoreExtensions/src/GPUArraysCoreExtensions.jl @@ -1,4 +1,12 @@ module GPUArraysCoreExtensions + +module Vendored + include(joinpath( + "..", "..", "..", "vendored", "TypeParameterAccessors", "src", + "TypeParameterAccessors.jl", + )) +end + include("gpuarrayscore.jl") end diff --git a/NDTensors/src/lib/GPUArraysCoreExtensions/src/gpuarrayscore.jl b/NDTensors/src/lib/GPUArraysCoreExtensions/src/gpuarrayscore.jl index 0042a05dab..1e37904042 100644 --- a/NDTensors/src/lib/GPUArraysCoreExtensions/src/gpuarrayscore.jl +++ b/NDTensors/src/lib/GPUArraysCoreExtensions/src/gpuarrayscore.jl @@ -1,15 +1,15 @@ using ..Expose: Exposed, unexpose -using TypeParameterAccessors: TypeParameterAccessors, type_parameters, set_type_parameters +using .Vendored.TypeParameterAccessors: TypeParameterAccessors, type_parameters, set_type_parameters function storagemode(object) - return storagemode(typeof(object)) + return storagemode(typeof(object)) end function storagemode(type::Type) - return type_parameters(type, storagemode) + return type_parameters(type, storagemode) end function set_storagemode(type::Type, param) - return set_type_parameters(type, storagemode, param) + return set_type_parameters(type, storagemode, param) end function cpu end diff --git a/NDTensors/src/lib/MetalExtensions/src/MetalExtensions.jl b/NDTensors/src/lib/MetalExtensions/src/MetalExtensions.jl index 504f13f7aa..d50f909aa0 100644 --- a/NDTensors/src/lib/MetalExtensions/src/MetalExtensions.jl +++ b/NDTensors/src/lib/MetalExtensions/src/MetalExtensions.jl @@ -1,4 +1,12 @@ module MetalExtensions + +module Vendored + include(joinpath( + "..", "..", "..", "vendored", "TypeParameterAccessors", "src", + "TypeParameterAccessors.jl", + )) +end + include("metal.jl") end diff --git a/NDTensors/src/lib/MetalExtensions/src/metal.jl b/NDTensors/src/lib/MetalExtensions/src/metal.jl index ab329d3249..3253c7ed66 100644 --- a/NDTensors/src/lib/MetalExtensions/src/metal.jl +++ b/NDTensors/src/lib/MetalExtensions/src/metal.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: TypeParameterAccessors, Position +using .Vendored.TypeParameterAccessors: TypeParameterAccessors, Position using ..GPUArraysCoreExtensions: storagemode # Implemented in NDTensorsMetalExt function mtl end @@ -11,5 +11,5 @@ function mtl end struct MtlArrayAdaptor{B} end function TypeParameterAccessors.position(::Type{<:MtlArrayAdaptor}, ::typeof(storagemode)) - return Position(1) + return Position(1) end diff --git a/NDTensors/src/lib/RankFactorization/src/RankFactorization.jl b/NDTensors/src/lib/RankFactorization/src/RankFactorization.jl index a2b950ac93..321a65112e 100644 --- a/NDTensors/src/lib/RankFactorization/src/RankFactorization.jl +++ b/NDTensors/src/lib/RankFactorization/src/RankFactorization.jl @@ -1,5 +1,14 @@ module RankFactorization + +module Vendored + include(joinpath( + "..", "..", "..", "vendored", "TypeParameterAccessors", "src", + "TypeParameterAccessors.jl", + )) +end + include("default_kwargs.jl") include("truncate_spectrum.jl") include("spectrum.jl") + end diff --git a/NDTensors/src/lib/RankFactorization/src/default_kwargs.jl b/NDTensors/src/lib/RankFactorization/src/default_kwargs.jl index fda4fe9786..9cb170b0e5 100644 --- a/NDTensors/src/lib/RankFactorization/src/default_kwargs.jl +++ b/NDTensors/src/lib/RankFactorization/src/default_kwargs.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type replace_nothing(::Nothing, replacement) = replacement replace_nothing(value, replacement) = value diff --git a/NDTensors/src/lib/RankFactorization/src/truncate_spectrum.jl b/NDTensors/src/lib/RankFactorization/src/truncate_spectrum.jl index d980623364..d05c2b2ebc 100644 --- a/NDTensors/src/lib/RankFactorization/src/truncate_spectrum.jl +++ b/NDTensors/src/lib/RankFactorization/src/truncate_spectrum.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type ## TODO write Exposed version of truncate function truncate!!(P::AbstractArray; kwargs...) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index 64135c6f72..bc128940e0 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -12,71 +12,71 @@ import .Expose: qr_positive, ql, ql_positive # ```julia # contract(T1, (1, -1), T2, (-1, 2)) # ``` -function Base.:*(T1::Tensor{<:Any,2,<:Dense}, T2::Tensor{<:Any,2,<:Dense}) - RM = matrix(T1) * matrix(T2) - indsR = (ind(T1, 1), ind(T2, 2)) - return tensor(Dense(vec(RM)), indsR) +function Base.:*(T1::Tensor{<:Any, 2, <:Dense}, T2::Tensor{<:Any, 2, <:Dense}) + RM = matrix(T1) * matrix(T2) + indsR = (ind(T1, 1), ind(T2, 2)) + return tensor(Dense(vec(RM)), indsR) end function LinearAlgebra.dot(x::Tensor, y::Tensor) - size(x) == size(y) || throw( - DimensionMismatch( - "dimensions must match in `dot(x::Tensor, y::Tensor)`: `x` has size `$(size(x))` while `y` has size `$(size(y))`.", - ), - ) - labels = ntuple(dim -> -dim, ndims(x)) - return contract(conj(x), labels, y, labels)[] + size(x) == size(y) || throw( + DimensionMismatch( + "dimensions must match in `dot(x::Tensor, y::Tensor)`: `x` has size `$(size(x))` while `y` has size `$(size(y))`.", + ), + ) + labels = ntuple(dim -> -dim, ndims(x)) + return contract(conj(x), labels, y, labels)[] end -function LinearAlgebra.exp(T::DenseTensor{ElT,2}) where {ElT<:Union{Real,Complex}} - expTM = exp(matrix(T)) - return tensor(Dense(vec(expTM)), inds(T)) +function LinearAlgebra.exp(T::DenseTensor{ElT, 2}) where {ElT <: Union{Real, Complex}} + expTM = exp(matrix(T)) + return tensor(Dense(vec(expTM)), inds(T)) end function LinearAlgebra.exp( - T::Hermitian{ElT,<:DenseTensor{ElT,2}} -) where {ElT<:Union{Real,Complex}} - # exp(::Hermitian/Symmetric) returns Hermitian/Symmetric, - # so extract the parent matrix - expTM = parent(exp(matrix(T))) - return tensor(Dense(vec(expTM)), inds(T)) + T::Hermitian{ElT, <:DenseTensor{ElT, 2}} + ) where {ElT <: Union{Real, Complex}} + # exp(::Hermitian/Symmetric) returns Hermitian/Symmetric, + # so extract the parent matrix + expTM = parent(exp(matrix(T))) + return tensor(Dense(vec(expTM)), inds(T)) end function svd_catch_error(A; kwargs...) - USV = try - svd(expose(A); kwargs...) - catch - return nothing - end - return USV + USV = try + svd(expose(A); kwargs...) + catch + return nothing + end + return USV end function lapack_svd_error_message(alg) - return "The SVD algorithm `\"$alg\"` has thrown an error,\n" * - "likely because of a convergance failure. You can try\n" * - "other SVD algorithms that may converge better using the\n" * - "`alg` (or `svd_alg` if called through `factorize` or MPS/MPO functionality) keyword argument:\n\n" * - " - \"divide_and_conquer\" is a divide-and-conquer algorithm\n" * - " (LAPACK's `gesdd`). It is fast, but may lead to some innacurate\n" * - " singular values for very ill-conditioned matrices.\n" * - " It also may sometimes fail to converge, leading to errors\n" * - " (in which case `\"qr_iteration\"` or `\"recursive\"` can be tried).\n\n" * - " - `\"qr_iteration\"` (LAPACK's `gesvd`) is typically slower \n" * - " than \"divide_and_conquer\", especially for large matrices,\n" * - " but is more accurate for very ill-conditioned matrices \n" * - " compared to `\"divide_and_conquer\"`.\n\n" * - " - `\"recursive\"` is ITensor's custom SVD algorithm. It is very\n" * - " reliable, but may be slow if high precision is needed.\n" * - " To get an `svd` of a matrix `A`, an eigendecomposition of\n" * - " ``A^{\\dagger} A`` is used to compute `U` and then a `qr` of\n" * - " ``A^{\\dagger} U`` is used to compute `V`. This is performed\n" * - " recursively to compute small singular values.\n" * - " - `\"qr_algorithm\"` is a CUDA.jl implemented SVD algorithm using QR.\n" * - " - `\"jacobi_algorithm\"` is a CUDA.jl implemented SVD algorithm.\n\n" * - "Returning `nothing`. For an output `F = svd(A, ...)` you can check if\n" * - "`isnothing(F)` in your code and try a different algorithm.\n\n" * - "To suppress this message in the future, you can wrap the `svd` call in the\n" * - "`@suppress` macro from the `Suppressor` package.\n" + return "The SVD algorithm `\"$alg\"` has thrown an error,\n" * + "likely because of a convergance failure. You can try\n" * + "other SVD algorithms that may converge better using the\n" * + "`alg` (or `svd_alg` if called through `factorize` or MPS/MPO functionality) keyword argument:\n\n" * + " - \"divide_and_conquer\" is a divide-and-conquer algorithm\n" * + " (LAPACK's `gesdd`). It is fast, but may lead to some innacurate\n" * + " singular values for very ill-conditioned matrices.\n" * + " It also may sometimes fail to converge, leading to errors\n" * + " (in which case `\"qr_iteration\"` or `\"recursive\"` can be tried).\n\n" * + " - `\"qr_iteration\"` (LAPACK's `gesvd`) is typically slower \n" * + " than \"divide_and_conquer\", especially for large matrices,\n" * + " but is more accurate for very ill-conditioned matrices \n" * + " compared to `\"divide_and_conquer\"`.\n\n" * + " - `\"recursive\"` is ITensor's custom SVD algorithm. It is very\n" * + " reliable, but may be slow if high precision is needed.\n" * + " To get an `svd` of a matrix `A`, an eigendecomposition of\n" * + " ``A^{\\dagger} A`` is used to compute `U` and then a `qr` of\n" * + " ``A^{\\dagger} U`` is used to compute `V`. This is performed\n" * + " recursively to compute small singular values.\n" * + " - `\"qr_algorithm\"` is a CUDA.jl implemented SVD algorithm using QR.\n" * + " - `\"jacobi_algorithm\"` is a CUDA.jl implemented SVD algorithm.\n\n" * + "Returning `nothing`. For an output `F = svd(A, ...)` you can check if\n" * + "`isnothing(F)` in your code and try a different algorithm.\n\n" * + "To suppress this message in the future, you can wrap the `svd` call in the\n" * + "`@suppress` macro from the `Suppressor` package.\n" end """ @@ -85,138 +85,138 @@ end svd of an order-2 DenseTensor """ function svd( - T::DenseTensor{ElT,2,IndsT}; - mindim=nothing, - maxdim=nothing, - cutoff=nothing, - use_absolute_cutoff=nothing, - use_relative_cutoff=nothing, - alg=nothing, - # Only used by BlockSparse svd - min_blockdim=nothing, -) where {ElT,IndsT} - alg = replace_nothing(alg, default_svd_alg(T)) - if alg == "divide_and_conquer" - MUSV = svd_catch_error(matrix(T); alg=LinearAlgebra.DivideAndConquer()) - if isnothing(MUSV) - # If "divide_and_conquer" fails, try "qr_iteration" - alg = "qr_iteration" - MUSV = svd_catch_error(matrix(T); alg=LinearAlgebra.QRIteration()) - if isnothing(MUSV) - # If "qr_iteration" fails, try "recursive" - alg = "recursive" + T::DenseTensor{ElT, 2, IndsT}; + mindim = nothing, + maxdim = nothing, + cutoff = nothing, + use_absolute_cutoff = nothing, + use_relative_cutoff = nothing, + alg = nothing, + # Only used by BlockSparse svd + min_blockdim = nothing, + ) where {ElT, IndsT} + alg = replace_nothing(alg, default_svd_alg(T)) + if alg == "divide_and_conquer" + MUSV = svd_catch_error(matrix(T); alg = LinearAlgebra.DivideAndConquer()) + if isnothing(MUSV) + # If "divide_and_conquer" fails, try "qr_iteration" + alg = "qr_iteration" + MUSV = svd_catch_error(matrix(T); alg = LinearAlgebra.QRIteration()) + if isnothing(MUSV) + # If "qr_iteration" fails, try "recursive" + alg = "recursive" + MUSV = svd_recursive(matrix(T)) + end + end + elseif alg == "qr_iteration" + MUSV = svd_catch_error(matrix(T); alg = LinearAlgebra.QRIteration()) + if isnothing(MUSV) + # If "qr_iteration" fails, try "recursive" + alg = "recursive" + MUSV = svd_recursive(matrix(T)) + end + elseif alg == "recursive" MUSV = svd_recursive(matrix(T)) - end + elseif alg == "qr_algorithm" || alg == "jacobi_algorithm" + MUSV = svd_catch_error(matrix(T); alg) + else + error( + "svd algorithm $alg is not currently supported. Please see the documentation for currently supported algorithms.", + ) end - elseif alg == "qr_iteration" - MUSV = svd_catch_error(matrix(T); alg=LinearAlgebra.QRIteration()) if isnothing(MUSV) - # If "qr_iteration" fails, try "recursive" - alg = "recursive" - MUSV = svd_recursive(matrix(T)) + if any(isnan, expose(T)) + println("SVD failed, the matrix you were trying to SVD contains NaNs.") + else + println(lapack_svd_error_message(alg)) + end + return nothing end - elseif alg == "recursive" - MUSV = svd_recursive(matrix(T)) - elseif alg == "qr_algorithm" || alg == "jacobi_algorithm" - MUSV = svd_catch_error(matrix(T); alg) - else - error( - "svd algorithm $alg is not currently supported. Please see the documentation for currently supported algorithms.", - ) - end - if isnothing(MUSV) - if any(isnan, expose(T)) - println("SVD failed, the matrix you were trying to SVD contains NaNs.") + MU, MS, MV = MUSV + conj!(MV) + #end # @timeit_debug + + P = MS .^ 2 + if any(!isnothing, (maxdim, cutoff)) + P, truncerr, _ = truncate!!( + P; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff + ) else - println(lapack_svd_error_message(alg)) + truncerr = 0.0 end - return nothing - end - MU, MS, MV = MUSV - conj!(MV) - #end # @timeit_debug - - P = MS .^ 2 - if any(!isnothing, (maxdim, cutoff)) - P, truncerr, _ = truncate!!( - P; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff - ) - else - truncerr = 0.0 - end - spec = Spectrum(P, truncerr) - dS = length(P) - if dS < length(MS) - MU = expose(MU)[:, 1:dS] - # Fails on some GPU backends like Metal. - # resize!(MS, dS) - MS = MS[1:dS] - MV = expose(MV)[:, 1:dS] - end - - # Make the new indices to go onto U and V - u = eltype(IndsT)(dS) - v = eltype(IndsT)(dS) - Uinds = IndsT((ind(T, 1), u)) - Sinds = IndsT((u, v)) - Vinds = IndsT((ind(T, 2), v)) - U = tensor(Dense(vec(MU)), Uinds) - S = tensor(Diag(MS), Sinds) - V = tensor(Dense(vec(MV)), Vinds) - return U, S, V, spec + spec = Spectrum(P, truncerr) + dS = length(P) + if dS < length(MS) + MU = expose(MU)[:, 1:dS] + # Fails on some GPU backends like Metal. + # resize!(MS, dS) + MS = MS[1:dS] + MV = expose(MV)[:, 1:dS] + end + + # Make the new indices to go onto U and V + u = eltype(IndsT)(dS) + v = eltype(IndsT)(dS) + Uinds = IndsT((ind(T, 1), u)) + Sinds = IndsT((u, v)) + Vinds = IndsT((ind(T, 2), v)) + U = tensor(Dense(vec(MU)), Uinds) + S = tensor(Diag(MS), Sinds) + V = tensor(Dense(vec(MV)), Vinds) + return U, S, V, spec end function LinearAlgebra.eigen( - T::Hermitian{ElT,<:DenseTensor{ElT,2,IndsT}}; - mindim=nothing, - maxdim=nothing, - cutoff=nothing, - use_absolute_cutoff=nothing, - use_relative_cutoff=nothing, -) where {ElT<:Union{Real,Complex},IndsT} - matrixT = matrix(T) - ## TODO Here I am calling parent to ensure that the correct `any` function - ## is envoked for non-cpu matrices - ## TODO use expose here - if any(!isfinite, parent(matrixT)) - throw( - ArgumentError( - "Trying to perform the eigendecomposition of a matrix containing NaNs or Infs" - ), - ) - end - - ### What do we do if DM is full of Nan or Inf? - DM, VM = eigen(expose(matrixT)) - - # Sort by largest to smallest eigenvalues - # TODO: Replace `cpu` with `unwrap_array_type` dispatch. - p = sortperm(cpu(DM); rev=true, by=abs) - DM = DM[p] - VM = VM[:, p] + T::Hermitian{ElT, <:DenseTensor{ElT, 2, IndsT}}; + mindim = nothing, + maxdim = nothing, + cutoff = nothing, + use_absolute_cutoff = nothing, + use_relative_cutoff = nothing, + ) where {ElT <: Union{Real, Complex}, IndsT} + matrixT = matrix(T) + ## TODO Here I am calling parent to ensure that the correct `any` function + ## is envoked for non-cpu matrices + ## TODO use expose here + if any(!isfinite, parent(matrixT)) + throw( + ArgumentError( + "Trying to perform the eigendecomposition of a matrix containing NaNs or Infs" + ), + ) + end - if any(!isnothing, (maxdim, cutoff)) - DM, truncerr, _ = truncate!!( - DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff - ) - dD = length(DM) - if dD < size(VM, 2) - VM = VM[:, 1:dD] + ### What do we do if DM is full of Nan or Inf? + DM, VM = eigen(expose(matrixT)) + + # Sort by largest to smallest eigenvalues + # TODO: Replace `cpu` with `unwrap_array_type` dispatch. + p = sortperm(cpu(DM); rev = true, by = abs) + DM = DM[p] + VM = VM[:, p] + + if any(!isnothing, (maxdim, cutoff)) + DM, truncerr, _ = truncate!!( + DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff + ) + dD = length(DM) + if dD < size(VM, 2) + VM = VM[:, 1:dD] + end + else + dD = length(DM) + truncerr = 0.0 end - else - dD = length(DM) - truncerr = 0.0 - end - spec = Spectrum(DM, truncerr) - - # Make the new indices to go onto V - l = eltype(IndsT)(dD) - r = eltype(IndsT)(dD) - Vinds = IndsT((dag(ind(T, 2)), dag(r))) - Dinds = IndsT((l, dag(r))) - V = tensor(Dense(vec(VM)), Vinds) - D = tensor(Diag(DM), Dinds) - return D, V, spec + spec = Spectrum(DM, truncerr) + + # Make the new indices to go onto V + l = eltype(IndsT)(dD) + r = eltype(IndsT)(dD) + Vinds = IndsT((dag(ind(T, 2)), dag(r))) + Dinds = IndsT((l, dag(r))) + V = tensor(Dense(vec(VM)), Vinds) + D = tensor(Diag(DM), Dinds) + return D, V, spec end """ @@ -232,30 +232,30 @@ Sampling is based on https://arxiv.org/abs/math-ph/0609050 such that in the case `n==m`, the unitary matrix will be sampled according to the Haar measure. """ -function random_unitary(::Type{ElT}, n::Int, m::Int) where {ElT<:Number} - return random_unitary(Random.default_rng(), ElT, n, m) +function random_unitary(::Type{ElT}, n::Int, m::Int) where {ElT <: Number} + return random_unitary(Random.default_rng(), ElT, n, m) end function random_unitary(rng::AbstractRNG, DataT::Type{<:AbstractArray}, n::Int, m::Int) - ElT = eltype(DataT) - if n < m - return DataT(random_unitary(rng, ElT, m, n)') - end - F = qr(randn(rng, ElT, n, m)) - Q = DataT(F.Q) - # The upper triangle of F.factors - # are the elements of R. - # Multiply cols of Q by the signs - # that would make diagonal of R - # non-negative: - for c in 1:size(Q, 2) - Q[:, c] .*= sign(F.factors[c, c]) - end - return Q + ElT = eltype(DataT) + if n < m + return DataT(random_unitary(rng, ElT, m, n)') + end + F = qr(randn(rng, ElT, n, m)) + Q = DataT(F.Q) + # The upper triangle of F.factors + # are the elements of R. + # Multiply cols of Q by the signs + # that would make diagonal of R + # non-negative: + for c in 1:size(Q, 2) + Q[:, c] .*= sign(F.factors[c, c]) + end + return Q end -function random_unitary(rng::AbstractRNG, ::Type{ElT}, n::Int, m::Int) where {ElT<:Number} - return random_unitary(rng, set_ndims(default_datatype(ElT), 2), n, m) +function random_unitary(rng::AbstractRNG, ::Type{ElT}, n::Int, m::Int) where {ElT <: Number} + return random_unitary(rng, set_ndims(default_datatype(ElT), 2), n, m) end random_unitary(n::Int, m::Int) = random_unitary(ComplexF64, n, m) @@ -270,99 +270,99 @@ identity, or if m > n O*transpose(O) is the identity. Optionally can pass a real number type as the first argument to obtain a matrix of that type. """ -random_orthog(::Type{ElT}, n::Int, m::Int) where {ElT<:Real} = random_unitary(ElT, n, m) +random_orthog(::Type{ElT}, n::Int, m::Int) where {ElT <: Real} = random_unitary(ElT, n, m) random_orthog(n::Int, m::Int) = random_orthog(Float64, n, m) function LinearAlgebra.eigen( - T::DenseTensor{ElT,2,IndsT}; - mindim=nothing, - maxdim=nothing, - cutoff=nothing, - use_absolute_cutoff=nothing, - use_relative_cutoff=nothing, -) where {ElT<:Union{Real,Complex},IndsT} - matrixT = matrix(T) - if any(!isfinite, matrixT) - throw( - ArgumentError( - "Trying to perform the eigendecomposition of a matrix containing NaNs or Infs" - ), - ) - end - - DM, VM = eigen(expose(matrixT)) - - # Sort by largest to smallest eigenvalues - #p = sortperm(DM; rev = true) - #DM = DM[p] - #VM = VM[:,p] + T::DenseTensor{ElT, 2, IndsT}; + mindim = nothing, + maxdim = nothing, + cutoff = nothing, + use_absolute_cutoff = nothing, + use_relative_cutoff = nothing, + ) where {ElT <: Union{Real, Complex}, IndsT} + matrixT = matrix(T) + if any(!isfinite, matrixT) + throw( + ArgumentError( + "Trying to perform the eigendecomposition of a matrix containing NaNs or Infs" + ), + ) + end - if any(!isnothing, (maxdim, cutoff)) - DM, truncerr, _ = truncate!!( - DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff - ) - dD = length(DM) - if dD < size(VM, 2) - VM = VM[:, 1:dD] + DM, VM = eigen(expose(matrixT)) + + # Sort by largest to smallest eigenvalues + #p = sortperm(DM; rev = true) + #DM = DM[p] + #VM = VM[:,p] + + if any(!isnothing, (maxdim, cutoff)) + DM, truncerr, _ = truncate!!( + DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff + ) + dD = length(DM) + if dD < size(VM, 2) + VM = VM[:, 1:dD] + end + else + dD = length(DM) + truncerr = 0.0 end - else - dD = length(DM) - truncerr = 0.0 - end - spec = Spectrum(abs.(DM), truncerr) - - i1, i2 = inds(T) - - # Make the new indices to go onto D and V - l = typeof(i1)(dD) - r = dag(sim(l)) - Dinds = (l, r) - Vinds = (dag(i2), r) - D = complex(tensor(Diag(DM), Dinds)) - V = complex(tensor(Dense(vec(VM)), Vinds)) - return D, V, spec + spec = Spectrum(abs.(DM), truncerr) + + i1, i2 = inds(T) + + # Make the new indices to go onto D and V + l = typeof(i1)(dD) + r = dag(sim(l)) + Dinds = (l, r) + Vinds = (dag(i2), r) + D = complex(tensor(Diag(DM), Dinds)) + V = complex(tensor(Dense(vec(VM)), Vinds)) + return D, V, spec end # LinearAlgebra.qr -function qr(T::DenseTensor{<:Any,2}; positive=false) - qxf = positive ? qr_positive : qr - return qx(qxf, T) +function qr(T::DenseTensor{<:Any, 2}; positive = false) + qxf = positive ? qr_positive : qr + return qx(qxf, T) end # NDTensors.Expose.ql -function ql(T::DenseTensor{<:Any,2}; positive=false) - qxf = positive ? ql_positive : ql - return qx(qxf, T) +function ql(T::DenseTensor{<:Any, 2}; positive = false) + qxf = positive ? ql_positive : ql + return qx(qxf, T) end # # Generic function for qr and ql decomposition of dense matrix. # The X tensor = R or L. # -function qx(qx::Function, T::DenseTensor{<:Any,2}) - QM, XM = qx(expose(matrix(T))) - # Be aware that if positive==false, then typeof(QM)=LinearAlgebra.QRCompactWYQ, not Matrix - # It gets converted to matrix below. - # Make the new indices to go onto Q and R - q, r = inds(T) - q = dim(q) < dim(r) ? sim(q) : sim(r) - IndsT = indstype(T) #get the index type - Qinds = IndsT((ind(T, 1), q)) - Xinds = IndsT((q, ind(T, 2))) - QM = convert(typeof(XM), QM) - ## Here I convert QM twice because of an issue in CUDA where convert does not take QM to be a UnifiedBuffer array - QM = convert(typeof(XM), QM) - Q = tensor(Dense(vec(QM)), Qinds) #Q was strided - X = tensor(Dense(vec(XM)), Xinds) - return Q, X +function qx(qx::Function, T::DenseTensor{<:Any, 2}) + QM, XM = qx(expose(matrix(T))) + # Be aware that if positive==false, then typeof(QM)=LinearAlgebra.QRCompactWYQ, not Matrix + # It gets converted to matrix below. + # Make the new indices to go onto Q and R + q, r = inds(T) + q = dim(q) < dim(r) ? sim(q) : sim(r) + IndsT = indstype(T) #get the index type + Qinds = IndsT((ind(T, 1), q)) + Xinds = IndsT((q, ind(T, 2))) + QM = convert(typeof(XM), QM) + ## Here I convert QM twice because of an issue in CUDA where convert does not take QM to be a UnifiedBuffer array + QM = convert(typeof(XM), QM) + Q = tensor(Dense(vec(QM)), Qinds) #Q was strided + X = tensor(Dense(vec(XM)), Xinds) + return Q, X end # Version of `sign` that returns one # if `x == 0`. function nonzero_sign(x) - iszero(x) && return one(x) - return sign(x) + iszero(x) && return one(x) + return sign(x) end # @@ -379,15 +379,15 @@ non-negative. Such a QR decomposition of a matrix is unique. Returns a tuple (Q,R). """ function qr_positive(M::AbstractMatrix) - sparseQ, R = qr(M) - Q = convert(typeof(R), sparseQ) - signs = nonzero_sign.(diag(R)) - Q = Q * Diagonal(signs) - R = Diagonal(conj.(signs)) * R - return (Q, R) + sparseQ, R = qr(M) + Q = convert(typeof(R), sparseQ) + signs = nonzero_sign.(diag(R)) + Q = Q * Diagonal(signs) + R = Diagonal(conj.(signs)) * R + return (Q, R) end -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type """ ql_positive(M::AbstractMatrix) @@ -397,23 +397,23 @@ non-negative. Such a QL decomposition of a matrix is unique. Returns a tuple (Q,L). """ function ql_positive(M::AbstractMatrix) - # TODO: Change to `isgpu`, or better yet rewrite - # in terms of broadcasting and linear algebra - # like `qr_positive`. - sparseQ, L = ql(M) - Q = convert(typeof(L), sparseQ) - nr, nc = size(L) - dc = nc > nr ? nc - nr : 0 #diag is shifted over by dc if nc>nr - for c in 1:(nc - dc) - if L[c, c + dc] != 0.0 #sign(0.0)==0.0 so we don't want to zero out a column of Q. - sign_Lc = sign(L[c, c + dc]) - if c <= nr && !isone(sign_Lc) - L[c, 1:(c + dc)] *= sign_Lc #only fip non-zero portion of the column. - Q[:, c] *= conj(sign_Lc) - end + # TODO: Change to `isgpu`, or better yet rewrite + # in terms of broadcasting and linear algebra + # like `qr_positive`. + sparseQ, L = ql(M) + Q = convert(typeof(L), sparseQ) + nr, nc = size(L) + dc = nc > nr ? nc - nr : 0 #diag is shifted over by dc if nc>nr + for c in 1:(nc - dc) + if L[c, c + dc] != 0.0 #sign(0.0)==0.0 so we don't want to zero out a column of Q. + sign_Lc = sign(L[c, c + dc]) + if c <= nr && !isone(sign_Lc) + L[c, 1:(c + dc)] *= sign_Lc #only fip non-zero portion of the column. + Q[:, c] *= conj(sign_Lc) + end + end end - end - return (Q, L) + return (Q, L) end # @@ -421,63 +421,63 @@ end # before letting lapack overwirte it. # function ql(A::AbstractMatrix) - Base.require_one_based_indexing(A) - T = eltype(A) - AA = similar(A, LinearAlgebra._qreltype(T), size(A)) - copyto!(expose(AA), expose(A)) - Q, L = ql!(AA) - return (Q, L) + Base.require_one_based_indexing(A) + T = eltype(A) + AA = similar(A, LinearAlgebra._qreltype(T), size(A)) + copyto!(expose(AA), expose(A)) + Q, L = ql!(AA) + return (Q, L) end # # This is where the low level call to lapack actually occurs. Most of the work is # about unpacking Q and L from the A matrix. # function ql!(A::StridedMatrix{<:LAPACK.BlasFloat}) - ## TODO is this really necessary here, we could create Expose function if - ## we need this function on CU/GPU - if iscu(A) - throw("Error: ql is not implemented in CUDA.jl") - end - tau = Base.similar(A, min(size(A)...)) - x = LAPACK.geqlf!(A, tau) - #save L from the lower portion of A, before orgql! mangles it! - nr, nc = size(A) - mn = min(nr, nc) - L = similar(A, (mn, nc)) - for r in 1:mn - for c in 1:(r + nc - mn) - L[r, c] = A[r + nr - mn, c] - end - for c in (r + 1 + nc - mn):nc - L[r, c] = 0.0 + ## TODO is this really necessary here, we could create Expose function if + ## we need this function on CU/GPU + if iscu(A) + throw("Error: ql is not implemented in CUDA.jl") end - end - # Now we need shift the orth vectors from the right side of Q over the left side, before - if (mn < nc) - for r in 1:nr - for c in 1:mn - A[r, c] = A[r, c + nc - mn] - end + tau = Base.similar(A, min(size(A)...)) + x = LAPACK.geqlf!(A, tau) + #save L from the lower portion of A, before orgql! mangles it! + nr, nc = size(A) + mn = min(nr, nc) + L = similar(A, (mn, nc)) + for r in 1:mn + for c in 1:(r + nc - mn) + L[r, c] = A[r + nr - mn, c] + end + for c in (r + 1 + nc - mn):nc + L[r, c] = 0.0 + end end - for r in 1:nr - A = A[:, 1:mn] #whack the extra columns in A. + # Now we need shift the orth vectors from the right side of Q over the left side, before + if (mn < nc) + for r in 1:nr + for c in 1:mn + A[r, c] = A[r, c + nc - mn] + end + end + for r in 1:nr + A = A[:, 1:mn] #whack the extra columns in A. + end end - end - LAPACK.orgql!(A, tau) - return A, L + LAPACK.orgql!(A, tau) + return A, L end # TODO: support alg keyword argument to choose the svd algorithm -function polar(T::DenseTensor{ElT,2,IndsT}) where {ElT,IndsT} - QM, RM = polar(matrix(T)) - dim = size(QM, 2) - # Make the new indices to go onto Q and R - q = eltype(IndsT)(dim) - # TODO: use push/pushfirst instead of a constructor - # call here - Qinds = IndsT((ind(T, 1), q)) - Rinds = IndsT((q, ind(T, 2))) - Q = tensor(Dense(vec(QM)), Qinds) - R = tensor(Dense(vec(RM)), Rinds) - return Q, R +function polar(T::DenseTensor{ElT, 2, IndsT}) where {ElT, IndsT} + QM, RM = polar(matrix(T)) + dim = size(QM, 2) + # Make the new indices to go onto Q and R + q = eltype(IndsT)(dim) + # TODO: use push/pushfirst instead of a constructor + # call here + Qinds = IndsT((ind(T, 1), q)) + Rinds = IndsT((q, ind(T, 2))) + Q = tensor(Dense(vec(QM)), Qinds) + R = tensor(Dense(vec(RM)), Rinds) + return Q, R end diff --git a/NDTensors/src/linearalgebra/svd.jl b/NDTensors/src/linearalgebra/svd.jl index bf7729e9ff..1dee094836 100644 --- a/NDTensors/src/linearalgebra/svd.jl +++ b/NDTensors/src/linearalgebra/svd.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type # The state of the `svd_recursive` algorithm. function svd_recursive_state(S::AbstractArray, thresh::Float64) return svd_recursive_state(unwrap_array_type(S), S, thresh) diff --git a/NDTensors/src/tensor/set_types.jl b/NDTensors/src/tensor/set_types.jl index 9830b81d93..66c3688f39 100644 --- a/NDTensors/src/tensor/set_types.jl +++ b/NDTensors/src/tensor/set_types.jl @@ -1,31 +1,31 @@ -using TypeParameterAccessors: TypeParameterAccessors, Position, parenttype +using .Vendored.TypeParameterAccessors: TypeParameterAccessors, Position, parenttype function TypeParameterAccessors.set_ndims(arraytype::Type{<:Tensor}, ndims) - # TODO: Implement something like: - # ```julia - # return set_storagetype(arraytype, set_ndims(storagetype(arraytype), ndims)) - # ``` - # However, we will also need to define `set_ndims(indstype(arraytype), ndims)` - # and use `set_indstype(arraytype, set_ndims(indstype(arraytype), ndims))`. - return error( - "Setting the number dimensions of the array type `$arraytype` (to `$ndims`) is not currently defined.", - ) + # TODO: Implement something like: + # ```julia + # return set_storagetype(arraytype, set_ndims(storagetype(arraytype), ndims)) + # ``` + # However, we will also need to define `set_ndims(indstype(arraytype), ndims)` + # and use `set_indstype(arraytype, set_ndims(indstype(arraytype), ndims))`. + return error( + "Setting the number dimensions of the array type `$arraytype` (to `$ndims`) is not currently defined.", + ) end function set_storagetype(tensortype::Type{<:Tensor}, storagetype) - return Tensor{eltype(tensortype),ndims(tensortype),storagetype,indstype(tensortype)} + return Tensor{eltype(tensortype), ndims(tensortype), storagetype, indstype(tensortype)} end # TODO: Modify the `storagetype` according to `inds`, such as the dimensions? # TODO: Make a version that accepts `indstype::Type`? function TypeParameterAccessors.set_indstype(tensortype::Type{<:Tensor}, inds::Tuple) - return Tensor{eltype(tensortype),length(inds),storagetype(tensortype),typeof(inds)} + return Tensor{eltype(tensortype), length(inds), storagetype(tensortype), typeof(inds)} end TypeParameterAccessors.parenttype(tensortype::Type{<:Tensor}) = storagetype(tensortype) function TypeParameterAccessors.parenttype(storagetype::Type{<:TensorStorage}) - return datatype(storagetype) + return datatype(storagetype) end function TypeParameterAccessors.position(::Type{<:Tensor}, ::typeof(parenttype)) - return Position(3) + return Position(3) end diff --git a/NDTensors/src/tensor/similar.jl b/NDTensors/src/tensor/similar.jl index 3facf61e36..1a64348cad 100644 --- a/NDTensors/src/tensor/similar.jl +++ b/NDTensors/src/tensor/similar.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: TypeParameterAccessors, set_indstype, similartype +using .Vendored.TypeParameterAccessors: TypeParameterAccessors, set_indstype, similartype # NDTensors.similar similar(tensor::Tensor) = setstorage(tensor, similar(storage(tensor))) @@ -8,42 +8,42 @@ similar(tensor::Tensor, eltype::Type) = setstorage(tensor, similar(storage(tenso # NDTensors.similar function similar(tensor::Tensor, dims::Tuple) - return setinds(setstorage(tensor, similar(storage(tensor), dims)), dims) + return setinds(setstorage(tensor, similar(storage(tensor), dims)), dims) end # NDTensors.similar function similar(tensor::Tensor, dims::Dims) - return setinds(setstorage(tensor, similar(storage(tensor), dims)), dims) + return setinds(setstorage(tensor, similar(storage(tensor), dims)), dims) end # NDTensors.similar function similar(tensortype::Type{<:Tensor}, dims::Tuple) - # TODO: Is there a better constructor pattern for this? - # Maybe use `setstorage(::Type{<:Tensor}, ...)` and - # `setinds(::Type{<:Tensor}, ...)`? - return similartype(tensortype, dims)( - AllowAlias(), similar(storagetype(tensortype), dims), dims - ) + # TODO: Is there a better constructor pattern for this? + # Maybe use `setstorage(::Type{<:Tensor}, ...)` and + # `setinds(::Type{<:Tensor}, ...)`? + return similartype(tensortype, dims)( + AllowAlias(), similar(storagetype(tensortype), dims), dims + ) end # NDTensors.similar function similar(tensortype::Type{<:Tensor}, dims::Dims) - # TODO: Is there a better constructor pattern for this? - # Maybe use `setstorage(::Type{<:Tensor}, ...)` and - # `setinds(::Type{<:Tensor}, ...)`? - return similartype(tensortype, dims)( - AllowAlias(), similar(storagetype(tensortype), dims), dims - ) + # TODO: Is there a better constructor pattern for this? + # Maybe use `setstorage(::Type{<:Tensor}, ...)` and + # `setinds(::Type{<:Tensor}, ...)`? + return similartype(tensortype, dims)( + AllowAlias(), similar(storagetype(tensortype), dims), dims + ) end # NDTensors.similar function similar(tensor::Tensor, eltype::Type, dims::Tuple) - return setinds(setstorage(tensor, similar(storage(tensor), eltype, dims)), dims) + return setinds(setstorage(tensor, similar(storage(tensor), eltype, dims)), dims) end # NDTensors.similar function similar(tensor::Tensor, eltype::Type, dims::Dims) - return setinds(setstorage(tensor, similar(storage(tensor), eltype, dims)), dims) + return setinds(setstorage(tensor, similar(storage(tensor), eltype, dims)), dims) end # Base overloads @@ -52,20 +52,20 @@ Base.similar(tensor::Tensor, eltype::Type) = NDTensors.similar(tensor, eltype) Base.similar(tensor::Tensor, dims::Tuple) = NDTensors.similar(tensor, dims) Base.similar(tensor::Tensor, dims::Dims) = NDTensors.similar(tensor, dims) function Base.similar(tensor::Tensor, eltype::Type, dims::Tuple) - return NDTensors.similar(tensor, eltype, dims) + return NDTensors.similar(tensor, eltype, dims) end function Base.similar(tensor::Tensor, eltype::Type, dims::Dims) - return NDTensors.similar(tensor, eltype, dims) + return NDTensors.similar(tensor, eltype, dims) end function TypeParameterAccessors.similartype(tensortype::Type{<:Tensor}, eltype::Type) - return set_storagetype(tensortype, similartype(storagetype(tensortype), eltype)) + return set_storagetype(tensortype, similartype(storagetype(tensortype), eltype)) end function TypeParameterAccessors.similartype(tensortype::Type{<:Tensor}, dims::Tuple) - tensortype_new_inds = set_indstype(tensortype, dims) - # Need to pass `dims` in case that information is needed to make a storage type, - # for example `BlockSparse` needs the number of dimensions. - storagetype_new_inds = similartype(storagetype(tensortype_new_inds), dims) - return set_storagetype(tensortype_new_inds, storagetype_new_inds) + tensortype_new_inds = set_indstype(tensortype, dims) + # Need to pass `dims` in case that information is needed to make a storage type, + # for example `BlockSparse` needs the number of dimensions. + storagetype_new_inds = similartype(storagetype(tensortype_new_inds), dims) + return set_storagetype(tensortype_new_inds, storagetype_new_inds) end diff --git a/NDTensors/src/tensorstorage/default_storage.jl b/NDTensors/src/tensorstorage/default_storage.jl index 872eefc233..c67a210f55 100644 --- a/NDTensors/src/tensorstorage/default_storage.jl +++ b/NDTensors/src/tensorstorage/default_storage.jl @@ -1,21 +1,21 @@ ## This is a fil which specifies the default storage type provided some set of parameters ## The parameters are the element type and storage type -default_datatype(eltype::Type=default_eltype()) = Vector{eltype} +default_datatype(eltype::Type = default_eltype()) = Vector{eltype} default_eltype() = Float64 -using TypeParameterAccessors: specify_default_type_parameters +using .Vendored.TypeParameterAccessors: specify_default_type_parameters ## TODO use multiple dispace to make this pick between dense and blocksparse function default_storagetype(datatype::Type{<:AbstractArray}, inds::Tuple) - datatype = specify_default_type_parameters(datatype) - return Dense{eltype(datatype),datatype} + datatype = specify_default_type_parameters(datatype) + return Dense{eltype(datatype), datatype} end function default_storagetype(datatype::Type{<:AbstractArray}) - return default_storagetype(datatype, ()) + return default_storagetype(datatype, ()) end default_storagetype(eltype::Type) = default_storagetype(default_datatype(eltype)) function default_storagetype(eltype::Type, inds::Tuple) - return default_storagetype(default_datatype(eltype), inds) + return default_storagetype(default_datatype(eltype), inds) end default_storagetype() = default_storagetype(default_eltype()) diff --git a/NDTensors/src/tensorstorage/set_types.jl b/NDTensors/src/tensorstorage/set_types.jl index 0d4156ed2c..87583be20b 100644 --- a/NDTensors/src/tensorstorage/set_types.jl +++ b/NDTensors/src/tensorstorage/set_types.jl @@ -1,7 +1,7 @@ -using TypeParameterAccessors: TypeParameterAccessors +using .Vendored.TypeParameterAccessors: TypeParameterAccessors function TypeParameterAccessors.set_ndims(arraytype::Type{<:TensorStorage}, ndims::Int) - # TODO: Change to this once `TensorStorage` types support wrapping - # non-AbstractVector types. - # return set_datatype(arraytype, set_ndims(datatype(arraytype), ndims)) - return arraytype + # TODO: Change to this once `TensorStorage` types support wrapping + # non-AbstractVector types. + # return set_datatype(arraytype, set_ndims(datatype(arraytype), ndims)) + return arraytype end diff --git a/NDTensors/src/tensorstorage/similar.jl b/NDTensors/src/tensorstorage/similar.jl index 97eba53089..5c12a3e321 100644 --- a/NDTensors/src/tensorstorage/similar.jl +++ b/NDTensors/src/tensorstorage/similar.jl @@ -1,58 +1,58 @@ -using TypeParameterAccessors: TypeParameterAccessors, set_ndims, similartype +using .Vendored.TypeParameterAccessors: TypeParameterAccessors, set_ndims, similartype # NDTensors.similar similar(storage::TensorStorage) = setdata(storage, NDTensors.similar(data(storage))) # NDTensors.similar function similar(storage::TensorStorage, eltype::Type) - return setdata(storage, NDTensors.similar(data(storage), eltype)) + return setdata(storage, NDTensors.similar(data(storage), eltype)) end # NDTensors.similar function similar(storage::TensorStorage, dims::Tuple) - # TODO: Don't convert to an `AbstractVector` with `vec`, once we support - # more general data types. - # return setdata(storage, NDTensors.similar(data(storage), dims)) - return setdata(storage, vec(NDTensors.similar(data(storage), dims))) + # TODO: Don't convert to an `AbstractVector` with `vec`, once we support + # more general data types. + # return setdata(storage, NDTensors.similar(data(storage), dims)) + return setdata(storage, vec(NDTensors.similar(data(storage), dims))) end # NDTensors.similar function similar(storage::TensorStorage, eltype::Type, dims::Tuple) - # TODO: Don't convert to an `AbstractVector` with `vec`, once we support - # more general data types. - # return setdata(storage, NDTensors.similar(data(storage), eltype, dims)) - return setdata(storage, vec(NDTensors.similar(data(storage), eltype, dims))) + # TODO: Don't convert to an `AbstractVector` with `vec`, once we support + # more general data types. + # return setdata(storage, NDTensors.similar(data(storage), eltype, dims)) + return setdata(storage, vec(NDTensors.similar(data(storage), eltype, dims))) end # NDTensors.similar function similar(storagetype::Type{<:TensorStorage}, eltype::Type, dims::Tuple) - return similar(similartype(storagetype, eltype), dims) + return similar(similartype(storagetype, eltype), dims) end # NDTensors.similar function similar(storagetype::Type{<:TensorStorage}, eltype::Type) - return error("Must specify dimensions.") + return error("Must specify dimensions.") end # NDTensors.similar function similar(storagetype::Type{<:TensorStorage}, dims::Tuple) - # TODO: Don't convert to an `AbstractVector` with `vec`, once we support - # more general data types. - # return setdata(storagetype, NDTensors.similar(datatype(storagetype), dims)) - return setdata(storagetype, vec(NDTensors.similar(datatype(storagetype), dims))) + # TODO: Don't convert to an `AbstractVector` with `vec`, once we support + # more general data types. + # return setdata(storagetype, NDTensors.similar(datatype(storagetype), dims)) + return setdata(storagetype, vec(NDTensors.similar(datatype(storagetype), dims))) end # NDTensors.similar function similar(storagetype::Type{<:TensorStorage}, dims::Dims) - # TODO: Don't convert to an `AbstractVector` with `prod`, once we support - # more general data types. - # return setdata(storagetype, NDTensors.similar(datatype(storagetype), dims)) - return setdata(storagetype, NDTensors.similar(datatype(storagetype), prod(dims))) + # TODO: Don't convert to an `AbstractVector` with `prod`, once we support + # more general data types. + # return setdata(storagetype, NDTensors.similar(datatype(storagetype), dims)) + return setdata(storagetype, NDTensors.similar(datatype(storagetype), prod(dims))) end # NDTensors.similar function similar(storagetype::Type{<:TensorStorage}, dims::DimOrInd...) - return similar(storagetype, NDTensors.to_shape(dims)) + return similar(storagetype, NDTensors.to_shape(dims)) end # Define Base.similar in terms of NDTensors.similar @@ -64,20 +64,20 @@ Base.similar(storage::TensorStorage, eltype::Type) = NDTensors.similar(storage, ## Base.similar(storage::TensorStorage, dims::DimOrInd...) = NDTensors.similar(storage, dims...) function TypeParameterAccessors.similartype( - storagetype::Type{<:TensorStorage}, eltype::Type -) - # TODO: Don't convert to an `AbstractVector` with `set_ndims(datatype, 1)`, once we support - # more general data types. - # return set_datatype(storagetype, NDTensors.similartype(datatype(storagetype), eltype)) - return set_datatype(storagetype, set_ndims(similartype(datatype(storagetype), eltype), 1)) + storagetype::Type{<:TensorStorage}, eltype::Type + ) + # TODO: Don't convert to an `AbstractVector` with `set_ndims(datatype, 1)`, once we support + # more general data types. + # return set_datatype(storagetype, NDTensors.similartype(datatype(storagetype), eltype)) + return set_datatype(storagetype, set_ndims(similartype(datatype(storagetype), eltype), 1)) end function TypeParameterAccessors.similartype(storagetype::Type{<:TensorStorage}, dims::Tuple) - # TODO: In the future, set the dimensions of the data type based on `dims`, once - # more general data types beyond `AbstractVector` are supported. - # `similartype` unwraps any wrapped data. - return set_ndims( - set_datatype(storagetype, set_ndims(similartype(datatype(storagetype)), 1)), - length(dims), - ) + # TODO: In the future, set the dimensions of the data type based on `dims`, once + # more general data types beyond `AbstractVector` are supported. + # `similartype` unwraps any wrapped data. + return set_ndims( + set_datatype(storagetype, set_ndims(similartype(datatype(storagetype)), 1)), + length(dims), + ) end diff --git a/NDTensors/src/truncate.jl b/NDTensors/src/truncate.jl index 5f46bb2655..73fefafdcf 100644 --- a/NDTensors/src/truncate.jl +++ b/NDTensors/src/truncate.jl @@ -1,4 +1,4 @@ -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type ## TODO write Exposed version of truncate function truncate!!(P::AbstractArray; kwargs...) return truncate!!(unwrap_array_type(P), P; kwargs...) @@ -10,7 +10,7 @@ function truncate!!(::Type{<:Array}, P::AbstractArray; kwargs...) return P, truncerr, docut end -using TypeParameterAccessors: unwrap_array_type +using .Vendored.TypeParameterAccessors: unwrap_array_type # GPU fallback version, convert to CPU. function truncate!!(::Type{<:AbstractArray}, P::AbstractArray; kwargs...) P_cpu = cpu(P) diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/TypeParameterAccessors.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/TypeParameterAccessors.jl new file mode 100755 index 0000000000..0a51102f6c --- /dev/null +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/TypeParameterAccessors.jl @@ -0,0 +1,23 @@ +module TypeParameterAccessors + +# Exports +export type_parameters, get_type_parameters +export nparameters, is_parameter_specified +export default_type_parameters +export set_type_parameters, set_default_type_parameters +export specify_type_parameters, specify_default_type_parameters +export unspecify_type_parameters + +# Imports +using SimpleTraits: SimpleTraits, @traitdef, @traitimpl + +include("type_utils.jl") +include("type_parameters.jl") + +# Implementations +include("ndims.jl") +include("base/abstractarray.jl") +include("base/similartype.jl") +include("base/linearalgebra.jl") + +end diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/base/abstractarray.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/base/abstractarray.jl new file mode 100755 index 0000000000..1d9e211ef9 --- /dev/null +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/base/abstractarray.jl @@ -0,0 +1,94 @@ +struct Self end +position(a, ::Self) = Position(0) +position(::Type, ::Self) = Position(0) +function set_type_parameters(type::Type, ::Self, param) + return error("Can't set the parent type of an unwrapped array type.") +end + +position(::Type{AbstractArray}, ::typeof(eltype)) = Position(1) +position(::Type{AbstractArray}, ::typeof(ndims)) = Position(2) +default_type_parameters(::Type{AbstractArray}) = (Float64, 1) + +position(::Type{<:Array}, ::typeof(eltype)) = Position(1) +position(::Type{<:Array}, ::typeof(ndims)) = Position(2) +default_type_parameters(::Type{<:Array}) = (Float64, 1) + +position(::Type{<:BitArray}, ::typeof(ndims)) = Position(1) +default_type_parameters(::Type{<:BitArray}) = (1,) + +function set_eltype(array::AbstractArray, param) + return convert(set_eltype(typeof(array), param), array) +end + +## This will fail if position of `ndims` is not defined for `type` +function set_ndims(type::Type{<:AbstractArray}, param) + return set_type_parameters(type, ndims, param) +end +function set_ndims(type::Type{<:AbstractArray}, param::NDims) + return set_type_parameters(type, ndims, ndims(param)) +end + +# Trait indicating if the AbstractArray type is an array wrapper. +# Assumes that it implements `NDTensors.parenttype`. +@traitdef IsWrappedArray{ArrayType} + +#! format: off +@traitimpl IsWrappedArray{ArrayType} <- is_wrapped_array(ArrayType) +#! format: on + +parenttype(type::Type{<:AbstractArray}) = type_parameters(type, parenttype) +parenttype(object::AbstractArray) = parenttype(typeof(object)) +position(::Type{<:AbstractArray}, ::typeof(parenttype)) = Self() + +is_wrapped_array(arraytype::Type{<:AbstractArray}) = (parenttype(arraytype) ≠ arraytype) +@inline is_wrapped_array(array::AbstractArray) = is_wrapped_array(typeof(array)) +@inline is_wrapped_array(object) = false + +using SimpleTraits: Not, @traitfn + +@traitfn function unwrap_array_type( + arraytype::Type{ArrayType} +) where {ArrayType;IsWrappedArray{ArrayType}} + return unwrap_array_type(parenttype(arraytype)) +end + +@traitfn function unwrap_array_type( + arraytype::Type{ArrayType} +) where {ArrayType;!IsWrappedArray{ArrayType}} + return arraytype +end + +# For working with instances. +unwrap_array_type(array::AbstractArray) = unwrap_array_type(typeof(array)) + +function set_parenttype(t::Type, param) + return set_type_parameters(t, parenttype, param) +end + +@traitfn function set_eltype( + type::Type{ArrayType}, param +) where {ArrayType<:AbstractArray;IsWrappedArray{ArrayType}} + new_parenttype = set_eltype(parenttype(type), param) + # Need to set both in one `set_type_parameters` call to avoid + # conflicts in type parameter constraints of certain wrapper types. + return set_type_parameters(type, (eltype, parenttype), (param, new_parenttype)) +end + +@traitfn function set_eltype( + type::Type{ArrayType}, param +) where {ArrayType<:AbstractArray;!IsWrappedArray{ArrayType}} + return set_type_parameters(type, eltype, param) +end + +for wrapper in [:PermutedDimsArray, :(Base.ReshapedArray), :SubArray] + @eval begin + position(type::Type{<:$wrapper}, ::typeof(eltype)) = Position(1) + position(type::Type{<:$wrapper}, ::typeof(ndims)) = Position(2) + end +end +for wrapper in [:(Base.ReshapedArray), :SubArray] + @eval position(type::Type{<:$wrapper}, ::typeof(parenttype)) = Position(3) +end +for wrapper in [:PermutedDimsArray] + @eval position(type::Type{<:$wrapper}, ::typeof(parenttype)) = Position(5) +end diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/base/linearalgebra.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/base/linearalgebra.jl new file mode 100755 index 0000000000..3d818cab16 --- /dev/null +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/base/linearalgebra.jl @@ -0,0 +1,25 @@ +using LinearAlgebra: + Adjoint, + Diagonal, + Hermitian, + LowerTriangular, + Symmetric, + Transpose, + UnitLowerTriangular, + UnitUpperTriangular, + UpperTriangular + +for wrapper in [ + :Transpose, + :Adjoint, + :Symmetric, + :Hermitian, + :UpperTriangular, + :LowerTriangular, + :UnitUpperTriangular, + :UnitLowerTriangular, + :Diagonal, +] + @eval position(::Type{<:$wrapper}, ::typeof(eltype)) = Position(1) + @eval position(::Type{<:$wrapper}, ::typeof(parenttype)) = Position(2) +end diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/base/similartype.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/base/similartype.jl new file mode 100755 index 0000000000..2d8175b012 --- /dev/null +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/base/similartype.jl @@ -0,0 +1,93 @@ +""" +`set_indstype` should be overloaded for +types with structured dimensions, +like `OffsetArrays` or named indices +(such as ITensors). +""" +function set_indstype(arraytype::Type{<:AbstractArray}, dims::Tuple) + return set_ndims(arraytype, NDims(length(dims))) +end + +function similartype(arraytype::Type{<:AbstractArray}, eltype::Type, ndims::NDims) + return similartype(similartype(arraytype, eltype), ndims) +end + +function similartype(arraytype::Type{<:AbstractArray}, eltype::Type, dims::Tuple) + return similartype(similartype(arraytype, eltype), dims) +end + +@traitfn function similartype( + arraytype::Type{ArrayT} +) where {ArrayT;!IsWrappedArray{ArrayT}} + return arraytype +end + +@traitfn function similartype( + arraytype::Type{ArrayT}, eltype::Type +) where {ArrayT;!IsWrappedArray{ArrayT}} + return set_eltype(arraytype, eltype) +end + +function similartype(arraytype::Type{<:BitArray}, elt::Type) + if is_parameter_specified(arraytype, ndims) + return Array{elt,ndims(arraytype)} + end + return Array{elt} +end + +@traitfn function similartype( + arraytype::Type{ArrayT}, dims::Tuple +) where {ArrayT;!IsWrappedArray{ArrayT}} + return set_indstype(arraytype, dims) +end + +@traitfn function similartype( + arraytype::Type{ArrayT}, ndims::NDims +) where {ArrayT;!IsWrappedArray{ArrayT}} + return set_ndims(arraytype, ndims) +end + +function similartype( + arraytype::Type{<:AbstractArray}, dim1::Base.DimOrInd, dim_rest::Base.DimOrInd... +) + return similartype(arraytype, (dim1, dim_rest...)) +end + +## Wrapped arrays +@traitfn function similartype(arraytype::Type{ArrayT}) where {ArrayT;IsWrappedArray{ArrayT}} + return similartype(unwrap_array_type(arraytype), NDims(arraytype)) +end + +@traitfn function similartype( + arraytype::Type{ArrayT}, eltype::Type +) where {ArrayT;IsWrappedArray{ArrayT}} + return similartype(unwrap_array_type(arraytype), eltype, NDims(arraytype)) +end + +@traitfn function similartype( + arraytype::Type{ArrayT}, dims::Tuple +) where {ArrayT;IsWrappedArray{ArrayT}} + return similartype(unwrap_array_type(arraytype), dims) +end + +@traitfn function similartype( + arraytype::Type{ArrayT}, ndims::NDims +) where {ArrayT;IsWrappedArray{ArrayT}} + return similartype(unwrap_array_type(arraytype), ndims) +end + +# This is for uniform `Diag` storage which uses +# a Number as the data type. +# TODO: Delete this when we change to using a +# `FillArray` instead. This is a stand-in +# to make things work with the current design. +function similartype(numbertype::Type{<:Number}) + return numbertype +end + +# Instances +function similartype(array::AbstractArray, eltype::Type, dims...) + return similartype(typeof(array), eltype, dims...) +end +similartype(array::AbstractArray, eltype::Type) = similartype(typeof(array), eltype) +similartype(array::AbstractArray, dims...) = similartype(typeof(array), dims...) diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/ndims.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/ndims.jl new file mode 100755 index 0000000000..f640763583 --- /dev/null +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/ndims.jl @@ -0,0 +1,6 @@ +struct NDims{ndims} end +Base.ndims(::NDims{ndims}) where {ndims} = ndims + +NDims(ndims::Integer) = NDims{ndims}() +NDims(arraytype::Type{<:AbstractArray}) = NDims(ndims(arraytype)) +NDims(array::AbstractArray) = NDims(typeof(array)) diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/type_parameters.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/type_parameters.jl new file mode 100755 index 0000000000..1dd31ce918 --- /dev/null +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/type_parameters.jl @@ -0,0 +1,331 @@ +""" + struct Position{P} end + +Singleton type to statically represent the type-parameter position. +This is meant for internal use as a `Val`-like structure to improve type-inference. +""" +struct Position{P} end +Position(pos::Int) = Position{pos}() + +Base.Int(pos::Position) = Int(typeof(pos)) +Base.Int(::Type{Position{P}}) where {P} = Int(P) +Base.to_index(pos::Position) = Base.to_index(typeof(pos)) +Base.to_index(::Type{P}) where {P<:Position} = Int(P) + +""" + position(type::Type, position_name)::Position + +An optional interface function. Defining this allows accessing a parameter +at the defined position using the `position_name`. + +For example, defining `TypeParameterAccessors.position(::Type{<:MyType}, ::typeof(eltype)) = Position(1)` +allows accessing the first type parameter with `type_parameters(MyType(...), eltype)`, +in addition to the standard `type_parameters(MyType(...), 1)` or `type_parameters(MyType(...), Position(1))`. +""" +function position end + +position(object, name) = position(typeof(object), name) +position(::Type, pos::Int) = Position(pos) +position(::Type, pos::Position) = pos + +function position(type::Type, name) + type′ = unspecify_type_parameters(type) + if type === type′ + # Fallback definition that determines the + # position automatically from the supertype of + # the type. + return position_from_supertype(type′, name) + end + return position(type′, name) +end + +# Automatically determine the position of a type parameter of a type given +# a supertype and the name of the parameter. +function position_from_supertype(type::Type, name) + type′ = unspecify_type_parameters(type) + supertype_pos = position(supertype(type′), name) + return position_from_supertype_position(type′, supertype_pos) +end + +# Automatically determine the position of a type parameter of a type given +# the supertype and the position of the corresponding parameter in the supertype. +@generated function position_from_supertype_position( + ::Type{T}, supertype_pos::Position +) where {T} + T′ = unspecify_type_parameters(T) + # The type parameters of the type as TypeVars. + # TODO: Ideally we would use `get_type_parameters` + # but that sometimes loses TypeVar names: + # https://github.com/ITensor/TypeParameterAccessors.jl/issues/30 + type_params = Base.unwrap_unionall(T′).parameters + # The type parameters of the immediate supertype as TypeVars. + # This has TypeVars with names that correspond to the names of + # the TypeVars of the type parameters of `T`, for example: + # ```julia + # julia> struct MyArray{B,A} <: AbstractArray{A,B} end + # + # julia> Base.unwrap_unionall(MyArray).parameters + # svec(B, A) + # + # julia> Base.unwrap_unionall(supertype(MyArray)).parameters + # svec(A, B) + # ``` + supertype_params = Base.unwrap_unionall(supertype(T′)).parameters + supertype_param = supertype_params[Int(supertype_pos)] + if !(supertype_param isa TypeVar) + error("Position not found.") + end + pos = findfirst(param -> (param.name == supertype_param.name), type_params) + if isnothing(pos) + return error("Position not found.") + end + return :(@inline; $(Position(pos))) +end + +function positions(::Type{T}, pos::Tuple) where {T} + return ntuple(length(pos)) do i + return position(T, pos[i]) + end +end + +""" + get_type_parameters(type_or_obj, [pos]) + +Return a tuple containing the type parameters of a given type or object. +Optionally you can specify a position to just get the parameter for that position, +or a tuple of positions to get a subset of parameters. + +If parameters are unspecified, returns a `TypeVar`. For a checked version, +see [`type_parameters`](@ref). +""" +function get_type_parameters end + +# This implementation is type-stable in 1.11, but not in 1.10. +# Attempts with `Base.@constprop :aggressive` failed, so generated function instead +# @inline type_parameters(::Type{T}) where {T} = Tuple(Base.unwrap_unionall(T).parameters) +@generated function get_type_parameters(::Type{T}) where {T} + params = wrap_symbol_quotenode.(Tuple(Base.unwrap_unionall(T).parameters)) + return :(@inline; ($(params...),)) +end +@inline get_type_parameters(::Type{T}, pos) where {T} = get_type_parameters( + T, position(T, pos) +) +@inline get_type_parameters(::Type{T}, ::Position{p}) where {T,p} = get_type_parameters(T)[p] +@inline get_type_parameters(::Type{T}, ::Position{0}) where {T} = T +@inline get_type_parameters(::Type{T}, pos::Tuple) where {T} = get_type_parameters.(T, pos) +@inline get_type_parameters(object, pos) = get_type_parameters(typeof(object), pos) +@inline get_type_parameters(object) = get_type_parameters(typeof(object)) + +""" + type_parameters(type_or_obj, [pos]) + +Return a tuple containing the type parameters of a given type or object. +Optionally you can specify a position to just get the parameter for that position, +or a tuple of positions to get a subset of parameters. + +Errors if parameters are unspecified. For an unchecked version, +see [`get_type_parameters`](@ref). +""" +function type_parameters end + +function type_parameters(::Type{T}) where {T} + params = get_type_parameters(T) + any(param -> param isa TypeVar, params) && + return error("One or more parameter is not specified.") + return params +end +@inline function type_parameters(::Type{T}, pos) where {T} + param = get_type_parameters(T, pos) + param isa TypeVar && return error("The parameter is not specified.") + return param +end +@inline type_parameters(::Type{T}, pos::Tuple) where {T} = type_parameters.(T, pos) +@inline type_parameters(object, pos) = type_parameters(typeof(object), pos) +@inline type_parameters(object) = type_parameters(typeof(object)) + +""" + nparameters(type_or_obj) + +Return the number of type parameters for a given type or object. +""" +nparameters(object) = nparameters(typeof(object)) +nparameters(::Type{T}) where {T} = length(get_type_parameters(T)) + +""" + is_parameter_specified(type::Type, pos) + +Return whether or not the type parameter at a given position is considered specified. +""" +function is_parameter_specified(::Type{T}, pos) where {T} + return !(get_type_parameters(T, pos) isa TypeVar) +end + +""" + unspecify_type_parameters(type::Type, [positions::Tuple]) + unspecify_type_parameters(type::Type, position) + +Return a new type where the type parameters at the given positions are unset. +""" +unspecify_type_parameters(::Type{T}) where {T} = Base.typename(T).wrapper +function unspecify_type_parameters(::Type{T}, pos::Tuple) where {T} + @inline + return unspecify_type_parameters(T, positions(T, pos)) +end +@generated function unspecify_type_parameters( + ::Type{T}, positions::Tuple{Vararg{Position}} +) where {T} + allparams = collect(Any, get_type_parameters(T)) + for pos in type_parameters(positions) + allparams[pos] = get_type_parameters(unspecify_type_parameters(T), Int(pos)) + end + type_expr = construct_type_expr(T, allparams) + return :(@inline; $type_expr) +end +unspecify_type_parameters(::Type{T}, pos) where {T} = unspecify_type_parameters(T, (pos,)) + +""" + set_type_parameters(type::Type, positions::Tuple, parameters::Tuple) + set_type_parameters(type::Type, position, parameter) + +Return a new type where the type parameters at the given positions are set to the provided values. +""" +function set_type_parameters( + ::Type{T}, pos::Tuple{Vararg{Any,N}}, parameters::Tuple{Vararg{Any,N}} +) where {T,N} + return set_type_parameters(T, positions(T, pos), parameters) +end +@generated function set_type_parameters( + ::Type{T}, positions::Tuple{Vararg{Position,N}}, params::Tuple{Vararg{Any,N}} +) where {T,N} + # collect parameters and change + allparams = collect(Any, get_type_parameters(T)) + for (i, pos) in enumerate(get_type_parameters(positions)) + allparams[pos] = :(params[$i]) + end + type_expr = construct_type_expr(T, allparams) + return :(@inline; $type_expr) +end +function set_type_parameters(::Type{T}, pos, param) where {T} + return set_type_parameters(T, (pos,), (param,)) +end + +""" + specify_type_parameters(type::Type, positions::Tuple, parameters::Tuple) + specify_type_parameters(type::Type, position, parameter) + +Return a new type where the type parameters at the given positions are set to the provided values, +only if they were previously unspecified. +""" +function specify_type_parameters( + ::Type{T}, pos::Tuple{Vararg{Any,N}}, parameters::Tuple{Vararg{Any,N}} +) where {T,N} + return specify_type_parameters(T, positions(T, pos), parameters) +end +function specify_type_parameters(::Type{T}, parameters::Tuple) where {T} + return specify_type_parameters(T, ntuple(identity, nparameters(T)), parameters) +end +@generated function specify_type_parameters( + ::Type{T}, positions::Tuple{Vararg{Position,N}}, params::Tuple{Vararg{Any,N}} +) where {T,N} + # collect parameters and change unspecified + allparams = collect(Any, get_type_parameters(T)) + for (i, pos) in enumerate(get_type_parameters(positions)) + if !is_parameter_specified(T, pos()) + allparams[pos] = :(params[$i]) + end + end + type_expr = construct_type_expr(T, allparams) + return :(@inline; $type_expr) +end +function specify_type_parameters(::Type{T}, pos, param) where {T} + return specify_type_parameters(T, (pos,), (param,)) +end + +""" + default_type_parameters(type::Type)::Tuple + +An optional interface function. Defining this allows filling type parameters +of the specified type with default values. + +This function should output a Tuple of the default values, with exactly +one for each type parameter slot of the type. +""" +function default_type_parameters(::Type{T}, pos) where {T} + return default_type_parameters(T, position(T, pos)) +end +function default_type_parameters(::Type{T}, ::Position{pos}) where {T,pos} + param = default_type_parameters(T)[pos] + if param isa UndefinedDefaultTypeParameter + return error("No default parameter defined at this position.") + end + return param +end +default_type_parameters(::Type{T}, pos::Tuple) where {T} = default_type_parameters.(T, pos) +default_type_parameters(t, pos) = default_type_parameters(typeof(t), pos) +default_type_parameters(t) = default_type_parameters(typeof(t)) +function default_type_parameters(type::Type) + type′ = unspecify_type_parameters(type) + if type === type′ + return default_type_parameters_from_supertype(type′) + end + return default_type_parameters(type′) +end + +struct UndefinedDefaultTypeParameter end + +# Determine the default type parameters of a type from the default type +# parameters of the supertype of the type. Uses similar logic as +# `position_from_supertype_position` for matching TypeVar names +# between the type and the supertype. Type parameters that exist +# in the type but not the supertype will have a default type parameter +# `UndefinedDefaultTypeParameter()`. Accessing those type parameters +# by name/position will throw an error. +@generated function default_type_parameters_from_supertype(::Type{T}) where {T} + T′ = unspecify_type_parameters(T) + supertype_default_type_params = default_type_parameters(supertype(T′)) + type_params = Base.unwrap_unionall(T′).parameters + supertype_params = Base.unwrap_unionall(supertype(T′)).parameters + defaults = Any[UndefinedDefaultTypeParameter() for _ in 1:nparameters(T′)] + for (supertype_param, supertype_default_type_param) in + zip(supertype_params, supertype_default_type_params) + if !(supertype_param isa TypeVar) + continue + end + param_position = findfirst(param -> (param.name == supertype_param.name), type_params) + defaults[param_position] = supertype_default_type_param + end + return :(@inline; $(Tuple(defaults))) +end + +""" + set_default_type_parameters(type::Type, [positions::Tuple]) + set_default_type_parameters(type::Type, position) + +Set the type parameters at the given positions to their default values. +""" +function set_default_type_parameters(::Type{T}, pos::Tuple) where {T} + return set_type_parameters(T, pos, default_type_parameters.(T, pos)) +end +function set_default_type_parameters(::Type{T}) where {T} + return set_default_type_parameters(T, ntuple(identity, nparameters(T))) +end +function set_default_type_parameters(::Type{T}, pos) where {T} + return set_default_type_parameters(T, (pos,)) +end + +""" + specify_default_type_parameters(type::Type, [positions::Tuple]) + specify_default_type_parameters(type::Type, position) + +Set the type parameters at the given positions to their default values, if they +had not been specified. +""" +function specify_default_type_parameters(::Type{T}, pos::Tuple) where {T} + return specify_type_parameters(T, pos, default_type_parameters.(T, pos)) +end +function specify_default_type_parameters(::Type{T}) where {T} + return specify_default_type_parameters(T, ntuple(identity, nparameters(T))) +end +function specify_default_type_parameters(::Type{T}, pos) where {T} + return specify_default_type_parameters(T, (pos,)) +end diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/type_utils.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/type_utils.jl new file mode 100755 index 0000000000..79a316dae7 --- /dev/null +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/type_utils.jl @@ -0,0 +1,14 @@ +"replace `Symbol`s with `QuoteNode`s to avoid expression interpolation" +wrap_symbol_quotenode(param) = param isa Symbol ? QuoteNode(param) : param + +"Construct the expression for qualifying a type with given parameters" +function construct_type_expr(type, parameters) + basetype = unspecify_type_parameters(type) + type_expr = Expr(:curly, basetype, wrap_symbol_quotenode.(parameters)...) + for parameter in reverse(parameters) + if parameter isa TypeVar + type_expr = Expr(:call, :UnionAll, parameter, type_expr) + end + end + return type_expr +end diff --git a/NDTensors/test/Project.toml b/NDTensors/test/Project.toml index 90742633a0..9193c2fefc 100644 --- a/NDTensors/test/Project.toml +++ b/NDTensors/test/Project.toml @@ -20,7 +20,6 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] diff --git a/NDTensors/vendor/Project.toml b/NDTensors/vendor/Project.toml new file mode 100644 index 0000000000..4c7dfde138 --- /dev/null +++ b/NDTensors/vendor/Project.toml @@ -0,0 +1,5 @@ +[deps] +PackageAnalyzer = "e713c705-17e4-4cec-abe0-95bf5bf3e10c" + +[compat] +PackageAnalyzer = "3.1" diff --git a/NDTensors/vendor/README.md b/NDTensors/vendor/README.md new file mode 100644 index 0000000000..7698ea761d --- /dev/null +++ b/NDTensors/vendor/README.md @@ -0,0 +1,13 @@ +# vendor + +Here we vendor some of our dependencies to load our own copies of them. This means: + +- the functions defined in them are not shared with the "real" package they come from, so we do not see method extensions of them +- we have our own private copy of the module, allowing us to load a different version than the user + +To update the dependencies, install the `vendor` environment and invoke `run.jl`: + +```sh +julia --project=vendor -e 'using Pkg; Pkg.instantiate()' +julia --project=vendor vendor/run.jl +``` \ No newline at end of file diff --git a/NDTensors/vendor/run.jl b/NDTensors/vendor/run.jl new file mode 100644 index 0000000000..0544fcddcc --- /dev/null +++ b/NDTensors/vendor/run.jl @@ -0,0 +1,18 @@ +using PackageAnalyzer: PackageAnalyzer + +deps = [("TypeParameterAccessors", v"0.3.10")] + +for (name, version) in deps + pkg = PackageAnalyzer.find_package(name; version) + local_path, reachable, _ = PackageAnalyzer.obtain_code(pkg) + @assert reachable + p = mkpath(joinpath(@__DIR__, "..", "src", "vendored", name)) + # remove any existing files + if isdir(p) + rm(p; recursive = true, force = true) + end + mkpath(joinpath(p, "src")) + cp(joinpath(local_path, "src"), joinpath(p, "src"); force = true) + # Make files writable so we can format. + chmod(p, 0o755; recursive = true) +end From e370eb0e648f1898d17864a14c47b6a12c6e11f2 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 10 Oct 2025 17:18:45 -0400 Subject: [PATCH 2/2] Runic format --- .../src/base/abstractarray.jl | 52 ++-- .../src/base/linearalgebra.jl | 42 +-- .../src/base/similartype.jl | 70 ++--- .../src/type_parameters.jl | 252 +++++++++--------- .../TypeParameterAccessors/src/type_utils.jl | 14 +- NDTensors/vendor/Project.toml | 2 + NDTensors/vendor/run.jl | 2 + 7 files changed, 219 insertions(+), 215 deletions(-) diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/base/abstractarray.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/base/abstractarray.jl index 1d9e211ef9..f45131e688 100755 --- a/NDTensors/src/vendored/TypeParameterAccessors/src/base/abstractarray.jl +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/base/abstractarray.jl @@ -2,7 +2,7 @@ struct Self end position(a, ::Self) = Position(0) position(::Type, ::Self) = Position(0) function set_type_parameters(type::Type, ::Self, param) - return error("Can't set the parent type of an unwrapped array type.") + return error("Can't set the parent type of an unwrapped array type.") end position(::Type{AbstractArray}, ::typeof(eltype)) = Position(1) @@ -17,15 +17,15 @@ position(::Type{<:BitArray}, ::typeof(ndims)) = Position(1) default_type_parameters(::Type{<:BitArray}) = (1,) function set_eltype(array::AbstractArray, param) - return convert(set_eltype(typeof(array), param), array) + return convert(set_eltype(typeof(array), param), array) end ## This will fail if position of `ndims` is not defined for `type` function set_ndims(type::Type{<:AbstractArray}, param) - return set_type_parameters(type, ndims, param) + return set_type_parameters(type, ndims, param) end function set_ndims(type::Type{<:AbstractArray}, param::NDims) - return set_type_parameters(type, ndims, ndims(param)) + return set_type_parameters(type, ndims, ndims(param)) end # Trait indicating if the AbstractArray type is an array wrapper. @@ -47,48 +47,48 @@ is_wrapped_array(arraytype::Type{<:AbstractArray}) = (parenttype(arraytype) ≠ using SimpleTraits: Not, @traitfn @traitfn function unwrap_array_type( - arraytype::Type{ArrayType} -) where {ArrayType;IsWrappedArray{ArrayType}} - return unwrap_array_type(parenttype(arraytype)) + arraytype::Type{ArrayType} + ) where {{ArrayType; IsWrappedArray{ArrayType}}} + return unwrap_array_type(parenttype(arraytype)) end @traitfn function unwrap_array_type( - arraytype::Type{ArrayType} -) where {ArrayType;!IsWrappedArray{ArrayType}} - return arraytype + arraytype::Type{ArrayType} + ) where {{ArrayType; !IsWrappedArray{ArrayType}}} + return arraytype end # For working with instances. unwrap_array_type(array::AbstractArray) = unwrap_array_type(typeof(array)) function set_parenttype(t::Type, param) - return set_type_parameters(t, parenttype, param) + return set_type_parameters(t, parenttype, param) end @traitfn function set_eltype( - type::Type{ArrayType}, param -) where {ArrayType<:AbstractArray;IsWrappedArray{ArrayType}} - new_parenttype = set_eltype(parenttype(type), param) - # Need to set both in one `set_type_parameters` call to avoid - # conflicts in type parameter constraints of certain wrapper types. - return set_type_parameters(type, (eltype, parenttype), (param, new_parenttype)) + type::Type{ArrayType}, param + ) where {{ArrayType <: AbstractArray; IsWrappedArray{ArrayType}}} + new_parenttype = set_eltype(parenttype(type), param) + # Need to set both in one `set_type_parameters` call to avoid + # conflicts in type parameter constraints of certain wrapper types. + return set_type_parameters(type, (eltype, parenttype), (param, new_parenttype)) end @traitfn function set_eltype( - type::Type{ArrayType}, param -) where {ArrayType<:AbstractArray;!IsWrappedArray{ArrayType}} - return set_type_parameters(type, eltype, param) + type::Type{ArrayType}, param + ) where {{ArrayType <: AbstractArray; !IsWrappedArray{ArrayType}}} + return set_type_parameters(type, eltype, param) end for wrapper in [:PermutedDimsArray, :(Base.ReshapedArray), :SubArray] - @eval begin - position(type::Type{<:$wrapper}, ::typeof(eltype)) = Position(1) - position(type::Type{<:$wrapper}, ::typeof(ndims)) = Position(2) - end + @eval begin + position(type::Type{<:$wrapper}, ::typeof(eltype)) = Position(1) + position(type::Type{<:$wrapper}, ::typeof(ndims)) = Position(2) + end end for wrapper in [:(Base.ReshapedArray), :SubArray] - @eval position(type::Type{<:$wrapper}, ::typeof(parenttype)) = Position(3) + @eval position(type::Type{<:$wrapper}, ::typeof(parenttype)) = Position(3) end for wrapper in [:PermutedDimsArray] - @eval position(type::Type{<:$wrapper}, ::typeof(parenttype)) = Position(5) + @eval position(type::Type{<:$wrapper}, ::typeof(parenttype)) = Position(5) end diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/base/linearalgebra.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/base/linearalgebra.jl index 3d818cab16..a7dfde2d35 100755 --- a/NDTensors/src/vendored/TypeParameterAccessors/src/base/linearalgebra.jl +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/base/linearalgebra.jl @@ -1,25 +1,25 @@ using LinearAlgebra: - Adjoint, - Diagonal, - Hermitian, - LowerTriangular, - Symmetric, - Transpose, - UnitLowerTriangular, - UnitUpperTriangular, - UpperTriangular + Adjoint, + Diagonal, + Hermitian, + LowerTriangular, + Symmetric, + Transpose, + UnitLowerTriangular, + UnitUpperTriangular, + UpperTriangular for wrapper in [ - :Transpose, - :Adjoint, - :Symmetric, - :Hermitian, - :UpperTriangular, - :LowerTriangular, - :UnitUpperTriangular, - :UnitLowerTriangular, - :Diagonal, -] - @eval position(::Type{<:$wrapper}, ::typeof(eltype)) = Position(1) - @eval position(::Type{<:$wrapper}, ::typeof(parenttype)) = Position(2) + :Transpose, + :Adjoint, + :Symmetric, + :Hermitian, + :UpperTriangular, + :LowerTriangular, + :UnitUpperTriangular, + :UnitLowerTriangular, + :Diagonal, + ] + @eval position(::Type{<:$wrapper}, ::typeof(eltype)) = Position(1) + @eval position(::Type{<:$wrapper}, ::typeof(parenttype)) = Position(2) end diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/base/similartype.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/base/similartype.jl index 2d8175b012..082c269bbb 100755 --- a/NDTensors/src/vendored/TypeParameterAccessors/src/base/similartype.jl +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/base/similartype.jl @@ -5,75 +5,75 @@ like `OffsetArrays` or named indices (such as ITensors). """ function set_indstype(arraytype::Type{<:AbstractArray}, dims::Tuple) - return set_ndims(arraytype, NDims(length(dims))) + return set_ndims(arraytype, NDims(length(dims))) end function similartype(arraytype::Type{<:AbstractArray}, eltype::Type, ndims::NDims) - return similartype(similartype(arraytype, eltype), ndims) + return similartype(similartype(arraytype, eltype), ndims) end function similartype(arraytype::Type{<:AbstractArray}, eltype::Type, dims::Tuple) - return similartype(similartype(arraytype, eltype), dims) + return similartype(similartype(arraytype, eltype), dims) end @traitfn function similartype( - arraytype::Type{ArrayT} -) where {ArrayT;!IsWrappedArray{ArrayT}} - return arraytype + arraytype::Type{ArrayT} + ) where {{ArrayT; !IsWrappedArray{ArrayT}}} + return arraytype end @traitfn function similartype( - arraytype::Type{ArrayT}, eltype::Type -) where {ArrayT;!IsWrappedArray{ArrayT}} - return set_eltype(arraytype, eltype) + arraytype::Type{ArrayT}, eltype::Type + ) where {{ArrayT; !IsWrappedArray{ArrayT}}} + return set_eltype(arraytype, eltype) end function similartype(arraytype::Type{<:BitArray}, elt::Type) - if is_parameter_specified(arraytype, ndims) - return Array{elt,ndims(arraytype)} - end - return Array{elt} + if is_parameter_specified(arraytype, ndims) + return Array{elt, ndims(arraytype)} + end + return Array{elt} end @traitfn function similartype( - arraytype::Type{ArrayT}, dims::Tuple -) where {ArrayT;!IsWrappedArray{ArrayT}} - return set_indstype(arraytype, dims) + arraytype::Type{ArrayT}, dims::Tuple + ) where {{ArrayT; !IsWrappedArray{ArrayT}}} + return set_indstype(arraytype, dims) end @traitfn function similartype( - arraytype::Type{ArrayT}, ndims::NDims -) where {ArrayT;!IsWrappedArray{ArrayT}} - return set_ndims(arraytype, ndims) + arraytype::Type{ArrayT}, ndims::NDims + ) where {{ArrayT; !IsWrappedArray{ArrayT}}} + return set_ndims(arraytype, ndims) end function similartype( - arraytype::Type{<:AbstractArray}, dim1::Base.DimOrInd, dim_rest::Base.DimOrInd... -) - return similartype(arraytype, (dim1, dim_rest...)) + arraytype::Type{<:AbstractArray}, dim1::Base.DimOrInd, dim_rest::Base.DimOrInd... + ) + return similartype(arraytype, (dim1, dim_rest...)) end ## Wrapped arrays -@traitfn function similartype(arraytype::Type{ArrayT}) where {ArrayT;IsWrappedArray{ArrayT}} - return similartype(unwrap_array_type(arraytype), NDims(arraytype)) +@traitfn function similartype(arraytype::Type{ArrayT}) where {{ArrayT; IsWrappedArray{ArrayT}}} + return similartype(unwrap_array_type(arraytype), NDims(arraytype)) end @traitfn function similartype( - arraytype::Type{ArrayT}, eltype::Type -) where {ArrayT;IsWrappedArray{ArrayT}} - return similartype(unwrap_array_type(arraytype), eltype, NDims(arraytype)) + arraytype::Type{ArrayT}, eltype::Type + ) where {{ArrayT; IsWrappedArray{ArrayT}}} + return similartype(unwrap_array_type(arraytype), eltype, NDims(arraytype)) end @traitfn function similartype( - arraytype::Type{ArrayT}, dims::Tuple -) where {ArrayT;IsWrappedArray{ArrayT}} - return similartype(unwrap_array_type(arraytype), dims) + arraytype::Type{ArrayT}, dims::Tuple + ) where {{ArrayT; IsWrappedArray{ArrayT}}} + return similartype(unwrap_array_type(arraytype), dims) end @traitfn function similartype( - arraytype::Type{ArrayT}, ndims::NDims -) where {ArrayT;IsWrappedArray{ArrayT}} - return similartype(unwrap_array_type(arraytype), ndims) + arraytype::Type{ArrayT}, ndims::NDims + ) where {{ArrayT; IsWrappedArray{ArrayT}}} + return similartype(unwrap_array_type(arraytype), ndims) end # This is for uniform `Diag` storage which uses @@ -82,12 +82,12 @@ end # `FillArray` instead. This is a stand-in # to make things work with the current design. function similartype(numbertype::Type{<:Number}) - return numbertype + return numbertype end # Instances function similartype(array::AbstractArray, eltype::Type, dims...) - return similartype(typeof(array), eltype, dims...) + return similartype(typeof(array), eltype, dims...) end similartype(array::AbstractArray, eltype::Type) = similartype(typeof(array), eltype) similartype(array::AbstractArray, dims...) = similartype(typeof(array), dims...) diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/type_parameters.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/type_parameters.jl index 1dd31ce918..301ad20abf 100755 --- a/NDTensors/src/vendored/TypeParameterAccessors/src/type_parameters.jl +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/type_parameters.jl @@ -10,7 +10,7 @@ Position(pos::Int) = Position{pos}() Base.Int(pos::Position) = Int(typeof(pos)) Base.Int(::Type{Position{P}}) where {P} = Int(P) Base.to_index(pos::Position) = Base.to_index(typeof(pos)) -Base.to_index(::Type{P}) where {P<:Position} = Int(P) +Base.to_index(::Type{P}) where {P <: Position} = Int(P) """ position(type::Type, position_name)::Position @@ -29,63 +29,63 @@ position(::Type, pos::Int) = Position(pos) position(::Type, pos::Position) = pos function position(type::Type, name) - type′ = unspecify_type_parameters(type) - if type === type′ - # Fallback definition that determines the - # position automatically from the supertype of - # the type. - return position_from_supertype(type′, name) - end - return position(type′, name) + type′ = unspecify_type_parameters(type) + if type === type′ + # Fallback definition that determines the + # position automatically from the supertype of + # the type. + return position_from_supertype(type′, name) + end + return position(type′, name) end # Automatically determine the position of a type parameter of a type given # a supertype and the name of the parameter. function position_from_supertype(type::Type, name) - type′ = unspecify_type_parameters(type) - supertype_pos = position(supertype(type′), name) - return position_from_supertype_position(type′, supertype_pos) + type′ = unspecify_type_parameters(type) + supertype_pos = position(supertype(type′), name) + return position_from_supertype_position(type′, supertype_pos) end # Automatically determine the position of a type parameter of a type given # the supertype and the position of the corresponding parameter in the supertype. @generated function position_from_supertype_position( - ::Type{T}, supertype_pos::Position -) where {T} - T′ = unspecify_type_parameters(T) - # The type parameters of the type as TypeVars. - # TODO: Ideally we would use `get_type_parameters` - # but that sometimes loses TypeVar names: - # https://github.com/ITensor/TypeParameterAccessors.jl/issues/30 - type_params = Base.unwrap_unionall(T′).parameters - # The type parameters of the immediate supertype as TypeVars. - # This has TypeVars with names that correspond to the names of - # the TypeVars of the type parameters of `T`, for example: - # ```julia - # julia> struct MyArray{B,A} <: AbstractArray{A,B} end - # - # julia> Base.unwrap_unionall(MyArray).parameters - # svec(B, A) - # - # julia> Base.unwrap_unionall(supertype(MyArray)).parameters - # svec(A, B) - # ``` - supertype_params = Base.unwrap_unionall(supertype(T′)).parameters - supertype_param = supertype_params[Int(supertype_pos)] - if !(supertype_param isa TypeVar) - error("Position not found.") - end - pos = findfirst(param -> (param.name == supertype_param.name), type_params) - if isnothing(pos) - return error("Position not found.") - end - return :(@inline; $(Position(pos))) + ::Type{T}, supertype_pos::Position + ) where {T} + T′ = unspecify_type_parameters(T) + # The type parameters of the type as TypeVars. + # TODO: Ideally we would use `get_type_parameters` + # but that sometimes loses TypeVar names: + # https://github.com/ITensor/TypeParameterAccessors.jl/issues/30 + type_params = Base.unwrap_unionall(T′).parameters + # The type parameters of the immediate supertype as TypeVars. + # This has TypeVars with names that correspond to the names of + # the TypeVars of the type parameters of `T`, for example: + # ```julia + # julia> struct MyArray{B,A} <: AbstractArray{A,B} end + # + # julia> Base.unwrap_unionall(MyArray).parameters + # svec(B, A) + # + # julia> Base.unwrap_unionall(supertype(MyArray)).parameters + # svec(A, B) + # ``` + supertype_params = Base.unwrap_unionall(supertype(T′)).parameters + supertype_param = supertype_params[Int(supertype_pos)] + if !(supertype_param isa TypeVar) + error("Position not found.") + end + pos = findfirst(param -> (param.name == supertype_param.name), type_params) + if isnothing(pos) + return error("Position not found.") + end + return :(@inline; $(Position(pos))) end function positions(::Type{T}, pos::Tuple) where {T} - return ntuple(length(pos)) do i - return position(T, pos[i]) - end + return ntuple(length(pos)) do i + return position(T, pos[i]) + end end """ @@ -104,13 +104,13 @@ function get_type_parameters end # Attempts with `Base.@constprop :aggressive` failed, so generated function instead # @inline type_parameters(::Type{T}) where {T} = Tuple(Base.unwrap_unionall(T).parameters) @generated function get_type_parameters(::Type{T}) where {T} - params = wrap_symbol_quotenode.(Tuple(Base.unwrap_unionall(T).parameters)) - return :(@inline; ($(params...),)) + params = wrap_symbol_quotenode.(Tuple(Base.unwrap_unionall(T).parameters)) + return :(@inline; ($(params...),)) end @inline get_type_parameters(::Type{T}, pos) where {T} = get_type_parameters( - T, position(T, pos) + T, position(T, pos) ) -@inline get_type_parameters(::Type{T}, ::Position{p}) where {T,p} = get_type_parameters(T)[p] +@inline get_type_parameters(::Type{T}, ::Position{p}) where {T, p} = get_type_parameters(T)[p] @inline get_type_parameters(::Type{T}, ::Position{0}) where {T} = T @inline get_type_parameters(::Type{T}, pos::Tuple) where {T} = get_type_parameters.(T, pos) @inline get_type_parameters(object, pos) = get_type_parameters(typeof(object), pos) @@ -129,15 +129,15 @@ see [`get_type_parameters`](@ref). function type_parameters end function type_parameters(::Type{T}) where {T} - params = get_type_parameters(T) - any(param -> param isa TypeVar, params) && - return error("One or more parameter is not specified.") - return params + params = get_type_parameters(T) + any(param -> param isa TypeVar, params) && + return error("One or more parameter is not specified.") + return params end @inline function type_parameters(::Type{T}, pos) where {T} - param = get_type_parameters(T, pos) - param isa TypeVar && return error("The parameter is not specified.") - return param + param = get_type_parameters(T, pos) + param isa TypeVar && return error("The parameter is not specified.") + return param end @inline type_parameters(::Type{T}, pos::Tuple) where {T} = type_parameters.(T, pos) @inline type_parameters(object, pos) = type_parameters(typeof(object), pos) @@ -157,7 +157,7 @@ nparameters(::Type{T}) where {T} = length(get_type_parameters(T)) Return whether or not the type parameter at a given position is considered specified. """ function is_parameter_specified(::Type{T}, pos) where {T} - return !(get_type_parameters(T, pos) isa TypeVar) + return !(get_type_parameters(T, pos) isa TypeVar) end """ @@ -168,18 +168,18 @@ Return a new type where the type parameters at the given positions are unset. """ unspecify_type_parameters(::Type{T}) where {T} = Base.typename(T).wrapper function unspecify_type_parameters(::Type{T}, pos::Tuple) where {T} - @inline - return unspecify_type_parameters(T, positions(T, pos)) + @inline + return unspecify_type_parameters(T, positions(T, pos)) end @generated function unspecify_type_parameters( - ::Type{T}, positions::Tuple{Vararg{Position}} -) where {T} - allparams = collect(Any, get_type_parameters(T)) - for pos in type_parameters(positions) - allparams[pos] = get_type_parameters(unspecify_type_parameters(T), Int(pos)) - end - type_expr = construct_type_expr(T, allparams) - return :(@inline; $type_expr) + ::Type{T}, positions::Tuple{Vararg{Position}} + ) where {T} + allparams = collect(Any, get_type_parameters(T)) + for pos in type_parameters(positions) + allparams[pos] = get_type_parameters(unspecify_type_parameters(T), Int(pos)) + end + type_expr = construct_type_expr(T, allparams) + return :(@inline; $type_expr) end unspecify_type_parameters(::Type{T}, pos) where {T} = unspecify_type_parameters(T, (pos,)) @@ -190,23 +190,23 @@ unspecify_type_parameters(::Type{T}, pos) where {T} = unspecify_type_parameters( Return a new type where the type parameters at the given positions are set to the provided values. """ function set_type_parameters( - ::Type{T}, pos::Tuple{Vararg{Any,N}}, parameters::Tuple{Vararg{Any,N}} -) where {T,N} - return set_type_parameters(T, positions(T, pos), parameters) + ::Type{T}, pos::Tuple{Vararg{Any, N}}, parameters::Tuple{Vararg{Any, N}} + ) where {T, N} + return set_type_parameters(T, positions(T, pos), parameters) end @generated function set_type_parameters( - ::Type{T}, positions::Tuple{Vararg{Position,N}}, params::Tuple{Vararg{Any,N}} -) where {T,N} - # collect parameters and change - allparams = collect(Any, get_type_parameters(T)) - for (i, pos) in enumerate(get_type_parameters(positions)) - allparams[pos] = :(params[$i]) - end - type_expr = construct_type_expr(T, allparams) - return :(@inline; $type_expr) + ::Type{T}, positions::Tuple{Vararg{Position, N}}, params::Tuple{Vararg{Any, N}} + ) where {T, N} + # collect parameters and change + allparams = collect(Any, get_type_parameters(T)) + for (i, pos) in enumerate(get_type_parameters(positions)) + allparams[pos] = :(params[$i]) + end + type_expr = construct_type_expr(T, allparams) + return :(@inline; $type_expr) end function set_type_parameters(::Type{T}, pos, param) where {T} - return set_type_parameters(T, (pos,), (param,)) + return set_type_parameters(T, (pos,), (param,)) end """ @@ -217,28 +217,28 @@ Return a new type where the type parameters at the given positions are set to th only if they were previously unspecified. """ function specify_type_parameters( - ::Type{T}, pos::Tuple{Vararg{Any,N}}, parameters::Tuple{Vararg{Any,N}} -) where {T,N} - return specify_type_parameters(T, positions(T, pos), parameters) + ::Type{T}, pos::Tuple{Vararg{Any, N}}, parameters::Tuple{Vararg{Any, N}} + ) where {T, N} + return specify_type_parameters(T, positions(T, pos), parameters) end function specify_type_parameters(::Type{T}, parameters::Tuple) where {T} - return specify_type_parameters(T, ntuple(identity, nparameters(T)), parameters) + return specify_type_parameters(T, ntuple(identity, nparameters(T)), parameters) end @generated function specify_type_parameters( - ::Type{T}, positions::Tuple{Vararg{Position,N}}, params::Tuple{Vararg{Any,N}} -) where {T,N} - # collect parameters and change unspecified - allparams = collect(Any, get_type_parameters(T)) - for (i, pos) in enumerate(get_type_parameters(positions)) - if !is_parameter_specified(T, pos()) - allparams[pos] = :(params[$i]) + ::Type{T}, positions::Tuple{Vararg{Position, N}}, params::Tuple{Vararg{Any, N}} + ) where {T, N} + # collect parameters and change unspecified + allparams = collect(Any, get_type_parameters(T)) + for (i, pos) in enumerate(get_type_parameters(positions)) + if !is_parameter_specified(T, pos()) + allparams[pos] = :(params[$i]) + end end - end - type_expr = construct_type_expr(T, allparams) - return :(@inline; $type_expr) + type_expr = construct_type_expr(T, allparams) + return :(@inline; $type_expr) end function specify_type_parameters(::Type{T}, pos, param) where {T} - return specify_type_parameters(T, (pos,), (param,)) + return specify_type_parameters(T, (pos,), (param,)) end """ @@ -251,24 +251,24 @@ This function should output a Tuple of the default values, with exactly one for each type parameter slot of the type. """ function default_type_parameters(::Type{T}, pos) where {T} - return default_type_parameters(T, position(T, pos)) + return default_type_parameters(T, position(T, pos)) end -function default_type_parameters(::Type{T}, ::Position{pos}) where {T,pos} - param = default_type_parameters(T)[pos] - if param isa UndefinedDefaultTypeParameter - return error("No default parameter defined at this position.") - end - return param +function default_type_parameters(::Type{T}, ::Position{pos}) where {T, pos} + param = default_type_parameters(T)[pos] + if param isa UndefinedDefaultTypeParameter + return error("No default parameter defined at this position.") + end + return param end default_type_parameters(::Type{T}, pos::Tuple) where {T} = default_type_parameters.(T, pos) default_type_parameters(t, pos) = default_type_parameters(typeof(t), pos) default_type_parameters(t) = default_type_parameters(typeof(t)) function default_type_parameters(type::Type) - type′ = unspecify_type_parameters(type) - if type === type′ - return default_type_parameters_from_supertype(type′) - end - return default_type_parameters(type′) + type′ = unspecify_type_parameters(type) + if type === type′ + return default_type_parameters_from_supertype(type′) + end + return default_type_parameters(type′) end struct UndefinedDefaultTypeParameter end @@ -281,20 +281,20 @@ struct UndefinedDefaultTypeParameter end # `UndefinedDefaultTypeParameter()`. Accessing those type parameters # by name/position will throw an error. @generated function default_type_parameters_from_supertype(::Type{T}) where {T} - T′ = unspecify_type_parameters(T) - supertype_default_type_params = default_type_parameters(supertype(T′)) - type_params = Base.unwrap_unionall(T′).parameters - supertype_params = Base.unwrap_unionall(supertype(T′)).parameters - defaults = Any[UndefinedDefaultTypeParameter() for _ in 1:nparameters(T′)] - for (supertype_param, supertype_default_type_param) in - zip(supertype_params, supertype_default_type_params) - if !(supertype_param isa TypeVar) - continue + T′ = unspecify_type_parameters(T) + supertype_default_type_params = default_type_parameters(supertype(T′)) + type_params = Base.unwrap_unionall(T′).parameters + supertype_params = Base.unwrap_unionall(supertype(T′)).parameters + defaults = Any[UndefinedDefaultTypeParameter() for _ in 1:nparameters(T′)] + for (supertype_param, supertype_default_type_param) in + zip(supertype_params, supertype_default_type_params) + if !(supertype_param isa TypeVar) + continue + end + param_position = findfirst(param -> (param.name == supertype_param.name), type_params) + defaults[param_position] = supertype_default_type_param end - param_position = findfirst(param -> (param.name == supertype_param.name), type_params) - defaults[param_position] = supertype_default_type_param - end - return :(@inline; $(Tuple(defaults))) + return :(@inline; $(Tuple(defaults))) end """ @@ -304,13 +304,13 @@ end Set the type parameters at the given positions to their default values. """ function set_default_type_parameters(::Type{T}, pos::Tuple) where {T} - return set_type_parameters(T, pos, default_type_parameters.(T, pos)) + return set_type_parameters(T, pos, default_type_parameters.(T, pos)) end function set_default_type_parameters(::Type{T}) where {T} - return set_default_type_parameters(T, ntuple(identity, nparameters(T))) + return set_default_type_parameters(T, ntuple(identity, nparameters(T))) end function set_default_type_parameters(::Type{T}, pos) where {T} - return set_default_type_parameters(T, (pos,)) + return set_default_type_parameters(T, (pos,)) end """ @@ -321,11 +321,11 @@ Set the type parameters at the given positions to their default values, if they had not been specified. """ function specify_default_type_parameters(::Type{T}, pos::Tuple) where {T} - return specify_type_parameters(T, pos, default_type_parameters.(T, pos)) + return specify_type_parameters(T, pos, default_type_parameters.(T, pos)) end function specify_default_type_parameters(::Type{T}) where {T} - return specify_default_type_parameters(T, ntuple(identity, nparameters(T))) + return specify_default_type_parameters(T, ntuple(identity, nparameters(T))) end function specify_default_type_parameters(::Type{T}, pos) where {T} - return specify_default_type_parameters(T, (pos,)) + return specify_default_type_parameters(T, (pos,)) end diff --git a/NDTensors/src/vendored/TypeParameterAccessors/src/type_utils.jl b/NDTensors/src/vendored/TypeParameterAccessors/src/type_utils.jl index 79a316dae7..a0779d162d 100755 --- a/NDTensors/src/vendored/TypeParameterAccessors/src/type_utils.jl +++ b/NDTensors/src/vendored/TypeParameterAccessors/src/type_utils.jl @@ -3,12 +3,12 @@ wrap_symbol_quotenode(param) = param isa Symbol ? QuoteNode(param) : param "Construct the expression for qualifying a type with given parameters" function construct_type_expr(type, parameters) - basetype = unspecify_type_parameters(type) - type_expr = Expr(:curly, basetype, wrap_symbol_quotenode.(parameters)...) - for parameter in reverse(parameters) - if parameter isa TypeVar - type_expr = Expr(:call, :UnionAll, parameter, type_expr) + basetype = unspecify_type_parameters(type) + type_expr = Expr(:curly, basetype, wrap_symbol_quotenode.(parameters)...) + for parameter in reverse(parameters) + if parameter isa TypeVar + type_expr = Expr(:call, :UnionAll, parameter, type_expr) + end end - end - return type_expr + return type_expr end diff --git a/NDTensors/vendor/Project.toml b/NDTensors/vendor/Project.toml index 4c7dfde138..10433d13c5 100644 --- a/NDTensors/vendor/Project.toml +++ b/NDTensors/vendor/Project.toml @@ -1,5 +1,7 @@ [deps] PackageAnalyzer = "e713c705-17e4-4cec-abe0-95bf5bf3e10c" +Runic = "62bfec6d-59d7-401d-8490-b29ee721c001" [compat] PackageAnalyzer = "3.1" +Runic = "1.5.1" diff --git a/NDTensors/vendor/run.jl b/NDTensors/vendor/run.jl index 0544fcddcc..5bf570c032 100644 --- a/NDTensors/vendor/run.jl +++ b/NDTensors/vendor/run.jl @@ -1,4 +1,5 @@ using PackageAnalyzer: PackageAnalyzer +using Runic: Runic deps = [("TypeParameterAccessors", v"0.3.10")] @@ -15,4 +16,5 @@ for (name, version) in deps cp(joinpath(local_path, "src"), joinpath(p, "src"); force = true) # Make files writable so we can format. chmod(p, 0o755; recursive = true) + Runic.main(["--inplace", p]) end