More missing ITensor broadcasting operations

[mtfishman]

The following broadcasting operation error or give incorrect results:

• A .= real.(A) .^ 2 .+ im .* imag.(A) .^2.
• real.(A) (for complex ITensor, returns a complex ITensor storing the real part, so not technically wrong but not desirable).
• A .+= c for c != 0 should error for QN ITensors (see comment below).
• A .= f.(B) should error for f(0) != 0 for QN ITensors (see comment below).

[mtfishman]

Look into Broadcast.flatten as a more generic way to forward the function call to NDTensors. Hopefully that wouldn't lead to type instabilities and performance issues.

[mtfishman]

Some notes on A .+= c and A .= f.(B) for QN ITensors, continued from a discussion here: #636 (comment). Right now there is this behavior:

i = Index([QN(0) => 1, QN(1) => 1])
A = emptyITensor(i', dag(i))
A[1, 1] = 2
A .= one.(A)
@show A

gives:

A = ITensor ord=2
Dim 1: (dim=2|id=210)' <Out>
 1: QN(0) => 1
 2: QN(1) => 1
Dim 2: (dim=2|id=210) <In>
 1: QN(0) => 1
 2: QN(1) => 1
NDTensors.BlockSparse{Float64, Vector{Float64}, 2}
 2×2
Block(1, 1)
 [1:1, 1:1]
 1.0

A question would be if this should fill in the other possible QN values, or just not be allowed. Personally I think it should just be ill-defined for QN ITensors. I think the "rule" would be that for A .= f.(B), for QN ITensors, if f(0) != 0 then the operation is ill-defined, and we should provide a different function like nonzeros(A) .= f.(nonzeros(B)) if you want the operation to only act on the nonzero values. This logic is a bit tricky to implement the way the ITensor broadcasting is written right now, it needs a rewriting.

For operations like A .+= c, currently this does the incorrect things for QN ITensors:

i = Index([QN(0) => 1, QN(1) => 1])
A = emptyITensor(i', dag(i))
A[1, 1] = 1
A .+= 1.5
@show A

Gives:

A = ITensor ord=2
Dim 1: (dim=2|id=571)' <Out>
 1: QN(0) => 1
 2: QN(1) => 1
Dim 2: (dim=2|id=571) <In>
 1: QN(0) => 1
 2: QN(1) => 1
NDTensors.BlockSparse{Float64, Vector{Float64}, 2}
 2×2
Block(1, 1)
 [1:1, 1:1]
 2.5

Block(2, 1)
 [2:2, 1:1]
 1.5

Block(1, 2)
 [1:1, 2:2]
 1.5

Block(2, 2)
 [2:2, 2:2]
 1.5

A question would be about how this should be defined for QN ITensors. Should it only add to the existing QN blocks, any possible nonzero flux block, etc? I think it should be disallowed, and adding to just the nonzero values or all of the possible flux-compatible nonzero values should be done with specialized syntax to make it clear that is what is going on. For example, Julia doesn't allow this operation for sparse arrays:

julia> using LinearAlgebra

julia> D = Diagonal(randn(4))
4×4 Diagonal{Float64, Vector{Float64}}:
 -1.28286   ⋅         ⋅         ⋅ 
   ⋅       0.15365    ⋅         ⋅ 
   ⋅        ⋅       -1.38377    ⋅ 
   ⋅        ⋅         ⋅       -1.87302

julia> D .+= 2.4
ERROR: ArgumentError: cannot set off-diagonal entry (2, 1) to a nonzero value (2.4)

I don't think that is a common enough operation to do to warrant making such a convenient and potentially confusing operation. I think for that usage we should provide separate functionality for that, like:

insert_zero_blocks!(A) # Insert blocks of zero consistent with the flux
nonzeros(A) .+= 1.5  # Now broadcast only on the nonzero values of A

Or alternative, you construct another QN ITensor B which is filled with the value 1.5 and add it in-place to A, which may be just as efficient as the above. Another reason for not making the above operations default behavior for broadcasting is that in general the functionality insert_zero_blocks!(A) can be very slow, so now we only do it when necessary, like in constructors. (Note that insert_zero_blocks!(A) doesn't actually exist in that form but that kind of operation is used in QN constructors in a different form.)