     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>());
-}
-
+void softmax(int row, int col, const float* x0, float* y)
+{
+    cudaStream_t t_stream = cudaStreamDefault;
+
+    int N = col;
+    int M = row;
+    cudaEvent_t done = nullptr;
+    cudaEventCreate(&done);
+    _assignColumnwiseSoftmaxOf<float> << <N, MAX_THREADS_PER_BLOCK, 0, t_stream >> >(x0, y, N, M);
+
+    cudaEventRecord(done);
+
+    cudaEventSynchronize(done);
+
+    cudaEventDestroy(done);
+}
+
+/// see http://research.microsoft.com/pubs/226641/CNTKBook-20160121.pdf
+/*
+void softmax_backward(const float  kSCALAR_MINUSONE, int row, int col, const float *fx, const float *dEdf, float *dEdx, float *tmp_one_row, float * gpu_gradient)
+{
+    InnerProduct(row, col, dEdf, row, col, fx, 1, col, tmp_one_row, true);
+    auto tb = SizeToBlockThreadPair(col * row);
+//    ScaleAndAdd<float>(kSCALAR_MINUSONE, 1, col, tmp_one_row, row, col, dEdf, row, col, gpu_gradient);
+    accBinaryExprKernel << <tb.first, tb.second >> > (row * col, fx, gpu_gradient, dEdx, FProduct());
+}
+*/
+
+void softmax_backward(int n, const float* fx, const float* dEdf, float* dEdx, float* gpu_ods) {
+    auto tb = SizeToBlockThreadPair(n);
+    float ods;
+    ker_dotproduct << <tb.first, tb.second >> >(n, fx, dEdf, gpu_ods);
+    cudaMemcpy(&ods, gpu_ods, sizeof(float), cudaMemcpyDeviceToHost);
+    accBinaryExprKernel << <tb.first, tb.second >> >(n, fx, dEdf, dEdx, FSoftmaxBackward(-ods));
+}
+
 // adapted from NVIDIA example
 __global__ void ker_pnlsoftmax(int n, int elem_idx, const float *x0, float* res, float* logz) {
   __shared__ float buf[256];

    cudaEventDestroy(done);
}

void logsoftmax_forward(int row, int col, const float* x0, float* y) 
{
    cudaStream_t t_stream = cudaStreamDefault;


    z[0] = buf[0];
}

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
