 
 
     template<typename m_t, typename v_t>  
-__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 *col, 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, const bool do_axpy, const bool do_axpby, const bool do_scale, const bool do_shift, const bool do_vshift, const 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, const bool do_axpy, const bool do_axpby, const bool do_scale, const bool do_shift, const bool do_vshift, const bool do_localdot)
 {
     UNUSED(T);
     int i = threadIdx.x+blockIdx.x*blockDim.x;
-
-    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*threadIdx.y+col[cs + tid + j*C]], val[cs+tid+j*C]);
-        }
-
-        if (do_shift) {
-            tmp = axpy<v_t,v_t>(tmp,rhs[rhs_lda*threadIdx.y+i],scale2<v_t,float>(shift[0],-1.f));
-        }
-        if (do_vshift) {
-            tmp = axpy<v_t,v_t>(tmp,rhs[rhs_lda*threadIdx.y+i],scale2<v_t,float>(shift[threadIdx.y],-1.f));
-        }
-        if (do_scale) {
-            tmp = scale<v_t>(alpha,tmp);
-        }
-        if (do_axpy) {
-            lhs[lhs_lda*threadIdx.y+i] = axpy<v_t,float>(lhs[lhs_lda*threadIdx.y+i],tmp,1.f);
-        } else if (do_axpby) {
-            lhs[lhs_lda*threadIdx.y+i] = axpy<v_t,float>(scale<v_t>(lhs[lhs_lda*threadIdx.y+i],beta),tmp,1.f);
-        } else {
-            lhs[lhs_lda*threadIdx.y+i] = tmp;
-        }
+    int col;
+
+    for (col=0; col<ncols; col++) {
+        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*blockDim.y+threadIdx.y)+mcol[cs + tid + j*C]], val[cs+tid+j*C]);
+            }
+
+            if (do_shift) {
+                tmp = axpy<v_t,v_t>(tmp,rhs[rhs_lda*(col*blockDim.y+threadIdx.y)+i],scale2<v_t,float>(shift[0],-1.f));
+            }
+            if (do_vshift) {
+                tmp = axpy<v_t,v_t>(tmp,rhs[rhs_lda*(col*blockDim.y+threadIdx.y)+i],scale2<v_t,float>(shift[col*blockDim.y+threadIdx.y],-1.f));
+            }
+            if (do_scale) {
+                tmp = scale<v_t>(alpha,tmp);
+            }
+            if (do_axpy) {
+                lhs[lhs_lda*(col*blockDim.y+threadIdx.y)+i] = axpy<v_t,float>(lhs[lhs_lda*(col*blockDim.y+threadIdx.y)+i],tmp,1.f);
+            } else if (do_axpby) {
+                lhs[lhs_lda*(col*blockDim.y+threadIdx.y)+i] = axpy<v_t,float>(scale<v_t>(lhs[lhs_lda*(col*blockDim.y+threadIdx.y)+i],beta),tmp,1.f);
+            } else {
+                lhs[lhs_lda*(col*blockDim.y+threadIdx.y)+i] = tmp;
+            }
+        }
+#ifdef LOCALDOT_ONTHEFLY 
+        if (do_localdot) {
+            v_t dot1, dot2, dot3;
+            zero<v_t>(dot1);
+            zero<v_t>(dot2);
+            zero<v_t>(dot3);
+
+            if (i<nrows) {
+                dot1 = axpy<v_t>(dot1,lhs[lhs_lda*(col*blockDim.y+threadIdx.y)+i],lhs[lhs_lda*(col*blockDim.y+threadIdx.y)+i]);
+                dot2 = axpy<v_t>(dot2,rhs[rhs_lda*(col*blockDim.y+threadIdx.y)+i],lhs[lhs_lda*(col*blockDim.y+threadIdx.y)+i]);
+                dot3 = axpy<v_t>(dot3,rhs[rhs_lda*(col*blockDim.y+threadIdx.y)+i],rhs[rhs_lda*(col*blockDim.y+threadIdx.y)+i]);
+            } else {
+                zero<v_t>(dot1);
+                zero<v_t>(dot2);
+                zero<v_t>(dot3);
+            }
+
+            dot1 = blockReduceSum(dot1);
+            dot2 = blockReduceSum(dot2);
+            dot3 = blockReduceSum(dot3);
+
+            if (threadIdx.x==0) {
+                localdot[0*ncols*blockDim.y*gridDim.x + (col*blockDim.y+threadIdx.y)*gridDim.x + blockIdx.x] = dot1;
+                localdot[1*ncols*blockDim.y*gridDim.x + (col*blockDim.y+threadIdx.y)*gridDim.x + blockIdx.x] = dot2;
+                localdot[2*ncols*blockDim.y*gridDim.x + (col*blockDim.y+threadIdx.y)*gridDim.x + blockIdx.x] = dot3;
+            }
+        }
+#endif
     }
+
 }
 
     template<typename m_t, typename v_t>

        }
        int j;
        v_t tmp;
        
        zero<v_t>(tmp);

        for (j=0; j<rowlen[i]; j++) {

                        lhs[i] = tmp;
                    }
                }

            }
        }
    }
