 extern "C" void dd_SELL_kernel_CU(ghost_mat_t *mat, ghost_vec_t *lhs, ghost_vec_t *rhs, int options)
 {
     CHOOSE_KERNEL(double,double);
-/*
-    SELL_kernel_CU_tmpl<double,double,options&GHOST_SPMVM_AXPY,GHOST_SPMVM_AXPBY,GHOST_SPMVM_APPLY_SCALE,GHOST_SPMVM_APPLY_SHIFT,GHOST_SPMVM_COMPUTE_LOCAL_DOTPRODUCT> 
-        <<< SELL_CUDA_NBLOCKS,SELL_CUDA_THREADSPERBLOCK >>> ((double *)lhs->CU_val,(double *)rhs->CU_val,options,SELL(mat)->cumat->nrows,SELL(mat)->cumat->nrowsPadded,SELL(mat)->cumat->rowLen,SELL(mat)->cumat->col,(double *)SELL(mat)->cumat->val,SELL(mat)->cumat->chunkStart,SELL(mat)->cumat->chunkLen,SELL(mat)->chunkHeight,SELL(mat)->T,*(double *)mat->traits->shift,*(double *)mat->traits->scale,*(double *)mat->traits->beta,(double *)lhs->traits->localdot);*/
 }
 
 extern "C" void ds_SELL_kernel_CU(ghost_mat_t *mat, ghost_vec_t *lhs, ghost_vec_t *rhs, int options)

        if (do_shift) {
            if (do_scale) {
                if (do_axpy) {
                    lhs[i] = axpy<v_t,float>(lhs[i],scale<v_t>(alpha,axpy<v_t,v_t>(rhs[i],tmp,scale2<v_t,float>(shift,-1.f))),1.f);
                } else if (do_axpby) {
                    lhs[i] = axpy<v_t,float>(scale<v_t>(lhs[i],beta),scale<v_t>(alpha,axpy<v_t,v_t>(rhs[i],tmp,scale2<v_t,float>(shift,-1.f))),1.f);
                } else {
                    lhs[i] = scale<v_t>(alpha,axpy<v_t,v_t>(rhs[i],tmp,scale2<v_t,float>(shift,-1.f)));
                }
            } else {
                if (do_axpy) {
                    lhs[i] = axpy<v_t,float>(lhs[i],axpy<v_t,v_t>(rhs[i],tmp,scale2<v_t,float>(shift,-1.f)),1.f);
                } else if (do_axpby) {
                    lhs[i] = axpy<v_t,float>(scale<v_t>(lhs[i],beta),axpy<v_t,v_t>(rhs[i],tmp,scale2<v_t,float>(shift,-1.f)),1.f);
                } else {
                    lhs[i] = axpy<v_t,v_t>(rhs[i],tmp,scale2<v_t,float>(shift,-1.f));
                }
            }
        } else {

        if (do_shift) {
            if (do_scale) {
                if (do_axpy) {
                    lhs[i] = axpy<v_t,float>(lhs[i],scale<v_t>(alpha,axpy<v_t,v_t>(tmp,rhs[i],scale2<v_t,float>(shift,-1.f))),1.f);
                } else if (do_axpby) {
                    lhs[i] = axpy<v_t,float>(scale<v_t>(lhs[i],beta),scale<v_t>(alpha,axpy<v_t,v_t>(tmp,rhs[i],scale2<v_t,float>(shift,-1.f))),1.f);
                } else {
                    lhs[i] = scale<v_t>(alpha,axpy<v_t,v_t>(tmp,rhs[i],scale2<v_t,float>(shift,-1.f)));
                }
            } else {
                if (do_axpy) {
                    lhs[i] = axpy<v_t,float>(lhs[i],axpy<v_t,v_t>(tmp,rhs[i],scale2<v_t,float>(shift,-1.f)),1.f);
                } else if (do_axpby) {
                    lhs[i] = axpy<v_t,float>(scale<v_t>(lhs[i],beta),axpy<v_t,v_t>(tmp,rhs[i],scale2<v_t,float>(shift,-1.f)),1.f);
                } else {
                    lhs[i] = axpy<v_t,v_t>(tmp,rhs[i],scale2<v_t,float>(shift,-1.f));
                }
            }
        } else {
