From 6438473f208ed8d568af2d155b63aa507838f6a6 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 13 Sep 2023 17:11:49 -0400 Subject: [PATCH 001/133] Create generic randn for CUDA and use randn to make on device --- NDTensors/ext/NDTensorsCUDAExt/fill.jl | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 NDTensors/ext/NDTensorsCUDAExt/fill.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/fill.jl b/NDTensors/ext/NDTensorsCUDAExt/fill.jl new file mode 100644 index 0000000000..a38e50b339 --- /dev/null +++ b/NDTensors/ext/NDTensorsCUDAExt/fill.jl @@ -0,0 +1,6 @@ +function NDTensors.generic_randn(arraytype::Type{<:CuArray}, dim::Integer=0) + arraytype_specified = NDTensors.SetParameters.set_unspecified_parameters( + NDTensors.leaf_parenttype(arraytype), NDTensors.SetParameters.DefaultParameters() + ) + return CUDA.randn(eltype(arraytype_specified), dim) +end \ No newline at end of file From 682bcf3c83f40701a1e9b6204e4c635cbbd5e6c7 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 13 Sep 2023 17:17:37 -0400 Subject: [PATCH 002/133] Make generic_zero for CUDA --- NDTensors/ext/NDTensorsCUDAExt/fill.jl | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/fill.jl b/NDTensors/ext/NDTensorsCUDAExt/fill.jl index a38e50b339..1a586ff220 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/fill.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/fill.jl @@ -1,6 +1,16 @@ +import NDTensors.SetParameters: set_unspecified_parameters, DefaultParameters +import NDTensors: leaf_parenttype + function NDTensors.generic_randn(arraytype::Type{<:CuArray}, dim::Integer=0) - arraytype_specified = NDTensors.SetParameters.set_unspecified_parameters( - NDTensors.leaf_parenttype(arraytype), NDTensors.SetParameters.DefaultParameters() + arraytype_specified = set_unspecified_parameters( + leaf_parenttype(arraytype), DefaultParameters() ) return CUDA.randn(eltype(arraytype_specified), dim) +end + +function NDTensors.generic_zeros(arraytype::Type{<:CuArray}, dim::Integer=0) + arraytype_specified = NDTensors.set_unspecified_parameters( + leaf_parenttype(arraytype), DefaultParameters() + ) + return CUDA.zeros(eltype(arraytype_specified), dim) end \ No newline at end of file From 4f879d5becd685de4dd3b2fa74b1d2632a1f57db Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 13 Sep 2023 17:22:36 -0400 Subject: [PATCH 003/133] Format --- NDTensors/ext/NDTensorsCUDAExt/fill.jl | 2 +- monorepo/Project.toml | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) create mode 100644 monorepo/Project.toml diff --git a/NDTensors/ext/NDTensorsCUDAExt/fill.jl b/NDTensors/ext/NDTensorsCUDAExt/fill.jl index 1a586ff220..984a7f8532 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/fill.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/fill.jl @@ -13,4 +13,4 @@ function NDTensors.generic_zeros(arraytype::Type{<:CuArray}, dim::Integer=0) leaf_parenttype(arraytype), DefaultParameters() ) return CUDA.zeros(eltype(arraytype_specified), dim) -end \ No newline at end of file +end diff --git a/monorepo/Project.toml b/monorepo/Project.toml new file mode 100644 index 0000000000..53e60981d4 --- /dev/null +++ b/monorepo/Project.toml @@ -0,0 +1,9 @@ +[deps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" +MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" +NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From b567ecd30f9523b945545b48373f018a0711b8ba Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 13 Sep 2023 17:29:03 -0400 Subject: [PATCH 004/133] remove monorepo --- monorepo/Project.toml | 9 --------- 1 file changed, 9 deletions(-) delete mode 100644 monorepo/Project.toml diff --git a/monorepo/Project.toml b/monorepo/Project.toml deleted file mode 100644 index 53e60981d4..0000000000 --- a/monorepo/Project.toml +++ /dev/null @@ -1,9 +0,0 @@ -[deps] -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" -Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" -ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" -JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" -MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" -NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf" -Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From a4baeb75bf4537248861678e41fa8ecefb9b9866 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 14 Sep 2023 10:35:44 -0400 Subject: [PATCH 005/133] Convert tests to CPU to not perform scalar operations --- NDTensors/test/combiner.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/NDTensors/test/combiner.jl b/NDTensors/test/combiner.jl index f8853ff8ff..b60c05bbd3 100644 --- a/NDTensors/test/combiner.jl +++ b/NDTensors/test/combiner.jl @@ -29,12 +29,12 @@ end @test output_tensor isa DenseTensor @test dims(output_tensor) == output_tensor_inds for i in 1:length(input_tensor) - @test input_tensor[i] == output_tensor[i] + @test NDTensors.cpu(input_tensor)[i] == NDTensors.cpu(output_tensor)[i] end # Test uncombining new_input_tensor = contract(output_tensor, (1, -1), combiner_tensor, (-1, 2, 3)) - @test new_input_tensor == input_tensor + @test NDTensors.cpu(new_input_tensor) == NDTensors.cpu(input_tensor) # Catch invalid combining input_tensor_inds = (d,) From 03a9b2764d18f5b505c75e35d610267664dfd7d6 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 14 Sep 2023 11:12:30 -0400 Subject: [PATCH 006/133] import -> using. Can use NDTenosrs randn! instead of another function --- NDTensors/ext/NDTensorsCUDAExt/fill.jl | 16 ++++++++-------- NDTensors/src/abstractarray/fill.jl | 8 ++------ 2 files changed, 10 insertions(+), 14 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/fill.jl b/NDTensors/ext/NDTensorsCUDAExt/fill.jl index 984a7f8532..dafa2fe694 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/fill.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/fill.jl @@ -1,12 +1,12 @@ -import NDTensors.SetParameters: set_unspecified_parameters, DefaultParameters -import NDTensors: leaf_parenttype +using NDTensors.SetParameters: set_unspecified_parameters, DefaultParameters +using NDTensors: leaf_parenttype -function NDTensors.generic_randn(arraytype::Type{<:CuArray}, dim::Integer=0) - arraytype_specified = set_unspecified_parameters( - leaf_parenttype(arraytype), DefaultParameters() - ) - return CUDA.randn(eltype(arraytype_specified), dim) -end +# function NDTensors.generic_randn(arraytype::Type{<:CuArray}, dim::Integer=0) +# arraytype_specified = set_unspecified_parameters( +# leaf_parenttype(arraytype), DefaultParameters() +# ) +# return CUDA.randn(eltype(arraytype_specified), dim) +# end function NDTensors.generic_zeros(arraytype::Type{<:CuArray}, dim::Integer=0) arraytype_specified = NDTensors.set_unspecified_parameters( diff --git a/NDTensors/src/abstractarray/fill.jl b/NDTensors/src/abstractarray/fill.jl index 85e8377d3a..439a115dd8 100644 --- a/NDTensors/src/abstractarray/fill.jl +++ b/NDTensors/src/abstractarray/fill.jl @@ -1,13 +1,9 @@ -function generic_randn(arraytype::Type{<:AbstractArray}, dim::Integer=0) +function generic_randn(arraytype::Type{<:AbstractArray}, dim::Integer=0; rng = Random.default_rng()) arraytype_specified = set_unspecified_parameters( leaf_parenttype(arraytype), DefaultParameters() ) data = similar(arraytype_specified, dim) - ElT = eltype(data) - for i in 1:length(data) - data[i] = randn(ElT) - end - return data + NDTensors.randn!(rng, data) end function generic_zeros(arraytype::Type{<:AbstractArray}, dim::Integer=0) From 154fac219b9a6c581fa99865e6e2815ab3754aac Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 21 Sep 2023 12:43:09 -0400 Subject: [PATCH 007/133] Temporarily get contract working by making a `mul!!` function --- NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl | 8 ++++++++ NDTensors/src/dense/tensoralgebra/contract.jl | 6 +++++- 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl index c144a4a802..7a363786ce 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl @@ -11,3 +11,11 @@ function NDTensors.svd_catch_error(A::CuMatrix; alg="JacobiAlgorithm") end return USV end + +using LinearAlgebra +function NDTensors.mul!!(CM::CuArray, AM::AbstractArray, BM::AbstractArray, α, β) + AM = AM isa CuArray ? AM : copy(AM) + BM = BM isa CuArray ? BM : copy(BM) + @show typeof(AM) + return mul!(CM, AM, BM, α, β) +end diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index b5faa91616..63a19bc098 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -393,7 +393,7 @@ function _contract!( #tC = similar(CM) #_gemm!(tA, tB, El(α), AM, BM, El(β), CM) - @strided mul!(CM, AM, BM, El(α), El(β)) + mul!!(CM, AM, BM, El(α), El(β)) if props.permuteC Cr = reshape(CM, props.newCrange) @@ -405,3 +405,7 @@ function _contract!( return CT end + +function mul!!(CM::AbstractArray, AM::AbstractArray, BM::AbstractArray, α, β) + @strided mul!(CM, AM, BM, α, β) +end From 3da5e96a6947d60ddd907a61e319f8d16c009578 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 21 Sep 2023 12:43:17 -0400 Subject: [PATCH 008/133] format --- NDTensors/src/abstractarray/fill.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/NDTensors/src/abstractarray/fill.jl b/NDTensors/src/abstractarray/fill.jl index 439a115dd8..016a839060 100644 --- a/NDTensors/src/abstractarray/fill.jl +++ b/NDTensors/src/abstractarray/fill.jl @@ -1,9 +1,11 @@ -function generic_randn(arraytype::Type{<:AbstractArray}, dim::Integer=0; rng = Random.default_rng()) +function generic_randn( + arraytype::Type{<:AbstractArray}, dim::Integer=0; rng=Random.default_rng() +) arraytype_specified = set_unspecified_parameters( leaf_parenttype(arraytype), DefaultParameters() ) data = similar(arraytype_specified, dim) - NDTensors.randn!(rng, data) + return NDTensors.randn!(rng, data) end function generic_zeros(arraytype::Type{<:AbstractArray}, dim::Integer=0) From cb36f23a9bc79a3d1288667d6e35e3465c422247 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 21 Sep 2023 12:47:27 -0400 Subject: [PATCH 009/133] remove change to combiner --- NDTensors/test/combiner.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/NDTensors/test/combiner.jl b/NDTensors/test/combiner.jl index b60c05bbd3..f8853ff8ff 100644 --- a/NDTensors/test/combiner.jl +++ b/NDTensors/test/combiner.jl @@ -29,12 +29,12 @@ end @test output_tensor isa DenseTensor @test dims(output_tensor) == output_tensor_inds for i in 1:length(input_tensor) - @test NDTensors.cpu(input_tensor)[i] == NDTensors.cpu(output_tensor)[i] + @test input_tensor[i] == output_tensor[i] end # Test uncombining new_input_tensor = contract(output_tensor, (1, -1), combiner_tensor, (-1, 2, 3)) - @test NDTensors.cpu(new_input_tensor) == NDTensors.cpu(input_tensor) + @test new_input_tensor == input_tensor # Catch invalid combining input_tensor_inds = (d,) From aef93a59967e75e775d39f37f8d152fb64a76907 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 21 Sep 2023 14:29:04 -0400 Subject: [PATCH 010/133] Use adapt to not copy elements --- NDTensors/ext/NDTensorsCUDAExt/adapt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/adapt.jl b/NDTensors/ext/NDTensorsCUDAExt/adapt.jl index 718901da4c..abc0c295c8 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/adapt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/adapt.jl @@ -16,7 +16,7 @@ buffertype(::NDTensorCuArrayAdaptor{B}) where {B} = B function Adapt.adapt_storage(adaptor::NDTensorCuArrayAdaptor, xs::AbstractArray) ElT = eltype(xs) BufT = buffertype(adaptor) - return isbits(xs) ? xs : CuArray{ElT,1,BufT}(xs) + return isbits(xs) ? xs : adapt(CuArray{ElT, 1, BufT},xs) end function NDTensors.adapt_storagetype( From 89d28375cf403dcd369fca2ed4000d08deb17f1a Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 21 Sep 2023 14:29:15 -0400 Subject: [PATCH 011/133] Remove commented code --- NDTensors/ext/NDTensorsCUDAExt/fill.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/fill.jl b/NDTensors/ext/NDTensorsCUDAExt/fill.jl index dafa2fe694..909cebf468 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/fill.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/fill.jl @@ -1,13 +1,6 @@ using NDTensors.SetParameters: set_unspecified_parameters, DefaultParameters using NDTensors: leaf_parenttype -# function NDTensors.generic_randn(arraytype::Type{<:CuArray}, dim::Integer=0) -# arraytype_specified = set_unspecified_parameters( -# leaf_parenttype(arraytype), DefaultParameters() -# ) -# return CUDA.randn(eltype(arraytype_specified), dim) -# end - function NDTensors.generic_zeros(arraytype::Type{<:CuArray}, dim::Integer=0) arraytype_specified = NDTensors.set_unspecified_parameters( leaf_parenttype(arraytype), DefaultParameters() From 3e93d61d84f404fafe37021a93dcbf9c6355d01f Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 21 Sep 2023 14:34:49 -0400 Subject: [PATCH 012/133] Remove CUDA fill functions and use more generic in abstractarray/fill.jl --- NDTensors/ext/NDTensorsCUDAExt/fill.jl | 9 --------- NDTensors/src/abstractarray/fill.jl | 2 +- 2 files changed, 1 insertion(+), 10 deletions(-) delete mode 100644 NDTensors/ext/NDTensorsCUDAExt/fill.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/fill.jl b/NDTensors/ext/NDTensorsCUDAExt/fill.jl deleted file mode 100644 index 909cebf468..0000000000 --- a/NDTensors/ext/NDTensorsCUDAExt/fill.jl +++ /dev/null @@ -1,9 +0,0 @@ -using NDTensors.SetParameters: set_unspecified_parameters, DefaultParameters -using NDTensors: leaf_parenttype - -function NDTensors.generic_zeros(arraytype::Type{<:CuArray}, dim::Integer=0) - arraytype_specified = NDTensors.set_unspecified_parameters( - leaf_parenttype(arraytype), DefaultParameters() - ) - return CUDA.zeros(eltype(arraytype_specified), dim) -end diff --git a/NDTensors/src/abstractarray/fill.jl b/NDTensors/src/abstractarray/fill.jl index 016a839060..7df74ad231 100644 --- a/NDTensors/src/abstractarray/fill.jl +++ b/NDTensors/src/abstractarray/fill.jl @@ -13,5 +13,5 @@ function generic_zeros(arraytype::Type{<:AbstractArray}, dim::Integer=0) leaf_parenttype(arraytype), DefaultParameters() ) ElT = eltype(arraytype_specified) - return fill!(similar(arraytype_specified, dim), zero(ElT)) + return arraytype_specified(zeros(ElT, dim)) end From 6b42956de81724525f169365f1fff50f85205b78 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 21 Sep 2023 14:35:55 -0400 Subject: [PATCH 013/133] Remove NDTensors. --- NDTensors/src/abstractarray/fill.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/abstractarray/fill.jl b/NDTensors/src/abstractarray/fill.jl index 7df74ad231..4d9bcf4721 100644 --- a/NDTensors/src/abstractarray/fill.jl +++ b/NDTensors/src/abstractarray/fill.jl @@ -5,7 +5,7 @@ function generic_randn( leaf_parenttype(arraytype), DefaultParameters() ) data = similar(arraytype_specified, dim) - return NDTensors.randn!(rng, data) + return randn!(rng, data) end function generic_zeros(arraytype::Type{<:AbstractArray}, dim::Integer=0) From fe470b743260e96e39db712e276aac5abeb0be2f Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 21 Sep 2023 14:37:25 -0400 Subject: [PATCH 014/133] format --- NDTensors/ext/NDTensorsCUDAExt/adapt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/adapt.jl b/NDTensors/ext/NDTensorsCUDAExt/adapt.jl index abc0c295c8..165b5405ce 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/adapt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/adapt.jl @@ -16,7 +16,7 @@ buffertype(::NDTensorCuArrayAdaptor{B}) where {B} = B function Adapt.adapt_storage(adaptor::NDTensorCuArrayAdaptor, xs::AbstractArray) ElT = eltype(xs) BufT = buffertype(adaptor) - return isbits(xs) ? xs : adapt(CuArray{ElT, 1, BufT},xs) + return isbits(xs) ? xs : adapt(CuArray{ElT,1,BufT}, xs) end function NDTensors.adapt_storagetype( From 65d71748eec06b0baae16a528c9b82428a1fe4bc Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 21 Sep 2023 16:22:25 -0400 Subject: [PATCH 015/133] Add comment about gpus --- NDTensors/src/linearalgebra/linearalgebra.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index b681cd4a9e..12062ea481 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -421,6 +421,7 @@ function qr_positive(M::AbstractMatrix) sparseQ, R = qr(M) Q = convert(typeof(R), sparseQ) nc = size(Q, 2) + ## TODO issue here for GPU because tying to access indices for c in 1:nc if R[c, c] != 0.0 #sign(0.0)==0.0 so we don't want to zero out a column of Q. sign_Rc = sign(R[c, c]) From 7232d31acc0361aefe8d0904def3bb12c2b396a7 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 25 Sep 2023 17:30:32 -0400 Subject: [PATCH 016/133] Do elementwise operations on data to avoid scalar indexing --- src/itensor.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/itensor.jl b/src/itensor.jl index 19286f065c..ce2e7fbff9 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1963,8 +1963,9 @@ map(f, x::ITensor) = itensor(map(f, tensor(x))) w .+= a .* v ``` """ -axpy!(a::Number, v::ITensor, w::ITensor) = (w .+= a .* v) - +function axpy!(a::Number, v::ITensor, w::ITensor) + itensor(data(w) .+= a .* data(v), inds(w)) +end """ axpby!(a,v,b,w) @@ -1982,11 +1983,11 @@ Scale the ITensor A by x in-place. May also be written `rmul!`. A .*= x ``` """ -scale!(T::ITensor, α::Number) = (T .*= α) +scale!(T::ITensor, α::Number) = itensor(data(T) .*= α, inds(T)) -rmul!(T::ITensor, α::Number) = (T .*= α) +rmul!(T::ITensor, α::Number) = itensor(data(T) .*= α, inds(T)) -lmul!(T::ITensor, α::Number) = (T .= α .* T) +lmul!(T::ITensor, α::Number) = itensor(data(T) .= α .* data(T), inds(T)) """ mul!(A::ITensor, x::Number, B::ITensor) @@ -1994,10 +1995,9 @@ lmul!(T::ITensor, α::Number) = (T .= α .* T) Scalar multiplication of ITensor B with x, and store the result in A. Like `A .= x .* B`. """ -mul!(R::ITensor, α::Number, T::ITensor) = (R .= α .* T) - -mul!(R::ITensor, T::ITensor, α::Number) = (R .= T .* α) +mul!(R::ITensor, α::Number, T::ITensor) = itensor(data(R) .= α * data(T), inds(T)) +mul!(R::ITensor, T::ITensor, α::Number) = itensor(data(R) .= data(T) * α, inds(T)) ######################### # End ITensor Operations # From 61972137f6faa73e3113c80d4cb755e568f4633b Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 25 Sep 2023 17:30:46 -0400 Subject: [PATCH 017/133] Force dot to return value on CPU --- src/itensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/itensor.jl b/src/itensor.jl index ce2e7fbff9..8dfc810e9e 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1892,7 +1892,7 @@ diag(T::ITensor) = diag(tensor(T)) mul!(C::ITensor, A::ITensor, B::ITensor, args...)::ITensor = contract!(C, A, B, args...) -dot(A::ITensor, B::ITensor) = (dag(A) * B)[] +dot(A::ITensor, B::ITensor) = NDTensors.cpu(dag(A) * B)[] inner(y::ITensor, A::ITensor, x::ITensor) = (dag(y) * A * x)[] inner(y::ITensor, x::ITensor) = (dag(y) * x)[] From 0ed96cb60821a90732c0a331b079c819edee4d13 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 25 Sep 2023 17:33:19 -0400 Subject: [PATCH 018/133] remove copy code --- NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl index 7a363786ce..ec7416c5f5 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl @@ -14,8 +14,5 @@ end using LinearAlgebra function NDTensors.mul!!(CM::CuArray, AM::AbstractArray, BM::AbstractArray, α, β) - AM = AM isa CuArray ? AM : copy(AM) - BM = BM isa CuArray ? BM : copy(BM) - @show typeof(AM) return mul!(CM, AM, BM, α, β) end From fa67f7c4e48336536eb21b8cc847027e0c2553f8 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 25 Sep 2023 17:33:53 -0400 Subject: [PATCH 019/133] format --- src/itensor.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/itensor.jl b/src/itensor.jl index 8dfc810e9e..fee56f5922 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1963,8 +1963,8 @@ map(f, x::ITensor) = itensor(map(f, tensor(x))) w .+= a .* v ``` """ -function axpy!(a::Number, v::ITensor, w::ITensor) - itensor(data(w) .+= a .* data(v), inds(w)) +function axpy!(a::Number, v::ITensor, w::ITensor) + return itensor(data(w) .+= a .* data(v), inds(w)) end """ axpby!(a,v,b,w) From 674c40c6ef7f45c6cf11c6dd8b094b2722ac8da9 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 26 Sep 2023 17:29:07 -0400 Subject: [PATCH 020/133] Bootleg fix for a conversion to UnifiedMemory issue --- NDTensors/src/linearalgebra/linearalgebra.jl | 1 + src/mps/dmrg.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index 12062ea481..8ad98525bc 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -400,6 +400,7 @@ function qx(qx::Function, T::DenseTensor{<:Any,2}; kwargs...) Qinds = IndsT((ind(T, 1), q)) Xinds = IndsT((q, ind(T, 2))) QM = convert(typeof(XM), QM) + QM = convert(typeof(XM), QM) Q = tensor(Dense(vec(QM)), Qinds) #Q was strided X = tensor(Dense(vec(XM)), Xinds) return Q, X diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index a6b78c6f8c..9774c38b82 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -278,7 +278,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) end energy = vals[1] - phi::ITensor = vecs[1] + phi::ITensor = itensor(storagetype(phi)(convert(typeof(data(phi)),data(vecs[1]))), inds(phi)) ortho = ha == 1 ? "left" : "right" From 79ec79f63c04fc735036d0548afe00ffdfcf6a48 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 26 Sep 2023 17:31:32 -0400 Subject: [PATCH 021/133] force itensor to use data to speed up computation --- src/tensor_operations/matrix_decomposition.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/tensor_operations/matrix_decomposition.jl b/src/tensor_operations/matrix_decomposition.jl index 03a622b692..9e93f3f03e 100644 --- a/src/tensor_operations/matrix_decomposition.jl +++ b/src/tensor_operations/matrix_decomposition.jl @@ -573,7 +573,8 @@ function factorize_eigen(A::ITensor, Linds...; kwargs...) # (Lis..., prime(Lis)...) delta_A2 = replaceinds(delta_A2, Lis, dag(simLis)) noprime!(delta_A2) - A2 += delta_A2 + itensor(data(A2) .+= data(delta_A2), inds(A2)) + #A2 += delta_A2 end F = eigen(A2, Lis, simLis; ishermitian=true, kwargs...) D, _, spec = F From 68aae9c45f6266e84f78c95d8ea421899c96d261 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 27 Sep 2023 10:02:28 -0400 Subject: [PATCH 022/133] remove unecessary code --- NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl index 7a363786ce..5796a9b521 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl @@ -12,10 +12,7 @@ function NDTensors.svd_catch_error(A::CuMatrix; alg="JacobiAlgorithm") return USV end -using LinearAlgebra function NDTensors.mul!!(CM::CuArray, AM::AbstractArray, BM::AbstractArray, α, β) - AM = AM isa CuArray ? AM : copy(AM) - BM = BM isa CuArray ? BM : copy(BM) - @show typeof(AM) + return mul!(CM, AM, BM, α, β) end From 3b808e1ca0a45562f085691b75f935c04232bbda Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 27 Sep 2023 10:02:43 -0400 Subject: [PATCH 023/133] Use all of linearalgebra --- NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl index 41f3e6ea57..f6eb791f93 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl @@ -4,7 +4,7 @@ using NDTensors using NDTensors.SetParameters using Adapt using Functors -using LinearAlgebra: BlasFloat +using LinearAlgebra if isdefined(Base, :get_extension) using CUDA From de4b1f9b3ee849baf33deb19945f796d5f8c2fc0 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 27 Sep 2023 16:03:57 -0400 Subject: [PATCH 024/133] grab data on cpu not gpu --- NDTensors/src/dense/tensoralgebra/contract.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 63a19bc098..8c0d4cee77 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -233,9 +233,11 @@ function _contract_scalar_maybe_perm!( β=zero(ElR), ) where {ElR,NR} if nnz(T₁) == 1 - _contract_scalar_maybe_perm!(R, labelsR, T₂, labelsT₂, α * T₁[1], β) + γ = α * NDTensors.cpu(T₁)[1] + _contract_scalar_maybe_perm!(R, labelsR, T₂, labelsT₂, γ, β) elseif nnz(T₂) == 1 - _contract_scalar_maybe_perm!(R, labelsR, T₁, labelsT₁, α * T₂[1], β) + γ = α * NDTensors.cpu(T₂)[1] + _contract_scalar_maybe_perm!(R, labelsR, T₁, labelsT₁, γ, β) else error("In _contract_scalar_perm!, one tensor must be a scalar") end @@ -345,7 +347,8 @@ function _contract!( tA = 'N' if props.permuteA #@timeit_debug timer "_contract!: permutedims A" begin - @strided Ap = permutedims(AT, props.PA) + #@strided + Ap = permutedims(AT, props.PA) #end # @timeit AM = transpose(reshape(Ap, (props.dmid, props.dleft))) else @@ -360,7 +363,8 @@ function _contract!( tB = 'N' if props.permuteB #@timeit_debug timer "_contract!: permutedims B" begin - @strided Bp = permutedims(BT, props.PB) + #@strided + Bp = permutedims(BT, props.PB) #end # @timeit BM = reshape(Bp, (props.dmid, props.dright)) else @@ -399,7 +403,8 @@ function _contract!( Cr = reshape(CM, props.newCrange) # TODO: use invperm(pC) here? #@timeit_debug timer "_contract!: permutedims C" begin - @strided CT .= permutedims(Cr, props.PC) + #@strided + CT .= permutedims(Cr, props.PC) #end # @timeit end From f69389f5e02783d8ade4ae786ef9219c8f1afb2a Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 27 Sep 2023 16:06:09 -0400 Subject: [PATCH 025/133] [experiment] remove strided. Theres an issue with broadcast here during apxy! --- NDTensors/src/dense/densetensor.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NDTensors/src/dense/densetensor.jl b/NDTensors/src/dense/densetensor.jl index 8a019e1f3a..769f47673d 100644 --- a/NDTensors/src/dense/densetensor.jl +++ b/NDTensors/src/dense/densetensor.jl @@ -243,7 +243,8 @@ function permutedims!( end RA = array(R) TA = array(T) - @strided RA .= f.(RA, permutedims(TA, perm)) + #@strided + RA .= f.(RA, permutedims(TA, perm)) return R end From 418043b5be2a3af88daa953e1556004f4cadc462 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 27 Sep 2023 16:06:24 -0400 Subject: [PATCH 026/133] [experiment] Add todo here --- NDTensors/src/linearalgebra/linearalgebra.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index b681cd4a9e..65c98b5437 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -218,6 +218,7 @@ function LinearAlgebra.eigen( use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff) matrixT = matrix(T) + ## TODO this doesn't work for GPU if any(!isfinite, matrixT) throw( ArgumentError( @@ -400,6 +401,7 @@ function qx(qx::Function, T::DenseTensor{<:Any,2}; kwargs...) Qinds = IndsT((ind(T, 1), q)) Xinds = IndsT((q, ind(T, 2))) QM = convert(typeof(XM), QM) + QM = convert(typeof(XM), QM) Q = tensor(Dense(vec(QM)), Qinds) #Q was strided X = tensor(Dense(vec(XM)), Xinds) return Q, X From b876aa1b976262317457aa4ffc5cf712cb951e9d Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 27 Sep 2023 16:07:17 -0400 Subject: [PATCH 027/133] [experiment] Trucate! does many scalar operations. Potentially could convert to CPU --- NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl | 1 + NDTensors/ext/NDTensorsCUDAExt/truncate.jl | 7 +++++++ 2 files changed, 8 insertions(+) create mode 100644 NDTensors/ext/NDTensorsCUDAExt/truncate.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl index f6eb791f93..383cd6b8ad 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl @@ -20,4 +20,5 @@ include("imports.jl") include("set_types.jl") include("adapt.jl") include("linearalgebra.jl") +include("truncate.jl") end diff --git a/NDTensors/ext/NDTensorsCUDAExt/truncate.jl b/NDTensors/ext/NDTensorsCUDAExt/truncate.jl new file mode 100644 index 0000000000..dda1059608 --- /dev/null +++ b/NDTensors/ext/NDTensorsCUDAExt/truncate.jl @@ -0,0 +1,7 @@ +# import NDTensors:truncate! +# function NDTensors.truncate!(P::CuArray{ElT}; kwargs...)::Tuple{ElT,ElT} where {ElT} +# p_cpu = NDTensors.cpu(P) +# truncerr, docut = NDTensors.truncate!(p_cpu) +# P = typeof(P)(p_cpu) +# return truncerr, docut +# end \ No newline at end of file From f03d5cf5cb502b670924c6ce7d8c5bff7167024d Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 27 Sep 2023 16:08:04 -0400 Subject: [PATCH 028/133] [experiment] Add comment for dot, its still broken for GPU --- src/itensor.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/itensor.jl b/src/itensor.jl index 19286f065c..648275e6f0 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1892,6 +1892,7 @@ diag(T::ITensor) = diag(tensor(T)) mul!(C::ITensor, A::ITensor, B::ITensor, args...)::ITensor = contract!(C, A, B, args...) +## TODO this operation does not work with GPU dot(A::ITensor, B::ITensor) = (dag(A) * B)[] inner(y::ITensor, A::ITensor, x::ITensor) = (dag(y) * A * x)[] @@ -1963,8 +1964,11 @@ map(f, x::ITensor) = itensor(map(f, tensor(x))) w .+= a .* v ``` """ -axpy!(a::Number, v::ITensor, w::ITensor) = (w .+= a .* v) - +function axpy!(a::Number, v::ITensor, w::ITensor) + #(w .+= a .* v) + w .+= a .* v + #axpy!(a, tensor(v), tensor(w)) +end """ axpby!(a,v,b,w) From c7732a3029f585052a7fad7eeaa0e62765b8c044 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 28 Sep 2023 09:29:50 -0400 Subject: [PATCH 029/133] Fix printing for GPU no scalar indexing --- NDTensors/src/tensor/tensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/tensor/tensor.jl b/NDTensors/src/tensor/tensor.jl index 84699222f1..548c8995ef 100644 --- a/NDTensors/src/tensor/tensor.jl +++ b/NDTensors/src/tensor/tensor.jl @@ -425,5 +425,5 @@ end # Printing # -print_tensor(io::IO, T::Tensor) = Base.print_array(io, T) +print_tensor(io::IO, T::Tensor) = Base.print_array(io, data(T)) print_tensor(io::IO, T::Tensor{<:Number,1}) = Base.print_array(io, reshape(T, (dim(T), 1))) From 2c486b312900934b1274a0fe512ffad3a34abd3a Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 28 Sep 2023 13:00:24 -0400 Subject: [PATCH 030/133] dot function calls permute and linearalgebra dot instead of gemm [experimental] --- NDTensors/src/tensoralgebra/generic_tensor_operations.jl | 4 ++++ src/itensor.jl | 6 +++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/NDTensors/src/tensoralgebra/generic_tensor_operations.jl b/NDTensors/src/tensoralgebra/generic_tensor_operations.jl index 170b472bd8..406a05260c 100644 --- a/NDTensors/src/tensoralgebra/generic_tensor_operations.jl +++ b/NDTensors/src/tensoralgebra/generic_tensor_operations.jl @@ -189,3 +189,7 @@ end function outer end const ⊗ = outer + +function LinearAlgebra.dot(A::Tensor, B::Tensor) + return LinearAlgebra.dot(data(A), data(B)) +end \ No newline at end of file diff --git a/src/itensor.jl b/src/itensor.jl index 648275e6f0..6c27746b58 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1893,8 +1893,12 @@ diag(T::ITensor) = diag(tensor(T)) mul!(C::ITensor, A::ITensor, B::ITensor, args...)::ITensor = contract!(C, A, B, args...) ## TODO this operation does not work with GPU -dot(A::ITensor, B::ITensor) = (dag(A) * B)[] +dot(A::ITensor, B::ITensor) = _dot(dag(A), B)#(dag(A) * B)[] +function _dot(A::ITensor, B::ITensor) + B = permute(B, inds(A)) + dot(tensor(A), tensor(B)) +end inner(y::ITensor, A::ITensor, x::ITensor) = (dag(y) * A * x)[] inner(y::ITensor, x::ITensor) = (dag(y) * x)[] From dcf247772ba3f8554a27214d7a1263df6a865e8a Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 29 Sep 2023 11:08:39 -0400 Subject: [PATCH 031/133] Revert axpy and working on broadcast --- src/itensor.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/itensor.jl b/src/itensor.jl index 6c27746b58..29855b0568 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1969,9 +1969,7 @@ w .+= a .* v ``` """ function axpy!(a::Number, v::ITensor, w::ITensor) - #(w .+= a .* v) w .+= a .* v - #axpy!(a, tensor(v), tensor(w)) end """ axpby!(a,v,b,w) From 72e0df905843cbc2b8dd63fd560ec49a122e8698 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 2 Oct 2023 12:53:56 -0400 Subject: [PATCH 032/133] add comments about broken code --- src/broadcast.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/broadcast.jl b/src/broadcast.jl index af3082ac1a..7f8067083a 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -378,6 +378,7 @@ end # C .= β .* C .+ α .* A .* B # +## TODO this code doesn't actually get called function Base.copyto!( T::ITensor, bc::Broadcasted{ITensorOpScalarStyle,<:Any,typeof(+),<:Tuple{Vararg{Broadcasted}}}, @@ -414,6 +415,7 @@ end # For A .+= α # +## TODO this code fails because of scalar indexing function Base.copyto!( T::ITensor, bc::Broadcasted{ @@ -495,6 +497,7 @@ end # For B .= f.(B) + g.(A) # +## TODO this code isn't called properly function Base.copyto!( R::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{Broadcasted}}} ) From 1cfaa2953d432357d8784c8465275b57f3bd1d0e Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 2 Oct 2023 16:30:48 -0400 Subject: [PATCH 033/133] Force truncate to be done on CPU [experimental] --- NDTensors/src/linearalgebra/linearalgebra.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index 65c98b5437..e66927ad79 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -235,9 +235,11 @@ function LinearAlgebra.eigen( VM = VM[:, p] if truncate + cpu_dm = NDTensors.cpu(DM) truncerr, _ = truncate!( - DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff, kwargs... + cpu_dm; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff, kwargs... ) + DM = adapt(typeof(DM), cpu_dm) dD = length(DM) if dD < size(VM, 2) VM = VM[:, 1:dD] From 442fe95afd247cf830c453db7559c69485f8e13a Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 2 Oct 2023 16:35:50 -0400 Subject: [PATCH 034/133] remove show --- src/broadcast.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index 7f8067083a..08a21a69f1 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -134,7 +134,6 @@ find_type(::Type{T}, ::Tuple{}) where {T} = nothing # function Base.copyto!(T::ITensor, bc::Broadcasted) - @show typeof(bc) return error( "The broadcasting operation you are attempting is not yet implemented for ITensors, please raise an issue if you would like it to be supported.", ) @@ -171,7 +170,7 @@ function Base.copyto!( ) α = find_type(Number, bc.args) A = find_type(ITensor, bc.args) - map!((t, a) -> bc.f(a, α), T, T, A) + map!((t, a) -> a / α, T, T, A) return T end @@ -223,7 +222,8 @@ function Base.copyto!( ) α = find_type(Number, bc.args) A = find_type(ITensor, bc.args) - map!((t, a) -> bc.f(α, a), T, T, A) + #map!((t, a) -> bc.f(α, a), T, T, A) + map!((t,a) -> /(α,a), T, T, A) return T end @@ -242,6 +242,7 @@ end # For B .= A .^ 2 # +## TODO this fails on GPU to compile function Base.copyto!( R::ITensor, bc::Broadcasted{ITensorOpScalarStyle,<:Any,typeof(Base.literal_pow)} ) @@ -281,11 +282,11 @@ end # function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{ITensor}}}) - return (r, t) -> bc.f(r, t) + return (r, t) -> +(r, t) end function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(-),<:Tuple{Vararg{ITensor}}}) - return (r, t) -> bc.f(r, t) + return (r, t) -> -(r, t) end function Base.copyto!( @@ -352,7 +353,7 @@ function Base.copyto!( powf = find_type(Base.RefValue{<:Function}, bc_bc.args) if !isnothing(α) && !isnothing(A) - map!((r, t) -> bc.f(r, bc_bc.f(t, α)), T, T, A) + map!((r, t) -> +(r, *(t, α)), T, T, A) elseif !isnothing(γ) && !isnothing(A) && !isnothing(powf) map!((r, t) -> bc.f(r, bc_bc.f(powf[], t, γ[])), T, T, A) else From 74791ac2639dca68d1bf7c1204f8df73b8417ca1 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 2 Oct 2023 16:39:29 -0400 Subject: [PATCH 035/133] unified memory hack [experimental] --- src/mps/dmrg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index a6b78c6f8c..4ef9976fdc 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -278,8 +278,8 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) end energy = vals[1] - phi::ITensor = vecs[1] - + phi::ITensor = itensor(convert(typeof(tensor(phi)), tensor(vecs[1]))) + ortho = ha == 1 ? "left" : "right" drho = nothing From ca023c6e84220947bc51b11cf0ccaee00f2fec46 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 2 Oct 2023 16:50:43 -0400 Subject: [PATCH 036/133] `any` performs scalar operations [experimental] --- NDTensors/src/linearalgebra/linearalgebra.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index e66927ad79..ebfabafed5 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -219,13 +219,13 @@ function LinearAlgebra.eigen( matrixT = matrix(T) ## TODO this doesn't work for GPU - if any(!isfinite, matrixT) - throw( - ArgumentError( - "Trying to perform the eigendecomposition of a matrix containing NaNs or Infs" - ), - ) - end + # if any(!isfinite, matrixT) + # throw( + # ArgumentError( + # "Trying to perform the eigendecomposition of a matrix containing NaNs or Infs" + # ), + # ) + # end DM, VM = eigen(matrixT) From e1ef74ec0d8fbeb297ac1f25b00204989351f0cb Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 2 Oct 2023 16:53:18 -0400 Subject: [PATCH 037/133] Make a sweepup function to deal with CUDA memory [experimental] --- NDTensors/ext/NDTensorsCUDAExt/factorize.jl | 4 ++++ src/mps/dmrg.jl | 5 +++++ 2 files changed, 9 insertions(+) create mode 100644 NDTensors/ext/NDTensorsCUDAExt/factorize.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/factorize.jl b/NDTensors/ext/NDTensorsCUDAExt/factorize.jl new file mode 100644 index 0000000000..212ada416a --- /dev/null +++ b/NDTensors/ext/NDTensorsCUDAExt/factorize.jl @@ -0,0 +1,4 @@ +function ITensors.sweepup(::Type{<:CuArray}) + GC.gc() + CUDA.reclaim() +end \ No newline at end of file diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index 4ef9976fdc..0ad8bf9166 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -346,6 +346,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) ) end end + sweepup(typeof(data(psi[1]))) if outputlevel >= 1 @printf( "After sweep %d energy=%s maxlinkdim=%d maxerr=%.2E time=%.3f\n", @@ -363,6 +364,10 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) return (energy, psi) end +function sweepup(::Type{<:AbstractArray}) + GC.gc() +end + function _dmrg_sweeps(; nsweeps, maxdim=typemax(Int), mindim=1, cutoff=1E-8, noise=0.0, kwargs... ) From 5569fcaf686e68e33aacbe714894db05c9771255 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 2 Oct 2023 17:27:48 -0400 Subject: [PATCH 038/133] format --- NDTensors/ext/NDTensorsCUDAExt/adapt.jl | 2 +- NDTensors/ext/NDTensorsCUDAExt/factorize.jl | 6 +++--- NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl | 1 - NDTensors/ext/NDTensorsCUDAExt/truncate.jl | 2 +- NDTensors/src/tensoralgebra/generic_tensor_operations.jl | 2 +- src/broadcast.jl | 2 +- src/itensor.jl | 4 ++-- src/mps/dmrg.jl | 4 ++-- 8 files changed, 11 insertions(+), 12 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/adapt.jl b/NDTensors/ext/NDTensorsCUDAExt/adapt.jl index abc0c295c8..165b5405ce 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/adapt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/adapt.jl @@ -16,7 +16,7 @@ buffertype(::NDTensorCuArrayAdaptor{B}) where {B} = B function Adapt.adapt_storage(adaptor::NDTensorCuArrayAdaptor, xs::AbstractArray) ElT = eltype(xs) BufT = buffertype(adaptor) - return isbits(xs) ? xs : adapt(CuArray{ElT, 1, BufT},xs) + return isbits(xs) ? xs : adapt(CuArray{ElT,1,BufT}, xs) end function NDTensors.adapt_storagetype( diff --git a/NDTensors/ext/NDTensorsCUDAExt/factorize.jl b/NDTensors/ext/NDTensorsCUDAExt/factorize.jl index 212ada416a..848461bdd2 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/factorize.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/factorize.jl @@ -1,4 +1,4 @@ function ITensors.sweepup(::Type{<:CuArray}) - GC.gc() - CUDA.reclaim() -end \ No newline at end of file + GC.gc() + return CUDA.reclaim() +end diff --git a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl index 5796a9b521..93c29ee5de 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl @@ -13,6 +13,5 @@ function NDTensors.svd_catch_error(A::CuMatrix; alg="JacobiAlgorithm") end function NDTensors.mul!!(CM::CuArray, AM::AbstractArray, BM::AbstractArray, α, β) - return mul!(CM, AM, BM, α, β) end diff --git a/NDTensors/ext/NDTensorsCUDAExt/truncate.jl b/NDTensors/ext/NDTensorsCUDAExt/truncate.jl index dda1059608..2f1ffa0a86 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/truncate.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/truncate.jl @@ -4,4 +4,4 @@ # truncerr, docut = NDTensors.truncate!(p_cpu) # P = typeof(P)(p_cpu) # return truncerr, docut -# end \ No newline at end of file +# end diff --git a/NDTensors/src/tensoralgebra/generic_tensor_operations.jl b/NDTensors/src/tensoralgebra/generic_tensor_operations.jl index 406a05260c..7e16da14bf 100644 --- a/NDTensors/src/tensoralgebra/generic_tensor_operations.jl +++ b/NDTensors/src/tensoralgebra/generic_tensor_operations.jl @@ -192,4 +192,4 @@ const ⊗ = outer function LinearAlgebra.dot(A::Tensor, B::Tensor) return LinearAlgebra.dot(data(A), data(B)) -end \ No newline at end of file +end diff --git a/src/broadcast.jl b/src/broadcast.jl index 08a21a69f1..6d8d83e6fe 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -223,7 +223,7 @@ function Base.copyto!( α = find_type(Number, bc.args) A = find_type(ITensor, bc.args) #map!((t, a) -> bc.f(α, a), T, T, A) - map!((t,a) -> /(α,a), T, T, A) + map!((t, a) -> /(α, a), T, T, A) return T end diff --git a/src/itensor.jl b/src/itensor.jl index 29855b0568..d62ca55d87 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1897,7 +1897,7 @@ dot(A::ITensor, B::ITensor) = _dot(dag(A), B)#(dag(A) * B)[] function _dot(A::ITensor, B::ITensor) B = permute(B, inds(A)) - dot(tensor(A), tensor(B)) + return dot(tensor(A), tensor(B)) end inner(y::ITensor, A::ITensor, x::ITensor) = (dag(y) * A * x)[] inner(y::ITensor, x::ITensor) = (dag(y) * x)[] @@ -1969,7 +1969,7 @@ w .+= a .* v ``` """ function axpy!(a::Number, v::ITensor, w::ITensor) - w .+= a .* v + return w .+= a .* v end """ axpby!(a,v,b,w) diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index 0ad8bf9166..4a69b54ddd 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -279,7 +279,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) energy = vals[1] phi::ITensor = itensor(convert(typeof(tensor(phi)), tensor(vecs[1]))) - + ortho = ha == 1 ? "left" : "right" drho = nothing @@ -365,7 +365,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) end function sweepup(::Type{<:AbstractArray}) - GC.gc() + return GC.gc() end function _dmrg_sweeps(; From 67d1a6b134ba15a919a7c13f4de0031be39a8158 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Tue, 3 Oct 2023 09:24:00 -0400 Subject: [PATCH 039/133] Revert to main --- src/itensor.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/itensor.jl b/src/itensor.jl index 94901f25f6..a1ce51097b 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1968,9 +1968,7 @@ map(f, x::ITensor) = itensor(map(f, tensor(x))) w .+= a .* v ``` """ -function axpy!(a::Number, v::ITensor, w::ITensor) - return w .+= a .* v -end +axpy!(a::Number, v::ITensor, w::ITensor) = (w .+= a .* v) """ axpby!(a,v,b,w) From 4ccc09273d845209451aaa526b409e36904fa6f5 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Tue, 3 Oct 2023 09:27:05 -0400 Subject: [PATCH 040/133] Remove factorize and truncate->sweepup [experimental] --- NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl | 2 +- .../ext/NDTensorsCUDAExt/{factorize.jl => sweeupup.jl} | 0 NDTensors/ext/NDTensorsCUDAExt/truncate.jl | 7 ------- 3 files changed, 1 insertion(+), 8 deletions(-) rename NDTensors/ext/NDTensorsCUDAExt/{factorize.jl => sweeupup.jl} (100%) delete mode 100644 NDTensors/ext/NDTensorsCUDAExt/truncate.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl index 383cd6b8ad..1d9658728f 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl @@ -20,5 +20,5 @@ include("imports.jl") include("set_types.jl") include("adapt.jl") include("linearalgebra.jl") -include("truncate.jl") +include("sweeupup.jl") end diff --git a/NDTensors/ext/NDTensorsCUDAExt/factorize.jl b/NDTensors/ext/NDTensorsCUDAExt/sweeupup.jl similarity index 100% rename from NDTensors/ext/NDTensorsCUDAExt/factorize.jl rename to NDTensors/ext/NDTensorsCUDAExt/sweeupup.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/truncate.jl b/NDTensors/ext/NDTensorsCUDAExt/truncate.jl deleted file mode 100644 index 2f1ffa0a86..0000000000 --- a/NDTensors/ext/NDTensorsCUDAExt/truncate.jl +++ /dev/null @@ -1,7 +0,0 @@ -# import NDTensors:truncate! -# function NDTensors.truncate!(P::CuArray{ElT}; kwargs...)::Tuple{ElT,ElT} where {ElT} -# p_cpu = NDTensors.cpu(P) -# truncerr, docut = NDTensors.truncate!(p_cpu) -# P = typeof(P)(p_cpu) -# return truncerr, docut -# end From c3cf6fbc4403c6942a73618614d2a54b3c15a94f Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Tue, 3 Oct 2023 12:00:34 -0400 Subject: [PATCH 041/133] remove sweepup function --- NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl | 1 - NDTensors/ext/NDTensorsCUDAExt/sweeupup.jl | 4 ---- src/mps/dmrg.jl | 5 ----- 3 files changed, 10 deletions(-) delete mode 100644 NDTensors/ext/NDTensorsCUDAExt/sweeupup.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl index 1d9658728f..f6eb791f93 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl @@ -20,5 +20,4 @@ include("imports.jl") include("set_types.jl") include("adapt.jl") include("linearalgebra.jl") -include("sweeupup.jl") end diff --git a/NDTensors/ext/NDTensorsCUDAExt/sweeupup.jl b/NDTensors/ext/NDTensorsCUDAExt/sweeupup.jl deleted file mode 100644 index 848461bdd2..0000000000 --- a/NDTensors/ext/NDTensorsCUDAExt/sweeupup.jl +++ /dev/null @@ -1,4 +0,0 @@ -function ITensors.sweepup(::Type{<:CuArray}) - GC.gc() - return CUDA.reclaim() -end diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index 8e0e45e2a7..88ff640bf0 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -347,7 +347,6 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) ) end end - sweepup(typeof(data(psi[1]))) if outputlevel >= 1 @printf( "After sweep %d energy=%s maxlinkdim=%d maxerr=%.2E time=%.3f\n", @@ -365,10 +364,6 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) return (energy, psi) end -function sweepup(::Type{<:AbstractArray}) - return GC.gc() -end - function _dmrg_sweeps(; nsweeps, maxdim=typemax(Int), mindim=1, cutoff=1E-8, noise=0.0, kwargs... ) From d3a7b24f4ab4aaacb481d2e05d02434e7e88e579 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Tue, 3 Oct 2023 12:23:21 -0400 Subject: [PATCH 042/133] Don't return anonymous function --- src/broadcast.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index 6d8d83e6fe..95cfd98daf 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -282,11 +282,11 @@ end # function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{ITensor}}}) - return (r, t) -> +(r, t) + return + end function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(-),<:Tuple{Vararg{ITensor}}}) - return (r, t) -> -(r, t) + return - end function Base.copyto!( From 7efb9dd26a6ccb6f901ef491f0cee6ad69ec5f47 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Tue, 3 Oct 2023 12:24:39 -0400 Subject: [PATCH 043/133] remove old line --- src/mps/dmrg.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index 88ff640bf0..fdf17bbefb 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -278,7 +278,6 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) end energy = vals[1] -# phi::ITensor = itensor(storagetype(phi)(convert(typeof(data(phi)),data(vecs[1]))), inds(phi)) phi::ITensor = itensor(convert(typeof(tensor(phi)), tensor(vecs[1]))) ortho = ha == 1 ? "left" : "right" From 5519053524967554ef816d01d7d24420fea944c3 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 3 Oct 2023 15:33:24 -0400 Subject: [PATCH 044/133] revert changes to itensors.jl --- src/itensor.jl | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/itensor.jl b/src/itensor.jl index dd44a79d01..779b8d29e4 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1893,12 +1893,8 @@ diag(T::ITensor) = diag(tensor(T)) mul!(C::ITensor, A::ITensor, B::ITensor, args...)::ITensor = contract!(C, A, B, args...) ## TODO this operation does not work with GPU -dot(A::ITensor, B::ITensor) = _dot(dag(A), B) +dot(A::ITensor, B::ITensor) = (dag(A) * B)[] -function _dot(A::ITensor, B::ITensor) - B = permute(B, inds(A)) - return dot(tensor(A), tensor(B)) -end inner(y::ITensor, A::ITensor, x::ITensor) = (dag(y) * A * x)[] inner(y::ITensor, x::ITensor) = (dag(y) * x)[] @@ -1986,11 +1982,11 @@ Scale the ITensor A by x in-place. May also be written `rmul!`. A .*= x ``` """ -scale!(T::ITensor, α::Number) = itensor(data(T) .*= α, inds(T)) +scale!(T::ITensor, α::Number) = (T .*= α) -rmul!(T::ITensor, α::Number) = itensor(data(T) .*= α, inds(T)) +rmul!(T::ITensor, α::Number) = (T .*= α) -lmul!(T::ITensor, α::Number) = itensor(data(T) .= α .* data(T), inds(T)) +lmul!(T::ITensor, α::Number) = (T .= α .* T) """ mul!(A::ITensor, x::Number, B::ITensor) @@ -1998,9 +1994,9 @@ lmul!(T::ITensor, α::Number) = itensor(data(T) .= α .* data(T), inds(T)) Scalar multiplication of ITensor B with x, and store the result in A. Like `A .= x .* B`. """ -mul!(R::ITensor, α::Number, T::ITensor) = itensor(data(R) .= α * data(T), inds(T)) +mul!(R::ITensor, α::Number, T::ITensor) = (R .= α .* T) -mul!(R::ITensor, T::ITensor, α::Number) = itensor(data(R) .= data(T) * α, inds(T)) +mul!(R::ITensor, T::ITensor, α::Number) = (R .= T .* α) ######################### # End ITensor Operations # From 69664cf549dca0dced3b9f6869b2c30feb8c7e64 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 3 Oct 2023 15:36:41 -0400 Subject: [PATCH 045/133] remove dot from NDTensors --- NDTensors/src/tensoralgebra/generic_tensor_operations.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/NDTensors/src/tensoralgebra/generic_tensor_operations.jl b/NDTensors/src/tensoralgebra/generic_tensor_operations.jl index 3088a65187..db0a740e22 100644 --- a/NDTensors/src/tensoralgebra/generic_tensor_operations.jl +++ b/NDTensors/src/tensoralgebra/generic_tensor_operations.jl @@ -189,7 +189,3 @@ end function outer end const ⊗ = outer - -function LinearAlgebra.dot(A::Tensor, B::Tensor) - return LinearAlgebra.dot(data(A), data(B)) -end From d5f5dad8516285438362f8334f6c2b7cfcf81eaa Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 11:53:47 -0400 Subject: [PATCH 046/133] Broadcast mostly working for GPU. Make more fmap functions and increase functionality --- src/broadcast.jl | 42 +++++++++++++++++++++++++++++------------- 1 file changed, 29 insertions(+), 13 deletions(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index 95cfd98daf..6259324f42 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -183,9 +183,9 @@ function Base.copyto!( ) T1, T2 = bc.args if R === T1 - map!((t1, t2) -> bc.f(t1, t2), R, T1, T2) + map!((t1, t2) -> /(t1, t2), R, T1, T2) elseif R === T2 - map!((t1, t2) -> bc.f(t2, t1), R, T2, T1) + map!((t1, t2) -> /(t2, t1), R, T2, T1) else error("When dividing two ITensors in-place, one must be the same as the output ITensor") end @@ -222,7 +222,6 @@ function Base.copyto!( ) α = find_type(Number, bc.args) A = find_type(ITensor, bc.args) - #map!((t, a) -> bc.f(α, a), T, T, A) map!((t, a) -> /(α, a), T, T, A) return T end @@ -242,7 +241,6 @@ end # For B .= A .^ 2 # -## TODO this fails on GPU to compile function Base.copyto!( R::ITensor, bc::Broadcasted{ITensorOpScalarStyle,<:Any,typeof(Base.literal_pow)} ) @@ -250,7 +248,7 @@ function Base.copyto!( powf = find_type(Base.RefValue{<:Function}, bc.args).x @assert !isnothing(powf) T = find_type(ITensor, bc.args) - map!((r, t) -> bc.f(^, t, α), R, R, T) + map!((r, t) -> Base.literal_pow(^, t, α), R, R, T) return R end @@ -277,18 +275,31 @@ function Base.copyto!( return T end -# -# B .+= A -# -function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{ITensor}}}) +function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(+)}) where {N} return + end -function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(-),<:Tuple{Vararg{ITensor}}}) +function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(-)}) where {N} return - end +function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(^)}) where {N} + return ^ +end + +function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(*)}) where {N} + return * +end + +function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(Base.literal_pow)}) where {N} + return Base.literal_pow +end + +# +# B .+= A +# + function Base.copyto!( T::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{ITensor}}} ) @@ -344,6 +355,8 @@ function Base.copyto!( ) C = find_type(ITensor, bc.args) bc_bc = find_type(Broadcasted, bc.args) + + bc = Broadcasted{ITensorStyle}(bc.f, (bc.args[1],), axes(bc)) if T === C A = find_type(ITensor, bc_bc.args) α = find_type(Number, bc_bc.args) @@ -351,11 +364,14 @@ function Base.copyto!( # Check if it is the case .^(::Int) γ = find_type(Base.RefValue{<:Val}, bc_bc.args) powf = find_type(Base.RefValue{<:Function}, bc_bc.args) + ## Putting fmap in the map call still doesn't actually grab the function and causes GPU to fail so just realize the function slightly earlier here + f1 = fmap(bc) + f2 = fmap(bc_bc) if !isnothing(α) && !isnothing(A) - map!((r, t) -> +(r, *(t, α)), T, T, A) + map!((r, t) -> f1(r, f2(t, α)), T, T, A) elseif !isnothing(γ) && !isnothing(A) && !isnothing(powf) - map!((r, t) -> bc.f(r, bc_bc.f(powf[], t, γ[])), T, T, A) + map!((r, t) -> f1(r, f2(powf[], t, γ[])), T, T, A) else # In-place contraction: # C .+= α .* A .* B @@ -498,7 +514,7 @@ end # For B .= f.(B) + g.(A) # -## TODO this code isn't called properly +## TODO check to see if this code is being called as expected function Base.copyto!( R::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{Broadcasted}}} ) From e607404b5869ad0e3b76c7888bc58abac8913653 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 11:53:54 -0400 Subject: [PATCH 047/133] Fix matrix decomp --- src/tensor_operations/matrix_decomposition.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/tensor_operations/matrix_decomposition.jl b/src/tensor_operations/matrix_decomposition.jl index 9e93f3f03e..03a622b692 100644 --- a/src/tensor_operations/matrix_decomposition.jl +++ b/src/tensor_operations/matrix_decomposition.jl @@ -573,8 +573,7 @@ function factorize_eigen(A::ITensor, Linds...; kwargs...) # (Lis..., prime(Lis)...) delta_A2 = replaceinds(delta_A2, Lis, dag(simLis)) noprime!(delta_A2) - itensor(data(A2) .+= data(delta_A2), inds(A2)) - #A2 += delta_A2 + A2 += delta_A2 end F = eigen(A2, Lis, simLis; ishermitian=true, kwargs...) D, _, spec = F From 7779d12ffb9233a17473a82afdbbb47aa51a6416 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 11:54:10 -0400 Subject: [PATCH 048/133] format --- src/broadcast.jl | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index 6259324f42..344a0878c7 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -275,24 +275,33 @@ function Base.copyto!( return T end - -function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(+)}) where {N} +function fmap( + bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(+)} +) where {N} return + end -function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(-)}) where {N} +function fmap( + bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(-)} +) where {N} return - end -function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(^)}) where {N} +function fmap( + bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(^)} +) where {N} return ^ end -function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(*)}) where {N} +function fmap( + bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(*)} +) where {N} return * end -function fmap(bc::Broadcasted{<:Union{ITensorStyle, ITensorOpScalarStyle},N,typeof(Base.literal_pow)}) where {N} +function fmap( + bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(Base.literal_pow)} +) where {N} return Base.literal_pow end From fdd8e1cbfd33dc414c49dc3ddf1785def8a0aec3 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 12:58:52 -0400 Subject: [PATCH 049/133] We don't actually need fmap just don't use bc.f directly --- src/broadcast.jl | 54 ++++++++++++++++++------------------------------ src/itensor.jl | 2 ++ 2 files changed, 22 insertions(+), 34 deletions(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index 344a0878c7..3dd8aaf715 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -170,7 +170,8 @@ function Base.copyto!( ) α = find_type(Number, bc.args) A = find_type(ITensor, bc.args) - map!((t, a) -> a / α, T, T, A) + f = bc.f + map!((t, a) -> f(a, α), T, T, A) return T end @@ -182,10 +183,13 @@ function Base.copyto!( R::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(/),<:Tuple{<:ITensor,<:ITensor}} ) T1, T2 = bc.args + f = bc.f if R === T1 - map!((t1, t2) -> /(t1, t2), R, T1, T2) + #map!((t1, t2) -> f(t1, t2), R, T1, T2) + map!(f,R, T1, T2) elseif R === T2 - map!((t1, t2) -> /(t2, t1), R, T2, T1) + #map!((t1, t2) -> f(t2, t1), R, T2, T1) + map!(f, R, T2, T1) else error("When dividing two ITensors in-place, one must be the same as the output ITensor") end @@ -222,7 +226,8 @@ function Base.copyto!( ) α = find_type(Number, bc.args) A = find_type(ITensor, bc.args) - map!((t, a) -> /(α, a), T, T, A) + f = bc.f + map!((t, a) -> f(α, a), T, T, A) return T end @@ -248,7 +253,8 @@ function Base.copyto!( powf = find_type(Base.RefValue{<:Function}, bc.args).x @assert !isnothing(powf) T = find_type(ITensor, bc.args) - map!((r, t) -> Base.literal_pow(^, t, α), R, R, T) + f = bc.f + map!((r, t) -> f(^, t, α), R, R, T) return R end @@ -275,34 +281,13 @@ function Base.copyto!( return T end -function fmap( - bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(+)} -) where {N} - return + +function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{ITensor}}}) + return (r, t) -> bc.f(r, t) end -function fmap( - bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(-)} -) where {N} - return - -end - -function fmap( - bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(^)} -) where {N} - return ^ -end - -function fmap( - bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(*)} -) where {N} - return * -end -function fmap( - bc::Broadcasted{<:Union{ITensorStyle,ITensorOpScalarStyle},N,typeof(Base.literal_pow)} -) where {N} - return Base.literal_pow +function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(-),<:Tuple{Vararg{ITensor}}}) + return (r, t) -> bc.f(r, t) end # @@ -374,8 +359,8 @@ function Base.copyto!( γ = find_type(Base.RefValue{<:Val}, bc_bc.args) powf = find_type(Base.RefValue{<:Function}, bc_bc.args) ## Putting fmap in the map call still doesn't actually grab the function and causes GPU to fail so just realize the function slightly earlier here - f1 = fmap(bc) - f2 = fmap(bc_bc) + f1 = bc.f + f2 = bc_bc.f if !isnothing(α) && !isnothing(A) map!((r, t) -> f1(r, f2(t, α)), T, T, A) @@ -391,7 +376,7 @@ function Base.copyto!( else A, B = bc_bc_bc.args end - mul!(T, A, B, β, bc.f(1)) + mul!(T, A, B, β, f1(1)) end else error("When adding two ITensors in-place, one must be the same as the output ITensor") @@ -511,8 +496,9 @@ function Base.copyto!( ) R̃ = find_type(ITensor, bc.args) bc2 = find_type(Broadcasted, bc.args) + f = bc2.f if R === R̃ - map!((r, t) -> r + bc2.f(t), R, R, bc2.args[1]) + map!((r, t) -> r + f(t), R, R, bc2.args[1]) else error("In C .= B .+ f.(A), C and B must be the same ITensor") end diff --git a/src/itensor.jl b/src/itensor.jl index 779b8d29e4..e95604db11 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1965,6 +1965,7 @@ w .+= a .* v ``` """ axpy!(a::Number, v::ITensor, w::ITensor) = (w .+= a .* v) + """ axpby!(a,v,b,w) @@ -1997,6 +1998,7 @@ Like `A .= x .* B`. mul!(R::ITensor, α::Number, T::ITensor) = (R .= α .* T) mul!(R::ITensor, T::ITensor, α::Number) = (R .= T .* α) + ######################### # End ITensor Operations # From 9316868cff33a754ef5bc22262efb43d76d0beaa Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 12:59:16 -0400 Subject: [PATCH 050/133] format --- src/broadcast.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index 3dd8aaf715..390fd1f872 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -186,7 +186,7 @@ function Base.copyto!( f = bc.f if R === T1 #map!((t1, t2) -> f(t1, t2), R, T1, T2) - map!(f,R, T1, T2) + map!(f, R, T1, T2) elseif R === T2 #map!((t1, t2) -> f(t2, t1), R, T2, T1) map!(f, R, T2, T1) @@ -285,7 +285,6 @@ function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{ITenso return (r, t) -> bc.f(r, t) end - function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(-),<:Tuple{Vararg{ITensor}}}) return (r, t) -> bc.f(r, t) end From 8bc1f3cb5786c588e9f25834d0a1a013f0a839ef Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 13:08:14 -0400 Subject: [PATCH 051/133] Fix broadcast --- src/broadcast.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index 390fd1f872..f746aae35e 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -185,11 +185,12 @@ function Base.copyto!( T1, T2 = bc.args f = bc.f if R === T1 - #map!((t1, t2) -> f(t1, t2), R, T1, T2) - map!(f, R, T1, T2) + map!((t1, t2) -> f(t1, t2), R, T1, T2) + ## I tried this and it is numberically wrong + #map!(f, R, T1, T2) elseif R === T2 - #map!((t1, t2) -> f(t2, t1), R, T2, T1) - map!(f, R, T2, T1) + map!((t1, t2) -> f(t2, t1), R, T2, T1) + #map!(f, R, T2, T1) else error("When dividing two ITensors in-place, one must be the same as the output ITensor") end From 4d4e10dbf74b1a87b510ed805ed6c92fc821d674 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 13:22:03 -0400 Subject: [PATCH 052/133] Don't need to remake bc --- src/broadcast.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index f746aae35e..209e38db0c 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -350,7 +350,6 @@ function Base.copyto!( C = find_type(ITensor, bc.args) bc_bc = find_type(Broadcasted, bc.args) - bc = Broadcasted{ITensorStyle}(bc.f, (bc.args[1],), axes(bc)) if T === C A = find_type(ITensor, bc_bc.args) α = find_type(Number, bc_bc.args) From 1b60b789172d1e80f65129d3e81cec4bd4d2273a Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 13:25:24 -0400 Subject: [PATCH 053/133] revert, dont use `data(T)` in print_tensor --- NDTensors/src/tensor/tensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/tensor/tensor.jl b/NDTensors/src/tensor/tensor.jl index f29cc3c4c7..7e4af3f403 100644 --- a/NDTensors/src/tensor/tensor.jl +++ b/NDTensors/src/tensor/tensor.jl @@ -424,5 +424,5 @@ end # Printing # -print_tensor(io::IO, T::Tensor) = Base.print_array(io, data(T)) +print_tensor(io::IO, T::Tensor) = Base.print_array(io, T) print_tensor(io::IO, T::Tensor{<:Number,1}) = Base.print_array(io, reshape(T, (dim(T), 1))) From d9fa5d621868391f3a3565d50545b4ec63879f56 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 14:23:24 -0400 Subject: [PATCH 054/133] Make sure to return the functions + and - --- src/broadcast.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index 209e38db0c..a75af92f35 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -283,11 +283,11 @@ function Base.copyto!( end function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{ITensor}}}) - return (r, t) -> bc.f(r, t) + return + end function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(-),<:Tuple{Vararg{ITensor}}}) - return (r, t) -> bc.f(r, t) + return - end # From cdf7731b07b32d37499fdf43c6d6e276d232ec78 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 14:26:00 -0400 Subject: [PATCH 055/133] Call any on the parent of matrixT (`hermiaitan{ElT, ParentT}`) --- NDTensors/src/linearalgebra/linearalgebra.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index e337a08a62..8abd81b5a3 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -218,14 +218,13 @@ function LinearAlgebra.eigen( use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff) matrixT = matrix(T) - ## TODO this doesn't work for GPU - # if any(!isfinite, matrixT) - # throw( - # ArgumentError( - # "Trying to perform the eigendecomposition of a matrix containing NaNs or Infs" - # ), - # ) - # end + if any(!isfinite, parent(matrixT)) + throw( + ArgumentError( + "Trying to perform the eigendecomposition of a matrix containing NaNs or Infs" + ), + ) + end DM, VM = eigen(matrixT) @@ -403,6 +402,7 @@ function qx(qx::Function, T::DenseTensor{<:Any,2}; kwargs...) 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) @@ -424,6 +424,7 @@ matrix is unique. Returns a tuple (Q,R). function qr_positive(M::AbstractMatrix) sparseQ, R = qr(M) Q = convert(typeof(R), sparseQ) + Q = convert(typeof(R), sparseQ) nc = size(Q, 2) ## TODO issue here for GPU because tying to access indices for c in 1:nc @@ -449,6 +450,7 @@ matrix is unique. Returns a tuple (Q,L). function ql_positive(M::AbstractMatrix) sparseQ, L = ql(M) Q = convert(typeof(L), sparseQ) + 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) From 48249843718d94cdcd87e4daae7e0bd9a70b7686 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 14:28:36 -0400 Subject: [PATCH 056/133] Add comment about `f = bc.f` --- src/broadcast.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/broadcast.jl b/src/broadcast.jl index a75af92f35..bbb4570790 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -170,6 +170,9 @@ function Base.copyto!( ) α = find_type(Number, bc.args) A = find_type(ITensor, bc.args) + ## GPU compilers can have a problem when map is + ## Given bc.f. map seems to make a closure with a + ## relatively complicated signature f = bc.f map!((t, a) -> f(a, α), T, T, A) return T From c02dad80e22fb855365837efbf5cbb6d90b5128e Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 15:26:45 -0400 Subject: [PATCH 057/133] Revert to using fill! --- NDTensors/src/abstractarray/fill.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/abstractarray/fill.jl b/NDTensors/src/abstractarray/fill.jl index 4d9bcf4721..146459a4ce 100644 --- a/NDTensors/src/abstractarray/fill.jl +++ b/NDTensors/src/abstractarray/fill.jl @@ -13,5 +13,5 @@ function generic_zeros(arraytype::Type{<:AbstractArray}, dim::Integer=0) leaf_parenttype(arraytype), DefaultParameters() ) ElT = eltype(arraytype_specified) - return arraytype_specified(zeros(ElT, dim)) + return fill!(similar(arraytype_specified, dim), zero(ElT)) end From d96f32ee8a5ce8d9855697af6fd5cce55aa2f25d Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 15:39:17 -0400 Subject: [PATCH 058/133] Don't need to convert positive twice because these functions don't work on gpu anyway --- NDTensors/src/linearalgebra/linearalgebra.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index 8abd81b5a3..bd90e9a24f 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -424,7 +424,6 @@ matrix is unique. Returns a tuple (Q,R). function qr_positive(M::AbstractMatrix) sparseQ, R = qr(M) Q = convert(typeof(R), sparseQ) - Q = convert(typeof(R), sparseQ) nc = size(Q, 2) ## TODO issue here for GPU because tying to access indices for c in 1:nc @@ -450,7 +449,6 @@ matrix is unique. Returns a tuple (Q,L). function ql_positive(M::AbstractMatrix) sparseQ, L = ql(M) Q = convert(typeof(L), sparseQ) - 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) From aac86810970523ea5d7e10efdff08223650b3819 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 15:40:42 -0400 Subject: [PATCH 059/133] Fix `[]` for GPU --- src/itensor.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/itensor.jl b/src/itensor.jl index e95604db11..c88aa18cf0 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1140,15 +1140,16 @@ A[i => 1, i' => 2] # 2.0, same as: A[i' => 2, i => 1] @propagate_inbounds (getindex(T::ITensor, ivs::Vararg{Any,N})::Any) where {N} = _getindex(tensor(T), ivs...) + ## Allowing one to get the first ITensor element if its an order 0 tensor or an order 1 tensor with a dimension of 1. Also convert GPU back to CPU @propagate_inbounds function getindex(T::ITensor)::Any - if order(T) != 0 - throw( + if order(T) != 0 && dim(T) != 1 + ITensors.throw( DimensionMismatch( "In scalar(T) or T[], ITensor T is not a scalar (it has indices $(inds(T)))." ), ) end - return tensor(T)[] + return NDTensors.cpu(tensor(T))[] end function _vals(is::Indices, I::String...) @@ -1893,7 +1894,7 @@ diag(T::ITensor) = diag(tensor(T)) mul!(C::ITensor, A::ITensor, B::ITensor, args...)::ITensor = contract!(C, A, B, args...) ## TODO this operation does not work with GPU -dot(A::ITensor, B::ITensor) = (dag(A) * B)[] +dot(A::ITensor, B::ITensor) = NDTensors.cpu(dag(A) * B)[] inner(y::ITensor, A::ITensor, x::ITensor) = (dag(y) * A * x)[] inner(y::ITensor, x::ITensor) = (dag(y) * x)[] From 4f941b797793a8c82c2be453f38e00883467e060 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 15:42:13 -0400 Subject: [PATCH 060/133] Add a comment --- src/mps/dmrg.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index fdf17bbefb..b149168c0e 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -278,6 +278,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) end energy = vals[1] + ## Right now this convert problem is prevelent in an eigen solver implemented in KrylovKit phi::ITensor = itensor(convert(typeof(tensor(phi)), tensor(vecs[1]))) ortho = ha == 1 ? "left" : "right" From d859624349c530b11a77896573ce301d8235bc86 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 6 Oct 2023 15:52:39 -0400 Subject: [PATCH 061/133] format --- src/itensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/itensor.jl b/src/itensor.jl index c88aa18cf0..69fc41c7da 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1140,7 +1140,7 @@ A[i => 1, i' => 2] # 2.0, same as: A[i' => 2, i => 1] @propagate_inbounds (getindex(T::ITensor, ivs::Vararg{Any,N})::Any) where {N} = _getindex(tensor(T), ivs...) - ## Allowing one to get the first ITensor element if its an order 0 tensor or an order 1 tensor with a dimension of 1. Also convert GPU back to CPU +## Allowing one to get the first ITensor element if its an order 0 tensor or an order 1 tensor with a dimension of 1. Also convert GPU back to CPU @propagate_inbounds function getindex(T::ITensor)::Any if order(T) != 0 && dim(T) != 1 ITensors.throw( From f19e52dda9acda7dea3ba5347572241002fe1cc6 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Sun, 8 Oct 2023 12:06:35 -0400 Subject: [PATCH 062/133] Make a better comment and add back the original line --- src/mps/dmrg.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index b149168c0e..57805db2cc 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -278,8 +278,10 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) end energy = vals[1] - ## Right now this convert problem is prevelent in an eigen solver implemented in KrylovKit + ## Right now there is a conversion problem in CUDA.jl where `UnifiedMemory` Arrays are being converted + ## into `DeviceMemory`. This conversion line is here temporarily to fix that problem when it arises phi::ITensor = itensor(convert(typeof(tensor(phi)), tensor(vecs[1]))) + #phi::ITensor = vecs[1] ortho = ha == 1 ? "left" : "right" From 0f310e9f30a5d2b643f2bae4afde1047f06ec93c Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Sun, 8 Oct 2023 12:08:32 -0400 Subject: [PATCH 063/133] itensor.throws -> throws --- src/itensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/itensor.jl b/src/itensor.jl index 69fc41c7da..6c41bfaf8e 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1143,7 +1143,7 @@ A[i => 1, i' => 2] # 2.0, same as: A[i' => 2, i => 1] ## Allowing one to get the first ITensor element if its an order 0 tensor or an order 1 tensor with a dimension of 1. Also convert GPU back to CPU @propagate_inbounds function getindex(T::ITensor)::Any if order(T) != 0 && dim(T) != 1 - ITensors.throw( + throw( DimensionMismatch( "In scalar(T) or T[], ITensor T is not a scalar (it has indices $(inds(T)))." ), From 46e00ce720495bb337912ba49a86cafdaebe24f5 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Sun, 8 Oct 2023 12:20:22 -0400 Subject: [PATCH 064/133] Add `iscu` function --- NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl | 2 ++ NDTensors/ext/NDTensorsCUDAExt/imports.jl | 2 +- NDTensors/src/abstractarray/similar.jl | 3 +++ 3 files changed, 6 insertions(+), 1 deletion(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl index f6eb791f93..24dc8d7ae0 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl @@ -17,6 +17,8 @@ else end include("imports.jl") + +iscu(::Type{<:CuArray}) = true include("set_types.jl") include("adapt.jl") include("linearalgebra.jl") diff --git a/NDTensors/ext/NDTensorsCUDAExt/imports.jl b/NDTensors/ext/NDTensorsCUDAExt/imports.jl index d8736ae0a2..387a5b91ba 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/imports.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/imports.jl @@ -1,6 +1,6 @@ import NDTensors: cu, set_ndims, set_eltype, set_eltype_if_unspecified, similartype import NDTensors: ContractionProperties, _contract!, GemmBackend, auto_select_backend, _gemm! -import NDTensors.SetParameters: nparameters, get_parameter, set_parameter, default_parameter +import NDTensors.SetParameters: nparameters, get_parameter, set_parameter, default_parameter, iscu import .CUDA: CuArrayAdaptor diff --git a/NDTensors/src/abstractarray/similar.jl b/NDTensors/src/abstractarray/similar.jl index 7e8313cd1c..73e4e0d820 100644 --- a/NDTensors/src/abstractarray/similar.jl +++ b/NDTensors/src/abstractarray/similar.jl @@ -56,6 +56,9 @@ function similar(arraytype::Type{<:AbstractArray}, dims::Tuple) return similartype(arraytype, shape)(undef, NDTensors.to_shape(arraytype, shape)) end +# For when there are CUArray specific issues inline +iscu(A::AbstractArray) = iscu(leaf_parenttype(A)) +iscu(::Type{<:AbstractArray}) = false # This function actually allocates the data. # Catches conversions of dimensions specified by ranges # dimensions specified by integers with `Base.to_shape`. From f212536b64bb9480d400d724884127d9f9035327 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Sun, 8 Oct 2023 12:29:06 -0400 Subject: [PATCH 065/133] import iscu from NDTensors not setparameters --- NDTensors/ext/NDTensorsCUDAExt/imports.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/imports.jl b/NDTensors/ext/NDTensorsCUDAExt/imports.jl index 387a5b91ba..96b8383c74 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/imports.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/imports.jl @@ -1,6 +1,6 @@ import NDTensors: cu, set_ndims, set_eltype, set_eltype_if_unspecified, similartype import NDTensors: - ContractionProperties, _contract!, GemmBackend, auto_select_backend, _gemm! -import NDTensors.SetParameters: nparameters, get_parameter, set_parameter, default_parameter, iscu + ContractionProperties, _contract!, GemmBackend, auto_select_backend, _gemm!, iscu +import NDTensors.SetParameters: nparameters, get_parameter, set_parameter, default_parameter import .CUDA: CuArrayAdaptor From 02d84c7df96551351c1db8a8cb714d66c26e5af2 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Sun, 8 Oct 2023 12:38:39 -0400 Subject: [PATCH 066/133] Throw an error if trying to compute qr or ql positive --- NDTensors/src/linearalgebra/linearalgebra.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index bd90e9a24f..c5ecab4baf 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -422,6 +422,9 @@ non-negative. Such a QR decomposition of a matrix is unique. Returns a tuple (Q,R). """ function qr_positive(M::AbstractMatrix) + if iscu(M) + throw("Currently qr positive methods do not work for CuArrays because they require scalar operations. Please convert to CPU Array or use generic qr") + end sparseQ, R = qr(M) Q = convert(typeof(R), sparseQ) nc = size(Q, 2) @@ -447,6 +450,9 @@ non-negative. Such a QL decomposition of a matrix is unique. Returns a tuple (Q,L). """ function ql_positive(M::AbstractMatrix) + if iscu(M) + throw("Currently ql positive methods do not work for CuArrays because they require scalar operations. Please convert to CPU Array or use generic ql") + end sparseQ, L = ql(M) Q = convert(typeof(L), sparseQ) nr, nc = size(L) From c770277490730afc3a70dcf8cb5418644abb3675 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Sun, 8 Oct 2023 13:02:18 -0400 Subject: [PATCH 067/133] Push CPU conversion to NDTensors and make iscu for Tensors --- NDTensors/src/dense/densetensor.jl | 4 ++++ NDTensors/src/tensor/tensor.jl | 2 +- src/itensor.jl | 2 +- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/NDTensors/src/dense/densetensor.jl b/NDTensors/src/dense/densetensor.jl index 769f47673d..fd0f72bbd0 100644 --- a/NDTensors/src/dense/densetensor.jl +++ b/NDTensors/src/dense/densetensor.jl @@ -80,6 +80,10 @@ end # Single index # +@propagate_inbounds function getindex(T::DenseTensor{<:Number}) + return (iscu(T) ? NDTensors.cpu(data(T))[] : data(T)[]) +end + @propagate_inbounds function getindex(T::DenseTensor{<:Number}, I::Integer...) Base.@_inline_meta return getindex(data(T), Base._sub2ind(T, I...)) diff --git a/NDTensors/src/tensor/tensor.jl b/NDTensors/src/tensor/tensor.jl index 7e4af3f403..87613bf760 100644 --- a/NDTensors/src/tensor/tensor.jl +++ b/NDTensors/src/tensor/tensor.jl @@ -278,7 +278,7 @@ matrix(T::Tensor{<:Number,2}) = array(T) vector(T::Tensor{<:Number,1}) = array(T) isempty(T::Tensor) = isempty(storage(T)) - +iscu(T::Tensor) = iscu(data(T)) # # Helper functions for BlockSparse-type storage # diff --git a/src/itensor.jl b/src/itensor.jl index 6c41bfaf8e..1fa64810d4 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -1149,7 +1149,7 @@ A[i => 1, i' => 2] # 2.0, same as: A[i' => 2, i => 1] ), ) end - return NDTensors.cpu(tensor(T))[] + return tensor(T)[] end function _vals(is::Indices, I::String...) From 3fadfcdf6c8151f9b1d3d67d75afc49f294a5fce Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Sun, 8 Oct 2023 13:02:42 -0400 Subject: [PATCH 068/133] Use the [] to convert `T1` and `T2` to CPU --- NDTensors/src/dense/tensoralgebra/contract.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 8c0d4cee77..009a2c5dd6 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -233,11 +233,9 @@ function _contract_scalar_maybe_perm!( β=zero(ElR), ) where {ElR,NR} if nnz(T₁) == 1 - γ = α * NDTensors.cpu(T₁)[1] - _contract_scalar_maybe_perm!(R, labelsR, T₂, labelsT₂, γ, β) + _contract_scalar_maybe_perm!(R, labelsR, T₂, labelsT₂, α * T₁[], β) elseif nnz(T₂) == 1 - γ = α * NDTensors.cpu(T₂)[1] - _contract_scalar_maybe_perm!(R, labelsR, T₁, labelsT₁, γ, β) + _contract_scalar_maybe_perm!(R, labelsR, T₁, labelsT₁, α * T₂[], β) else error("In _contract_scalar_perm!, one tensor must be a scalar") end From 2b9ba994787a8d39ab9eeb8a77cbdc24465bef83 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Sun, 8 Oct 2023 13:03:03 -0400 Subject: [PATCH 069/133] format --- NDTensors/src/linearalgebra/linearalgebra.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index c5ecab4baf..ab52e658f5 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -423,7 +423,9 @@ matrix is unique. Returns a tuple (Q,R). """ function qr_positive(M::AbstractMatrix) if iscu(M) - throw("Currently qr positive methods do not work for CuArrays because they require scalar operations. Please convert to CPU Array or use generic qr") + throw( + "Currently qr positive methods do not work for CuArrays because they require scalar operations. Please convert to CPU Array or use generic qr", + ) end sparseQ, R = qr(M) Q = convert(typeof(R), sparseQ) @@ -451,7 +453,9 @@ matrix is unique. Returns a tuple (Q,L). """ function ql_positive(M::AbstractMatrix) if iscu(M) - throw("Currently ql positive methods do not work for CuArrays because they require scalar operations. Please convert to CPU Array or use generic ql") + throw( + "Currently ql positive methods do not work for CuArrays because they require scalar operations. Please convert to CPU Array or use generic ql", + ) end sparseQ, L = ql(M) Q = convert(typeof(L), sparseQ) From 499a616c076659c28704e774faf9345b4c08a18a Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Sun, 8 Oct 2023 13:36:15 -0400 Subject: [PATCH 070/133] Print warning because they do work when scalar allowed --- NDTensors/src/linearalgebra/linearalgebra.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index ab52e658f5..000cc8f7bd 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -423,8 +423,8 @@ matrix is unique. Returns a tuple (Q,R). """ function qr_positive(M::AbstractMatrix) if iscu(M) - throw( - "Currently qr positive methods do not work for CuArrays because they require scalar operations. Please convert to CPU Array or use generic qr", + println( + "WARNING!!: Currently qr positive methods are not efficient for CuArrays because they require scalar operations. Please convert to CPU Array or use generic qr", ) end sparseQ, R = qr(M) @@ -453,8 +453,8 @@ matrix is unique. Returns a tuple (Q,L). """ function ql_positive(M::AbstractMatrix) if iscu(M) - throw( - "Currently ql positive methods do not work for CuArrays because they require scalar operations. Please convert to CPU Array or use generic ql", + println( + "WARNING!! Currently the ql positive methods are not efficient for CuArrays because they require scalar operations. Please convert to CPU Array or use generic ql", ) end sparseQ, L = ql(M) From 86c8c3297987fa299b5da22c1701e943abe34471 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Mon, 9 Oct 2023 09:43:10 -0400 Subject: [PATCH 071/133] Use recursion to unwrap types --- NDTensors/src/abstractarray/similar.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/abstractarray/similar.jl b/NDTensors/src/abstractarray/similar.jl index 73e4e0d820..2b65a611c2 100644 --- a/NDTensors/src/abstractarray/similar.jl +++ b/NDTensors/src/abstractarray/similar.jl @@ -58,7 +58,7 @@ end # For when there are CUArray specific issues inline iscu(A::AbstractArray) = iscu(leaf_parenttype(A)) -iscu(::Type{<:AbstractArray}) = false +iscu(A::Type{<:AbstractArray}) = (leaf_parenttype(A) == A ? false : iscu(leaf_parenttype(A))) # This function actually allocates the data. # Catches conversions of dimensions specified by ranges # dimensions specified by integers with `Base.to_shape`. From 2122ae3253b13d9cc853102d981402e2ddfc92a8 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 09:46:05 -0400 Subject: [PATCH 072/133] Only call leaf_parenttype once --- NDTensors/src/abstractarray/similar.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/abstractarray/similar.jl b/NDTensors/src/abstractarray/similar.jl index 2b65a611c2..e4906b3221 100644 --- a/NDTensors/src/abstractarray/similar.jl +++ b/NDTensors/src/abstractarray/similar.jl @@ -57,7 +57,7 @@ function similar(arraytype::Type{<:AbstractArray}, dims::Tuple) end # For when there are CUArray specific issues inline -iscu(A::AbstractArray) = iscu(leaf_parenttype(A)) +iscu(A::AbstractArray) = iscu(typof(A)) iscu(A::Type{<:AbstractArray}) = (leaf_parenttype(A) == A ? false : iscu(leaf_parenttype(A))) # This function actually allocates the data. # Catches conversions of dimensions specified by ranges From 078877092006e61fa8e476bb96a32c6e98b7e467 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 09:51:45 -0400 Subject: [PATCH 073/133] Fix spelling error --- NDTensors/src/abstractarray/similar.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/abstractarray/similar.jl b/NDTensors/src/abstractarray/similar.jl index e4906b3221..ace2ff2a4d 100644 --- a/NDTensors/src/abstractarray/similar.jl +++ b/NDTensors/src/abstractarray/similar.jl @@ -57,7 +57,7 @@ function similar(arraytype::Type{<:AbstractArray}, dims::Tuple) end # For when there are CUArray specific issues inline -iscu(A::AbstractArray) = iscu(typof(A)) +iscu(A::AbstractArray) = iscu(typeof(A)) iscu(A::Type{<:AbstractArray}) = (leaf_parenttype(A) == A ? false : iscu(leaf_parenttype(A))) # This function actually allocates the data. # Catches conversions of dimensions specified by ranges From 78ecb1b51ff6f19c71064c4ae56c1886165f0ab1 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 10:18:14 -0400 Subject: [PATCH 074/133] Use iscu to just adapt CUDA computations. skip on CPU --- src/mps/dmrg.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index 57805db2cc..619f240506 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -280,7 +280,8 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) energy = vals[1] ## Right now there is a conversion problem in CUDA.jl where `UnifiedMemory` Arrays are being converted ## into `DeviceMemory`. This conversion line is here temporarily to fix that problem when it arises - phi::ITensor = itensor(convert(typeof(tensor(phi)), tensor(vecs[1]))) + ## Adapt is only called when using CUDA backend. CPU will work as implemented previously. + phi::ITensor = NDTensors.iscu(tensor(phi)) ? itensor(adapt(typeof(data(phi)), tensor(vecs[1]))) : vecs[1] #phi::ITensor = vecs[1] ortho = ha == 1 ? "left" : "right" From 2d46f5a63aaa5cedd1352a99156077563066acd3 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 10:18:50 -0400 Subject: [PATCH 075/133] format --- NDTensors/src/abstractarray/similar.jl | 4 +++- src/mps/dmrg.jl | 6 +++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/NDTensors/src/abstractarray/similar.jl b/NDTensors/src/abstractarray/similar.jl index ace2ff2a4d..873e22e44e 100644 --- a/NDTensors/src/abstractarray/similar.jl +++ b/NDTensors/src/abstractarray/similar.jl @@ -58,7 +58,9 @@ end # For when there are CUArray specific issues inline iscu(A::AbstractArray) = iscu(typeof(A)) -iscu(A::Type{<:AbstractArray}) = (leaf_parenttype(A) == A ? false : iscu(leaf_parenttype(A))) +function iscu(A::Type{<:AbstractArray}) + return (leaf_parenttype(A) == A ? false : iscu(leaf_parenttype(A))) +end # This function actually allocates the data. # Catches conversions of dimensions specified by ranges # dimensions specified by integers with `Base.to_shape`. diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index 619f240506..d1d43903d1 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -281,7 +281,11 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) ## Right now there is a conversion problem in CUDA.jl where `UnifiedMemory` Arrays are being converted ## into `DeviceMemory`. This conversion line is here temporarily to fix that problem when it arises ## Adapt is only called when using CUDA backend. CPU will work as implemented previously. - phi::ITensor = NDTensors.iscu(tensor(phi)) ? itensor(adapt(typeof(data(phi)), tensor(vecs[1]))) : vecs[1] + phi::ITensor = if NDTensors.iscu(tensor(phi)) + itensor(adapt(typeof(data(phi)), tensor(vecs[1]))) + else + vecs[1] + end #phi::ITensor = vecs[1] ortho = ha == 1 ? "left" : "right" From 3e86fe95d53f9cab149a401106a8ce79468b98eb Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 10:50:27 -0400 Subject: [PATCH 076/133] Move iscu to its own file --- NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl | 2 +- NDTensors/ext/NDTensorsCUDAExt/iscu.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) create mode 100644 NDTensors/ext/NDTensorsCUDAExt/iscu.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl index 24dc8d7ae0..5ac53aa674 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl @@ -18,8 +18,8 @@ end include("imports.jl") -iscu(::Type{<:CuArray}) = true include("set_types.jl") +include("iscu.jl") include("adapt.jl") include("linearalgebra.jl") end diff --git a/NDTensors/ext/NDTensorsCUDAExt/iscu.jl b/NDTensors/ext/NDTensorsCUDAExt/iscu.jl new file mode 100644 index 0000000000..522cddd311 --- /dev/null +++ b/NDTensors/ext/NDTensorsCUDAExt/iscu.jl @@ -0,0 +1 @@ +iscu(::Type{<:CuArray}) = true \ No newline at end of file From 10e26cc4b07abd15bab55ff9c2010a8ac34bb96a Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 11:00:22 -0400 Subject: [PATCH 077/133] Update mul!! system --- NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl | 2 +- NDTensors/src/dense/tensoralgebra/contract.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl index 93c29ee5de..767fff4135 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl @@ -12,6 +12,6 @@ function NDTensors.svd_catch_error(A::CuMatrix; alg="JacobiAlgorithm") return USV end -function NDTensors.mul!!(CM::CuArray, AM::AbstractArray, BM::AbstractArray, α, β) +function NDTensors.mul!!(::Type{<:CuArray}, CM, AM, BM, α, β) return mul!(CM, AM, BM, α, β) end diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 009a2c5dd6..dbb4598814 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -395,7 +395,7 @@ function _contract!( #tC = similar(CM) #_gemm!(tA, tB, El(α), AM, BM, El(β), CM) - mul!!(CM, AM, BM, El(α), El(β)) + CM = mul!!(leaf_parenttype(CM), CM, AM, BM, El(α), El(β)) if props.permuteC Cr = reshape(CM, props.newCrange) @@ -409,6 +409,6 @@ function _contract!( return CT end -function mul!!(CM::AbstractArray, AM::AbstractArray, BM::AbstractArray, α, β) - @strided mul!(CM, AM, BM, α, β) +function mul!!(::Type{<:AbstractArray}, CM, AM, BM, α, β) + return @strided mul!(CM, AM, BM, α, β) end From 64c2c9ab33ffc780687097f653a6abdc1a5c6d93 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 11:00:42 -0400 Subject: [PATCH 078/133] format --- NDTensors/ext/NDTensorsCUDAExt/iscu.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/iscu.jl b/NDTensors/ext/NDTensorsCUDAExt/iscu.jl index 522cddd311..5c3fc95a25 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/iscu.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/iscu.jl @@ -1 +1 @@ -iscu(::Type{<:CuArray}) = true \ No newline at end of file +iscu(::Type{<:CuArray}) = true From 394d9d88ee0b85ca1b47eb89c39ef66af5e076dc Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 15:17:47 -0400 Subject: [PATCH 079/133] Remove space --- NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl index 5ac53aa674..369c398d98 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl @@ -17,7 +17,6 @@ else end include("imports.jl") - include("set_types.jl") include("iscu.jl") include("adapt.jl") From 8668e0f979eb8a8e7d4de35a042355a725668e10 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 16:32:22 -0400 Subject: [PATCH 080/133] convert to CPU for qr_positive. And throw error if using ql with CUDA --- NDTensors/src/linearalgebra/linearalgebra.jl | 21 +++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index 000cc8f7bd..0c832bc15d 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -422,10 +422,10 @@ non-negative. Such a QR decomposition of a matrix is unique. Returns a tuple (Q,R). """ function qr_positive(M::AbstractMatrix) - if iscu(M) - println( - "WARNING!!: Currently qr positive methods are not efficient for CuArrays because they require scalar operations. Please convert to CPU Array or use generic qr", - ) + iscuda = iscu(M) + if iscuda + cutype = typeof(M) + M = NDTensors.cpu(M) end sparseQ, R = qr(M) Q = convert(typeof(R), sparseQ) @@ -440,6 +440,11 @@ function qr_positive(M::AbstractMatrix) end end end + if iscuda + M = cutype(M) + Q = cutype(Q) + R = cutype(R) + end return (Q, R) end @@ -452,11 +457,6 @@ non-negative. Such a QL decomposition of a matrix is unique. Returns a tuple (Q,L). """ function ql_positive(M::AbstractMatrix) - if iscu(M) - println( - "WARNING!! Currently the ql positive methods are not efficient for CuArrays because they require scalar operations. Please convert to CPU Array or use generic ql", - ) - end sparseQ, L = ql(M) Q = convert(typeof(L), sparseQ) nr, nc = size(L) @@ -489,6 +489,9 @@ end # about unpacking Q and L from the A matrix. # function ql!(A::StridedMatrix{<:LAPACK.BlasFloat}) + 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! From 2e9ecac9056bf18b7652edb40b77de7485880688 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 16:35:00 -0400 Subject: [PATCH 081/133] add iscu for itensors --- src/itensor.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/itensor.jl b/src/itensor.jl index 1fa64810d4..e9b980c12f 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -843,7 +843,7 @@ size(A::ITensor, d::Int) = size(tensor(A), d) _isemptyscalar(A::ITensor) = _isemptyscalar(tensor(A)) _isemptyscalar(A::Tensor) = ndims(A) == 0 && isemptystorage(A) && eltype(A) === EmptyNumber - +iscu(A::ITensor) = NDTensors.iscu(tensor(A)) """ dir(A::ITensor, i::Index) @@ -1893,7 +1893,6 @@ diag(T::ITensor) = diag(tensor(T)) mul!(C::ITensor, A::ITensor, B::ITensor, args...)::ITensor = contract!(C, A, B, args...) -## TODO this operation does not work with GPU dot(A::ITensor, B::ITensor) = NDTensors.cpu(dag(A) * B)[] inner(y::ITensor, A::ITensor, x::ITensor) = (dag(y) * A * x)[] From 8d89f7e4a7c000e0120586a34d826825ec155a87 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 16:40:59 -0400 Subject: [PATCH 082/133] define parenttype for Tensor and TensorStorage --- NDTensors/src/tensor/set_types.jl | 3 +++ NDTensors/src/tensor/tensor.jl | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/NDTensors/src/tensor/set_types.jl b/NDTensors/src/tensor/set_types.jl index 77f7130d6b..1d5df2af93 100644 --- a/NDTensors/src/tensor/set_types.jl +++ b/NDTensors/src/tensor/set_types.jl @@ -23,3 +23,6 @@ end function set_indstype(tensortype::Type{<:Tensor}, inds::Tuple) return Tensor{eltype(tensortype),length(inds),storagetype(tensortype),typeof(inds)} end + +parenttype(tensortype::Type{<:Tensor}) = storagetype(tensortype) +parenttype(storagetype::Type{<:TensorStorage}) = datatype(storagetype) \ No newline at end of file diff --git a/NDTensors/src/tensor/tensor.jl b/NDTensors/src/tensor/tensor.jl index 87613bf760..7e4af3f403 100644 --- a/NDTensors/src/tensor/tensor.jl +++ b/NDTensors/src/tensor/tensor.jl @@ -278,7 +278,7 @@ matrix(T::Tensor{<:Number,2}) = array(T) vector(T::Tensor{<:Number,1}) = array(T) isempty(T::Tensor) = isempty(storage(T)) -iscu(T::Tensor) = iscu(data(T)) + # # Helper functions for BlockSparse-type storage # From 9b1dee05f54220dc29b3bc91af46c4226a67eac7 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 17:11:14 -0400 Subject: [PATCH 083/133] call ndtensors.iscu --- src/itensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/itensor.jl b/src/itensor.jl index e9b980c12f..a56c219aad 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -843,7 +843,7 @@ size(A::ITensor, d::Int) = size(tensor(A), d) _isemptyscalar(A::ITensor) = _isemptyscalar(tensor(A)) _isemptyscalar(A::Tensor) = ndims(A) == 0 && isemptystorage(A) && eltype(A) === EmptyNumber -iscu(A::ITensor) = NDTensors.iscu(tensor(A)) +NDTensors.iscu(A::ITensor) = NDTensors.iscu(tensor(A)) """ dir(A::ITensor, i::Index) From d0557265732941a4df7e698f77b1184d743e9f15 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 17:11:29 -0400 Subject: [PATCH 084/133] Create set_types.jl for leaf_parenttype --- src/ITensors.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ITensors.jl b/src/ITensors.jl index c1effe5ebe..d571e436bf 100644 --- a/src/ITensors.jl +++ b/src/ITensors.jl @@ -138,6 +138,7 @@ include("broadcast.jl") include("tensor_operations/matrix_decomposition.jl") include("iterativesolvers.jl") include("adapt.jl") +include("set_types.jl") ##################################### # Experimental ITensor Functions From 4bb0444d89c4ea82a451d281688534275a099cab Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 17:11:39 -0400 Subject: [PATCH 085/133] create leaf_parenttype for itensor --- src/set_types.jl | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 src/set_types.jl diff --git a/src/set_types.jl b/src/set_types.jl new file mode 100644 index 0000000000..f07a8d2e4a --- /dev/null +++ b/src/set_types.jl @@ -0,0 +1,2 @@ +import NDTensors:leaf_parenttype +NDTensors.leaf_parenttype(T::ITensor) = leaf_parenttype(NDTensors.datatype(T)) From 227f7ea1e289c92dd9a623173c0e2a629bf56a67 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 17:11:59 -0400 Subject: [PATCH 086/133] Update this code per matts comments --- src/mps/dmrg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index d1d43903d1..f81d12b288 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -281,8 +281,8 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...) ## Right now there is a conversion problem in CUDA.jl where `UnifiedMemory` Arrays are being converted ## into `DeviceMemory`. This conversion line is here temporarily to fix that problem when it arises ## Adapt is only called when using CUDA backend. CPU will work as implemented previously. - phi::ITensor = if NDTensors.iscu(tensor(phi)) - itensor(adapt(typeof(data(phi)), tensor(vecs[1]))) + phi::ITensor = if NDTensors.iscu(phi) && NDTensors.iscu(vecs[1]) + adapt(set_eltype(leaf_parenttype(phi), eltype(vecs[1])), vecs[1]) else vecs[1] end From 595d5337e4f22a56058c450a107c408a007ff493 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 9 Oct 2023 17:12:30 -0400 Subject: [PATCH 087/133] format --- NDTensors/src/tensor/set_types.jl | 2 +- src/set_types.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/NDTensors/src/tensor/set_types.jl b/NDTensors/src/tensor/set_types.jl index 1d5df2af93..eb1fc548d0 100644 --- a/NDTensors/src/tensor/set_types.jl +++ b/NDTensors/src/tensor/set_types.jl @@ -25,4 +25,4 @@ function set_indstype(tensortype::Type{<:Tensor}, inds::Tuple) end parenttype(tensortype::Type{<:Tensor}) = storagetype(tensortype) -parenttype(storagetype::Type{<:TensorStorage}) = datatype(storagetype) \ No newline at end of file +parenttype(storagetype::Type{<:TensorStorage}) = datatype(storagetype) diff --git a/src/set_types.jl b/src/set_types.jl index f07a8d2e4a..a43f84d15e 100644 --- a/src/set_types.jl +++ b/src/set_types.jl @@ -1,2 +1,2 @@ -import NDTensors:leaf_parenttype +import NDTensors: leaf_parenttype NDTensors.leaf_parenttype(T::ITensor) = leaf_parenttype(NDTensors.datatype(T)) From 1631b414f95c4570eb1c4a5e7f7bdd666b340b10 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 10 Oct 2023 09:23:48 -0400 Subject: [PATCH 088/133] Launch mul!! on all 3 tensor types --- NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl | 2 +- NDTensors/src/dense/tensoralgebra/contract.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl index 767fff4135..c70f5536c1 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl @@ -12,6 +12,6 @@ function NDTensors.svd_catch_error(A::CuMatrix; alg="JacobiAlgorithm") return USV end -function NDTensors.mul!!(::Type{<:CuArray}, CM, AM, BM, α, β) +function NDTensors.mul!!(::Type{<:CuArray}, CM, ::Type{<:CuArray}, AM, ::Type{<:CuArray}, BM, α, β) return mul!(CM, AM, BM, α, β) end diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index dbb4598814..386223ac25 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -395,7 +395,7 @@ function _contract!( #tC = similar(CM) #_gemm!(tA, tB, El(α), AM, BM, El(β), CM) - CM = mul!!(leaf_parenttype(CM), CM, AM, BM, El(α), El(β)) + CM = mul!!(leaf_parenttype(CM), CM, leaf_parenttype(AM), AM, leaf_parenttype(BM), BM, El(α), El(β)) if props.permuteC Cr = reshape(CM, props.newCrange) @@ -409,6 +409,6 @@ function _contract!( return CT end -function mul!!(::Type{<:AbstractArray}, CM, AM, BM, α, β) +function mul!!(::Type{<:AbstractArray}, CM, ::Type{<:AbstractArray}, AM, ::Type{<:AbstractArray}, BM, α, β) return @strided mul!(CM, AM, BM, α, β) end From f7d5e4456301ca28e19f3cf22a6cac8d2940d22b Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 10 Oct 2023 09:30:55 -0400 Subject: [PATCH 089/133] update itensor leaf_parenttype --- src/imports.jl | 1 + src/set_types.jl | 3 +-- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/imports.jl b/src/imports.jl index 1df7f62a18..7bac316e50 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -126,6 +126,7 @@ using ITensors.NDTensors: eachdiagblock, enable_auto_fermion, fill!!, + leaf_parenttype, randn!!, set_eltype, single_precision, diff --git a/src/set_types.jl b/src/set_types.jl index a43f84d15e..8d866d7c27 100644 --- a/src/set_types.jl +++ b/src/set_types.jl @@ -1,2 +1 @@ -import NDTensors: leaf_parenttype -NDTensors.leaf_parenttype(T::ITensor) = leaf_parenttype(NDTensors.datatype(T)) +leaf_parenttype(T::ITensor) = leaf_parenttype(tensor(T)) From 015686dcb0e1b2acf7f307049a7340eb242f7d0e Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 10 Oct 2023 09:57:34 -0400 Subject: [PATCH 090/133] move leaf to import not using --- src/imports.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/imports.jl b/src/imports.jl index 7bac316e50..cbdb3839d2 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -126,7 +126,6 @@ using ITensors.NDTensors: eachdiagblock, enable_auto_fermion, fill!!, - leaf_parenttype, randn!!, set_eltype, single_precision, @@ -159,6 +158,7 @@ import ITensors.NDTensors: inds, insertblock!, insert_diag_blocks!, + leaf_parenttype, matrix, maxdim, mindim, From 379a6c5de6cb482eec2aa405ef0e5f67f499b82b Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 10 Oct 2023 10:09:38 -0400 Subject: [PATCH 091/133] Fix qr_positive conversion and convert cu to cpu for ql --- NDTensors/src/linearalgebra/linearalgebra.jl | 28 ++++++++++++++++---- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index 0c832bc15d..89a390c73a 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -424,7 +424,7 @@ matrix is unique. Returns a tuple (Q,R). function qr_positive(M::AbstractMatrix) iscuda = iscu(M) if iscuda - cutype = typeof(M) + cutype = leaf_parenttype(M) M = NDTensors.cpu(M) end sparseQ, R = qr(M) @@ -441,9 +441,8 @@ function qr_positive(M::AbstractMatrix) end end if iscuda - M = cutype(M) - Q = cutype(Q) - R = cutype(R) + Q = adapt(cutype,Q) + R = adapt(cutype,R) end return (Q, R) end @@ -457,6 +456,11 @@ non-negative. Such a QL decomposition of a matrix is unique. Returns a tuple (Q,L). """ function ql_positive(M::AbstractMatrix) + iscuda = iscu(M) + if iscuda + cutype = leaf_parenttype(M) + M = NDTensors.cpu(M) + end sparseQ, L = ql(M) Q = convert(typeof(L), sparseQ) nr, nc = size(L) @@ -470,6 +474,10 @@ function ql_positive(M::AbstractMatrix) end end end + if iscuda + Q = adapt(cutype,Q) + L = adapt(cutype,L) + end return (Q, L) end @@ -482,7 +490,17 @@ function ql(A::AbstractMatrix; kwargs...) T = eltype(A) AA = similar(A, LinearAlgebra._qreltype(T), size(A)) copyto!(AA, A) - return ql!(AA; kwargs...) + iscuda = iscu(AA) + if iscuda + cutype = leaf_parenttype(AA) + AA = NDTensors.cpu(AA) + end + Q,L = ql!(AA; kwargs...) + if iscuda + Q = adapt(cutype,Q) + L = adapt(cutype,L) + end + return (Q,L) end # # This is where the low level call to lapack actually occurs. Most of the work is From 93096dd0703706c3d6652349b8768d0580638880 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 10 Oct 2023 10:14:20 -0400 Subject: [PATCH 092/133] ql test no longer broken --- NDTensors/test/linearalgebra.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/NDTensors/test/linearalgebra.jl b/NDTensors/test/linearalgebra.jl index a643837e2d..1ffda564ab 100644 --- a/NDTensors/test/linearalgebra.jl +++ b/NDTensors/test/linearalgebra.jl @@ -51,10 +51,6 @@ devs = devices_list(copy(ARGS)) A[i, :] = A[1, :] end end - if qx == ql && dev != NDTensors.cpu - @test_broken qx(A; positive=positive) - continue - end Q, X = qx(A; positive=positive) #X is R or L. @test A ≈ Q * X atol = eps @test array(Q)' * array(Q) ≈ Id atol = eps From 8da4495d47d0553b5adb56be096f3ac0d815821d Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 10 Oct 2023 10:15:43 -0400 Subject: [PATCH 093/133] format --- NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl | 4 +++- NDTensors/src/dense/tensoralgebra/contract.jl | 15 +++++++++++++-- NDTensors/src/linearalgebra/linearalgebra.jl | 16 ++++++++-------- 3 files changed, 24 insertions(+), 11 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl index c70f5536c1..cb669a0841 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl @@ -12,6 +12,8 @@ function NDTensors.svd_catch_error(A::CuMatrix; alg="JacobiAlgorithm") return USV end -function NDTensors.mul!!(::Type{<:CuArray}, CM, ::Type{<:CuArray}, AM, ::Type{<:CuArray}, BM, α, β) +function NDTensors.mul!!( + ::Type{<:CuArray}, CM, ::Type{<:CuArray}, AM, ::Type{<:CuArray}, BM, α, β +) return mul!(CM, AM, BM, α, β) end diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 386223ac25..3a427d3be2 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -395,7 +395,9 @@ function _contract!( #tC = similar(CM) #_gemm!(tA, tB, El(α), AM, BM, El(β), CM) - CM = mul!!(leaf_parenttype(CM), CM, leaf_parenttype(AM), AM, leaf_parenttype(BM), BM, El(α), El(β)) + CM = mul!!( + leaf_parenttype(CM), CM, leaf_parenttype(AM), AM, leaf_parenttype(BM), BM, El(α), El(β) + ) if props.permuteC Cr = reshape(CM, props.newCrange) @@ -409,6 +411,15 @@ function _contract!( return CT end -function mul!!(::Type{<:AbstractArray}, CM, ::Type{<:AbstractArray}, AM, ::Type{<:AbstractArray}, BM, α, β) +function mul!!( + ::Type{<:AbstractArray}, + CM, + ::Type{<:AbstractArray}, + AM, + ::Type{<:AbstractArray}, + BM, + α, + β, +) return @strided mul!(CM, AM, BM, α, β) end diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index 89a390c73a..578d8d27b2 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -441,8 +441,8 @@ function qr_positive(M::AbstractMatrix) end end if iscuda - Q = adapt(cutype,Q) - R = adapt(cutype,R) + Q = adapt(cutype, Q) + R = adapt(cutype, R) end return (Q, R) end @@ -475,8 +475,8 @@ function ql_positive(M::AbstractMatrix) end end if iscuda - Q = adapt(cutype,Q) - L = adapt(cutype,L) + Q = adapt(cutype, Q) + L = adapt(cutype, L) end return (Q, L) end @@ -495,12 +495,12 @@ function ql(A::AbstractMatrix; kwargs...) cutype = leaf_parenttype(AA) AA = NDTensors.cpu(AA) end - Q,L = ql!(AA; kwargs...) + Q, L = ql!(AA; kwargs...) if iscuda - Q = adapt(cutype,Q) - L = adapt(cutype,L) + Q = adapt(cutype, Q) + L = adapt(cutype, L) end - return (Q,L) + return (Q, L) end # # This is where the low level call to lapack actually occurs. Most of the work is From e2f83e838950da54c285c5333f5c0b535cbb77d3 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 10 Oct 2023 11:37:42 -0400 Subject: [PATCH 094/133] create a permutedims!! function and implement in code --- .../ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl | 1 + NDTensors/ext/NDTensorsCUDAExt/permutedims.jl | 3 +++ NDTensors/src/NDTensors.jl | 1 + NDTensors/src/abstractarray/permutedims.jl | 3 +++ NDTensors/src/arraytensor/array.jl | 2 +- NDTensors/src/dense/densetensor.jl | 5 ++--- NDTensors/src/dense/tensoralgebra/contract.jl | 18 ++++++++---------- 7 files changed, 19 insertions(+), 14 deletions(-) create mode 100644 NDTensors/ext/NDTensorsCUDAExt/permutedims.jl create mode 100644 NDTensors/src/abstractarray/permutedims.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl index 369c398d98..cb9b1d0bbc 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl @@ -21,4 +21,5 @@ include("set_types.jl") include("iscu.jl") include("adapt.jl") include("linearalgebra.jl") +include("permutedims.jl") end diff --git a/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl b/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl new file mode 100644 index 0000000000..eeca2441f0 --- /dev/null +++ b/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl @@ -0,0 +1,3 @@ +function NDTensors.permutedims!!(::Type{<:CuArray}, M, perm) + return permutedims(M, perm) +end \ No newline at end of file diff --git a/NDTensors/src/NDTensors.jl b/NDTensors/src/NDTensors.jl index 8445b925cb..f742e331a1 100644 --- a/NDTensors/src/NDTensors.jl +++ b/NDTensors/src/NDTensors.jl @@ -50,6 +50,7 @@ include("abstractarray/set_types.jl") include("abstractarray/to_shape.jl") include("abstractarray/similar.jl") include("abstractarray/ndims.jl") +include("abstractarray/permutedims.jl") include("abstractarray/fill.jl") include("array/set_types.jl") include("tupletools.jl") diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl new file mode 100644 index 0000000000..6fdd9b322b --- /dev/null +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -0,0 +1,3 @@ +function permutedims!!(::Type{<:AbstractArray}, M, perm) + return @strided Mdest = permutedims(M, perm) +end \ No newline at end of file diff --git a/NDTensors/src/arraytensor/array.jl b/NDTensors/src/arraytensor/array.jl index 40e89c5181..c6e1b2b82f 100644 --- a/NDTensors/src/arraytensor/array.jl +++ b/NDTensors/src/arraytensor/array.jl @@ -61,6 +61,6 @@ end function permutedims!( output_array::MatrixOrArrayStorage, array::MatrixOrArrayStorage, perm, f::Function ) - @strided output_array .= f.(output_array, permutedims(array, perm)) + output_array .= f.(output_array, permutedims!!(leaf_parenttype(array),array, perm)) return output_array end diff --git a/NDTensors/src/dense/densetensor.jl b/NDTensors/src/dense/densetensor.jl index fd0f72bbd0..0743ca4ae0 100644 --- a/NDTensors/src/dense/densetensor.jl +++ b/NDTensors/src/dense/densetensor.jl @@ -199,7 +199,7 @@ function permutedims!( ) where {N,StoreT<:StridedArray} RA = array(R) TA = array(T) - @strided RA .= permutedims(TA, perm) + RA .= permutedims!!(leaf_parenttype(TA), TA, perm) return R end @@ -247,8 +247,7 @@ function permutedims!( end RA = array(R) TA = array(T) - #@strided - RA .= f.(RA, permutedims(TA, perm)) + RA .= f.(RA, permutedims!!(leaf_parenttype(TA), TA, perm)) return R end diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 3a427d3be2..44b6787992 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -175,14 +175,14 @@ function _contract_scalar_perm!( if iszero(α) fill!(Rᵃ, 0) else - @strided Rᵃ .= α .* permutedims(Tᵃ, perm) + Rᵃ .= α .* permutedims!!(leaf_parenttype(Tᵃ), Tᵃ, perm) end elseif isone(β) if iszero(α) # Rᵃ .= Rᵃ # No-op else - @strided Rᵃ .= α .* permutedims(Tᵃ, perm) .+ Rᵃ + Rᵃ .= α .* permutedims!!(leaf_parenttype(Tᵃ), Tᵃ, perm) .+ Rᵃ end else if iszero(α) @@ -343,10 +343,10 @@ function _contract!( β::Number=zero(El), ) where {El,NC,NA,NB} tA = 'N' + t = TimerOutput() if props.permuteA #@timeit_debug timer "_contract!: permutedims A" begin - #@strided - Ap = permutedims(AT, props.PA) + Ap = permutedims!!(leaf_parenttype(AT), AT, props.PA) #end # @timeit AM = transpose(reshape(Ap, (props.dmid, props.dleft))) else @@ -361,8 +361,7 @@ function _contract!( tB = 'N' if props.permuteB #@timeit_debug timer "_contract!: permutedims B" begin - #@strided - Bp = permutedims(BT, props.PB) + Bp = permutedims!!(leaf_parenttype(BT), BT, props.PB) #end # @timeit BM = reshape(Bp, (props.dmid, props.dright)) else @@ -379,7 +378,7 @@ function _contract!( # we need to make sure C is permuted to the same # ordering as A B which is the inverse of props.PC if β ≠ 0 - CM = reshape(permutedims(CT, invperm(props.PC)), (props.dleft, props.dright)) + CM = reshape(permutedims!!(leaf_parenttype(CT), CT, invperm(props.PC)), (props.dleft, props.dright)) else # Need to copy here since we will be permuting # into C later @@ -402,9 +401,8 @@ function _contract!( if props.permuteC Cr = reshape(CM, props.newCrange) # TODO: use invperm(pC) here? - #@timeit_debug timer "_contract!: permutedims C" begin - #@strided - CT .= permutedims(Cr, props.PC) + #@timeit_debug timer "_contract!: permutedims C" begin + CT .= permutedims!!(leaf_parenttype(Cr), Cr, props.PC) #end # @timeit end From 3c6f073726493964921771be53a95a2887be92b1 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Tue, 10 Oct 2023 11:42:03 -0400 Subject: [PATCH 095/133] format --- NDTensors/ext/NDTensorsCUDAExt/permutedims.jl | 2 +- NDTensors/src/abstractarray/permutedims.jl | 2 +- NDTensors/src/arraytensor/array.jl | 2 +- NDTensors/src/dense/tensoralgebra/contract.jl | 5 ++++- 4 files changed, 7 insertions(+), 4 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl b/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl index eeca2441f0..b7830cc576 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl @@ -1,3 +1,3 @@ function NDTensors.permutedims!!(::Type{<:CuArray}, M, perm) return permutedims(M, perm) -end \ No newline at end of file +end diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index 6fdd9b322b..4bf1e39dd8 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -1,3 +1,3 @@ function permutedims!!(::Type{<:AbstractArray}, M, perm) return @strided Mdest = permutedims(M, perm) -end \ No newline at end of file +end diff --git a/NDTensors/src/arraytensor/array.jl b/NDTensors/src/arraytensor/array.jl index c6e1b2b82f..c2950b2c20 100644 --- a/NDTensors/src/arraytensor/array.jl +++ b/NDTensors/src/arraytensor/array.jl @@ -61,6 +61,6 @@ end function permutedims!( output_array::MatrixOrArrayStorage, array::MatrixOrArrayStorage, perm, f::Function ) - output_array .= f.(output_array, permutedims!!(leaf_parenttype(array),array, perm)) + output_array .= f.(output_array, permutedims!!(leaf_parenttype(array), array, perm)) return output_array end diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 44b6787992..8b764c1913 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -378,7 +378,10 @@ function _contract!( # we need to make sure C is permuted to the same # ordering as A B which is the inverse of props.PC if β ≠ 0 - CM = reshape(permutedims!!(leaf_parenttype(CT), CT, invperm(props.PC)), (props.dleft, props.dright)) + CM = reshape( + permutedims!!(leaf_parenttype(CT), CT, invperm(props.PC)), + (props.dleft, props.dright), + ) else # Need to copy here since we will be permuting # into C later From 3aecf80442273d90a88b1289ec5293bec5b83107 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 11 Oct 2023 15:14:46 -0400 Subject: [PATCH 096/133] Remove cuda permutedims --- NDTensors/ext/NDTensorsCUDAExt/permutedims.jl | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 NDTensors/ext/NDTensorsCUDAExt/permutedims.jl diff --git a/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl b/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl deleted file mode 100644 index b7830cc576..0000000000 --- a/NDTensors/ext/NDTensorsCUDAExt/permutedims.jl +++ /dev/null @@ -1,3 +0,0 @@ -function NDTensors.permutedims!!(::Type{<:CuArray}, M, perm) - return permutedims(M, perm) -end From 5634b1bfc79b1d329eb22e8ec1759468a73bc1ff Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 11 Oct 2023 15:15:11 -0400 Subject: [PATCH 097/133] More robust permutedims functions --- NDTensors/src/abstractarray/permutedims.jl | 24 ++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index 4bf1e39dd8..dbee912311 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -1,3 +1,23 @@ -function permutedims!!(::Type{<:AbstractArray}, M, perm) - return @strided Mdest = permutedims(M, perm) +function permutedims(::Type{<:Array}, M, perm) + return @strided Mdest = Base.permutedims(M, perm) end + +function permutedims!(::Type{<:Array}, Mdest, ::Type{<:Array}, M, perm) + return @strided Mdest .= Base.permutedims(M, perm) +end + +function permutedims!!(::Type{<:Array}, B, ::Type{<:Array}, A, perm, f) + @strided B .= f.(B, Base.permutedims(A, perm)) +end + +function NDTensors.permutedims(::Type{<:AbstractArray}, M, perm) + return Base.permutedims(M, perm) +end + +function permutedims!(::Type{<:AbstractArray}, Mdest, ::Type{<:AbstractArray}, M, perm) + return Mdest .= Base.permutedims(M, perm) +end + +function permutedims!!(::Type{<:AbstractArray}, B, ::Type{<:AbstractArray}, A, perm, f) + B .= f.(B, Base.permutedims(A, perm)) +end \ No newline at end of file From be69b1308530ab0980f613a556e8407b4b375ae7 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 11 Oct 2023 15:19:19 -0400 Subject: [PATCH 098/133] Add a comment/idea --- NDTensors/src/abstractarray/permutedims.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index dbee912311..9407ec02e5 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -2,6 +2,7 @@ function permutedims(::Type{<:Array}, M, perm) return @strided Mdest = Base.permutedims(M, perm) end +## TODO is it possible to do something like f = identity and remove the permutedims!! funciton? function permutedims!(::Type{<:Array}, Mdest, ::Type{<:Array}, M, perm) return @strided Mdest .= Base.permutedims(M, perm) end From 86bd2943f7ef318684aa558fc5d0b05d459a9cd4 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 11 Oct 2023 15:21:20 -0400 Subject: [PATCH 099/133] remove include --- NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl index cb9b1d0bbc..369c398d98 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/NDTensorsCUDAExt.jl @@ -21,5 +21,4 @@ include("set_types.jl") include("iscu.jl") include("adapt.jl") include("linearalgebra.jl") -include("permutedims.jl") end From 0dd361c5eea135eed3405042aa799b4812e00b84 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 11 Oct 2023 15:42:24 -0400 Subject: [PATCH 100/133] update permutedim calls --- NDTensors/src/arraytensor/array.jl | 2 +- NDTensors/src/dense/densetensor.jl | 5 ++--- NDTensors/src/dense/tensoralgebra/contract.jl | 12 ++++++------ 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/NDTensors/src/arraytensor/array.jl b/NDTensors/src/arraytensor/array.jl index c2950b2c20..8abe19d74c 100644 --- a/NDTensors/src/arraytensor/array.jl +++ b/NDTensors/src/arraytensor/array.jl @@ -61,6 +61,6 @@ end function permutedims!( output_array::MatrixOrArrayStorage, array::MatrixOrArrayStorage, perm, f::Function ) - output_array .= f.(output_array, permutedims!!(leaf_parenttype(array), array, perm)) + output_array = permutedims!!(leaf_parenttype(output_array), output_array, leaf_parenttype(array), array, perm, f) return output_array end diff --git a/NDTensors/src/dense/densetensor.jl b/NDTensors/src/dense/densetensor.jl index 0743ca4ae0..d3d5582fcb 100644 --- a/NDTensors/src/dense/densetensor.jl +++ b/NDTensors/src/dense/densetensor.jl @@ -199,7 +199,7 @@ function permutedims!( ) where {N,StoreT<:StridedArray} RA = array(R) TA = array(T) - RA .= permutedims!!(leaf_parenttype(TA), TA, perm) + RA .= permutedims!(leaf_parenttype(RA), RA, leaf_parenttype(TA), TA, perm) return R end @@ -247,8 +247,7 @@ function permutedims!( end RA = array(R) TA = array(T) - RA .= f.(RA, permutedims!!(leaf_parenttype(TA), TA, perm)) - return R + return permutedims!!(leaf_parenttype(RA), RA, leaf_parenttype(TA), TA, perm, f) end """ diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 8b764c1913..6b1ceef1e8 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -175,14 +175,14 @@ function _contract_scalar_perm!( if iszero(α) fill!(Rᵃ, 0) else - Rᵃ .= α .* permutedims!!(leaf_parenttype(Tᵃ), Tᵃ, perm) + Rᵃ = permutedims!!(leaf_parenttype(Rᵃ), Rᵃ, leaf_parenttype(Tᵃ), Tᵃ, perm, (r,t) -> *(α, t)) end elseif isone(β) if iszero(α) # Rᵃ .= Rᵃ # No-op else - Rᵃ .= α .* permutedims!!(leaf_parenttype(Tᵃ), Tᵃ, perm) .+ Rᵃ + Rᵃ = permutedims!!(leaf_parenttype(Rᵃ), Rᵃ, leaf_parenttype(Tᵃ), Tᵃ, perm, (r,t) -> r + α * t) end else if iszero(α) @@ -346,7 +346,7 @@ function _contract!( t = TimerOutput() if props.permuteA #@timeit_debug timer "_contract!: permutedims A" begin - Ap = permutedims!!(leaf_parenttype(AT), AT, props.PA) + Ap = permutedims(leaf_parenttype(AT), AT, props.PA) #end # @timeit AM = transpose(reshape(Ap, (props.dmid, props.dleft))) else @@ -361,7 +361,7 @@ function _contract!( tB = 'N' if props.permuteB #@timeit_debug timer "_contract!: permutedims B" begin - Bp = permutedims!!(leaf_parenttype(BT), BT, props.PB) + Bp = permutedims(leaf_parenttype(BT), BT, props.PB) #end # @timeit BM = reshape(Bp, (props.dmid, props.dright)) else @@ -379,7 +379,7 @@ function _contract!( # ordering as A B which is the inverse of props.PC if β ≠ 0 CM = reshape( - permutedims!!(leaf_parenttype(CT), CT, invperm(props.PC)), + permutedims(leaf_parenttype(CT), CT, invperm(props.PC)), (props.dleft, props.dright), ) else @@ -405,7 +405,7 @@ function _contract!( Cr = reshape(CM, props.newCrange) # TODO: use invperm(pC) here? #@timeit_debug timer "_contract!: permutedims C" begin - CT .= permutedims!!(leaf_parenttype(Cr), Cr, props.PC) + CT .= permutedims(leaf_parenttype(Cr), Cr, props.PC) #end # @timeit end From b08342d75128d055c6f6f94657940a5023ccc27e Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 11 Oct 2023 15:43:37 -0400 Subject: [PATCH 101/133] Make abstract mul!! call, not CuArray mul!! --- NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl | 6 ------ NDTensors/src/dense/tensoralgebra/contract.jl | 12 +++++++++--- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl index cb669a0841..c144a4a802 100644 --- a/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl +++ b/NDTensors/ext/NDTensorsCUDAExt/linearalgebra.jl @@ -11,9 +11,3 @@ function NDTensors.svd_catch_error(A::CuMatrix; alg="JacobiAlgorithm") end return USV end - -function NDTensors.mul!!( - ::Type{<:CuArray}, CM, ::Type{<:CuArray}, AM, ::Type{<:CuArray}, BM, α, β -) - return mul!(CM, AM, BM, α, β) -end diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 6b1ceef1e8..9a1982aae5 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -413,14 +413,20 @@ function _contract!( end function mul!!( - ::Type{<:AbstractArray}, + ::Type{<:Array}, CM, - ::Type{<:AbstractArray}, + ::Type{<:Array}, AM, - ::Type{<:AbstractArray}, + ::Type{<:Array}, BM, α, β, ) return @strided mul!(CM, AM, BM, α, β) end + +function NDTensors.mul!!( + ::Type{<:AbstractArray}, CM, ::Type{<:AbstractArray}, AM, ::Type{<:AbstractArray}, BM, α, β +) + return mul!(CM, AM, BM, α, β) +end From 383e51f5d9cf6a688865beb3648385c6d8b15dac Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Wed, 11 Oct 2023 15:44:02 -0400 Subject: [PATCH 102/133] format --- NDTensors/src/abstractarray/permutedims.jl | 4 +-- NDTensors/src/arraytensor/array.jl | 4 ++- NDTensors/src/dense/tensoralgebra/contract.jl | 29 ++++++++++--------- 3 files changed, 20 insertions(+), 17 deletions(-) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index 9407ec02e5..a23d6f89ac 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -20,5 +20,5 @@ function permutedims!(::Type{<:AbstractArray}, Mdest, ::Type{<:AbstractArray}, M end function permutedims!!(::Type{<:AbstractArray}, B, ::Type{<:AbstractArray}, A, perm, f) - B .= f.(B, Base.permutedims(A, perm)) -end \ No newline at end of file + return B .= f.(B, Base.permutedims(A, perm)) +end diff --git a/NDTensors/src/arraytensor/array.jl b/NDTensors/src/arraytensor/array.jl index 8abe19d74c..997204394d 100644 --- a/NDTensors/src/arraytensor/array.jl +++ b/NDTensors/src/arraytensor/array.jl @@ -61,6 +61,8 @@ end function permutedims!( output_array::MatrixOrArrayStorage, array::MatrixOrArrayStorage, perm, f::Function ) - output_array = permutedims!!(leaf_parenttype(output_array), output_array, leaf_parenttype(array), array, perm, f) + output_array = permutedims!!( + leaf_parenttype(output_array), output_array, leaf_parenttype(array), array, perm, f + ) return output_array end diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 9a1982aae5..9d2e012d52 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -175,14 +175,18 @@ function _contract_scalar_perm!( if iszero(α) fill!(Rᵃ, 0) else - Rᵃ = permutedims!!(leaf_parenttype(Rᵃ), Rᵃ, leaf_parenttype(Tᵃ), Tᵃ, perm, (r,t) -> *(α, t)) + Rᵃ = permutedims!!( + leaf_parenttype(Rᵃ), Rᵃ, leaf_parenttype(Tᵃ), Tᵃ, perm, (r, t) -> *(α, t) + ) end elseif isone(β) if iszero(α) # Rᵃ .= Rᵃ # No-op else - Rᵃ = permutedims!!(leaf_parenttype(Rᵃ), Rᵃ, leaf_parenttype(Tᵃ), Tᵃ, perm, (r,t) -> r + α * t) + Rᵃ = permutedims!!( + leaf_parenttype(Rᵃ), Rᵃ, leaf_parenttype(Tᵃ), Tᵃ, perm, (r, t) -> r + α * t + ) end else if iszero(α) @@ -379,8 +383,7 @@ function _contract!( # ordering as A B which is the inverse of props.PC if β ≠ 0 CM = reshape( - permutedims(leaf_parenttype(CT), CT, invperm(props.PC)), - (props.dleft, props.dright), + permutedims(leaf_parenttype(CT), CT, invperm(props.PC)), (props.dleft, props.dright) ) else # Need to copy here since we will be permuting @@ -412,21 +415,19 @@ function _contract!( return CT end -function mul!!( - ::Type{<:Array}, +function mul!!(::Type{<:Array}, CM, ::Type{<:Array}, AM, ::Type{<:Array}, BM, α, β) + return @strided mul!(CM, AM, BM, α, β) +end + +function NDTensors.mul!!( + ::Type{<:AbstractArray}, CM, - ::Type{<:Array}, + ::Type{<:AbstractArray}, AM, - ::Type{<:Array}, + ::Type{<:AbstractArray}, BM, α, β, -) - return @strided mul!(CM, AM, BM, α, β) -end - -function NDTensors.mul!!( - ::Type{<:AbstractArray}, CM, ::Type{<:AbstractArray}, AM, ::Type{<:AbstractArray}, BM, α, β ) return mul!(CM, AM, BM, α, β) end From 7ce845bb7687fdb112e66cfbc1696cf51e4b8de3 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 09:36:23 -0400 Subject: [PATCH 103/133] Move array specif functions to permutedims --- NDTensors/src/abstractarray/permutedims.jl | 13 ------------- NDTensors/src/array/permutedims.jl | 12 ++++++++++++ 2 files changed, 12 insertions(+), 13 deletions(-) create mode 100644 NDTensors/src/array/permutedims.jl diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index a23d6f89ac..eec998ca4c 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -1,16 +1,3 @@ -function permutedims(::Type{<:Array}, M, perm) - return @strided Mdest = Base.permutedims(M, perm) -end - -## TODO is it possible to do something like f = identity and remove the permutedims!! funciton? -function permutedims!(::Type{<:Array}, Mdest, ::Type{<:Array}, M, perm) - return @strided Mdest .= Base.permutedims(M, perm) -end - -function permutedims!!(::Type{<:Array}, B, ::Type{<:Array}, A, perm, f) - @strided B .= f.(B, Base.permutedims(A, perm)) -end - function NDTensors.permutedims(::Type{<:AbstractArray}, M, perm) return Base.permutedims(M, perm) end diff --git a/NDTensors/src/array/permutedims.jl b/NDTensors/src/array/permutedims.jl new file mode 100644 index 0000000000..0c381fdda5 --- /dev/null +++ b/NDTensors/src/array/permutedims.jl @@ -0,0 +1,12 @@ +function permutedims(::Type{<:Array}, M, perm) + return @strided Mdest = Base.permutedims(M, perm) +end + +## TODO is it possible to do something like f = identity and remove the permutedims!! funciton? +function permutedims!(::Type{<:Array}, Mdest, ::Type{<:Array}, M, perm) + return @strided Mdest .= Base.permutedims(M, perm) +end + +function permutedims!!(::Type{<:Array}, B, ::Type{<:Array}, A, perm, f) + @strided B .= f.(B, Base.permutedims(A, perm)) +end \ No newline at end of file From e4afbadd32afaa1b9b2567027377ed147171dc64 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 09:41:10 -0400 Subject: [PATCH 104/133] Make a choke point function for permutedims --- NDTensors/src/abstractarray/permutedims.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index eec998ca4c..cb7c0052ae 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -6,6 +6,10 @@ function permutedims!(::Type{<:AbstractArray}, Mdest, ::Type{<:AbstractArray}, M return Mdest .= Base.permutedims(M, perm) end +function permutedims!!(B::AbstractArray, A::AbstractArray, perm, f) + return permutedims!!(leaf_parenttype(B), B, leaf_parenttype(A), A, perm, f) +end + function permutedims!!(::Type{<:AbstractArray}, B, ::Type{<:AbstractArray}, A, perm, f) return B .= f.(B, Base.permutedims(A, perm)) end From b06606511ce38ff40b1d949295a82bdbde8f9452 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 09:53:45 -0400 Subject: [PATCH 105/133] Add array/permutedims.jl --- NDTensors/src/NDTensors.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/NDTensors/src/NDTensors.jl b/NDTensors/src/NDTensors.jl index f742e331a1..c3eb488da7 100644 --- a/NDTensors/src/NDTensors.jl +++ b/NDTensors/src/NDTensors.jl @@ -53,6 +53,7 @@ include("abstractarray/ndims.jl") include("abstractarray/permutedims.jl") include("abstractarray/fill.jl") include("array/set_types.jl") +include("array/permutedims.jl") include("tupletools.jl") include("emptynumber.jl") include("nodata.jl") From d9caa442dd6423de5f6a14711a19a4eaf95d09ae Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 09:54:02 -0400 Subject: [PATCH 106/133] create mul.jl functions --- NDTensors/src/NDTensors.jl | 2 ++ NDTensors/src/abstractarray/mul.jl | 12 ++++++++++++ NDTensors/src/array/mul.jl | 3 +++ 3 files changed, 17 insertions(+) create mode 100644 NDTensors/src/abstractarray/mul.jl create mode 100644 NDTensors/src/array/mul.jl diff --git a/NDTensors/src/NDTensors.jl b/NDTensors/src/NDTensors.jl index c3eb488da7..9e10d33711 100644 --- a/NDTensors/src/NDTensors.jl +++ b/NDTensors/src/NDTensors.jl @@ -52,8 +52,10 @@ include("abstractarray/similar.jl") include("abstractarray/ndims.jl") include("abstractarray/permutedims.jl") include("abstractarray/fill.jl") +include("abstractarray/mul.jl") include("array/set_types.jl") include("array/permutedims.jl") +include("array/mul.jl") include("tupletools.jl") include("emptynumber.jl") include("nodata.jl") diff --git a/NDTensors/src/abstractarray/mul.jl b/NDTensors/src/abstractarray/mul.jl new file mode 100644 index 0000000000..52b2652182 --- /dev/null +++ b/NDTensors/src/abstractarray/mul.jl @@ -0,0 +1,12 @@ +function NDTensors.mul!!( + ::Type{<:AbstractArray}, + CM, + ::Type{<:AbstractArray}, + AM, + ::Type{<:AbstractArray}, + BM, + α, + β, +) + return mul!(CM, AM, BM, α, β) +end \ No newline at end of file diff --git a/NDTensors/src/array/mul.jl b/NDTensors/src/array/mul.jl new file mode 100644 index 0000000000..a8c90e5b94 --- /dev/null +++ b/NDTensors/src/array/mul.jl @@ -0,0 +1,3 @@ +function mul!!(::Type{<:Array}, CM, ::Type{<:Array}, AM, ::Type{<:Array}, BM, α, β) + return @strided mul!(CM, AM, BM, α, β) +end \ No newline at end of file From 573500ccb7e506d65ec0aea04b85ef972f0eb397 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 10:01:29 -0400 Subject: [PATCH 107/133] Move abstract array functions out of densetensor --- NDTensors/src/abstractarray/mul.jl | 6 +- .../abstractarray/tensoralgebra/contract.jl | 122 +++++++++++++++ NDTensors/src/dense/tensoralgebra/contract.jl | 146 ------------------ 3 files changed, 127 insertions(+), 147 deletions(-) create mode 100644 NDTensors/src/abstractarray/tensoralgebra/contract.jl diff --git a/NDTensors/src/abstractarray/mul.jl b/NDTensors/src/abstractarray/mul.jl index 52b2652182..203e69ff30 100644 --- a/NDTensors/src/abstractarray/mul.jl +++ b/NDTensors/src/abstractarray/mul.jl @@ -1,4 +1,8 @@ -function NDTensors.mul!!( +function mul!!(CM::AbstractArray, AM::AbstractArray, BM::AbstractArray, α, β) + return mul!!(leaf_parenttype(CM), CM, leaf_parenttype(AM), AM, leaf_parenttype(BM), BM, α, β) +end + +function mul!!( ::Type{<:AbstractArray}, CM, ::Type{<:AbstractArray}, diff --git a/NDTensors/src/abstractarray/tensoralgebra/contract.jl b/NDTensors/src/abstractarray/tensoralgebra/contract.jl new file mode 100644 index 0000000000..64fd9d1fa4 --- /dev/null +++ b/NDTensors/src/abstractarray/tensoralgebra/contract.jl @@ -0,0 +1,122 @@ +# generic matmul +function _gemm!( + ::GemmBackend{:Generic}, + tA, + tB, + alpha::AT, + A::AbstractVecOrMat, + B::AbstractVecOrMat, + beta::BT, + C::AbstractVecOrMat, +) where {AT,BT} + mul!(C, tA == 'T' ? transpose(A) : A, tB == 'T' ? transpose(B) : B, alpha, beta) + return C +end + +function contraction_output(tensor1::DenseTensor, tensor2::DenseTensor, indsR) + tensortypeR = contraction_output_type(typeof(tensor1), typeof(tensor2), indsR) + return NDTensors.similar(tensortypeR, indsR) +end + +# Non-trivial permutation +function _contract_scalar_perm!( + Rᵃ::AbstractArray{ElR}, Tᵃ::AbstractArray, perm, α, β=zero(ElR) +) where {ElR} + if iszero(β) + if iszero(α) + fill!(Rᵃ, 0) + else + Rᵃ = permutedims!!(Rᵃ, Tᵃ, perm, (r, t) -> α * t) + end + elseif isone(β) + if iszero(α) + # Rᵃ .= Rᵃ + # No-op + else + Rᵃ = permutedims!!(Rᵃ, Tᵃ, perm, (r, t) -> r + α * t) + end + else + if iszero(α) + # Rᵃ .= β .* Rᵃ + LinearAlgebra.scal!(length(Rᵃ), β, Rᵃ, 1) + else + Rᵃ .= α .* permutedims(Tᵃ, perm) .+ β .* Rᵃ + end + end + return Rᵃ +end + +function _contract!( + CT::AbstractArray{El,NC}, + AT::AbstractArray{El,NA}, + BT::AbstractArray{El,NB}, + props::ContractionProperties, + α::Number=one(El), + β::Number=zero(El), +) where {El,NC,NA,NB} + tA = 'N' + if props.permuteA + #@timeit_debug timer "_contract!: permutedims A" begin + Ap = permutedims(leaf_parenttype(AT), AT, props.PA) + #end # @timeit + AM = transpose(reshape(Ap, (props.dmid, props.dleft))) + else + #A doesn't have to be permuted + if Atrans(props) + AM = transpose(reshape(AT, (props.dmid, props.dleft))) + else + AM = reshape(AT, (props.dleft, props.dmid)) + end + end + + tB = 'N' + if props.permuteB + #@timeit_debug timer "_contract!: permutedims B" begin + Bp = permutedims(leaf_parenttype(BT), BT, props.PB) + #end # @timeit + BM = reshape(Bp, (props.dmid, props.dright)) + else + if Btrans(props) + BM = transpose(reshape(BT, (props.dright, props.dmid))) + else + BM = reshape(BT, (props.dmid, props.dright)) + end + end + + # TODO: this logic may be wrong + if props.permuteC + # if we are computing C = α * A B + β * C + # we need to make sure C is permuted to the same + # ordering as A B which is the inverse of props.PC + if β ≠ 0 + CM = reshape( + permutedims(leaf_parenttype(CT), CT, invperm(props.PC)), (props.dleft, props.dright) + ) + else + # Need to copy here since we will be permuting + # into C later + CM = reshape(copy(CT), (props.dleft, props.dright)) + end + else + if Ctrans(props) + CM = transpose(reshape(CT, (props.dright, props.dleft))) + else + CM = reshape(CT, (props.dleft, props.dright)) + end + end + + #tC = similar(CM) + #_gemm!(tA, tB, El(α), AM, BM, El(β), CM) + CM = mul!!(CM, AM, BM, El(α), El(β) + ) + + if props.permuteC + Cr = reshape(CM, props.newCrange) + # TODO: use invperm(pC) here? + #@timeit_debug timer "_contract!: permutedims C" begin + CT .= permutedims(leaf_parenttype(Cr), Cr, props.PC) + #end # @timeit + end + + return CT +end \ No newline at end of file diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 9d2e012d52..6857828705 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -59,26 +59,6 @@ function _gemm!( #end # @timeit end -# generic matmul -function _gemm!( - ::GemmBackend{:Generic}, - tA, - tB, - alpha::AT, - A::AbstractVecOrMat, - B::AbstractVecOrMat, - beta::BT, - C::AbstractVecOrMat, -) where {AT,BT} - mul!(C, tA == 'T' ? transpose(A) : A, tB == 'T' ? transpose(B) : B, alpha, beta) - return C -end - -function contraction_output(tensor1::DenseTensor, tensor2::DenseTensor, indsR) - tensortypeR = contraction_output_type(typeof(tensor1), typeof(tensor2), indsR) - return NDTensors.similar(tensortypeR, indsR) -end - # Both are scalar-like tensors function _contract_scalar!( R::DenseTensor{ElR}, @@ -167,38 +147,6 @@ function _contract_scalar_noperm!( return R end -# Non-trivial permutation -function _contract_scalar_perm!( - Rᵃ::AbstractArray{ElR}, Tᵃ::AbstractArray, perm, α, β=zero(ElR) -) where {ElR} - if iszero(β) - if iszero(α) - fill!(Rᵃ, 0) - else - Rᵃ = permutedims!!( - leaf_parenttype(Rᵃ), Rᵃ, leaf_parenttype(Tᵃ), Tᵃ, perm, (r, t) -> *(α, t) - ) - end - elseif isone(β) - if iszero(α) - # Rᵃ .= Rᵃ - # No-op - else - Rᵃ = permutedims!!( - leaf_parenttype(Rᵃ), Rᵃ, leaf_parenttype(Tᵃ), Tᵃ, perm, (r, t) -> r + α * t - ) - end - else - if iszero(α) - # Rᵃ .= β .* Rᵃ - LinearAlgebra.scal!(length(Rᵃ), β, Rᵃ, 1) - else - Rᵃ .= α .* permutedims(Tᵃ, perm) .+ β .* Rᵃ - end - end - return Rᵃ -end - function _contract_scalar_maybe_perm!( ::Order{N}, R::DenseTensor{ElR,NR}, labelsR, T::DenseTensor, labelsT, α, β=zero(ElR) ) where {ElR,NR,N} @@ -337,97 +285,3 @@ function _contract!( return _contract!(C, A, B, props, α, β) end - -function _contract!( - CT::AbstractArray{El,NC}, - AT::AbstractArray{El,NA}, - BT::AbstractArray{El,NB}, - props::ContractionProperties, - α::Number=one(El), - β::Number=zero(El), -) where {El,NC,NA,NB} - tA = 'N' - t = TimerOutput() - if props.permuteA - #@timeit_debug timer "_contract!: permutedims A" begin - Ap = permutedims(leaf_parenttype(AT), AT, props.PA) - #end # @timeit - AM = transpose(reshape(Ap, (props.dmid, props.dleft))) - else - #A doesn't have to be permuted - if Atrans(props) - AM = transpose(reshape(AT, (props.dmid, props.dleft))) - else - AM = reshape(AT, (props.dleft, props.dmid)) - end - end - - tB = 'N' - if props.permuteB - #@timeit_debug timer "_contract!: permutedims B" begin - Bp = permutedims(leaf_parenttype(BT), BT, props.PB) - #end # @timeit - BM = reshape(Bp, (props.dmid, props.dright)) - else - if Btrans(props) - BM = transpose(reshape(BT, (props.dright, props.dmid))) - else - BM = reshape(BT, (props.dmid, props.dright)) - end - end - - # TODO: this logic may be wrong - if props.permuteC - # if we are computing C = α * A B + β * C - # we need to make sure C is permuted to the same - # ordering as A B which is the inverse of props.PC - if β ≠ 0 - CM = reshape( - permutedims(leaf_parenttype(CT), CT, invperm(props.PC)), (props.dleft, props.dright) - ) - else - # Need to copy here since we will be permuting - # into C later - CM = reshape(copy(CT), (props.dleft, props.dright)) - end - else - if Ctrans(props) - CM = transpose(reshape(CT, (props.dright, props.dleft))) - else - CM = reshape(CT, (props.dleft, props.dright)) - end - end - - #tC = similar(CM) - #_gemm!(tA, tB, El(α), AM, BM, El(β), CM) - CM = mul!!( - leaf_parenttype(CM), CM, leaf_parenttype(AM), AM, leaf_parenttype(BM), BM, El(α), El(β) - ) - - if props.permuteC - Cr = reshape(CM, props.newCrange) - # TODO: use invperm(pC) here? - #@timeit_debug timer "_contract!: permutedims C" begin - CT .= permutedims(leaf_parenttype(Cr), Cr, props.PC) - #end # @timeit - end - - return CT -end - -function mul!!(::Type{<:Array}, CM, ::Type{<:Array}, AM, ::Type{<:Array}, BM, α, β) - return @strided mul!(CM, AM, BM, α, β) -end - -function NDTensors.mul!!( - ::Type{<:AbstractArray}, - CM, - ::Type{<:AbstractArray}, - AM, - ::Type{<:AbstractArray}, - BM, - α, - β, -) - return mul!(CM, AM, BM, α, β) -end From 0a8e8953edfff0b4e3f0d928373ad4b9edeec017 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 10:03:27 -0400 Subject: [PATCH 108/133] Create choke point for permutedims! function --- NDTensors/src/abstractarray/permutedims.jl | 4 ++++ NDTensors/src/dense/densetensor.jl | 4 ++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index cb7c0052ae..108aa037d8 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -2,6 +2,10 @@ function NDTensors.permutedims(::Type{<:AbstractArray}, M, perm) return Base.permutedims(M, perm) end +function permutedims!(Mdest::AbstractArray, M::AbstractArray, perm) + return permutedims!(leaf_parenttype(Mdest), Mdest, leaf_parenttype(M), M, perm) +end + function permutedims!(::Type{<:AbstractArray}, Mdest, ::Type{<:AbstractArray}, M, perm) return Mdest .= Base.permutedims(M, perm) end diff --git a/NDTensors/src/dense/densetensor.jl b/NDTensors/src/dense/densetensor.jl index d3d5582fcb..bd149920e8 100644 --- a/NDTensors/src/dense/densetensor.jl +++ b/NDTensors/src/dense/densetensor.jl @@ -199,7 +199,7 @@ function permutedims!( ) where {N,StoreT<:StridedArray} RA = array(R) TA = array(T) - RA .= permutedims!(leaf_parenttype(RA), RA, leaf_parenttype(TA), TA, perm) + RA = permutedims!(RA, TA, perm) return R end @@ -247,7 +247,7 @@ function permutedims!( end RA = array(R) TA = array(T) - return permutedims!!(leaf_parenttype(RA), RA, leaf_parenttype(TA), TA, perm, f) + return permutedims!!(RA, TA, perm, f) end """ From 09e83dd8010ef89c34e912fdb4ffc15fd8a44920 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 10:24:54 -0400 Subject: [PATCH 109/133] Remove this function as it created an infinite loop --- NDTensors/src/abstractarray/permutedims.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index 108aa037d8..cb7c0052ae 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -2,10 +2,6 @@ function NDTensors.permutedims(::Type{<:AbstractArray}, M, perm) return Base.permutedims(M, perm) end -function permutedims!(Mdest::AbstractArray, M::AbstractArray, perm) - return permutedims!(leaf_parenttype(Mdest), Mdest, leaf_parenttype(M), M, perm) -end - function permutedims!(::Type{<:AbstractArray}, Mdest, ::Type{<:AbstractArray}, M, perm) return Mdest .= Base.permutedims(M, perm) end From 06d6a046b6b1d80e5ade186c1cd6cbbf5fd5d46e Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 10:25:23 -0400 Subject: [PATCH 110/133] rewrite this permutedims function --- NDTensors/src/dense/densetensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/dense/densetensor.jl b/NDTensors/src/dense/densetensor.jl index bd149920e8..d14e4f1d3c 100644 --- a/NDTensors/src/dense/densetensor.jl +++ b/NDTensors/src/dense/densetensor.jl @@ -199,7 +199,7 @@ function permutedims!( ) where {N,StoreT<:StridedArray} RA = array(R) TA = array(T) - RA = permutedims!(RA, TA, perm) + RA = permutedims!(leaf_parenttype(RA), RA, leaf_parenttype(TA), TA, perm) return R end From 062efebf381c4b85d26fa519971af6e18edfe7d7 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 10:26:37 -0400 Subject: [PATCH 111/133] Move abstract array contract functions to different file. --- NDTensors/src/NDTensors.jl | 1 + .../abstractarray/tensoralgebra/contract.jl | 67 +++++++++++++++++-- NDTensors/src/dense/tensoralgebra/contract.jl | 62 +---------------- 3 files changed, 66 insertions(+), 64 deletions(-) diff --git a/NDTensors/src/NDTensors.jl b/NDTensors/src/NDTensors.jl index 9e10d33711..4d62501244 100644 --- a/NDTensors/src/NDTensors.jl +++ b/NDTensors/src/NDTensors.jl @@ -70,6 +70,7 @@ include("tensor/similar.jl") include("adapt.jl") include("tensoralgebra/generic_tensor_operations.jl") include("tensoralgebra/contraction_logic.jl") +include("abstractarray/tensoralgebra/contract.jl") ##################################### # DenseTensor and DiagTensor diff --git a/NDTensors/src/abstractarray/tensoralgebra/contract.jl b/NDTensors/src/abstractarray/tensoralgebra/contract.jl index 64fd9d1fa4..86412bf88d 100644 --- a/NDTensors/src/abstractarray/tensoralgebra/contract.jl +++ b/NDTensors/src/abstractarray/tensoralgebra/contract.jl @@ -1,3 +1,65 @@ +using LinearAlgebra: BlasFloat +export backend_auto, backend_blas, backend_generic + +@eval struct GemmBackend{T} + (f::Type{<:GemmBackend})() = $(Expr(:new, :f)) +end +GemmBackend(s) = GemmBackend{Symbol(s)}() +macro GemmBackend_str(s) + return :(GemmBackend{$(Expr(:quote, Symbol(s)))}) +end + +const gemm_backend = Ref(:Auto) +function backend_auto() + return gemm_backend[] = :Auto +end +function backend_blas() + return gemm_backend[] = :BLAS +end +function backend_generic() + return gemm_backend[] = :Generic +end + +@inline function auto_select_backend( + ::Type{<:StridedVecOrMat{<:BlasFloat}}, + ::Type{<:StridedVecOrMat{<:BlasFloat}}, + ::Type{<:StridedVecOrMat{<:BlasFloat}}, +) + return GemmBackend(:BLAS) +end + +@inline function auto_select_backend( + ::Type{<:AbstractVecOrMat}, ::Type{<:AbstractVecOrMat}, ::Type{<:AbstractVecOrMat} +) + return GemmBackend(:Generic) +end + +function _gemm!( + tA, tB, alpha, A::TA, B::TB, beta, C::TC +) where {TA<:AbstractVecOrMat,TB<:AbstractVecOrMat,TC<:AbstractVecOrMat} + if gemm_backend[] == :Auto + _gemm!(auto_select_backend(TA, TB, TC), tA, tB, alpha, A, B, beta, C) + else + _gemm!(GemmBackend(gemm_backend[]), tA, tB, alpha, A, B, beta, C) + end +end + +# BLAS matmul +function _gemm!( + ::GemmBackend{:BLAS}, + tA, + tB, + alpha, + A::AbstractVecOrMat, + B::AbstractVecOrMat, + beta, + C::AbstractVecOrMat, +) + #@timeit_debug timer "BLAS.gemm!" begin + return BLAS.gemm!(tA, tB, alpha, A, B, beta, C) + #end # @timeit +end + # generic matmul function _gemm!( ::GemmBackend{:Generic}, @@ -13,11 +75,6 @@ function _gemm!( return C end -function contraction_output(tensor1::DenseTensor, tensor2::DenseTensor, indsR) - tensortypeR = contraction_output_type(typeof(tensor1), typeof(tensor2), indsR) - return NDTensors.similar(tensortypeR, indsR) -end - # Non-trivial permutation function _contract_scalar_perm!( Rᵃ::AbstractArray{ElR}, Tᵃ::AbstractArray, perm, α, β=zero(ElR) diff --git a/NDTensors/src/dense/tensoralgebra/contract.jl b/NDTensors/src/dense/tensoralgebra/contract.jl index 6857828705..e1d99ce2eb 100644 --- a/NDTensors/src/dense/tensoralgebra/contract.jl +++ b/NDTensors/src/dense/tensoralgebra/contract.jl @@ -1,62 +1,6 @@ -export backend_auto, backend_blas, backend_generic - -@eval struct GemmBackend{T} - (f::Type{<:GemmBackend})() = $(Expr(:new, :f)) -end -GemmBackend(s) = GemmBackend{Symbol(s)}() -macro GemmBackend_str(s) - return :(GemmBackend{$(Expr(:quote, Symbol(s)))}) -end - -const gemm_backend = Ref(:Auto) -function backend_auto() - return gemm_backend[] = :Auto -end -function backend_blas() - return gemm_backend[] = :BLAS -end -function backend_generic() - return gemm_backend[] = :Generic -end - -@inline function auto_select_backend( - ::Type{<:StridedVecOrMat{<:BlasFloat}}, - ::Type{<:StridedVecOrMat{<:BlasFloat}}, - ::Type{<:StridedVecOrMat{<:BlasFloat}}, -) - return GemmBackend(:BLAS) -end - -@inline function auto_select_backend( - ::Type{<:AbstractVecOrMat}, ::Type{<:AbstractVecOrMat}, ::Type{<:AbstractVecOrMat} -) - return GemmBackend(:Generic) -end - -function _gemm!( - tA, tB, alpha, A::TA, B::TB, beta, C::TC -) where {TA<:AbstractVecOrMat,TB<:AbstractVecOrMat,TC<:AbstractVecOrMat} - if gemm_backend[] == :Auto - _gemm!(auto_select_backend(TA, TB, TC), tA, tB, alpha, A, B, beta, C) - else - _gemm!(GemmBackend(gemm_backend[]), tA, tB, alpha, A, B, beta, C) - end -end - -# BLAS matmul -function _gemm!( - ::GemmBackend{:BLAS}, - tA, - tB, - alpha, - A::AbstractVecOrMat, - B::AbstractVecOrMat, - beta, - C::AbstractVecOrMat, -) - #@timeit_debug timer "BLAS.gemm!" begin - return BLAS.gemm!(tA, tB, alpha, A, B, beta, C) - #end # @timeit +function contraction_output(tensor1::DenseTensor, tensor2::DenseTensor, indsR) + tensortypeR = contraction_output_type(typeof(tensor1), typeof(tensor2), indsR) + return NDTensors.similar(tensortypeR, indsR) end # Both are scalar-like tensors From d579f5e8b65d80cafb113e12e3fc8a1979fc5cf9 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 10:26:58 -0400 Subject: [PATCH 112/133] remove unecessary call --- NDTensors/src/dense/dense.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/NDTensors/src/dense/dense.jl b/NDTensors/src/dense/dense.jl index 12e9b05882..96afd695b5 100644 --- a/NDTensors/src/dense/dense.jl +++ b/NDTensors/src/dense/dense.jl @@ -1,7 +1,6 @@ # # Dense storage # -using LinearAlgebra: BlasFloat struct Dense{ElT,DataT<:AbstractArray} <: TensorStorage{ElT} data::DataT From 851614a1a88058586d4e5936d3b52b2c0ae85f27 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 13:58:25 -0400 Subject: [PATCH 113/133] format --- NDTensors/src/abstractarray/mul.jl | 6 ++++-- NDTensors/src/abstractarray/tensoralgebra/contract.jl | 5 ++--- NDTensors/src/array/mul.jl | 2 +- NDTensors/src/array/permutedims.jl | 2 +- 4 files changed, 8 insertions(+), 7 deletions(-) diff --git a/NDTensors/src/abstractarray/mul.jl b/NDTensors/src/abstractarray/mul.jl index 203e69ff30..5234afc03e 100644 --- a/NDTensors/src/abstractarray/mul.jl +++ b/NDTensors/src/abstractarray/mul.jl @@ -1,5 +1,7 @@ function mul!!(CM::AbstractArray, AM::AbstractArray, BM::AbstractArray, α, β) - return mul!!(leaf_parenttype(CM), CM, leaf_parenttype(AM), AM, leaf_parenttype(BM), BM, α, β) + return mul!!( + leaf_parenttype(CM), CM, leaf_parenttype(AM), AM, leaf_parenttype(BM), BM, α, β + ) end function mul!!( @@ -13,4 +15,4 @@ function mul!!( β, ) return mul!(CM, AM, BM, α, β) -end \ No newline at end of file +end diff --git a/NDTensors/src/abstractarray/tensoralgebra/contract.jl b/NDTensors/src/abstractarray/tensoralgebra/contract.jl index 86412bf88d..8ca7ba2139 100644 --- a/NDTensors/src/abstractarray/tensoralgebra/contract.jl +++ b/NDTensors/src/abstractarray/tensoralgebra/contract.jl @@ -164,8 +164,7 @@ function _contract!( #tC = similar(CM) #_gemm!(tA, tB, El(α), AM, BM, El(β), CM) - CM = mul!!(CM, AM, BM, El(α), El(β) - ) + CM = mul!!(CM, AM, BM, El(α), El(β)) if props.permuteC Cr = reshape(CM, props.newCrange) @@ -176,4 +175,4 @@ function _contract!( end return CT -end \ No newline at end of file +end diff --git a/NDTensors/src/array/mul.jl b/NDTensors/src/array/mul.jl index a8c90e5b94..e8746319cc 100644 --- a/NDTensors/src/array/mul.jl +++ b/NDTensors/src/array/mul.jl @@ -1,3 +1,3 @@ function mul!!(::Type{<:Array}, CM, ::Type{<:Array}, AM, ::Type{<:Array}, BM, α, β) return @strided mul!(CM, AM, BM, α, β) -end \ No newline at end of file +end diff --git a/NDTensors/src/array/permutedims.jl b/NDTensors/src/array/permutedims.jl index 0c381fdda5..a04854f315 100644 --- a/NDTensors/src/array/permutedims.jl +++ b/NDTensors/src/array/permutedims.jl @@ -9,4 +9,4 @@ end function permutedims!!(::Type{<:Array}, B, ::Type{<:Array}, A, perm, f) @strided B .= f.(B, Base.permutedims(A, perm)) -end \ No newline at end of file +end From cb962438c919542b10a269473465ca4e677198da Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 16:01:43 -0400 Subject: [PATCH 114/133] These functions could cause conflicts ``` WARNING: Method definition permutedims(AbstractArray{T, N} where N where T, Any) in module PermutedDimsArrays at permuteddimsarray.jl:137 overwritten in module NDTensors at /mnt/home/kpierce/.julia/dev/ITensors/NDTensors/src/abstractarray/permutedims.jl:1. ** incremental compilation may be fatally broken for this module ** ``` --- NDTensors/src/abstractarray/permutedims.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index cb7c0052ae..7c7e8a9bbd 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -1,7 +1,15 @@ +function NDTensors.permutedims(M::AbstractArray, perm) + permutedims(leaf_parenttype(M), M, perm) +end + function NDTensors.permutedims(::Type{<:AbstractArray}, M, perm) return Base.permutedims(M, perm) end +function NDTensors.permutedims!(Mdest::AbstractArray, M::AbstractArray, perm) + permutedims!(leaf_parenttype(Mdest), Mdest, leaf_parenttype(M), M, perm) +end + function permutedims!(::Type{<:AbstractArray}, Mdest, ::Type{<:AbstractArray}, M, perm) return Mdest .= Base.permutedims(M, perm) end From 6c902bd956082cf4c061c00e3766725e7120b323 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Thu, 12 Oct 2023 16:02:13 -0400 Subject: [PATCH 115/133] format --- NDTensors/src/abstractarray/permutedims.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index 7c7e8a9bbd..a1b0f44bb3 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -1,5 +1,5 @@ function NDTensors.permutedims(M::AbstractArray, perm) - permutedims(leaf_parenttype(M), M, perm) + return permutedims(leaf_parenttype(M), M, perm) end function NDTensors.permutedims(::Type{<:AbstractArray}, M, perm) @@ -7,7 +7,7 @@ function NDTensors.permutedims(::Type{<:AbstractArray}, M, perm) end function NDTensors.permutedims!(Mdest::AbstractArray, M::AbstractArray, perm) - permutedims!(leaf_parenttype(Mdest), Mdest, leaf_parenttype(M), M, perm) + return permutedims!(leaf_parenttype(Mdest), Mdest, leaf_parenttype(M), M, perm) end function permutedims!(::Type{<:AbstractArray}, Mdest, ::Type{<:AbstractArray}, M, perm) From 5c74d2a06054d69383f487a7e526461f8ecd8f79 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 13 Oct 2023 14:40:12 -0400 Subject: [PATCH 116/133] Updated and working NDTensors.permutedims --- NDTensors/src/abstractarray/permutedims.jl | 12 +++++++++--- NDTensors/src/array/permutedims.jl | 3 ++- NDTensors/src/imports.jl | 2 -- src/imports.jl | 3 ++- src/tensor_operations/permutations.jl | 4 ++-- 5 files changed, 15 insertions(+), 9 deletions(-) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index a1b0f44bb3..9494ff2804 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -1,15 +1,21 @@ -function NDTensors.permutedims(M::AbstractArray, perm) +## NOTICE!!: Here we are not importing Base.permutedims or Base.permutedims! but +## are writing our own implementation. This allows us to +# NDTensors.permutedims +function permutedims(M::AbstractArray, perm) return permutedims(leaf_parenttype(M), M, perm) end -function NDTensors.permutedims(::Type{<:AbstractArray}, M, perm) +# NDTensors.permutedims +function permutedims(::Type{<:AbstractArray}, M, perm) return Base.permutedims(M, perm) end -function NDTensors.permutedims!(Mdest::AbstractArray, M::AbstractArray, perm) +# NDTensors.permutedims! +function permutedims!(Mdest::AbstractArray, M::AbstractArray, perm) return permutedims!(leaf_parenttype(Mdest), Mdest, leaf_parenttype(M), M, perm) end +# NDTensors.permutedims! function permutedims!(::Type{<:AbstractArray}, Mdest, ::Type{<:AbstractArray}, M, perm) return Mdest .= Base.permutedims(M, perm) end diff --git a/NDTensors/src/array/permutedims.jl b/NDTensors/src/array/permutedims.jl index a04854f315..63bac52118 100644 --- a/NDTensors/src/array/permutedims.jl +++ b/NDTensors/src/array/permutedims.jl @@ -1,8 +1,9 @@ +# NDTensors.permutedims function permutedims(::Type{<:Array}, M, perm) return @strided Mdest = Base.permutedims(M, perm) end -## TODO is it possible to do something like f = identity and remove the permutedims!! funciton? +# NDTensors.permutedims! function permutedims!(::Type{<:Array}, Mdest, ::Type{<:Array}, M, perm) return @strided Mdest .= Base.permutedims(M, perm) end diff --git a/NDTensors/src/imports.jl b/NDTensors/src/imports.jl index 13b466fca4..4b58569f21 100644 --- a/NDTensors/src/imports.jl +++ b/NDTensors/src/imports.jl @@ -31,8 +31,6 @@ import Base: iterate, length, map, - permutedims, - permutedims!, print, promote_rule, randn, diff --git a/src/imports.jl b/src/imports.jl index cbdb3839d2..02e5830eb5 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -51,7 +51,6 @@ import Base: map, map!, ndims, - permutedims, print, promote_rule, push!, @@ -127,6 +126,8 @@ using ITensors.NDTensors: enable_auto_fermion, fill!!, randn!!, + permutedims, + permutedims!, set_eltype, single_precision, timer, diff --git a/src/tensor_operations/permutations.jl b/src/tensor_operations/permutations.jl index 5f7bb174e0..b85e4f5345 100644 --- a/src/tensor_operations/permutations.jl +++ b/src/tensor_operations/permutations.jl @@ -55,12 +55,12 @@ function permute(T::ITensor, new_inds...; kwargs...) end # TODO: move to NDTensors -function permutedims(::AllowAlias, T::Tensor, perm) +function NDTensors.permutedims(::AllowAlias, T::Tensor, perm) return NDTensors.is_trivial_permutation(perm) ? T : permutedims(NeverAlias(), T, perm) end # TODO: move to NDTensors, define `permutedims` in terms of `NeverAlias` -function permutedims(::NeverAlias, T::Tensor, perm) +function NDTensors.permutedims(::NeverAlias, T::Tensor, perm) return permutedims(T, perm) end From 211eb8446d5be792a7654fe28c17deb727f4fc30 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 13 Oct 2023 17:02:40 -0400 Subject: [PATCH 117/133] Permute in place --- NDTensors/src/abstractarray/permutedims.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index 9494ff2804..3e3571954f 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -17,7 +17,7 @@ end # NDTensors.permutedims! function permutedims!(::Type{<:AbstractArray}, Mdest, ::Type{<:AbstractArray}, M, perm) - return Mdest .= Base.permutedims(M, perm) + return Base.permutedims!(Mdest, M, perm) end function permutedims!!(B::AbstractArray, A::AbstractArray, perm, f) From 23e76599fddec7259f56c3707682aecc2c604dee Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 13 Oct 2023 17:03:40 -0400 Subject: [PATCH 118/133] Call NDTensors.permute in tests --- NDTensors/test/arraytensor/array.jl | 2 +- NDTensors/test/arraytensor/blocksparsearray.jl | 2 +- NDTensors/test/blocksparse.jl | 4 ++-- NDTensors/test/combiner.jl | 4 ++-- NDTensors/test/dense.jl | 4 ++-- NDTensors/test/diag.jl | 2 +- test/ITensorNetworkMaps/utils/utils.jl | 2 +- test/base/test_itensor_scalar_contract.jl | 8 ++++---- 8 files changed, 14 insertions(+), 14 deletions(-) diff --git a/NDTensors/test/arraytensor/array.jl b/NDTensors/test/arraytensor/array.jl index 39809d0db3..a8542e2f93 100644 --- a/NDTensors/test/arraytensor/array.jl +++ b/NDTensors/test/arraytensor/array.jl @@ -32,7 +32,7 @@ using NDTensors: storage, storagetype T1r = randn!(similar(T1)) @test Array(T1r + T1) ≈ Array(T1r) + Array(T1) - @test Array(permutedims(T1, (2, 1))) ≈ permutedims(Array(T1), (2, 1)) + @test Array(NDTensors.permutedims(T1, (2, 1))) ≈ NDTensors.permutedims(Array(T1), (2, 1)) U, S, V = svd(T1) @test U * S * V ≈ T1 diff --git a/NDTensors/test/arraytensor/blocksparsearray.jl b/NDTensors/test/arraytensor/blocksparsearray.jl index 5d508608a5..9c7173b943 100644 --- a/NDTensors/test/arraytensor/blocksparsearray.jl +++ b/NDTensors/test/arraytensor/blocksparsearray.jl @@ -38,7 +38,7 @@ using NDTensors: storage, storagetype @test_broken T1r = randn!(similar(T1)) @test_broken Array(T1r + T1) ≈ Array(T1r) + Array(T1) - @test_broken Array(permutedims(T1, (2, 1))) ≈ permutedims(Array(T1), (2, 1)) + @test_broken Array(NDTensors.permutedims(T1, (2, 1))) ≈ NDTensors.permutedims(Array(T1), (2, 1)) # TODO: Not implemented yet. ## U, S, V = svd(T1) diff --git a/NDTensors/test/blocksparse.jl b/NDTensors/test/blocksparse.jl index 6a541f15a1..2e552551dc 100644 --- a/NDTensors/test/blocksparse.jl +++ b/NDTensors/test/blocksparse.jl @@ -82,7 +82,7 @@ end @test C[I] == A[I] + B[I] end - Ap = permutedims(A, (2, 1)) + Ap = NDTensors.permutedims(A, (2, 1)) @test blockdims(Ap, (1, 2)) == (4, 3) @test blockdims(Ap, (2, 1)) == (5, 2) @@ -162,7 +162,7 @@ end @test nnzblocks(A) == nnzblocks(B) @test nnz(A) == nnz(B) - Ap = permutedims(A, (3, 2, 1)) + Ap = NDTensors.(A, (3, 2, 1)) for (bAp, bB) in zip(eachnzblock(Ap), eachnzblock(B)) blockAp = blockview(Ap, bAp) diff --git a/NDTensors/test/combiner.jl b/NDTensors/test/combiner.jl index f8853ff8ff..ace35975df 100644 --- a/NDTensors/test/combiner.jl +++ b/NDTensors/test/combiner.jl @@ -67,14 +67,14 @@ end output_tensor = contract(input_tensor, (1, -1, -2), combiner_tensor, (2, -1, -2)) @test output_tensor isa BlockSparseTensor @test dims(output_tensor) == dims(output_tensor_inds) - output_tensor = permutedims(output_tensor, (2, 1)) + output_tensor = NDTensors.permutedims(output_tensor, (2, 1)) for i in 1:length(input_tensor) @test input_tensor[i] == output_tensor[i] end # Test uncombining. Broken for inds that are not `Index`. new_input_tensor = contract(output_tensor, (1, -1), combiner_tensor, (-1, 2, 3)) - new_input_tensor = permutedims(new_input_tensor, (3, 1, 2)) + new_input_tensor = NDTensors.permutedims(new_input_tensor, (3, 1, 2)) @test new_input_tensor == input_tensor # Catch invalid combining diff --git a/NDTensors/test/dense.jl b/NDTensors/test/dense.jl index 80be64c0ed..0ba8ba8432 100644 --- a/NDTensors/test/dense.jl +++ b/NDTensors/test/dense.jl @@ -67,7 +67,7 @@ end @test C[I] == A[I] + B[I] end - Ap = permutedims(A, (2, 1)) + Ap = NDTensors.permutedims(A, (2, 1)) for I in eachindex(A) @test A[I] == Ap[NDTensors.permute(I, (2, 1))] @@ -222,7 +222,7 @@ end T2 = dev(randomTensor(complex(elt), (1, 1))) NDTensors.contract!(R, (2, 1, 3), T1, (1, 2, -1), T2, (-1, 1)) @test !any(isnan, R) - @test convert(Array, R) ≈ permutedims(convert(Array, T1), (2, 1, 3)) * T2[1] + @test convert(Array, R) ≈ NDTensors.permutedims(convert(Array, T1), (2, 1, 3)) * T2[1] end end end diff --git a/NDTensors/test/diag.jl b/NDTensors/test/diag.jl index 8422082e92..4e3464fc72 100644 --- a/NDTensors/test/diag.jl +++ b/NDTensors/test/diag.jl @@ -28,7 +28,7 @@ using Test D = dev(tensor(Diag(vr), (d, d))) @test Array(D) == NDTensors.LinearAlgebra.diagm(0 => vr) @test matrix(D) == NDTensors.LinearAlgebra.diagm(0 => vr) - @test permutedims(D, (2, 1)) == D + @test NDTensors.permutedims(D, (2, 1)) == D end end @testset "DiagTensor contractions" begin diff --git a/test/ITensorNetworkMaps/utils/utils.jl b/test/ITensorNetworkMaps/utils/utils.jl index 02275aee89..31a62f4e8c 100644 --- a/test/ITensorNetworkMaps/utils/utils.jl +++ b/test/ITensorNetworkMaps/utils/utils.jl @@ -33,7 +33,7 @@ function ITensors.prime(::typeof(linkinds), tn::InfTN) return tn_p end -interleave(x::Vector, y::Vector) = permutedims([x y])[:] +interleave(x::Vector, y::Vector) = NDTensors.permutedims([x y])[:] function ITensors.linkind(tn::InfTN, e) return commonind(tn[e[1]], tn[e[2]]) diff --git a/test/base/test_itensor_scalar_contract.jl b/test/base/test_itensor_scalar_contract.jl index c322139bba..63e311e59e 100644 --- a/test/base/test_itensor_scalar_contract.jl +++ b/test/base/test_itensor_scalar_contract.jl @@ -65,14 +65,14 @@ end R .= A .* B @test !any(isnan, R) - @test array(R) ≈ permutedims(array(A), (2, 1, 3)) * array(B)[] + @test array(R) ≈ NDTensors.permutedims(array(A), (2, 1, 3)) * array(B)[] R .= NaN @test any(isnan, R) R .= B .* A @test !any(isnan, R) - @test array(R) ≈ permutedims(array(A), (2, 1, 3)) * array(B)[] + @test array(R) ≈ NDTensors.permutedims(array(A), (2, 1, 3)) * array(B)[] end @testset "General contraction, permutation" for ElA in BlasFloats, ElB in BlasFloats @@ -89,7 +89,7 @@ end R .= A .* B @test !any(isnan, R) @test reshape(array(R), 6, 2) ≈ - reshape(permutedims(array(A), (2, 1, 3)), 6, 2) * array(B) + reshape(NDTensors.permutedims(array(A), (2, 1, 3)), 6, 2) * array(B) R .= NaN @test any(isnan, R) @@ -97,6 +97,6 @@ end R .= B .* A @test !any(isnan, R) @test reshape(array(R), 6, 2) ≈ - reshape(permutedims(array(A), (2, 1, 3)), 6, 2) * array(B) + reshape(NDTensors.permutedims(array(A), (2, 1, 3)), 6, 2) * array(B) end end From bf82369bdea4c0f1136ba1148440b9d058e84a3e Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Fri, 13 Oct 2023 17:04:09 -0400 Subject: [PATCH 119/133] Format --- NDTensors/test/arraytensor/blocksparsearray.jl | 3 ++- NDTensors/test/dense.jl | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/NDTensors/test/arraytensor/blocksparsearray.jl b/NDTensors/test/arraytensor/blocksparsearray.jl index 9c7173b943..8abe6de266 100644 --- a/NDTensors/test/arraytensor/blocksparsearray.jl +++ b/NDTensors/test/arraytensor/blocksparsearray.jl @@ -38,7 +38,8 @@ using NDTensors: storage, storagetype @test_broken T1r = randn!(similar(T1)) @test_broken Array(T1r + T1) ≈ Array(T1r) + Array(T1) - @test_broken Array(NDTensors.permutedims(T1, (2, 1))) ≈ NDTensors.permutedims(Array(T1), (2, 1)) + @test_broken Array(NDTensors.permutedims(T1, (2, 1))) ≈ + NDTensors.permutedims(Array(T1), (2, 1)) # TODO: Not implemented yet. ## U, S, V = svd(T1) diff --git a/NDTensors/test/dense.jl b/NDTensors/test/dense.jl index 0ba8ba8432..46d85e9a07 100644 --- a/NDTensors/test/dense.jl +++ b/NDTensors/test/dense.jl @@ -222,7 +222,8 @@ end T2 = dev(randomTensor(complex(elt), (1, 1))) NDTensors.contract!(R, (2, 1, 3), T1, (1, 2, -1), T2, (-1, 1)) @test !any(isnan, R) - @test convert(Array, R) ≈ NDTensors.permutedims(convert(Array, T1), (2, 1, 3)) * T2[1] + @test convert(Array, R) ≈ + NDTensors.permutedims(convert(Array, T1), (2, 1, 3)) * T2[1] end end end From 6445af676e7ca68376274ffb3fd66c93830af9b6 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Sat, 14 Oct 2023 18:53:31 -0400 Subject: [PATCH 120/133] Make `base.permutedims(Tensor) =NDTensors.permutedims` --- NDTensors/src/NDTensors.jl | 1 + NDTensors/src/tensor/permutedims.jl | 1 + NDTensors/test/arraytensor/array.jl | 2 +- NDTensors/test/arraytensor/blocksparsearray.jl | 4 ++-- NDTensors/test/blocksparse.jl | 4 ++-- NDTensors/test/combiner.jl | 4 ++-- NDTensors/test/dense.jl | 4 ++-- NDTensors/test/diag.jl | 2 +- 8 files changed, 12 insertions(+), 10 deletions(-) create mode 100644 NDTensors/src/tensor/permutedims.jl diff --git a/NDTensors/src/NDTensors.jl b/NDTensors/src/NDTensors.jl index 4d62501244..641253dec0 100644 --- a/NDTensors/src/NDTensors.jl +++ b/NDTensors/src/NDTensors.jl @@ -67,6 +67,7 @@ include("tensor/tensor.jl") include("dims.jl") include("tensor/set_types.jl") include("tensor/similar.jl") +include("tensor/permutedims.jl") include("adapt.jl") include("tensoralgebra/generic_tensor_operations.jl") include("tensoralgebra/contraction_logic.jl") diff --git a/NDTensors/src/tensor/permutedims.jl b/NDTensors/src/tensor/permutedims.jl new file mode 100644 index 0000000000..516da64e35 --- /dev/null +++ b/NDTensors/src/tensor/permutedims.jl @@ -0,0 +1 @@ +Base.permutedims(A::Tensor, perm) = NDTensors.permutedims(A, perm) \ No newline at end of file diff --git a/NDTensors/test/arraytensor/array.jl b/NDTensors/test/arraytensor/array.jl index a8542e2f93..39809d0db3 100644 --- a/NDTensors/test/arraytensor/array.jl +++ b/NDTensors/test/arraytensor/array.jl @@ -32,7 +32,7 @@ using NDTensors: storage, storagetype T1r = randn!(similar(T1)) @test Array(T1r + T1) ≈ Array(T1r) + Array(T1) - @test Array(NDTensors.permutedims(T1, (2, 1))) ≈ NDTensors.permutedims(Array(T1), (2, 1)) + @test Array(permutedims(T1, (2, 1))) ≈ permutedims(Array(T1), (2, 1)) U, S, V = svd(T1) @test U * S * V ≈ T1 diff --git a/NDTensors/test/arraytensor/blocksparsearray.jl b/NDTensors/test/arraytensor/blocksparsearray.jl index 8abe6de266..526122b2e1 100644 --- a/NDTensors/test/arraytensor/blocksparsearray.jl +++ b/NDTensors/test/arraytensor/blocksparsearray.jl @@ -38,8 +38,8 @@ using NDTensors: storage, storagetype @test_broken T1r = randn!(similar(T1)) @test_broken Array(T1r + T1) ≈ Array(T1r) + Array(T1) - @test_broken Array(NDTensors.permutedims(T1, (2, 1))) ≈ - NDTensors.permutedims(Array(T1), (2, 1)) + @test_broken Array(permutedims(T1, (2, 1))) ≈ + permutedims(Array(T1), (2, 1)) # TODO: Not implemented yet. ## U, S, V = svd(T1) diff --git a/NDTensors/test/blocksparse.jl b/NDTensors/test/blocksparse.jl index 2e552551dc..8f52f51f14 100644 --- a/NDTensors/test/blocksparse.jl +++ b/NDTensors/test/blocksparse.jl @@ -82,7 +82,7 @@ end @test C[I] == A[I] + B[I] end - Ap = NDTensors.permutedims(A, (2, 1)) + Ap = permutedims(A, (2, 1)) @test blockdims(Ap, (1, 2)) == (4, 3) @test blockdims(Ap, (2, 1)) == (5, 2) @@ -162,7 +162,7 @@ end @test nnzblocks(A) == nnzblocks(B) @test nnz(A) == nnz(B) - Ap = NDTensors.(A, (3, 2, 1)) + Ap = NDTensors.permutedims(A, (3, 2, 1)) for (bAp, bB) in zip(eachnzblock(Ap), eachnzblock(B)) blockAp = blockview(Ap, bAp) diff --git a/NDTensors/test/combiner.jl b/NDTensors/test/combiner.jl index ace35975df..f8853ff8ff 100644 --- a/NDTensors/test/combiner.jl +++ b/NDTensors/test/combiner.jl @@ -67,14 +67,14 @@ end output_tensor = contract(input_tensor, (1, -1, -2), combiner_tensor, (2, -1, -2)) @test output_tensor isa BlockSparseTensor @test dims(output_tensor) == dims(output_tensor_inds) - output_tensor = NDTensors.permutedims(output_tensor, (2, 1)) + output_tensor = permutedims(output_tensor, (2, 1)) for i in 1:length(input_tensor) @test input_tensor[i] == output_tensor[i] end # Test uncombining. Broken for inds that are not `Index`. new_input_tensor = contract(output_tensor, (1, -1), combiner_tensor, (-1, 2, 3)) - new_input_tensor = NDTensors.permutedims(new_input_tensor, (3, 1, 2)) + new_input_tensor = permutedims(new_input_tensor, (3, 1, 2)) @test new_input_tensor == input_tensor # Catch invalid combining diff --git a/NDTensors/test/dense.jl b/NDTensors/test/dense.jl index 46d85e9a07..a38bb00eee 100644 --- a/NDTensors/test/dense.jl +++ b/NDTensors/test/dense.jl @@ -67,7 +67,7 @@ end @test C[I] == A[I] + B[I] end - Ap = NDTensors.permutedims(A, (2, 1)) + Ap = permutedims(A, (2, 1)) for I in eachindex(A) @test A[I] == Ap[NDTensors.permute(I, (2, 1))] @@ -223,7 +223,7 @@ end NDTensors.contract!(R, (2, 1, 3), T1, (1, 2, -1), T2, (-1, 1)) @test !any(isnan, R) @test convert(Array, R) ≈ - NDTensors.permutedims(convert(Array, T1), (2, 1, 3)) * T2[1] + permutedims(convert(Array, T1), (2, 1, 3)) * T2[1] end end end diff --git a/NDTensors/test/diag.jl b/NDTensors/test/diag.jl index 4e3464fc72..8422082e92 100644 --- a/NDTensors/test/diag.jl +++ b/NDTensors/test/diag.jl @@ -28,7 +28,7 @@ using Test D = dev(tensor(Diag(vr), (d, d))) @test Array(D) == NDTensors.LinearAlgebra.diagm(0 => vr) @test matrix(D) == NDTensors.LinearAlgebra.diagm(0 => vr) - @test NDTensors.permutedims(D, (2, 1)) == D + @test permutedims(D, (2, 1)) == D end end @testset "DiagTensor contractions" begin From 2c18c13ab564f342290c627241f7cd971f0c840e Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Sat, 14 Oct 2023 18:58:58 -0400 Subject: [PATCH 121/133] Make permutedims match bangbang code flow --- NDTensors/src/abstractarray/permutedims.jl | 10 ++++++++-- NDTensors/src/array/permutedims.jl | 3 ++- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index 3e3571954f..d9afb74e25 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -24,6 +24,12 @@ function permutedims!!(B::AbstractArray, A::AbstractArray, perm, f) return permutedims!!(leaf_parenttype(B), B, leaf_parenttype(A), A, perm, f) end -function permutedims!!(::Type{<:AbstractArray}, B, ::Type{<:AbstractArray}, A, perm, f) - return B .= f.(B, Base.permutedims(A, perm)) +function permutedims!!(Bleaftype::Type{<:AbstractArray}, B, Aleaftype::Type{<:AbstractArray}, A, perm, f) + permutedims!(Bleaftype, B, Aleaftype, A, perm, f) + return B +end + +function permutedims!(::Type{<:AbstractArray}, B, ::Type{<:AbstractArray}, A, perm, f) + B .= f.(B, Base.permutedims(A, perm)) + return B end diff --git a/NDTensors/src/array/permutedims.jl b/NDTensors/src/array/permutedims.jl index 63bac52118..73d7e5afe4 100644 --- a/NDTensors/src/array/permutedims.jl +++ b/NDTensors/src/array/permutedims.jl @@ -8,6 +8,7 @@ function permutedims!(::Type{<:Array}, Mdest, ::Type{<:Array}, M, perm) return @strided Mdest .= Base.permutedims(M, perm) end -function permutedims!!(::Type{<:Array}, B, ::Type{<:Array}, A, perm, f) +function permutedims!(::Type{<:Array}, B, ::Type{<:Array}, A, perm, f) @strided B .= f.(B, Base.permutedims(A, perm)) + return B end From f1a1b6117909a253497ab892926fc2feb929632a Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Sat, 14 Oct 2023 18:59:15 -0400 Subject: [PATCH 122/133] Have mul call mul! then return --- NDTensors/src/abstractarray/mul.jl | 4 +++- NDTensors/src/array/mul.jl | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/NDTensors/src/abstractarray/mul.jl b/NDTensors/src/abstractarray/mul.jl index 5234afc03e..e5f6bdffc2 100644 --- a/NDTensors/src/abstractarray/mul.jl +++ b/NDTensors/src/abstractarray/mul.jl @@ -2,6 +2,7 @@ function mul!!(CM::AbstractArray, AM::AbstractArray, BM::AbstractArray, α, β) return mul!!( leaf_parenttype(CM), CM, leaf_parenttype(AM), AM, leaf_parenttype(BM), BM, α, β ) + return CM; end function mul!!( @@ -14,5 +15,6 @@ function mul!!( α, β, ) - return mul!(CM, AM, BM, α, β) + mul!(CM, AM, BM, α, β) + return CM end diff --git a/NDTensors/src/array/mul.jl b/NDTensors/src/array/mul.jl index e8746319cc..f5ca2ec67e 100644 --- a/NDTensors/src/array/mul.jl +++ b/NDTensors/src/array/mul.jl @@ -1,3 +1,4 @@ function mul!!(::Type{<:Array}, CM, ::Type{<:Array}, AM, ::Type{<:Array}, BM, α, β) - return @strided mul!(CM, AM, BM, α, β) + @strided CM = mul!(CM, AM, BM, α, β) + return CM end From b5fd9d1d7ff258f325ea36fbcff02778bcd12c73 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Sat, 14 Oct 2023 18:59:36 -0400 Subject: [PATCH 123/133] format --- NDTensors/src/abstractarray/mul.jl | 2 +- NDTensors/src/abstractarray/permutedims.jl | 4 +++- NDTensors/src/array/mul.jl | 2 +- NDTensors/src/tensor/permutedims.jl | 2 +- NDTensors/test/arraytensor/blocksparsearray.jl | 3 +-- NDTensors/test/dense.jl | 3 +-- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/NDTensors/src/abstractarray/mul.jl b/NDTensors/src/abstractarray/mul.jl index e5f6bdffc2..954ecca3e1 100644 --- a/NDTensors/src/abstractarray/mul.jl +++ b/NDTensors/src/abstractarray/mul.jl @@ -2,7 +2,7 @@ function mul!!(CM::AbstractArray, AM::AbstractArray, BM::AbstractArray, α, β) return mul!!( leaf_parenttype(CM), CM, leaf_parenttype(AM), AM, leaf_parenttype(BM), BM, α, β ) - return CM; + return CM end function mul!!( diff --git a/NDTensors/src/abstractarray/permutedims.jl b/NDTensors/src/abstractarray/permutedims.jl index d9afb74e25..f23db4ea73 100644 --- a/NDTensors/src/abstractarray/permutedims.jl +++ b/NDTensors/src/abstractarray/permutedims.jl @@ -24,7 +24,9 @@ function permutedims!!(B::AbstractArray, A::AbstractArray, perm, f) return permutedims!!(leaf_parenttype(B), B, leaf_parenttype(A), A, perm, f) end -function permutedims!!(Bleaftype::Type{<:AbstractArray}, B, Aleaftype::Type{<:AbstractArray}, A, perm, f) +function permutedims!!( + Bleaftype::Type{<:AbstractArray}, B, Aleaftype::Type{<:AbstractArray}, A, perm, f +) permutedims!(Bleaftype, B, Aleaftype, A, perm, f) return B end diff --git a/NDTensors/src/array/mul.jl b/NDTensors/src/array/mul.jl index f5ca2ec67e..30c8338513 100644 --- a/NDTensors/src/array/mul.jl +++ b/NDTensors/src/array/mul.jl @@ -1,4 +1,4 @@ function mul!!(::Type{<:Array}, CM, ::Type{<:Array}, AM, ::Type{<:Array}, BM, α, β) - @strided CM = mul!(CM, AM, BM, α, β) + @strided CM = mul!(CM, AM, BM, α, β) return CM end diff --git a/NDTensors/src/tensor/permutedims.jl b/NDTensors/src/tensor/permutedims.jl index 516da64e35..f79e229504 100644 --- a/NDTensors/src/tensor/permutedims.jl +++ b/NDTensors/src/tensor/permutedims.jl @@ -1 +1 @@ -Base.permutedims(A::Tensor, perm) = NDTensors.permutedims(A, perm) \ No newline at end of file +Base.permutedims(A::Tensor, perm) = NDTensors.permutedims(A, perm) diff --git a/NDTensors/test/arraytensor/blocksparsearray.jl b/NDTensors/test/arraytensor/blocksparsearray.jl index 526122b2e1..5d508608a5 100644 --- a/NDTensors/test/arraytensor/blocksparsearray.jl +++ b/NDTensors/test/arraytensor/blocksparsearray.jl @@ -38,8 +38,7 @@ using NDTensors: storage, storagetype @test_broken T1r = randn!(similar(T1)) @test_broken Array(T1r + T1) ≈ Array(T1r) + Array(T1) - @test_broken Array(permutedims(T1, (2, 1))) ≈ - permutedims(Array(T1), (2, 1)) + @test_broken Array(permutedims(T1, (2, 1))) ≈ permutedims(Array(T1), (2, 1)) # TODO: Not implemented yet. ## U, S, V = svd(T1) diff --git a/NDTensors/test/dense.jl b/NDTensors/test/dense.jl index a38bb00eee..80be64c0ed 100644 --- a/NDTensors/test/dense.jl +++ b/NDTensors/test/dense.jl @@ -222,8 +222,7 @@ end T2 = dev(randomTensor(complex(elt), (1, 1))) NDTensors.contract!(R, (2, 1, 3), T1, (1, 2, -1), T2, (-1, 1)) @test !any(isnan, R) - @test convert(Array, R) ≈ - permutedims(convert(Array, T1), (2, 1, 3)) * T2[1] + @test convert(Array, R) ≈ permutedims(convert(Array, T1), (2, 1, 3)) * T2[1] end end end From 50eb9a87b157e30e5a041534770ac847fab56d06 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Sat, 14 Oct 2023 19:22:35 -0400 Subject: [PATCH 124/133] Use base.permutedims not NDTensors --- test/ITensorNetworkMaps/utils/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/ITensorNetworkMaps/utils/utils.jl b/test/ITensorNetworkMaps/utils/utils.jl index 31a62f4e8c..7b4254c4a5 100644 --- a/test/ITensorNetworkMaps/utils/utils.jl +++ b/test/ITensorNetworkMaps/utils/utils.jl @@ -33,7 +33,7 @@ function ITensors.prime(::typeof(linkinds), tn::InfTN) return tn_p end -interleave(x::Vector, y::Vector) = NDTensors.permutedims([x y])[:] +interleave(x::Vector, y::Vector) = Base.permutedims([x y])[:] function ITensors.linkind(tn::InfTN, e) return commonind(tn[e[1]], tn[e[2]]) From bf86c839a93961d37280030820172c8d52b26cac Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 16 Oct 2023 16:14:34 -0400 Subject: [PATCH 125/133] Use simplified functions to dispatch later on parenttype --- NDTensors/src/abstractarray/tensoralgebra/contract.jl | 8 ++++---- NDTensors/src/dense/densetensor.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/NDTensors/src/abstractarray/tensoralgebra/contract.jl b/NDTensors/src/abstractarray/tensoralgebra/contract.jl index 8ca7ba2139..52be037edb 100644 --- a/NDTensors/src/abstractarray/tensoralgebra/contract.jl +++ b/NDTensors/src/abstractarray/tensoralgebra/contract.jl @@ -114,7 +114,7 @@ function _contract!( tA = 'N' if props.permuteA #@timeit_debug timer "_contract!: permutedims A" begin - Ap = permutedims(leaf_parenttype(AT), AT, props.PA) + Ap = permutedims(AT, props.PA) #end # @timeit AM = transpose(reshape(Ap, (props.dmid, props.dleft))) else @@ -129,7 +129,7 @@ function _contract!( tB = 'N' if props.permuteB #@timeit_debug timer "_contract!: permutedims B" begin - Bp = permutedims(leaf_parenttype(BT), BT, props.PB) + Bp = permutedims(BT, props.PB) #end # @timeit BM = reshape(Bp, (props.dmid, props.dright)) else @@ -147,7 +147,7 @@ function _contract!( # ordering as A B which is the inverse of props.PC if β ≠ 0 CM = reshape( - permutedims(leaf_parenttype(CT), CT, invperm(props.PC)), (props.dleft, props.dright) + permutedims(CT, invperm(props.PC)), (props.dleft, props.dright) ) else # Need to copy here since we will be permuting @@ -170,7 +170,7 @@ function _contract!( Cr = reshape(CM, props.newCrange) # TODO: use invperm(pC) here? #@timeit_debug timer "_contract!: permutedims C" begin - CT .= permutedims(leaf_parenttype(Cr), Cr, props.PC) + CT .= permutedims(Cr, props.PC) #end # @timeit end diff --git a/NDTensors/src/dense/densetensor.jl b/NDTensors/src/dense/densetensor.jl index d14e4f1d3c..bd149920e8 100644 --- a/NDTensors/src/dense/densetensor.jl +++ b/NDTensors/src/dense/densetensor.jl @@ -199,7 +199,7 @@ function permutedims!( ) where {N,StoreT<:StridedArray} RA = array(R) TA = array(T) - RA = permutedims!(leaf_parenttype(RA), RA, leaf_parenttype(TA), TA, perm) + RA = permutedims!(RA, TA, perm) return R end From a836b89d4cf0cd42a21995cb77457739a6172655 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 16 Oct 2023 16:15:20 -0400 Subject: [PATCH 126/133] format --- NDTensors/src/abstractarray/tensoralgebra/contract.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/NDTensors/src/abstractarray/tensoralgebra/contract.jl b/NDTensors/src/abstractarray/tensoralgebra/contract.jl index 52be037edb..5db9bf5370 100644 --- a/NDTensors/src/abstractarray/tensoralgebra/contract.jl +++ b/NDTensors/src/abstractarray/tensoralgebra/contract.jl @@ -146,9 +146,7 @@ function _contract!( # we need to make sure C is permuted to the same # ordering as A B which is the inverse of props.PC if β ≠ 0 - CM = reshape( - permutedims(CT, invperm(props.PC)), (props.dleft, props.dright) - ) + CM = reshape(permutedims(CT, invperm(props.PC)), (props.dleft, props.dright)) else # Need to copy here since we will be permuting # into C later From 709a83b1f27f6a83e3e4e20c8dbae9dddf1f4889 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Mon, 16 Oct 2023 16:38:49 -0400 Subject: [PATCH 127/133] Update NDTensors/src/dense/densetensor.jl Co-authored-by: Matt Fishman --- NDTensors/src/dense/densetensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/dense/densetensor.jl b/NDTensors/src/dense/densetensor.jl index bd149920e8..14b8e571d2 100644 --- a/NDTensors/src/dense/densetensor.jl +++ b/NDTensors/src/dense/densetensor.jl @@ -199,7 +199,7 @@ function permutedims!( ) where {N,StoreT<:StridedArray} RA = array(R) TA = array(T) - RA = permutedims!(RA, TA, perm) +permutedims!(RA, TA, perm) return R end From eb26ef1436876e0d7d0296ea7f6533deed628a17 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Mon, 16 Oct 2023 16:39:01 -0400 Subject: [PATCH 128/133] Update NDTensors/src/array/mul.jl Co-authored-by: Matt Fishman --- NDTensors/src/array/mul.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/array/mul.jl b/NDTensors/src/array/mul.jl index 30c8338513..54104ce502 100644 --- a/NDTensors/src/array/mul.jl +++ b/NDTensors/src/array/mul.jl @@ -1,4 +1,4 @@ function mul!!(::Type{<:Array}, CM, ::Type{<:Array}, AM, ::Type{<:Array}, BM, α, β) - @strided CM = mul!(CM, AM, BM, α, β) + @strided mul!(CM, AM, BM, α, β) return CM end From e032199ac90d102c02a0a933e7ccf70198c4d269 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Mon, 16 Oct 2023 16:39:11 -0400 Subject: [PATCH 129/133] Update NDTensors/src/array/permutedims.jl Co-authored-by: Matt Fishman --- NDTensors/src/array/permutedims.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/array/permutedims.jl b/NDTensors/src/array/permutedims.jl index 73d7e5afe4..60a0ecfbf9 100644 --- a/NDTensors/src/array/permutedims.jl +++ b/NDTensors/src/array/permutedims.jl @@ -1,6 +1,6 @@ # NDTensors.permutedims function permutedims(::Type{<:Array}, M, perm) - return @strided Mdest = Base.permutedims(M, perm) + return @strided Base.permutedims(M, perm) end # NDTensors.permutedims! From 7d93e0cfe4cf23b9c82fd21486bf739729097dc0 Mon Sep 17 00:00:00 2001 From: Karl Pierce Date: Mon, 16 Oct 2023 16:39:27 -0400 Subject: [PATCH 130/133] Update NDTensors/src/array/permutedims.jl Co-authored-by: Matt Fishman --- NDTensors/src/array/permutedims.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NDTensors/src/array/permutedims.jl b/NDTensors/src/array/permutedims.jl index 60a0ecfbf9..7a033201fb 100644 --- a/NDTensors/src/array/permutedims.jl +++ b/NDTensors/src/array/permutedims.jl @@ -5,7 +5,8 @@ end # NDTensors.permutedims! function permutedims!(::Type{<:Array}, Mdest, ::Type{<:Array}, M, perm) - return @strided Mdest .= Base.permutedims(M, perm) + @strided Mdest .= Base.permutedims(M, perm) + return Mdest end function permutedims!(::Type{<:Array}, B, ::Type{<:Array}, A, perm, f) From cc9276729c1b3e31b5ac69702cba092f09fe79b6 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 16 Oct 2023 16:42:19 -0400 Subject: [PATCH 131/133] Add --- NDTensors/src/linearalgebra/linearalgebra.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index 578d8d27b2..b998348d98 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -218,6 +218,8 @@ function LinearAlgebra.eigen( use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff) matrixT = matrix(T) + ## Here calling parent to ensure that the correct `any` function + ## is envoked for non-cpu matrices if any(!isfinite, parent(matrixT)) throw( ArgumentError( From ebcea17f1f0bded67fad62c2344dea8130d801c2 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 16 Oct 2023 16:42:35 -0400 Subject: [PATCH 132/133] appended to previous commit --- NDTensors/src/linearalgebra/linearalgebra.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl index b998348d98..e2c8133eb9 100644 --- a/NDTensors/src/linearalgebra/linearalgebra.jl +++ b/NDTensors/src/linearalgebra/linearalgebra.jl @@ -218,7 +218,7 @@ function LinearAlgebra.eigen( use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff) matrixT = matrix(T) - ## Here calling parent to ensure that the correct `any` function + ## TODO Here I am calling parent to ensure that the correct `any` function ## is envoked for non-cpu matrices if any(!isfinite, parent(matrixT)) throw( From 262e92432153bba7819c931d8aa8c01f6b051d18 Mon Sep 17 00:00:00 2001 From: kmp5VT Date: Mon, 16 Oct 2023 16:42:55 -0400 Subject: [PATCH 133/133] format --- NDTensors/src/dense/densetensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NDTensors/src/dense/densetensor.jl b/NDTensors/src/dense/densetensor.jl index 14b8e571d2..ed773d6f0e 100644 --- a/NDTensors/src/dense/densetensor.jl +++ b/NDTensors/src/dense/densetensor.jl @@ -199,7 +199,7 @@ function permutedims!( ) where {N,StoreT<:StridedArray} RA = array(R) TA = array(T) -permutedims!(RA, TA, perm) + permutedims!(RA, TA, perm) return R end