         if (indims == 1) {
             threads = dim3(TPB, 1, 1);
             blocks = dim3(divup(odims[0], TPB));
-            diff2_1D<<<blocks, threads>>>(out, in, dim, oElem);
-        } else if (indims == 2) {
-            diff2_2D<<<blocks, threads>>>(out, in, dim, oElem, odims[0], ostrides[1], istrides[1]);
-        } else if (indims == 3) {
-            diff2_3D<<<blocks, threads>>>(out, in, dim, oElem, odims[0], odims[1],
-                ostrides[1], ostrides[2], istrides[1], istrides[2], blocksPerMatX);
-        } else if (indims == 4) {
-            diff1_4D<<<blocks, threads>>>(out, in, dim, oElem, odims[0], odims[1], odims[2],
+            diff_1D<T, true><<<blocks, threads>>>(out, in, dim, oElem);
+        } else if (indims < 4) {
+            diff_4D<T, true><<<blocks, threads>>>(out, in, dim, oElem, odims[0], odims[1], odims[2],
                                           ostrides[1], ostrides[2], ostrides[3],
                                           istrides[1], istrides[2], istrides[3],
                                           blocksPerMatX, blocksPerMatY);

        if (indims == 1) {
            threads = dim3(TPB, 1, 1);
            blocks = dim3(divup(odims[0], TPB));
            diff1_1D<<<blocks, threads>>>(out, in, dim, oElem);
        } else if (indims == 2) {
            diff1_2D<<<blocks, threads>>>(out, in, dim, oElem, odims[0], ostrides[1], istrides[1]);
        } else if (indims == 3) {
            diff1_3D<<<blocks, threads>>>(out, in, dim, oElem, odims[0], odims[1],
                    ostrides[1], ostrides[2], istrides[1], istrides[2], blocksPerMatX);
        } else if (indims == 4) {
            diff1_4D<<<blocks, threads>>>(out, in, dim, oElem, odims[0], odims[1], odims[2],
                                          ostrides[1], ostrides[2], ostrides[3],
                                          istrides[1], istrides[2], istrides[3],
                                          blocksPerMatX, blocksPerMatY);

                         (idy + 2 * isDim1) * istrides1 + (idx + 2 * isDim0);
        unsigned oMem = idw * ostrides3 + idz * ostrides2 + idy * ostrides1 + idx;

        iMem2 *= D;

        diff_this(out, in, oMem, iMem0, iMem1, iMem2);
    }

    ///////////////////////////////////////////////////////////////////////////
    // Wrapper functions
