[NDTensors] [BUG] Outer fails for two GPU diag tensors

[kmp5VT]

Description of bug

The outer! function fails when contracting two diag tensors.

Minimal runnable code

 using ITensors, CUDA
 i,j = Index.((10,11))
 A, B = cu(diag_itensor(1, i)), cu(diag_itensor(1, j))
 A * B

Actual output or behavior
Scalar indexing error

Output of minimal runnable code

Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore should be avoided.

If you want to allow scalar iteration, use `allowscalar` or `@allowscalar`
to enable scalar iteration globally or for the operations in question.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:44
  [2] errorscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex
    @ ~/.julia/packages/GPUArrays/w335n/src/host/indexing.jl:50 [inlined]
  [6] getindex
    @ ~/.julia/packages/NDTensors/iVL17/src/tensorstorage/tensorstorage.jl:28 [inlined]
  [7] getdiagindex
    @ ~/.julia/packages/NDTensors/iVL17/src/diag/diagtensor.jl:57 [inlined]
  [8] outer!(R::DenseTensor{Float32, 2, Tuple{…}, Dense{…}}, T1::DiagTensor{Float32, 1, Tuple{…}, Diag{…}}, T2::DiagTensor{Float32, 1, Tuple{…}, Diag{…}})
    @ NDTensors ~/.julia/packages/NDTensors/iVL17/src/diag/tensoralgebra/outer.jl:6
  [9] outer!!
    @ ~/.julia/packages/NDTensors/iVL17/src/tensoroperations/generic_tensor_operations.jl:228 [inlined]
 [10] contract!!(output_tensor::DenseTensor{…}, labelsoutput_tensor::Tuple{…}, tensor1::DiagTensor{…}, labelstensor1::Tuple{…}, tensor2::DiagTensor{…}, labelstensor2::Tuple{…}, α::Int64, β::Int64)
    @ NDTensors ~/.julia/packages/NDTensors/iVL17/src/tensoroperations/generic_tensor_operations.jl:199
 [11] contract!!
    @ ~/.julia/packages/NDTensors/iVL17/src/tensoroperations/generic_tensor_operations.jl:188 [inlined]
 [12] contract(tensor1::DiagTensor{…}, labelstensor1::Tuple{…}, tensor2::DiagTensor{…}, labelstensor2::Tuple{…}, labelsoutput_tensor::Tuple{…})
    @ NDTensors ~/.julia/packages/NDTensors/iVL17/src/tensoroperations/generic_tensor_operations.jl:113
 [13] contract(::Type{…}, tensor1::DiagTensor{…}, labels_tensor1::Tuple{…}, tensor2::DiagTensor{…}, labels_tensor2::Tuple{…})
    @ NDTensors ~/.julia/packages/NDTensors/iVL17/src/tensoroperations/generic_tensor_operations.jl:91
 [14] contract
    @ ~/.julia/packages/SimpleTraits/7VJph/src/SimpleTraits.jl:332 [inlined]
 [15] _contract(A::DiagTensor{Float32, 1, Tuple{Index{…}}, Diag{Float32, CuArray{…}}}, B::DiagTensor{Float32, 1, Tuple{Index{…}}, Diag{Float32, CuArray{…}}})
    @ ITensors ~/.julia/packages/ITensors/iPhbw/src/tensor_operations/tensor_algebra.jl:3
 [16] _contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/iPhbw/src/tensor_operations/tensor_algebra.jl:9
 [17] contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/iPhbw/src/tensor_operations/tensor_algebra.jl:74

Version information

• Output from versioninfo():

julia> versioninfo()
Julia Version 1.12.1
Commit ba1e628ee49 (2025-10-17 13:02 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 64 × Intel(R) Xeon(R) Platinum 8358 CPU @ 2.60GHz
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, icelake-server)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 64 virtual cores)
Environment:
  LD_LIBRARY_PATH = /mnt/sw/nix/store/1fjr9h1k55mcxz6kvyszh7qb4cvnl4wp-openblas-0.3.29/lib

• Output from using Pkg; Pkg.status("ITensors"):

  [052768ef] CUDA v5.9.3
  [9136182c] ITensors v0.9.14

[mtfishman]

Thanks for the report, that's probably a bit tricky to get working on GPU in general, a workaround could be to call dense on the inputs. Do you have a use case for this kind of operation?

[kmp5VT]

I am using this specific operation in the ITensorCPD.jl code but the performance of this operation is not critical so it should be okay to switch to Dense for now. Do you think it would be reasonable to add an expose definition for GPU array types to convert to dense? I also tried this operation on Metal and it fails for the same reason..

[mtfishman]

Good question, probably that is worth doing so that it "just works".

[kmp5VT]

I can submit a PR later today

[kmp5VT]

@mtfishman It looks like a more generic fix (which works for gpu arrays) is to utilize the outer product of two vectors, i.e.

function outer!(
        R::DenseTensor{<:Number, NR}, T1::DiagTensor{<:Number, N1}, T2::DiagTensor{<:Number, N2}
    ) where {NR, N1, N2}
    array(R) .= array(T1) .* array(T2')
    return R
end

Is this an okay fox for you? Also does the ' of a tensor compute the complex conjugate or should I put dag(T2)' ?

Also, in trying to write a unit test I found another potential bug.
In this function in NDTensors/diag/tensoralgebra/outer.jl

function outer(T1::DiagTensor{ElT1, N1}, T2::DiagTensor{ElT2, N2}) where {ElT1, ElT2, N1, N2}
    indsR = unioninds(inds(T1), inds(T2))
    R = tensor(Dense(generic_zeros(promote_type(ElT1, ElT2), dim(indsR))), indsR)
    outer!(R, T1, T2)
    return R
end

the function unioninds does not exist in the NDTensors namespace and is only a function in ITensors.

[kmp5VT]

See #1684