Support for eigen of dual operator

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensors"
 uuid = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>", "Miles Stoudenmire <mstoudenmire@flatironinstitute.org>"]
-version = "0.9.10"
+version = "0.9.11"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/tensor_operations/matrix_decomposition.jl

@@ -353,17 +353,22 @@ function eigen(
   end
 
   # <fermions>
+  L_arrow_dir = Out
   if hasqns(A) && using_auto_fermion()
-    if !all(i -> dir(i) == Out, Lis)
-      error("With auto_fermion enabled, left inds in eigen must have Out arrows")
-    end
-    if !all(i -> dir(i) == In, Ris)
-      error("With auto_fermion enabled, right inds in eigen must have Out arrows")
+    # Make arrows of combined ITensor match those of index sets
+    if all(i -> dir(i) == Out, Lis) && all(i -> dir(i) == In, Ris)
+      L_arrow_dir = Out
+    elseif all(i -> dir(i) == In, Lis) && all(i -> dir(i) == Out, Ris)
+      L_arrow_dir = In
+    else
+      error(
+        "With auto_fermion enabled, index sets in eigen must have all arrows the same, and opposite between the sets",
+      )
     end
   end
 
-  CL = combiner(Lis...; dir=Out, tags="CMB,left")
-  CR = combiner(dag(Ris)...; dir=Out, tags="CMB,right")
+  CL = combiner(Lis...; dir=L_arrow_dir, tags="CMB,left")
+  CR = combiner(dag(Ris)...; dir=L_arrow_dir, tags="CMB,right")
 
   AC = A * dag(CR) * CL
 

test/base/test_decomp.jl

@@ -500,4 +500,44 @@ end
     @test norm(U * B - phi) < 1E-5
     @test dim(commonind(U, B)) <= 4
   end
+
+  @testset "Eigen of Fermionic Matrices" begin
+    ITensors.enable_auto_fermion()
+    s = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "s,Site,Fermion")
+    t = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "t,Site,Fermion")
+
+    #
+    # HPSD Operator (Out,In) case
+    #
+    M = random_itensor(s, dag(t))
+    O = prime(M, s)*dag(M)
+
+    @test dir(inds(O)[1]) == ITensors.Out
+    @test dir(inds(O)[2]) == ITensors.In
+    linds = [s']
+    rinds = [dag(s)]
+    D_O, U = eigen(O, linds, rinds; ishermitian=true)
+    @test norm(prime(U)*D_O*dag(U)-O) < 1E-10
+    @test all(>=(0.0), diag(array(D_O)))
+
+    #
+    # HPSD Dual operator (In,Out) case
+    #
+    # Make ρ out of two squared states
+    # to populate both blocks: (0,0) and (1,1)
+    ψ0 = random_itensor(t, s)
+    ρ0 = prime(dag(ψ0), s)*ψ0
+
+    ψ2 = random_itensor(QN("Nf", 2, -1), t, s)
+    ρ2 = prime(dag(ψ2), s)*ψ2
+
+    ρ = ρ0/2 + ρ2/2
+    @test dir(inds(ρ)[1]) == ITensors.In
+    @test dir(inds(ρ)[2]) == ITensors.Out
+
+    D_ρ, U = eigen(ρ, [dag(s)'], [s]; ishermitian=true)
+    @test all(>=(0.0), diag(array(D_ρ)))
+    @test norm(prime(U)*D_ρ*dag(U)-ρ) < 1E-10
+    ITensors.disable_auto_fermion()
+  end
 end