[ITensors] Fix issue with directsum involving EmptyStorage (#1443)

src/tensor_operations/tensor_algebra.jl

@@ -254,12 +254,23 @@ function directsum_projectors(
   # Or with new notation:
   # D1 = zeros(elt1, dag(i), ij)
   # D2 = zeros(elt1, dag(j), ij)
-  D1 = zeros_itensor(elt1, dag(i), ij)
-  D2 = zeros_itensor(elt1, dag(j), ij)
+  elt = promote_type(elt1, elt2)
+  D1 = zeros_itensor(elt, dag(i), ij)
+  D2 = zeros_itensor(elt, dag(j), ij)
   directsum_projectors!(tensor(D1), tensor(D2))
   return D1, D2
 end
 
+function directsum_projectors(
+  ::Type{<:EmptyNumber}, ::Type{<:EmptyNumber}, ::Index, ::Index, ::Index
+)
+  return error(
+    "It is not possible to call directsum on two tensors with element type EmptyNumber.
+If you are inputting ITensors constructed like `ITensor(i, j)`, try specifying the element type, 
+e.g. `ITensor(Float64, i, j)`, or fill them with zero values, e.g. `ITensor(zero(Float64), i, j)`.",
+  )
+end
+
 function check_directsum_inds(A::ITensor, I, B::ITensor, J)
   a = uniqueinds(A, I)
   b = uniqueinds(B, J)

test/base/test_itensor.jl

@@ -1805,6 +1805,8 @@ end
     A = randomITensor(i1, i2, j)
     B = randomITensor(i1, i2, k)
     C = randomITensor(i1, i2, l)
+    D = ITensor(i1, i2, k)
+    F = ITensor(i1, i2, j)
 
     S, s = directsum(A => index_op(j), B => index_op(k))
     @test dim(s) == dim(j) + dim(k)
@@ -1814,6 +1816,16 @@ end
     @test dim(s) == dim(j) + dim(k)
     @test hassameinds(S, (i1, i2, s))
 
+    S, s = (A => index_op(j)) ⊕ (D => index_op(k))
+    @test dim(s) == dim(j) + dim(k)
+    @test hassameinds(S, (i1, i2, s))
+
+    @test_throws ErrorException (F => index_op(j)) ⊕ (D => index_op(k))
+
+    S, s = (D => index_op(k)) ⊕ (A => index_op(j))
+    @test dim(s) == dim(j) + dim(k)
+    @test hassameinds(S, (i1, i2, s))
+
     S, s = directsum(A => index_op(j), B => index_op(k), C => index_op(l))
     @test dim(s) == dim(j) + dim(k) + dim(l)
     @test hassameinds(S, (i1, i2, s))