     ghost_lidx K, ghost_lidx ldv, ghost_lidx ldw)
 {
     __shared__ oT blockStorage[BLOCKSIZE];
+    zero(blockStorage[threadIdx.x]);
     int tidx = threadIdx.x + blockIdx.x * blockDim.x;
     int n = tidx % N;
+    int m = (tidx / N) % M;
 
-    if (blockDim.x * gridDim.x / N == tidx / N) return;
+    if (blockDim.x * gridDim.x / N / M == tidx / N / M) return;
 
 
-    oT threadSum[M];
-    for (int chunk = tidx / N; chunk < nchunks; chunk += blockDim.x * gridDim.x / N) {
+    oT threadSum;
+    zero(threadSum);
+    for (int chunk = tidx / N / M; chunk < nchunks; chunk += blockDim.x * gridDim.x / N / M) {
         for (int c = 0; c < C; c++) {
             oT wsum;
             zero(wsum);
 
             for (int j = 0; j < chunkLen[chunk]; j++) {
                 ghost_gidx idx = chunkStart[chunk] + j * C + c;
-                wsum += (oT)Aval[idx] * (oT)wval[A->col[idx] * ldw + n];
+                wsum = axpy(wsum, (oT)Aval[idx], (oT)wval[col[idx] * ldw + n]);
             }
             if (chunk * C + c < K) {
-                for (int m = 0; m < M; m++) {
-                    threadSum[m] += wsum * (oT)vval[(chunk * C + c) * ldv + m];
-                }
+                threadSum = axpy(threadSum, wsum, (oT)vval[(chunk * C + c) * ldv + m]);
             }
         }
     }
 
-    for (int m = 0; m < M; m++) {
-        __syncthreads();
-        blockStorage[threadIdx.x] = threadSum[m];
-        __syncthreads();
+    __syncthreads();
+    blockStorage[threadIdx.x] = threadSum;
+    __syncthreads();
 
-        if (threadIdx.x < N) {
-            oT blockSum;
-            zero(blockSum);
-            for (int i = threadIdx.x; i < BLOCKSIZE; i += N) {
-                blockSum += blockStorage[i];
-            }
-            out[blockIdx.x * N * M + n * M + m] = blockSum;
+    if (threadIdx.x < M * N) {
+        oT blockSum;
+        zero(blockSum);
+        for (int i = threadIdx.x; i < BLOCKSIZE; i += M * N) {
+          blockSum = accu(blockSum, blockStorage[i]);
         }
+
+        out[blockIdx.x * N * M + n * M + m] = blockSum;
     }
 }
 
 template<typename T, typename oT>
-ghost_error typed_tsmtspmtsm(ghost_densemat *v, ghost_densemat *w, ghost_densemat *x,
+ghost_error typed_tsmtspmtsm(ghost_densemat *x, ghost_densemat *v, ghost_densemat *w,
     ghost_sparsemat *A, void *pAlpha, void *pBeta)
 {
-
-    T *wval = (T *)w->val;
-    T *vval = (T *)v->val;
-    oT *xval = (oT *)x->val;
+    ghost_error ret;
     oT alpha = *(oT *)pAlpha;
     oT beta = *(oT *)pBeta;
-    ghost_lidx C = A->traits.C;
+
     const ghost_lidx ldv = v->stride;
     const ghost_lidx ldw = w->stride;
     const ghost_lidx ldx = x->stride;
 
     ghost_lidx N = w->traits.ncols;
     ghost_lidx M = v->traits.ncols;
-    int nchunks = SPM_NCHUNKS(A);
+
+
+    int deviceUsed;
+    cudaGetDevice(&deviceUsed);
+    cudaDeviceProp prop;
+    cudaGetDeviceProperties(&prop, deviceUsed);
+    int numBlocks;
+
+    int const blockSize = 256;
+    CUDA_CALL(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
+                  &numBlocks, tsmtspmtsm_kernel<T, oT, blockSize>, blockSize, 0),
+        ret);
+    int blockCount = prop.multiProcessorCount * numBlocks;
 
     void *d_temp_storage = NULL;
     size_t required_temp_storage_bytes = M * N * blockCount * sizeof(oT);
-    GPU_ERROR(cudaMalloc(&d_temp_storage, required_temp_storage_bytes));
+    CUDA_CALL(cudaMalloc(&d_temp_storage, required_temp_storage_bytes), ret);
+    CUDA_CALL(cudaMemset(d_temp_storage, 0, required_temp_storage_bytes), ret);
 
+    tsmtspmtsm_kernel<T, oT, blockSize><<<blockCount, blockSize>>>((oT *)d_temp_storage,
+        (T *)w->cu_val, (T *)v->cu_val, (T *)A->cu_val, A->cu_chunkStart, A->cu_chunkLen, A->cu_col,
+        SPM_NCHUNKS(A), A->traits.C, M, N, SPM_NROWS(A), ldv, ldw);
 
-    tsmtspmtsm_kernel<T, oT>((oT *)d_temp_storage, wval, vval, (T *)A->cu_val, A->cu_chunkStart,
-        A->cu_chunkLen, A->cu_col, A->nchunks, A->traits.C, M, N, SPM_NROWS(A), ldv, ldw);
+    deviceReduce<oT><<<M * N / 256 + 1, 256>>>(
+        (oT *)d_temp_storage, (oT *)x->cu_val, alpha, beta, blockCount, M, N, ldx);
 
-    deviceReduce<oT>
-        <<<M * N / 256 + 1, 256>>>((oT *)d_temp_storage, xval, alpha, beta, blockCount, M, N, ldx);
-
-    GPU_ERROR(cudaFree(d_temp_storage));
-    return GHOST_SUCCESS;
+    CUDA_CALL(cudaFree(d_temp_storage), ret);
+    return ret;
 }
 }
 
 
-ghost_error tsmtspmtsm_var2_cuda(ghost_densemat *v, ghost_densemat *w, ghost_densemat *x,
+ghost_error tsmtspmtsm_var2_cuda(ghost_densemat *x, ghost_densemat *v, ghost_densemat *w,
     ghost_sparsemat *A, void *pAlpha, void *pBeta)
 {
     if (x->traits.datatype & GHOST_DT_REAL && w->traits.datatype & GHOST_DT_REAL) {

            return typed_tsmtspmtsm<float, double>(x, v, w, A, pAlpha, pBeta);
        }
    }
    return GHOST_SUCCESS;
}


namespace {


template<typename oT>
__global__ void deviceReduce(oT *blockResults, oT *result, oT alpha, oT beta, int blockCount,
    ghost_lidx M, ghost_lidx N, ghost_lidx ldx)
