             tmp = scale<v_t>(alpha,tmp);
         }
         if (do_axpby) {
-            lhs[lhs_lda*i+col] = axpy<v_t,float>(lhs[lhs_lda*i+col],tmp,1.f);
+            lhs[lhs_lda*i+col] = axpy<v_t,v_t>(tmp,lhs[lhs_lda*i+col],beta);
         } else {
             lhs[lhs_lda*i+col] = tmp;
-        }
-    }
-#ifdef LOCALDOT_ONTHEFLY 
-#endif
-
-}
-
-    template<typename m_t, typename v_t, bool do_axpby, bool do_scale, bool do_vshift, bool do_localdot>  
-__global__ void SELL_kernel_CU_tmpl(v_t *lhs, int lhs_lda, v_t *rhs, int rhs_lda, ghost_spmv_flags_t flags, int nrows, int nrowspadded, int ncols, ghost_lidx_t *rowlen, ghost_lidx_t *mcol, m_t *val, ghost_lidx_t *chunkstart, ghost_lidx_t *chunklen, int C, int T, v_t *shift, v_t alpha, v_t beta, v_t *localdot)
-{
-    UNUSED(T);
-    int i = threadIdx.x+blockIdx.x*blockDim.x;
-    int col = blockDim.y*blockIdx.y+threadIdx.y;
-
-    if (i<nrows) {
-        int cs, tid;
-        if (C == blockDim.x) {
-            cs = chunkstart[blockIdx.x];
-            tid = threadIdx.x;
-        } else {
-            cs = chunkstart[i/C];
-            tid = threadIdx.x%C;
-        }
-        int j;
-        v_t tmp;
-
-        zero<v_t>(tmp);
-
-        for (j=0; j<rowlen[i]; j++) {
-            tmp = axpy<v_t,m_t>(tmp, rhs[rhs_lda*col+mcol[cs + tid + j*C]], val[cs+tid+j*C]);
-        }
-
-        if (do_vshift) {
-            tmp = axpy<v_t,v_t>(tmp,rhs[rhs_lda*col+i],scale2<v_t,float>(shift[col*blockDim.y+threadIdx.y],-1.f));
-        }
-        if (do_scale) {
-            tmp = scale<v_t>(alpha,tmp);
-        }
-        if (do_axpby) {
-            lhs[lhs_lda*col+i] = axpy<v_t,float>(scale<v_t>(lhs[lhs_lda*col+i],beta),tmp,1.f);
-        } else {
-            lhs[lhs_lda*col+i] = tmp;
         }
     }
 #ifdef LOCALDOT_ONTHEFLY 

    dt2 *cu_shift = NULL;\
    dt2_host *localdot = NULL;\
    dt2 *shift, scale, beta;\
    GHOST_SPMV_PARSE_ARGS(flags,argp,scale,beta,shift,localdot,dt2_host,dt2);\
    if (flags & GHOST_SPMV_AXPY) {\
        dt2_host hbeta = 1.;\
        beta = *((dt2 *)&hbeta);\
    }\
    if (flags & GHOST_SPMV_SHIFT) {\
        ERROR_LOG("Currently not working!");\
    }\
    GHOST_CALL_RETURN(ghost_cu_malloc((void **)&cu_localdot,sizeof(dt2)*rhs->traits.ncols*3*SELL_CUDA_NBLOCKS));\
    size_t shiftsize = sizeof(dt2)*(flags & GHOST_SPMV_SHIFT?1:(flags & GHOST_SPMV_VSHIFT?rhs->traits.ncols:0));\
    GHOST_CALL_RETURN(ghost_cu_malloc((void **)&cu_shift,shiftsize));\
    GHOST_CALL_RETURN(ghost_cu_upload(cu_shift,shift,shiftsize));\
    struct cudaDeviceProp prop;\
    CUDA_CALL_RETURN(cudaGetDeviceProperties(&prop,cu_device));\
    GHOST_INSTR_START(spmv_cuda)\
    if (rhs->traits.ncols > 1) {\
        if (SELL(mat)->T > 1) {\
            WARNING_LOG("SELL-T kernel for multiple vectors nor implemented, falling back to SELL-1!");\
        }\
        int blockheight = PAD((int)ceil((double)SELL_CUDA_THREADSPERBLOCK/MIN(rhs->traits.ncols,MAX_COLS_PER_BLOCK)),SELL(mat)->chunkHeight);\
        if (blockheight*MIN(rhs->traits.ncols,MAX_COLS_PER_BLOCK) > 1024) {\
            WARNING_LOG("Too many threads! (FIXME)");\
        }\
        if (rhs->traits.ncols > MAX_COLS_PER_BLOCK) {\
            WARNING_LOG("Will have a loop over the vectors!");\
        }\
        if (rhs->traits.storage == GHOST_DENSEMAT_COLMAJOR) {\
            dim3 block(SELL_CUDA_THREADSPERBLOCK/MIN(MAX_COLS_PER_BLOCK,rhs->traits.ncols),MIN(MAX_COLS_PER_BLOCK,rhs->traits.ncols));\
            dim3 grid((int)ceil(mat->nrowsPadded/(double)block.x),(int)(ceil(rhs->traits.ncols/(double)MAX_COLS_PER_BLOCK)));\
            size_t reqSmem = 0;\
            if (flags & GHOST_SPMV_DOT) {\
                reqSmem = sizeof(dt2)*32*block.y;\
            }\
            if (prop.sharedMemPerBlock < reqSmem) {\
                WARNING_LOG("Not enough shared memory available! CUDA kernel will not execute!");\
            }\
            INFO_LOG("grid %dx%d block %dx%d shmem %zu",grid.x,grid.y,block.x,block.y,reqSmem);\
            SWITCH_BOOLS(SELL_kernel_CU_tmpl,dt1,dt2,grid,block,reqSmem)\
            /*SELL_kernel_CU_tmpl<dt1,dt2><<<(int)ceil(mat->nrowsPadded/(double)blockheight),block,reqSmem>>>((dt2 *)lhs->cu_val,lhs->traits.nrowspadded,(dt2 *)rhs->cu_val,rhs->traits.nrowspadded,flags,mat->nrows,mat->nrowsPadded,rhs->traits.ncols/block.y,SELL(mat)->cumat->rowLen,SELL(mat)->cumat->col,(dt1 *)SELL(mat)->cumat->val,SELL(mat)->cumat->chunkStart,SELL(mat)->cumat->chunkLen,SELL(mat)->chunkHeight,SELL(mat)->T,cu_shift,scale,beta,cu_localdot,flags&GHOST_SPMV_AXPY,flags&GHOST_SPMV_AXPBY,flags&GHOST_SPMV_SCALE,flags&GHOST_SPMV_SHIFT,flags&GHOST_SPMV_VSHIFT,flags&GHOST_SPMV_DOT);*/\
        } else {\
            INFO_LOG("Experimental row-major CUDA SELL-SpMMV");\
            dim3 block(MIN(MAX_COLS_PER_BLOCK,rhs->traits.ncols),SELL_CUDA_THREADSPERBLOCK/MIN(MAX_COLS_PER_BLOCK,rhs->traits.ncols));\
            dim3 grid((int)ceil(mat->nrowsPadded/(double)block.y),(int)(ceil(rhs->traits.ncols/(double)MAX_COLS_PER_BLOCK)));\
            INFO_LOG("grid %dx%d block %dx%d",grid.x,grid.y,block.x,block.y);\
            SWITCH_BOOLS(SELL_kernel_CU_rm_tmpl,dt1,dt2,grid,block)\
            /*SELL_kernel_CU_rm_tmpl<dt1,dt2><<<(int)ceil(mat->nrowsPadded/(double)32),newblock,reqSmem>>>((dt2 *)lhs->cu_val,lhs->traits.ncolspadded,(dt2 *)rhs->cu_val,rhs->traits.ncolspadded,flags,mat->nrows,mat->nrowsPadded,rhs->traits.ncols/newblock.x,SELL(mat)->cumat->rowLen,SELL(mat)->cumat->col,(dt1 *)SELL(mat)->cumat->val,SELL(mat)->cumat->chunkStart,SELL(mat)->cumat->chunkLen,SELL(mat)->chunkHeight,SELL(mat)->T,cu_shift,scale,beta,cu_localdot,flags&GHOST_SPMV_AXPY,flags&GHOST_SPMV_AXPBY,flags&GHOST_SPMV_SCALE,flags&GHOST_SPMV_SHIFT,flags&GHOST_SPMV_VSHIFT,flags&GHOST_SPMV_DOT);*/\
        }\
    } else {\
        int blockheight = PAD((int)ceil((double)SELL_CUDA_THREADSPERBLOCK/rhs->traits.ncols),SELL(mat)->chunkHeight);\
        if (blockheight*rhs->traits.ncols > 1024) {\
            WARNING_LOG("Too many threads! (FIXME)");\
        }\
        size_t reqSmem;\
        ghost_datatype_size(&reqSmem,lhs->traits.datatype);\
        reqSmem *= blockheight*rhs->traits.ncols;\
        if (prop.sharedMemPerBlock < reqSmem) {\
            WARNING_LOG("Not enough shared memory available! CUDA kernel will not execute!");\
        }\
        dim3 block(blockheight,rhs->traits.ncols);\
        SWITCH_BOOLS(SELL_kernel_CU_tmpl,dt1,dt2,SELL_CUDA_NBLOCKS,SELL_CUDA_THREADSPERBLOCK)\
    }\
    cudaDeviceSynchronize();\
    GHOST_INSTR_STOP(spmv_cuda)\


    if (i<nrows) {
        int cs, tid;
        if (C == blockDim.y) {
            cs = chunkstart[blockIdx.x];
            tid = threadIdx.y;
        } else {
