             tid = threadIdx.x%C;
         }
         int j;
-        v_t *tmp = ((v_t *)shared)+ncols*threadIdx.x;
-        m_t matval;
+        v_t tmp;
         
-        int c;
-        for (c=0; c<ncols; c++) {
-            zero<v_t>(tmp[c]);
-        }
+        zero<v_t>(tmp);
 
         for (j=0; j<rowlen[i]; j++) {
-            matval = val[cs+tid+j*C];
-            for (c=0; c<ncols; c++) {
-                tmp[c] = axpy<v_t,m_t>(tmp[c], rhs[rhs_lda*c+col[cs + tid + j*C]], matval);
-            }
-        }
-
-        for (c=0; c<ncols; c++) {
-            if (!do_shift) {
-                if (!do_scale) {
-                    if (do_axpy) {
-                        lhs[lhs_lda*c+i] = axpy<v_t,float>(lhs[lhs_lda*c+i],tmp[c],1.f);
-                    } else if (do_axpby) {
-                        lhs[i] = axpy<v_t,float>(scale<v_t>(lhs[i],beta),tmp[c],1.f);
-                    } else {
-                        lhs[i] = tmp[c];
-                    }
-                } else {
-                    if (do_axpy) {
-                        lhs[i] = axpy<v_t,float>(lhs[i],scale<v_t>(alpha,tmp[c]),1.f);
-                    } else if (do_axpby) {
-                        lhs[i] = axpy<v_t,float>(scale<v_t>(lhs[i],beta),scale<v_t>(alpha,tmp[c]),1.f);
-                    } else {
-                        lhs[i] = scale<v_t>(alpha,tmp[c]);
-                    }
+            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_scale) {
-                    if (do_axpy) {
-                        lhs[i] = axpy<v_t,float>(lhs[i],axpy<v_t,v_t>(tmp[c],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[c],rhs[i],scale2<v_t,float>(shift,-1.f)),1.f);
-                    } else {
-                        lhs[i] = axpy<v_t,v_t>(tmp[c],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[c],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[c],rhs[i],scale2<v_t,float>(shift,-1.f))),1.f);
-                    } else {
-                        lhs[i] = scale<v_t>(alpha,axpy<v_t,v_t>(tmp[c],rhs[i],scale2<v_t,float>(shift,-1.f)));
-                    }
-                }
-                
-            }
+                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)));
+                }
+            }
+            
         }
 
     }


#include "ghost/cu_complex.h"

#define SELL_CUDA_NBLOCKS (int)ceil(mat->nrowsPadded/(double)(SELL_CUDA_THREADSPERBLOCK/SELL(mat)->T)) 
//#define SELLT_STRIDE_ONE

extern __shared__ char shared[];

#include "ghost/log.h"
#include "ghost/error.h"
#include "ghost/util.h"
#include "ghost/math.h"

#include <cuComplex.h>
#include <stdio.h>
#include <cuda_runtime.h>
