diff --git a/src/contractMPO.jl b/src/contractMPO.jl
index c471b1a..801eddd 100644
--- a/src/contractMPO.jl
+++ b/src/contractMPO.jl
@@ -3,6 +3,8 @@ function contract_mpo_mpo(M1::MPO, M2::MPO; alg::String = "densitymatrix", kwarg
         return contract_densitymatrix(M1, M2; kwargs...)
     elseif alg == "fit"
         return contract_fit(M1, M2; kwargs...)
+    elseif alg == "zipup"
+        return ITensors.contract(M1, M2; alg = "zipup", kwargs...)
     elseif alg == "naive"
         return ITensors.contract(M1, M2; alg = "naive", kwargs...)
     else
diff --git a/test/densitymatrix.jl b/test/contractMPO.jl
similarity index 80%
rename from test/densitymatrix.jl
rename to test/contractMPO.jl
index 0d97f55..0d40a96 100644
--- a/test/densitymatrix.jl
+++ b/test/contractMPO.jl
@@ -2,10 +2,13 @@ using Test
 
 import FastMPOContractions as FMPOC
 using ITensors
+using Random
 
+algs = ["densitymatrix", "fit", "zipup"]
 
+@testset "MPO-MPO contraction (x-y-z)" for alg in algs
+    Random.seed!(1234)
 
-@testset "densitymatrix (x-y-z)" begin
     R = 3
     sitesx = [Index(2, "Qubit,x=$n") for n = 1:R]
     sitesy = [Index(2, "Qubit,y=$n") for n = 1:R]
@@ -16,11 +19,11 @@ using ITensors
     a = _random_mpo(sitesa)
     b = _random_mpo(sitesb)
     ab_ref = contract(a, b; alg = "naive")
-    ab = FMPOC.contract_densitymatrix(a, b)
+    ab = FMPOC.contract_mpo_mpo(a, b; alg=alg)
     @test ab_ref ≈ ab
 end
 
-@testset "densitymatrix (xk-y-z)" begin
+@testset "MPO-MPO contraction (xk-y-z)" for alg in algs
     R = 3
     sitesx = [Index(2, "Qubit,x=$n") for n = 1:R]
     sitesk = [Index(2, "Qubit,k=$n") for n = 1:R]
@@ -32,11 +35,11 @@ end
     a = _random_mpo(sitesa)
     b = _random_mpo(sitesb)
     ab_ref = contract(a, b; alg = "naive")
-    ab = FMPOC.contract_densitymatrix(a, b)
+    ab = FMPOC.contract_mpo_mpo(a, b; alg=alg)
     @test ab_ref ≈ ab
 end
 
-@testset "densitymatrix (xk-y-zl)" begin
+@testset "MPO-MPO contraction (xk-y-zl)" for alg in algs
     R = 3
     sitesx = [Index(2, "Qubit,x=$n") for n = 1:R]
     sitesk = [Index(2, "Qubit,k=$n") for n = 1:R]
@@ -49,11 +52,11 @@ end
     a = _random_mpo(sitesa)
     b = _random_mpo(sitesb)
     ab_ref = contract(a, b; alg = "naive")
-    ab = FMPOC.contract_densitymatrix(a, b)
+    ab = FMPOC.contract_mpo_mpo(a, b; alg=alg)
     @test ab_ref ≈ ab
 end
 
-@testset "densitymatrix (xk-ym-zl)" begin
+@testset "MPO-MPO contraction (xk-ym-zl)" for alg in algs
     R = 3
     sitesx = [Index(2, "Qubit,x=$n") for n = 1:R]
     sitesk = [Index(2, "Qubit,k=$n") for n = 1:R]
@@ -67,6 +70,6 @@ end
     a = _random_mpo(sitesa)
     b = _random_mpo(sitesb)
     ab_ref = contract(a, b; alg = "naive")
-    ab = FMPOC.contract_densitymatrix(a, b)
+    ab = FMPOC.contract_mpo_mpo(a, b; alg=alg)
     @test ab_ref ≈ ab
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 952c3d9..51f4158 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,7 +6,7 @@ include("test_with_jet.jl")
 
 include("_util.jl")
 
-include("densitymatrix.jl")
+include("contractMPO.jl")
 include("fitalgorithm.jl")
 include("util.jl")
 include("fitalgorithm_sum.jl")