         for (j=0; j<rowlen[i]; j++) {
             tmp = axpy<v_t,m_t>(tmp, rhs[col[i + j*nrowspadded]], val[i + j*nrowspadded]);
         }
-        
+
         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 (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/rhs->traits.ncols),SELL(mat)->chunkHeight);\
            if (blockheight*rhs->traits.ncols > 1024) {\
                WARNING_LOG("Too many threads! (FIXME)");\
            }\
            dim3 block(blockheight,rhs->traits.ncols);\
            SELL_kernel_CU_tmpl<dt1,dt2><<<(int)ceil(mat->nrowsPadded/(double)blockheight),block>>>((dt2 *)lhs->cu_val,lhs->traits.nrowspadded,(dt2 *)rhs->cu_val,rhs->traits.nrowspadded,flags,mat->nrows,mat->nrowsPadded,rhs->traits.ncols,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 {\
            if (SELL(mat)->chunkHeight == mat->nrowsPadded) {\
                if (SELL(mat)->T > 1) {\

#include <stdio.h>
#include <cuda_runtime.h>
#include <cuda.h>
#include <complex.h>

#include "ghost/cu_complex.h"

#define MAX_COLS_PER_BLOCK 2
#define SELL_CUDA_NBLOCKS (int)ceil(mat->nrowsPadded/ceil((double)(SELL_CUDA_THREADSPERBLOCK/((double)SELL(mat)->T*(double)(MIN(rhs->traits.ncols,MAX_COLS_PER_BLOCK))))))
//#define SELLT_STRIDE_ONE
#define LOCALDOT_ONTHEFLY

extern __shared__ char shared[];

#define CALL(func,dt1,dt2,dt2_host,b1,b2,b3,b4,b5,...){\
    dt2 shift, scale, beta;\
    GHOST_SPMV_PARSE_ARGS(flags,argp,scale,beta,shift,localdot,dt2_host,dt2);\
    func<dt1,dt2,b1,b2,b3,b4,b5><<<__VA_ARGS__>>>((dt2 *)lhs->cu_val,(dt2 *)rhs->cu_val,flags,mat->nrows,mat->nrowsPadded,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,shift,scale,beta,(dt2 *)cu_localdot);\
}\

#define SWITCH_BOOLS(func,dt1,dt2,dt2_host,...)\
    if (flags & GHOST_SPMV_AXPY) {\
        if (flags & GHOST_SPMV_AXPBY) {\
            if (flags & GHOST_SPMV_SCALE) {\
                if (flags & GHOST_SPMV_SHIFT) {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,true,true,true,true,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,true,true,true,true,false,__VA_ARGS__)\
                    }\
                } else {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,true,true,true,false,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,true,true,true,false,false,__VA_ARGS__)\
                    }\
                }\
            } else {\
                if (flags & GHOST_SPMV_SHIFT) {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,true,true,false,true,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,true,true,false,true,false,__VA_ARGS__)\
                    }\
                } else {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,true,true,false,false,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,true,true,false,false,false,__VA_ARGS__)\
                    }\
                }\
            }\
        } else {\
            if (flags & GHOST_SPMV_SCALE) {\
                if (flags & GHOST_SPMV_SHIFT) {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,true,false,true,true,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,true,false,true,true,false,__VA_ARGS__)\
                    }\
                } else {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,true,false,true,false,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,true,false,true,false,false,__VA_ARGS__)\
                    }\
                }\
            } else {\
                if (flags & GHOST_SPMV_SHIFT) {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,true,false,false,true,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,true,false,false,true,false,__VA_ARGS__)\
                    }\
                } else {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,true,false,false,false,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,true,false,false,false,false,__VA_ARGS__)\
                    }\
                }\
            }\
        }\
    } else {\
        if (flags & GHOST_SPMV_AXPBY) {\
            if (flags & GHOST_SPMV_SCALE) {\
                if (flags & GHOST_SPMV_SHIFT) {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,false,true,true,true,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,false,true,true,true,false,__VA_ARGS__)\
                    }\
                } else {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,false,true,true,false,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,false,true,true,false,false,__VA_ARGS__)\
                    }\
                }\
            } else {\
                if (flags & GHOST_SPMV_SHIFT) {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,false,true,false,true,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,false,true,false,true,false,__VA_ARGS__)\
                    }\
                } else {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,false,true,false,false,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,false,true,false,false,false,__VA_ARGS__)\
                    }\
                }\
            }\
        } else {\
            if (flags & GHOST_SPMV_SCALE) {\
                if (flags & GHOST_SPMV_SHIFT) {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,false,false,true,true,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,false,false,true,true,false,__VA_ARGS__)\
                    }\
                } else {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,false,false,true,false,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,false,false,true,false,false,__VA_ARGS__)\
                    }\
                }\
            } else {\
                if (flags & GHOST_SPMV_SHIFT) {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,false,false,false,true,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,false,false,false,true,false,__VA_ARGS__)\
                    }\
                } else {\
                    if (flags & GHOST_SPMV_DOT) {\
                        CALL(func,dt1,dt2,dt2_host,false,false,false,false,true,__VA_ARGS__)\
                    } else {\
                        CALL(func,dt1,dt2,dt2_host,false,false,false,false,false,__VA_ARGS__)\
                    }\
                }\
            }\
        }\
    }\


#ifdef LOCALDOT_ONTHEFLY
#define PROCESS_LOCALDOT(dt2_host)\
        GHOST_INSTR_START(spmv_cuda_dot_reduction)\
        int block, col;\
        INFO_LOG("Experimental local dot product with final reduction over %d blocks!",SELL_CUDA_NBLOCKS);\
        dt2_host *localdot_blocks;\
        GHOST_CALL_RETURN(ghost_malloc((void **)&localdot_blocks,sizeof(dt2_host)*rhs->traits.ncols*3*SELL_CUDA_NBLOCKS));\
        GHOST_CALL_RETURN(ghost_cu_download(localdot_blocks,cu_localdot,sizeof(dt2_host)*rhs->traits.ncols*3*SELL_CUDA_NBLOCKS));\
        _Pragma("omp parallel for private(block)")\
        for (col=0; col<rhs->traits.ncols; col++) {\
            localdot[col                      ] = 0;\
            localdot[col + 1*rhs->traits.ncols] = 0;\
            localdot[col + 2*rhs->traits.ncols] = 0;\
            for (block=0; block<SELL_CUDA_NBLOCKS; block++) {\
                localdot[col                      ] += localdot_blocks[                                        col*SELL_CUDA_NBLOCKS + block];\
                localdot[col + 1*rhs->traits.ncols] += localdot_blocks[1*SELL_CUDA_NBLOCKS*rhs->traits.ncols + col*SELL_CUDA_NBLOCKS + block];\
                localdot[col + 2*rhs->traits.ncols] += localdot_blocks[2*SELL_CUDA_NBLOCKS*rhs->traits.ncols + col*SELL_CUDA_NBLOCKS + block];\
            }\
        }\
        free(localdot_blocks);\
        GHOST_INSTR_STOP(spmv_cuda_dot_reduction)
#else
#define PROCESS_LOCALDOT(dt2_host)\
        GHOST_INSTR_START(spmv_cuda_dot)\
          INFO_LOG("Not doing the local dot product on-the-fly!");\
          memset(localdot,0,rhs->traits.ncols*3*sizeof(dt2_host));\
          lhs->dot(lhs,&localdot[0],lhs);\
          lhs->dot(lhs,&localdot[rhs->traits.ncols],rhs);\
          lhs->dot(rhs,&localdot[2*rhs->traits.ncols],rhs);\
          GHOST_INSTR_STOP(spmv_cuda_dot)
#endif

#define CHOOSE_KERNEL(dt1,dt2,dt2_host) {\
    ghost_error_t ret = GHOST_SUCCESS;\
    int cu_device;\
    GHOST_CALL_RETURN(ghost_cu_device(&cu_device));\
    dt2 *cu_localdot = NULL;\
    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);\
    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));\
    if ((SELL(mat)->T > 128) || (SELL(mat)->T == 0) || (SELL(mat)->T & (SELL(mat)->T-1)))\
    WARNING_LOG("Invalid T: %d (must be power of two and T <= 128)",SELL(mat)->T);\
    struct cudaDeviceProp prop;\
    CUDA_CALL_RETURN(cudaGetDeviceProperties(&prop,cu_device));\
    GHOST_INSTR_START(spmv_cuda)\
    if (rhs->traits.flags & (GHOST_DENSEMAT_VIEW | GHOST_DENSEMAT_SCATTERED)) {\
        if (!ghost_bitmap_iscompact(rhs->ldmask)) {\
