                             const OpContext& ctx, const std::vector<TBlob>& inputs,
                             const std::vector<OpReqType>& req,
                             const std::vector<TBlob>& outputs) {
-
+  using namespace mshadow;
+  CHECK_EQ(inputs.size(), 5U);
+  const TBlob ograd = inputs[0];
+  const TBlob in_data = inputs[1];
+  const TBlob gamma = inputs[2];
+  const TBlob mean_data = inputs[3];
+  const TBlob std_data = inputs[4];
+  const TBlob data_grad = outputs[0];
+  const TBlob gamma_grad = outputs[1];
+  const TBlob beta_grad = outputs[2];
+
+  // Make sure the inputs are contiguous
+  CHECK_EQ(ograd.CheckContiguous(), true);
+  CHECK_EQ(in_data.CheckContiguous(), true);
+  CHECK_EQ(gamma.CheckContiguous(), true);
+  CHECK_EQ(mean_data.CheckContiguous(), true);
+  CHECK_EQ(std_data.CheckContiguous(), true);
+  int nbatch = in_data.shape_.ProdShape(0, in_data.ndim() - 1);
+  int nchannel = in_data.shape_[in_data.ndim() - 1];
+  int data_req = req[0];
+  int gamma_req = req[1];
+  int beta_req = req[2];
+  CHECK_NE(data_req, kWriteInplace);
+  CHECK_NE(gamma_req, kWriteInplace);
+  CHECK_NE(beta_req, kWriteInplace);
+  cudaStream_t stream = Stream<gpu>::GetStream(ctx.get_stream<gpu>());
+
+  // Calculate the gradient for gamma/beta
+  CHECK_EQ(gamma_grad.CheckContiguous(), true);
+  CHECK_EQ(beta_grad.CheckContiguous(), true);
+  const dim3 dimBlock(32, 4);
+  const int LOAD_UNROLL = 4;
+  if(gamma_req != kNullOp || beta_req != kNullOp) {
+    MSHADOW_REAL_TYPE_SWITCH(in_data.type_flag_, DType, {
+      DType* gamma_grad_ptr = (gamma_req != kNullOp) ? gamma.dptr<DType>() : nullptr;
+      DType* beta_grad_ptr = (beta_req != kNullOp) ? beta.dptr<DType>() : nullptr;
+      if(gamma_req == kAddTo && beta_req != kAddTo) {
+        LayerNormFusedBackwardKernel_GammaBeta<LOAD_UNROLL, true, false> <<<1, dimBlock, stream>>>
+          (nbatch, nchannel, in_data.dptr<DType>(), out_grad.dptr<DType>(), mean_data.dptr<DType>(),
+           std_data.dptr<DType>(), gamma_grad_ptr, beta_grad_ptr);
+      } else if(gamma_req != kAddTo && beta_req == kAddTo) {
+        LayerNormFusedBackwardKernel_GammaBeta<LOAD_UNROLL, false, true> <<<1, dimBlock, stream>>>
+          (nbatch, nchannel, in_data.dptr<DType>(), out_grad.dptr<DType>(), mean_data.dptr<DType>(),
+           std_data.dptr<DType>(), gamma_grad_ptr, beta_grad_ptr);
+      } else if(gamma_req == kAddTo && beta_req == kAddTo) {
+        LayerNormFusedBackwardKernel_GammaBeta<LOAD_UNROLL, true, true> <<<1, dimBlock, stream>>>
+          (nbatch, nchannel, in_data.dptr<DType>(), out_grad.dptr<DType>(), mean_data.dptr<DType>(),
+           std_data.dptr<DType>(), gamma_grad_ptr, beta_grad_ptr);
+      } else {
+        LayerNormFusedBackwardKernel_GammaBeta<LOAD_UNROLL, false, false> <<<1, dimBlock, stream>>>
+          (nbatch, nchannel, in_data.dptr<DType>(), out_grad.dptr<DType>(), mean_data.dptr<DType>(),
+            std_data.dptr<DType>(), gamma_grad_ptr, beta_grad_ptr);
+      }
+    });
+  }
+
+  // Calculate the gradient for data
+  CHECK_EQ(data_grad.CheckContiguous(), true);
+  int ngrid_x = (nbatch > kMaxGridDim) ? (nbatch + kBaseGridNum - 1) / kBaseGridNum : nbatch;
+  int ngrid_y = (nbatch > kMaxGridDim) ? kBaseGridNum : 1;
+  const dim3 dimGrid(ngrid_x, ngrid_y);
+  if(nchannel <= 32) {
+    dimBlock.y = 1;
+  } else {
+    dimBlock.y = 2;
+  }
+  if(data_req != kNullOp) {
+    MSHADOW_REAL_TYPE_SWITCH(in_data.type_flag_, DType, {
+      if(data_req == kAddTo) {
+        LayerNormFusedBackwardKernel_Data<LOAD_UNROLL, true> <<<dimGrid, dimBlock, stream>>>
+          (nbatch, nchannel, in_data.dptr<DType>(), out_grad.dptr<DType>(), mean_data.dptr<DType>(),
+           std_data.dptr<DType>(), gamma.dptr<DType>(), data_grad.dptr<DType>());
+      } else {
+        LayerNormFusedBackwardKernel_Data<LOAD_UNROLL, false> <<<dimGrid, dimBlock, stream>>>
+          (nbatch, nchannel, in_data.dptr<DType>(), out_grad.dptr<DType>(), mean_data.dptr<DType>(),
+           std_data.dptr<DType>(), gamma.dptr<DType>(), data_grad.dptr<DType>());
+      }
+    });
+  }
 }
 
 template<>

      for(int i = 0; i < LOAD_UNROLL; ++i) {
        DType ele_og = out_grad[b * nchannel + l + i];
        DType ele_x = in_data[b * nchannel + l + i];
        if(!no_gamma){
          local_gamma_grad[i] += (ele_x - local_mean[i]) / local_std[i] * ele_og;
        }
        if(!no_beta) {
          local_beta_grad[i] += ele_og;
        }
      }
    }
    // Write to destination
    for(int i = 0; i < LOAD_UNROLL; ++i) {
      if(!no_gamma) {
        gamma_grad[l + i] += local_gamma_grad[i];
      }
      if(!no_beta) {
        beta_grad[l + i] += local_beta_grad[i];
      }
    }
  }

    for(int l = tid; l < nchannel; l += nthread) {
      DType ele_out_grad = out_grad[bid * nchannel + l];
      DType ele_x = in_data[bid * nchannel + l];
      if(data_addto) {
        data_grad[bid * nchannel + l] +=
          ele_out_grad * gamma[l] * invstd_eps - sum_val0
                                               - (ele_x - mean) * invstd_eps * sum_val1;
      } else {
        data_grad[bid * nchannel + l] =
          ele_out_grad * gamma[l] * invstd_eps - sum_val0
                                               - (ele_x - mean) * invstd_eps * sum_val1;
      }
    }
  }
