diff --git a/modules/dnn/perf/perf_slice.cpp b/modules/dnn/perf/perf_slice.cpp
new file mode 100644
index 0000000000..06d81cc3dc
--- /dev/null
+++ b/modules/dnn/perf/perf_slice.cpp
@@ -0,0 +1,93 @@
+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html.
+
+#include "perf_precomp.hpp"
+
+namespace opencv_test {
+
+struct Layer_Slice : public TestBaseWithParam<tuple<Backend, Target>>
+{
+    void test_slice(const std::vector<int>& input_shape, int axis, int begin, int end, int step = 1)
+    {
+        int backendId = get<0>(GetParam());
+        int targetId = get<1>(GetParam());
+
+        Mat data(input_shape, CV_32FC1);
+        randu(data, 0.f, 1.f);
+
+        Net net;
+        LayerParams lp;
+        lp.type = "Slice";
+        lp.name = "testLayer";
+        lp.set("axis", axis);
+
+        std::vector<int> begins(input_shape.size(), 0);
+        std::vector<int> ends = input_shape;
+        std::vector<int> steps(input_shape.size(), 1);
+
+        begins[axis] = begin;
+        ends[axis] = end;
+        steps[axis] = step;
+
+        lp.set("begin", DictValue::arrayInt(&begins[0], begins.size()));
+        lp.set("end", DictValue::arrayInt(&ends[0], ends.size()));
+        if (step != 1) {
+            lp.set("steps", DictValue::arrayInt(&steps[0], steps.size()));
+        }
+
+        int id = net.addLayerToPrev(lp.name, lp.type, lp);
+        net.connect(0, 0, id, 0);
+
+        net.setInputsNames({"data"});
+
+        // warmup
+        {
+            net.setInput(data, "data");
+            net.setPreferableBackend(backendId);
+            net.setPreferableTarget(targetId);
+            Mat out = net.forward();
+        }
+
+        TEST_CYCLE()
+        {
+            Mat res = net.forward();
+        }
+
+        SANITY_CHECK_NOTHING();
+    }
+};
+
+PERF_TEST_P_(Layer_Slice, Slice_Contiguous_Axis0)
+{
+    test_slice({64, 128, 128}, 0, 10, 54);
+}
+
+PERF_TEST_P_(Layer_Slice, Slice_Contiguous_Axis2)
+{
+    test_slice({64, 128, 128}, 2, 10, 118);
+}
+
+PERF_TEST_P_(Layer_Slice, Slice_Small_Middle)
+{
+    test_slice({32, 64, 32}, 1, 20, 40);
+}
+
+PERF_TEST_P_(Layer_Slice, Slice_Strided_Axis0_Step2)
+{
+    // Strided slice on outer axis.
+    // [64, 128, 128] -> [0:64:2, ...]
+    test_slice({64, 128, 128}, 0, 0, 64, 2);
+}
+
+PERF_TEST_P_(Layer_Slice, Slice_Strided_Axis2_Step2)
+{
+    // Strided slice on inner axis.
+    // [64, 128, 128] -> [..., 0:128:2]
+    test_slice({64, 128, 128}, 2, 0, 128, 2);
+}
+
+
+INSTANTIATE_TEST_CASE_P(/**/, Layer_Slice, dnnBackendsAndTargets(false, false, true, false, false, false, false, false));
+
+} // namespace opencv_test
diff --git a/modules/dnn/src/layers/slice2_layer.cpp b/modules/dnn/src/layers/slice2_layer.cpp
index cd3a460bef..d7edd7148b 100644
--- a/modules/dnn/src/layers/slice2_layer.cpp
+++ b/modules/dnn/src/layers/slice2_layer.cpp
@@ -26,109 +26,7 @@ namespace dnn
     Opset's 1 to 13 are covered.
 */
 
-/* Slice op for CPU.
-   starts_, ends_ and steps_ must contain as many elements as
-   the dimensionality in inp and out.
-*/
-static void slice(const Mat& inp, const int* starts_,
-                  const int*, const int* steps_,
-                  Mat& out)
-{
-    /// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    /// in this function steps can be negative, so
-    /// please don't replace int64_t's with size_t's
-    /// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    enum {SLICE_MAX_DIMS=7};
 
-    CV_Assert_N(inp.isContinuous(), out.isContinuous());
-    CV_Assert(inp.type() == out.type());
-    CV_Assert_N(inp.dims <= SLICE_MAX_DIMS, inp.dims == out.dims);
-
-    MatShape inpShape = inp.shape();
-    MatShape outShape = out.shape();
-    int64_t esz = (int64_t)inp.elemSize();
-
-    int ndims = inpShape.dims;
-    int starts[SLICE_MAX_DIMS], steps[SLICE_MAX_DIMS];
-    int inpsz[SLICE_MAX_DIMS], outsz[SLICE_MAX_DIMS];
-    int64_t inpstep[SLICE_MAX_DIMS];
-
-    int delta = SLICE_MAX_DIMS - ndims;
-    bool emptyOut = false;
-
-    for (int i = 0; i < SLICE_MAX_DIMS; i++) {
-        inpsz[i] = outsz[i] = steps[i] = 1;
-        starts[i] = 0;
-    }
-
-    for (int i = 0; i < ndims; i++) {
-        inpsz[delta + i] = inpShape[i];
-        outsz[delta + i] = outShape[i];
-        starts[delta + i] = starts_[i];
-        steps[delta + i] = steps_[i];
-        if (outShape[i] == 0)
-            emptyOut = true;
-    }
-
-    for (int i = SLICE_MAX_DIMS-1; i >= 0; i--)
-        inpstep[i] = i == SLICE_MAX_DIMS-1 ? 1 : inpstep[i+1]*inpsz[i+1];
-
-    const uchar* inptr0 = inp.data;
-
-    for (int i = 0; i < SLICE_MAX_DIMS; i++) {
-        inptr0 += starts[i]*inpstep[i]*esz;
-        inpstep[i] *= steps[i];
-    }
-
-    int sz0 = outsz[6], sz1 = outsz[5];
-    int sz2 = outsz[4], sz3 = outsz[3];
-    int sz4 = outsz[2], sz5 = outsz[1], sz6 = outsz[0];
-    int64_t p0 = inpstep[6], p1 = inpstep[5];
-    int64_t p2 = inpstep[4], p3 = inpstep[3];
-    int64_t p4 = inpstep[2], p5 = inpstep[1], p6 = inpstep[0];
-
-    #undef CV_IMPLEMENT_SLICE
-    #define CV_IMPLEMENT_SLICE(typ) \
-        typ* outptr = (typ*)(out.data); \
-        for(int i6 = 0; i6 < sz6; i6++) { \
-        for(int i5 = 0; i5 < sz5; i5++) { \
-        for(int i4 = 0; i4 < sz4; i4++) { \
-        for(int i3 = 0; i3 < sz3; i3++) { \
-        for(int i2 = 0; i2 < sz2; i2++) { \
-        for(int i1 = 0; i1 < sz1; i1++, outptr += sz0) { \
-            const typ* inptr = (const typ*)inptr0 + i6*p6 + \
-                    i5*p5 + i4*p4 + i3*p3 + i2*p2 + i1*p1; \
-            int i0 = 0; \
-            if (p0 == 1) { \
-                for (; i0 < sz0; i0++) \
-                    outptr[i0] = inptr[i0]; \
-            } \
-            else { \
-                for (; i0 <= sz0 - 4; i0 += 4) { \
-                    int64_t ip0 = i0*p0; \
-                    typ t0 = inptr[ip0], t1 = inptr[ip0 + p0]; \
-                    typ t2 = inptr[ip0 + p0*2], t3 = inptr[ip0 + p0*3]; \
-                    outptr[i0] = t0; outptr[i0+1] = t1; \
-                    outptr[i0+2] = t2; outptr[i0+3] = t3; \
-                } \
-                for (; i0 < sz0; i0++) \
-                    outptr[i0] = inptr[i0*p0]; \
-            } \
-        }}}}}}
-
-    if (emptyOut) return;
-    if (esz == 4) {
-        CV_IMPLEMENT_SLICE(int)
-    } else if (esz == 2) {
-        CV_IMPLEMENT_SLICE(int16_t)
-    } else if (esz == 1) {
-        CV_IMPLEMENT_SLICE(int8_t)
-    } else if (esz == 8) {
-        CV_IMPLEMENT_SLICE(int64_t)
-    } else {
-        CV_Error(Error::StsNotImplemented, "");
-    }
-}
 
 class Slice2LayerImpl CV_FINAL : public Slice2Layer
 {
@@ -281,6 +179,160 @@ public:
     {
     }
 
+
+private:
+    template <typename T>
+    class ParallelSlice : public cv::ParallelLoopBody
+    {
+    public:
+        ParallelSlice(const Mat& inp, Mat& out,
+                      const std::vector<Range>& ranges,
+                      const std::vector<int>& steps)
+            : inp_(inp), out_(out), ranges_(ranges), steps_(steps)
+        {
+            dims_ = inp.dims;
+            es_ = inp.elemSize();
+
+            inp_strides_.resize(dims_);
+            out_strides_.resize(dims_);
+            for(int i=0; i<dims_; ++i) {
+                inp_strides_[i] = inp.step.p[i];
+                out_strides_[i] = out.step.p[i];
+            }
+        }
+
+        void operator()(const Range& range) const CV_OVERRIDE
+        {
+            int b = ranges_[0].start;
+            int s = steps_[0];
+
+            const uchar* src_base = inp_.ptr();
+            uchar* dst_base = out_.ptr();
+
+            for (int i = range.start; i < range.end; ++i)
+            {
+                int k = b + i * s;
+                size_t src_offset = k * inp_strides_[0];
+                size_t dst_offset = i * out_strides_[0];
+                if (dims_ == 1)
+                    std::memcpy(dst_base + dst_offset, src_base + src_offset, es_);
+                else
+                    recursive_copy(1, src_base + src_offset, dst_base + dst_offset);
+            }
+        }
+
+        void recursive_copy(int dim, const uchar* src_ptr, uchar* dst_ptr) const
+        {
+            if (dim >= dims_) return;
+
+            int begin = ranges_[dim].start;
+            int end = ranges_[dim].end;
+            int step = steps_[dim];
+
+            if (dim == dims_ - 1)
+            {
+                if (step == 1)
+                {
+                    size_t count = end - begin;
+                    std::memcpy(dst_ptr, src_ptr + begin * inp_strides_[dim], count * es_);
+                }
+                else
+                {
+                    const uchar* s_ptr = src_ptr + begin * inp_strides_[dim];
+                    uchar* d_ptr = dst_ptr;
+                    size_t s_stride = step * inp_strides_[dim];
+                    size_t d_stride = out_strides_[dim];
+
+                    if (step > 0)
+                    {
+                        for (int k = begin; k < end; k += step)
+                        {
+                            *(T*)d_ptr = *(const T*)s_ptr;
+                            s_ptr += s_stride;
+                            d_ptr += d_stride;
+                        }
+                    }
+                    else
+                    {
+                        for (int k = begin; k > end; k += step)
+                        {
+                            *(T*)d_ptr = *(const T*)s_ptr;
+                            s_ptr += s_stride;
+                            d_ptr += d_stride;
+                        }
+                    }
+                }
+                return;
+            }
+
+            if (step == 1 && is_fully_contiguous(dim))
+            {
+               size_t count = end - begin;
+               size_t bytes = count * inp_strides_[dim];
+               std::memcpy(dst_ptr, src_ptr + begin * inp_strides_[dim], bytes);
+               return;
+            }
+
+            size_t src_stride = step * inp_strides_[dim];
+            size_t dst_stride = out_strides_[dim];
+
+            const uchar* s_ptr = src_ptr + begin * inp_strides_[dim];
+            uchar* d_ptr = dst_ptr;
+
+            if (step > 0)
+            {
+                for (int k = begin; k < end; k += step)
+                {
+                    recursive_copy(dim + 1, s_ptr, d_ptr);
+                    s_ptr += src_stride;
+                    d_ptr += dst_stride;
+                }
+            }
+            else
+            {
+                for (int k = begin; k > end; k += step)
+                {
+                    recursive_copy(dim + 1, s_ptr, d_ptr);
+                    s_ptr += src_stride;
+                    d_ptr += dst_stride;
+                }
+            }
+        }
+
+        bool is_fully_contiguous(int dim) const
+        {
+            size_t expected_step = es_;
+            for (int d = dims_ - 1; d >= dim; --d)
+            {
+                 if (inp_.step[d] != expected_step) return false;
+                 if (d > dim) {
+                     if (steps_[d] != 1) return false;
+                     if (ranges_[d].start != 0 || ranges_[d].end != inp_.size[d]) return false;
+                     expected_step *= inp_.size[d];
+                 }
+            }
+            return true;
+        }
+
+    private:
+        const Mat& inp_;
+        Mat& out_;
+        const std::vector<Range>& ranges_;
+        const std::vector<int>& steps_;
+        int dims_;
+        size_t es_;
+        std::vector<size_t> inp_strides_;
+        std::vector<size_t> out_strides_;
+    };
+
+    template <typename T>
+    void run_parallel(const Mat& inp, Mat& out, const std::vector<Range>& ranges, const std::vector<int>& steps)
+    {
+        ParallelSlice<T> body(inp, out, ranges, steps);
+        int dim0_size = out.size[0];
+        parallel_for_(Range(0, dim0_size), body);
+    }
+
     void forward(InputArrayOfArrays inputs_arr,
                  OutputArrayOfArrays outputs_arr,
                  OutputArrayOfArrays) CV_OVERRIDE
@@ -320,8 +372,14 @@ public:
         MatShape outShape = getOutShape(inpShape, *starts_, *ends_, *axes_, steps,
                                         allStarts, allEnds, allSteps);
 
-        int outKind = outputs_arr.kind();
+        std::vector<Range> ranges;
+        std::vector<int> steps_vec;
+        for (int i = 0; i < inpShape.dims; ++i) {
+            ranges.push_back(Range(allStarts[i], allEnds[i]));
+            steps_vec.push_back(allSteps[i]);
+        }
 
+        int outKind = outputs_arr.kind();
         CV_Assert(outKind == _InputArray::STD_VECTOR_MAT ||
                   outKind == _InputArray::STD_VECTOR_UMAT);
 
@@ -330,23 +388,31 @@ public:
             std::vector<Mat>& outs = outputs_arr.getMatVecRef();
             outs.resize(1);
             outs[0].fit(outShape, inpType);
-            runOp(inp, allStarts, allEnds, allSteps, outs[0]);
-        } else {
-            // [TODO] more efficient OpenCL implementation
-            Mat inp = inputs_arr.getMat(0);
-            std::vector<UMat>& outs = outputs_arr.getUMatVecRef();
-            outs.resize(1);
-            outs[0].fit(outShape, inpType);
-            Mat temp(outShape, inpType);
-            runOp(inp, allStarts, allEnds, allSteps, temp);
-            temp.copyTo(outs[0]);
-        }
-    }
 
-    void runOp(const Mat& inp, const int* starts_,
-               const int* ends_, const int* steps_, Mat& out)
-    {
-        slice(inp, starts_, ends_, steps_, out);
+            if (inp.depth() == CV_32S) run_parallel<int32_t>(inp, outs[0], ranges, steps_vec);
+            else if (inp.depth() == CV_64S) run_parallel<int64_t>(inp, outs[0], ranges, steps_vec);
+            else if (inp.depth() == CV_16F) run_parallel<int16_t>(inp, outs[0], ranges, steps_vec);
+            else if (inp.depth() == CV_8S) run_parallel<int8_t>(inp, outs[0], ranges, steps_vec);
+            else if (inp.depth() == CV_8U) run_parallel<uint8_t>(inp, outs[0], ranges, steps_vec);
+            else if (inp.depth() == CV_Bool) run_parallel<uint8_t>(inp, outs[0], ranges, steps_vec);
+            else run_parallel<float>(inp, outs[0], ranges, steps_vec);
+        } else {
+             Mat inp = inputs_arr.getMat(0);
+             std::vector<UMat>& outs = outputs_arr.getUMatVecRef();
+             outs.resize(1);
+             outs[0].fit(outShape, inpType);
+             Mat temp(outShape, inpType);
+
+             if (inp.depth() == CV_32S) run_parallel<int32_t>(inp, temp, ranges, steps_vec);
+             else if (inp.depth() == CV_64S) run_parallel<int64_t>(inp, temp, ranges, steps_vec);
+             else if (inp.depth() == CV_16F) run_parallel<int16_t>(inp, temp, ranges, steps_vec);
+             else if (inp.depth() == CV_8S) run_parallel<int8_t>(inp, temp, ranges, steps_vec);
+             else if (inp.depth() == CV_8U) run_parallel<uint8_t>(inp, temp, ranges, steps_vec);
+             else if (inp.depth() == CV_Bool) run_parallel<uint8_t>(inp, temp, ranges, steps_vec);
+             else run_parallel<float>(inp, temp, ranges, steps_vec);
+
+             temp.copyTo(outs[0]);
+        }
     }
 };
 
diff --git a/modules/dnn/src/layers/slice_layer.cpp b/modules/dnn/src/layers/slice_layer.cpp
index 91264d2677..d86d66f6e7 100644
--- a/modules/dnn/src/layers/slice_layer.cpp
+++ b/modules/dnn/src/layers/slice_layer.cpp
@@ -191,7 +191,7 @@ public:
                     CV_Assert(step != 0);
                     if (step < 0)
                         neg_step_dims.push_back(i);
-                    if (std::abs(step) > 1)
+                    if (std::abs(step) > 1 || step < 0)
                         hasSteps = true;
                     sliceSteps[0][i] = step;
                 }
@@ -636,31 +636,18 @@ public:
                 }
             }
         }
-        else
+        if (hasSteps)
         {
-            int dimsNum = inpMat.dims;
+            if (inpMat.depth() == CV_32S) run_parallel<int32_t>(inpMat, outputs);
+            else if (inpMat.depth() == CV_64S) run_parallel<int64_t>(inpMat, outputs);
+            else if (inpMat.depth() == CV_16F) run_parallel<int16_t>(inpMat, outputs);
+            else if (inpMat.depth() == CV_8S) run_parallel<int8_t>(inpMat, outputs);
+            else if (inpMat.depth() == CV_8U) run_parallel<uint8_t>(inpMat, outputs);
+            else if (inpMat.depth() == CV_Bool) run_parallel<uint8_t>(inpMat, outputs);
+            else run_parallel<float>(inpMat, outputs);
 
             for (size_t i = 0; i < outputs.size(); i++)
-            {
-                std::vector<int> inpIdx(dimsNum, 0);
-                std::vector<int> outIdx(dimsNum, 0);
-                if (inpMat.type() == CV_32S)
-                    getSliceRecursive<int32_t>(inpMat, inpIdx, finalSliceRanges[i], sliceSteps[i], 0, dimsNum, outputs[i], outIdx);
-                else if (inpMat.type() == CV_64S)
-                    getSliceRecursive<int64_t>(inpMat, inpIdx, finalSliceRanges[i], sliceSteps[i], 0, dimsNum, outputs[i], outIdx);
-                else if (inpMat.type() == CV_16F)
-                    getSliceRecursive<int16_t>(inpMat, inpIdx, finalSliceRanges[i], sliceSteps[i], 0, dimsNum, outputs[i], outIdx);
-                else if (inpMat.type() == CV_8S)
-                    getSliceRecursive<int8_t>(inpMat, inpIdx, finalSliceRanges[i], sliceSteps[i], 0, dimsNum, outputs[i], outIdx);
-                else if (inpMat.type() == CV_8U)
-                    getSliceRecursive<uint8_t>(inpMat, inpIdx, finalSliceRanges[i], sliceSteps[i], 0, dimsNum, outputs[i], outIdx);
-                else if (inpMat.type() == CV_Bool)
-                    getSliceRecursive<bool>(inpMat, inpIdx, finalSliceRanges[i], sliceSteps[i], 0, dimsNum, outputs[i], outIdx);
-                else
-                    getSliceRecursive<float>(inpMat, inpIdx, finalSliceRanges[i], sliceSteps[i], 0, dimsNum, outputs[i], outIdx);
-                // flip for negative steps
                 flip(outputs[i]);
-            }
         }
     }
 
@@ -844,28 +831,141 @@ public:
 
 private:
     template <typename T>
-    void getSliceRecursive(const Mat &inpMat, std::vector<int> &inpIdx,
-                           const std::vector<Range> &sliceRanges,
-                           const std::vector<int> &sliceSteps, int dim, int dimsNum,
-                           Mat &outputs, std::vector<int> &outIdx)
+    class ParallelSlice : public cv::ParallelLoopBody
     {
-        int begin = sliceRanges[dim].start;
-        int end = sliceRanges[dim].end;
-        int step = !sliceSteps.empty() ? sliceSteps[dim] : 1;
-
-        // TODO optimization is required (for 2D tail case at least)
-        for (int k = begin, j = 0; k < end; k += step, j++)
+    public:
+        ParallelSlice(const Mat& inp, Mat& out,
+                      const std::vector<Range>& ranges,
+                      const std::vector<int>& steps)
+            : inp_(inp), out_(out), ranges_(ranges), steps_(steps)
         {
-            inpIdx[dim] = k;
-            outIdx[dim] = j;
+            dims_ = inp.dims;
+            es_ = inp.elemSize();
 
-            if (dim + 1 < dimsNum)
-                getSliceRecursive<T>(inpMat, inpIdx, sliceRanges, sliceSteps, dim + 1, dimsNum, outputs, outIdx);
-            else
-                outputs.at<T>(outIdx.data()) = inpMat.at<T>(inpIdx.data());
+            inp_strides_.resize(dims_);
+            out_strides_.resize(dims_);
+            for(int i=0; i<dims_; ++i) {
+                inp_strides_[i] = inp.step.p[i];
+                out_strides_[i] = out.step.p[i];
+            }
+        }
+
+        void operator()(const Range& range) const CV_OVERRIDE
+        {
+            int b = ranges_[0].start;
+            int s = steps_[0];
+
+            const uchar* src_base = inp_.ptr();
+            uchar* dst_base = out_.ptr();
+
+            for (int i = range.start; i < range.end; ++i)
+            {
+                int k = b + i * s;
+                size_t src_offset = k * inp_strides_[0];
+                size_t dst_offset = i * out_strides_[0];
+                if (dims_ == 1)
+                    std::memcpy(dst_base + dst_offset, src_base + src_offset, es_);
+                else
+                    recursive_copy(1, src_base + src_offset, dst_base + dst_offset);
+            }
+        }
+
+        void recursive_copy(int dim, const uchar* src_ptr, uchar* dst_ptr) const
+        {
+            if (dim >= dims_) return;
+
+            int begin = ranges_[dim].start;
+            int end = ranges_[dim].end;
+            int step = steps_[dim];
+
+            if (dim == dims_ - 1)
+            {
+                if (step == 1)
+                {
+                    size_t count = end - begin;
+                    std::memcpy(dst_ptr, src_ptr + begin * inp_strides_[dim], count * es_);
+                }
+                else
+                {
+                    const uchar* s_ptr = src_ptr + begin * inp_strides_[dim];
+                    uchar* d_ptr = dst_ptr;
+                    size_t s_stride = step * inp_strides_[dim];
+                    size_t d_stride = out_strides_[dim];
+
+                    for (int k = begin; k < end; k += step)
+                    {
+                        *(T*)d_ptr = *(const T*)s_ptr;
+                        s_ptr += s_stride;
+                        d_ptr += d_stride;
+                    }
+                }
+                return;
+            }
+
+            if (step == 1 && is_fully_contiguous(dim))
+            {
+               size_t count = end - begin;
+               size_t bytes = count * inp_strides_[dim];
+               std::memcpy(dst_ptr, src_ptr + begin * inp_strides_[dim], bytes);
+               return;
+            }
+
+            size_t src_stride = step * inp_strides_[dim];
+            size_t dst_stride = out_strides_[dim];
+
+            const uchar* s_ptr = src_ptr + begin * inp_strides_[dim];
+            uchar* d_ptr = dst_ptr;
+
+            for (int k = begin; k < end; k += step)
+            {
+                recursive_copy(dim + 1, s_ptr, d_ptr);
+                s_ptr += src_stride;
+                d_ptr += dst_stride;
+            }
+        }
+
+        bool is_fully_contiguous(int dim) const
+        {
+            size_t expected_step = es_;
+            for (int d = dims_ - 1; d >= dim; --d)
+            {
+                 if (inp_.step[d] != expected_step) return false;
+                 if (d > dim) {
+                     if (steps_[d] != 1) return false;
+                     if (ranges_[d].start != 0 || ranges_[d].end != inp_.size[d]) return false;
+                     expected_step *= inp_.size[d];
+                 }
+            }
+            return true;
+        }
+
+    private:
+        const Mat& inp_;
+        Mat& out_;
+        const std::vector<Range>& ranges_;
+        const std::vector<int>& steps_;
+        int dims_;
+        size_t es_;
+        std::vector<size_t> inp_strides_;
+        std::vector<size_t> out_strides_;
+    };
+
+    template <typename T>
+    void run_parallel(const Mat& inp, std::vector<Mat>& outputs)
+    {
+        int dimsNum = inp.dims;
+        for(size_t i=0; i<outputs.size(); ++i) {
+             std::vector<int> steps;
+             if (!sliceSteps.empty() && i < sliceSteps.size() && !sliceSteps[i].empty())
+                 steps = sliceSteps[i];
+             else
+                 steps.assign(dimsNum, 1);
+
+             ParallelSlice<T> body(inp, outputs[i], finalSliceRanges[i], steps);
+             int dim0_size = outputs[i].size[0];
+             parallel_for_(Range(0, dim0_size), body);
         }
     }
-
     void flip(Mat& output) // break if 1d tensor?
     {
         for (int i = 0; i < neg_step_dims.size(); ++i)