 
 #include "ghost/cu_complex.h"
 
-#define SELL_CUDA_NBLOCKS (int)ceil(mat->nrowsPadded/(double)(SELL_CUDA_THREADSPERBLOCK/SELL(mat)->T)) 
+#define SELL_CUDA_NBLOCKS (int)ceil(mat->nrowsPadded/ceil((double)(SELL_CUDA_THREADSPERBLOCK/((double)SELL(mat)->T*(double)rhs->traits.ncols))))
 //#define SELLT_STRIDE_ONE
 
 extern __shared__ char shared[];

#include "ghost/log.h"
#include "ghost/error.h"
#include "ghost/util.h"
#include <cuComplex.h>
#include <stdio.h>
#include <cuda_runtime.h>

            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) {
            if (!do_scale) {
                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[i] = axpy<v_t,float>(scale<v_t>(lhs[i],beta),tmp,1.f);
                } else {
                    lhs[i] = tmp;
                }
            } else {
                if (do_axpy) {
                    lhs[i] = axpy<v_t,float>(lhs[i],scale<v_t>(alpha,tmp),1.f);
                } else if (do_axpby) {
                    lhs[i] = axpy<v_t,float>(scale<v_t>(lhs[i],beta),scale<v_t>(alpha,tmp),1.f);
                } else {
                    lhs[i] = scale<v_t>(alpha,tmp);
                }
            }
        } else {
            if (!do_scale) {
                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 {
                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)));
                }
            }
            
        }

    }
