diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py
index e9f5dbd5..075b9605 100644
--- a/pyerrors/linalg.py
+++ b/pyerrors/linalg.py
@@ -187,40 +187,46 @@ def jack_matmul(*operands):
     """
 
     if any(isinstance(o[0, 0], CObs) for o in operands):
+        name = operands[0][0, 0].real.names[0]
+        idl = operands[0][0, 0].real.idl[name]
+
         def _exp_to_jack(matrix):
             base_matrix = np.empty_like(matrix)
             for (n, m), entry in np.ndenumerate(matrix):
                 base_matrix[n, m] = entry.real.export_jackknife() + 1j * entry.imag.export_jackknife()
             return base_matrix
 
-        def _imp_from_jack(matrix, name):
+        def _imp_from_jack(matrix):
             base_matrix = np.empty_like(matrix)
             for (n, m), entry in np.ndenumerate(matrix):
-                base_matrix[n, m] = CObs(import_jackknife(entry.real, name),
-                                         import_jackknife(entry.imag, name))
+                base_matrix[n, m] = CObs(import_jackknife(entry.real, name, [idl]),
+                                         import_jackknife(entry.imag, name, [idl]))
             return base_matrix
 
         r = _exp_to_jack(operands[0])
         for op in operands[1:]:
             r = r @ _exp_to_jack(op)
-        return _imp_from_jack(r, op.ravel()[0].real.names[0])
+        return _imp_from_jack(r)
     else:
+        name = operands[0][0, 0].names[0]
+        idl = operands[0][0, 0].idl[name]
+
         def _exp_to_jack(matrix):
             base_matrix = np.empty_like(matrix)
             for (n, m), entry in np.ndenumerate(matrix):
                 base_matrix[n, m] = entry.export_jackknife()
             return base_matrix
 
-        def _imp_from_jack(matrix, name):
+        def _imp_from_jack(matrix):
             base_matrix = np.empty_like(matrix)
             for (n, m), entry in np.ndenumerate(matrix):
-                base_matrix[n, m] = import_jackknife(entry, name)
+                base_matrix[n, m] = import_jackknife(entry, name, [idl])
             return base_matrix
 
         r = _exp_to_jack(operands[0])
         for op in operands[1:]:
             r = r @ _exp_to_jack(op)
-        return _imp_from_jack(r, op.ravel()[0].names[0])
+        return _imp_from_jack(r)
 
 
 def inv(x):
diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 1e594ea4..56a67e57 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -1559,7 +1559,7 @@ def load_object(path):
         return pickle.load(file)
 
 
-def import_jackknife(jacks, name):
+def import_jackknife(jacks, name, idl=None):
     """Imports jackknife samples and returns an Obs
 
     Parameters
@@ -1573,7 +1573,7 @@ def import_jackknife(jacks, name):
     length = len(jacks) - 1
     prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
     samples = jacks[1:] @ prj
-    new_obs = Obs([samples], [name])
+    new_obs = Obs([samples], [name], idl=idl)
     new_obs._value = jacks[0]
     return new_obs
 
diff --git a/tests/linalg_test.py b/tests/linalg_test.py
index 4cfa9107..5af84b52 100644
--- a/tests/linalg_test.py
+++ b/tests/linalg_test.py
@@ -92,6 +92,19 @@ def test_multi_dot():
             assert e.is_zero(), t
 
 
+def test_jack_multi_dot():
+    for dim in [2, 4, 8]:
+        my_array = get_real_matrix(dim)
+
+        tt = pe.linalg.jack_matmul(my_array, my_array, my_array) - pe.linalg.matmul(my_array, my_array, my_array)
+
+        for t, e in np.ndenumerate(tt):
+            e.gamma_method()
+            assert e.is_zero_within_error(0.01)
+            assert e.is_zero(atol=1e-1), t
+            assert np.isclose(e.value, 0.0)
+
+
 def test_matmul_irregular_histories():
     dim = 2
     length = 500