Avoid copy in matrix contract

This reverts commit c9db523.

Author: lkdvos

src/algorithms/derivatives/hamiltonian_derivatives.jl

@@ -113,10 +113,8 @@ function prepare_operator!!(
         H.E
     end
 
-    # O3′ = Core.Compiler.return_type(prepare_operator!!, Tuple{O3, typeof(backend), typeof(allocator)})
-    # A = ismissing(H.A) ? H.A : prepare_operator!!(H.A, backend, allocator)
-    O3′ = O3
-    A = H.A
+    O3′ = Core.Compiler.return_type(prepare_operator!!, Tuple{O3, typeof(backend), typeof(allocator)})
+    A = ismissing(H.A) ? H.A : prepare_operator!!(H.A, backend, allocator)
 
     return JordanMPO_AC_Hamiltonian{O1, O2, O3′}(D, I, E, C, B, A)
 end
@@ -330,10 +328,8 @@ function prepare_operator!!(
         H.EE
     end
 
-    # O4′ = Core.Compiler.return_type(prepare_operator!!, Tuple{O4, typeof(backend), typeof(allocator)})
-    # AA = ismissing(H.AA) ? H.AA : prepare_operator!!(H.AA, backend, allocator)
-    O4′ = O4
-    AA = H.AA
+    O4′ = Core.Compiler.return_type(prepare_operator!!, Tuple{O4, typeof(backend), typeof(allocator)})
+    AA = prepare_operator!!(H.AA, backend, allocator)
 
     return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4′}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE)
 end

src/algorithms/derivatives/mpo_derivatives.jl

@@ -183,24 +183,24 @@ function (H::PrecomputedDerivative)(x::AbstractTensorMap)
             @inbounds sz, str, offset = structure_R.fusiontreestructure[i]
             r = TensorKit.Strided.StridedView(R_fused.data, sz, str, offset)
 
-            # if sz[2] < sz[3]
-            #     for k in axes(r, 2)
-            #         C = xr[:, k, :]
-            #         B = r[:, k, :]
-            #         LinearAlgebra.BLAS.gemm!('N', 'N', one(TC), x, B, zero(TC), C)
-            #     end
-            # else
-            #     for k in axes(r, 3)
-            #         C = xr[:, :, k]
-            #         B = r[:, :, k]
-            #         LinearAlgebra.BLAS.gemm!('N', 'N', one(TC), x, B, zero(TC), C)
-            #     end
-            # end
-
-            TensorOperations.tensorcontract!(
-                xr, x, ((1,), (2,)), false,
-                r, ((1,), (2, 3)), false, ((1, 2), (3,)), One(), Zero(), H.backend, H.allocator
-            )
+            if TensorOperations.isblascontractable(r, ((1,), (2, 3))) &&
+                    TensorOperations.isblasdestination(xr, ((1,), (2, 3)))
+                C = TensorKit.Strided.sreshape(xr, size(xr, 1), size(xr, 2) * size(xr, 3))
+                B = TensorKit.Strided.sreshape(r, size(r, 1), size(r, 2) * size(r, 3))
+                LinearAlgebra.BLAS.gemm!('N', 'N', one(TC), x, B, zero(TC), C)
+            elseif sz[2] < sz[3]
+                for k in axes(r, 2)
+                    C = xr[:, k, :]
+                    B = r[:, k, :]
+                    LinearAlgebra.BLAS.gemm!('N', 'N', one(TC), x, B, zero(TC), C)
+                end
+            else
+                for k in axes(r, 3)
+                    C = xr[:, :, k]
+                    B = r[:, :, k]
+                    LinearAlgebra.BLAS.gemm!('N', 'N', one(TC), x, B, zero(TC), C)
+                end
+            end
         else
             zerovector!(xr)
         end