       <<<DIV_UP(H * W * N, BLOCK_SIZE1D), BLOCK_SIZE1D>>>(H * W * N, vOnes, 1);
   CHECK_LAST_CUDA_CALL();
   for (int z = 0; z < 4; ++z) {
-    CHECK_CUBLAS_CALL(cublasSetStream(handle, stream[z * 4 + 0]));
+    CHECK_CUBLAS_CALL(cublasSetStream(handle, stream[z]));
     CHECK_CUBLAS_CALL(gemv_gpu<T>(
         handle, CUBLAS_OP_T, H * W * N, 5 * D,
         scale, dQ_ptr(z, 0, 0, 0, 0, 0), 6 * D,

        1.0, dW_ptr(z, 0, 0, 0), 5 * D));
  }

  // dJ/dRy
  for (int z = 0; z < 4; ++z) {
    CHECK_CUBLAS_CALL(cublasSetStream(handle, stream[z * 4 + 2]));
    for (int y = 0; y < H; ++y) {
      for (int x = 0; x < W; ++x) {
        const int yp = (z == 0 || z == 1) ? y - 1 : y + 1;  // previous y
        if (yp >= 0 && yp < H) {
          CHECK_CUBLAS_CALL(gemm_gpu<T>(
              handle, CUBLAS_OP_T, CUBLAS_OP_N, D, 5 * D, N,
              scale, O_ptr(yp, x, 0, z, 0), 4 * D,
              dQ_ptr(z, y, x, 0, 0, 0), 6 * D,
              1.0, dRy_ptr(z, 0, 0, 0), 5 * D));
        }
      }
    }
  }

  // dJ/dRx
  for (int z = 0; z < 4; ++z) {
    CHECK_CUBLAS_CALL(cublasSetStream(handle, stream[z * 4 + 3]));
    for (int y = 0; y < H; ++y) {
      for (int x = 0; x < W; ++x) {
        const int xp = (z == 0 || z == 2) ? x - 1 : x + 1;  // previous x
        if (xp >= 0 && xp < W) {
          CHECK_CUBLAS_CALL(gemm_gpu<T>(
              handle, CUBLAS_OP_T, CUBLAS_OP_N, D, 5 * D, N,
              scale, O_ptr(y, xp, 0, z, 0), 4 * D,
              dQ_ptr(z, y, x, 0, 0, 0), 6 * D,
              1.0, dRx_ptr(z, 0, 0, 0), 5 * D));
        }
      }
    }
  }

  STREAMS_SYNCHRONIZE(4 * 4);


  // dJ/dW
  for (int z = 0; z < 4; ++z) {
    CHECK_CUBLAS_CALL(cublasSetStream(handle, stream[4 + z]));
    CHECK_CUBLAS_CALL(gemm_gpu<T>(
        handle, CUBLAS_OP_T, CUBLAS_OP_N, K, 5 * D, H * W * N,
        scale, I, K,
