 }
 
     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)
+__global__ void SELL_kernel_CU_tmpl(v_t * const __restrict__ lhs, const int lhs_lda, const v_t * const __restrict__ rhs, const int rhs_lda, const ghost_spmv_flags_t flags, const int nrows, const int nrowspadded, const int ncols, const ghost_lidx_t * const __restrict__ rowlen, const ghost_lidx_t * const __restrict__ mcol, const m_t * const __restrict__ val, const ghost_lidx_t * const __restrict__ chunkstart, const ghost_lidx_t * const __restrict__ chunklen, const int C, const int T, const v_t * const __restrict__ shift, const v_t alpha, const v_t beta, v_t * const __restrict__ localdot)
 {
     UNUSED(T);
     int i = threadIdx.x+blockIdx.x*blockDim.x;

#ifdef LOCALDOT_ONTHEFLY 

    double3 dot;

    i = threadIdx.x+blockIdx.x*blockDim.y;
    col = blockDim.x*blockIdx.y+threadIdx.y;

    __syncthreads();
    if (i<nrows) {
        dot.x = lhs[lhs_lda*i+col] * lhs[lhs_lda*i+col];
        dot.y = rhs[rhs_lda*i+col] * lhs[lhs_lda*i+col];
        dot.z = rhs[rhs_lda*i+col] * rhs[rhs_lda*i+col];


    if (i<nrows) {
        int cs, tid;
        if (32 == blockDim.y) {
            cs = chunkstart[blockIdx.x];
            tid = threadIdx.y;
        } else {
            cs = chunkstart[i/32];
            tid = threadIdx.y%32;
        }
        int j;
        double tmp = 0.;

        for (j=0; j<rowlen[i]; j++) {
            tmp += rhs[rhs_lda*mcol[cs + tid + j*32]+col] * val[cs+tid+j*32];
        }

        lhs[lhs_lda*i+col] = lhs[lhs_lda*i+col]*beta + alpha*(tmp-rhs[rhs_lda*i+col]*shift[col]);
