     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);

    accBinaryExprKernel << <tb.first, tb.second >> >(col * row, dEdf, gpu_softmax, dEdx, FSubtract());
}

void logsoftmax_backward(int n, const float* fx, const float* dEdf, float* dEdx) 
{
    /*
    float off_diag_sum = 0;
    for (auto p : as_vector(dEdf))
    off_diag_sum += p;
    off_diag_sum *= -1;
    *dEdxi += (*fx).binaryExpr(*dEdf, FLogSoftmaxBackward(off_diag_sum));
    */
    thrust::device_ptr<float> dp = thrust::device_pointer_cast((float*)fx);
    thrust::device_ptr<float> de = thrust::device_pointer_cast((float*)dEdf);
    thrust::device_ptr<float> dr = thrust::device_pointer_cast(dEdx);
    thrust::device_vector<float> dtemp(n);
//    thrust::transform(dp, dp + n, de, dtemp.begin(), FWeightedError());
    float off_diag_sum  = - thrust::reduce(de, de + n);
    thrust::transform(dp, dp + n, de, dtemp.begin(), FLogSoftmaxBackward(off_diag_sum)); 
    thrust::transform(dtemp.begin(), dtemp.end(), dr, dr, thrust::plus<float>());
}

// adapted from NVIDIA example
__global__ void ker_pnlsoftmax(int n, int elem_idx, const float *x0, float* res, float* logz) {
  __shared__ float buf[256];

  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];
