   ker_l2_norm_reducer<<<tb.first,tb.second>>>(n, x0, y, square, accumulate);
 }
 
-// adapted from NVIDIA example
-__global__ void ker_softmax(int n, const float *x0, float* res) {
-  __shared__ float buf[256];
-  for (int i = threadIdx.x; i < 256; i += blockDim.x) {
-    float me = __int_as_float(0xff800000);
-    for (int pos = i; pos < n; pos += 256) {
-      const float d = x0[pos];
-      me = d > me ? d : me;
-    }
-    buf[i] = me;
-  }
-  for (int stride = 128; stride > 0; stride >>= 1) {
-    __syncthreads();
-    for (int i = threadIdx.x; i < stride; i += blockDim.x)
-        buf[i] = buf[i] > buf[stride + i] ? buf[i] : buf[stride + i];
-  }
-  __syncthreads();
-  const float max_elem = buf[0];
-  for (int i = threadIdx.x; i < 256; i += blockDim.x) {
-    float sum = 0;
-    for (int pos = i; pos < n; pos += 256)
-      sum += expf(x0[pos] - max_elem);
-    buf[i] = sum;
-  }
-  for (int stride = 128; stride > 0; stride >>= 1) {
-    __syncthreads();
-    for (int i = threadIdx.x; i < stride; i += blockDim.x)
-        buf[i] += buf[stride + i];
-  }
-  __syncthreads();
-  float lz = log(buf[0]) + max_elem;
-  int i = threadIdx.x + blockIdx.x * blockDim.x;
-  while (i < n) {
-    res[i] = exp(x0[i] - lz);
-    i += gridDim.x * blockDim.x;
-  }
-}
-
-void softmax(int n, const float* x0, float* y) {
-  auto tb = SizeToBlockThreadPair(n);
-  ker_softmax<<<tb.first,tb.second>>>(n, x0, y);
-}
-
 // 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];

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

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

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