     z[0] = buf[0];
 }
 
-void softmax_backward(int n, const float* fx, const float* dEdf, float* dEdx) {
-  auto tb = SizeToBlockThreadPair(n);
-  float* gpu_ods;
-  float ods;
-  cudaMalloc((void **)&gpu_ods, sizeof(float));
-  ker_dotproduct<<<tb.first, tb.second>>>(n, fx, dEdf, gpu_ods);
-  cudaMemcpy(&ods, gpu_ods, sizeof(float), cudaMemcpyDeviceToHost);
-  cudaFree(gpu_ods);
-  accBinaryExprKernel<<<tb.first, tb.second>>>(n, fx, dEdf, dEdx, FSoftmaxBackward(-ods));
-}
-
 void VectorSum(int rows, int cols, const float * a, float* c, const bool isColWise)
 {
     assert(rows > 0 && cols > 0); // converting from size_t to int may cause overflow

    thrust::transform(X.begin(), X.end(), Y.begin(), Y.begin(), saxpy_functor(A));
}

void add_to(int n, const float* x, float *y)
{
    thrust::device_ptr<float> src_ptr = thrust::device_pointer_cast((float*)x);
    thrust::device_ptr<float> tgt_ptr = thrust::device_pointer_cast(y);
    // Y <- X + Y
    thrust::transform(src_ptr, src_ptr + n, tgt_ptr, tgt_ptr, thrust::plus<float>());
}

void set_to_value_of(int n, float* x0, float val)
{
    thrust::device_ptr<float> dev_ptr = thrust::device_pointer_cast(x0);

  ker_l2_norm_reducer<<<tb.first,tb.second>>>(n, x0, y, square, accumulate);
}

// A kernel to calculate the dot product between two arrays
__global__ void ker_dotproduct(int n, const float* x, const float* y, float* z) {
  __shared__ float buf[256];
