 {
     int idx = blockIdx.x*blockDim.x+threadIdx.x;
 
-    for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
-        if (rowmask[idx]) {
+    if (rowmask || colmask) {
+        for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
+            if (rowmask[idx]) {
+                ghost_idx_t v;
+                for (v=0; v<ncols; v++) {
+                    if (colmask[v]) {
+                        vec[v*nrowspadded+idx] = a;
+                    }
+                }
+            }
+        }
+    } else {
+        for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
             ghost_idx_t v;
             for (v=0; v<ncols; v++) {
-                if (colmask[v]) {
-                    vec[v*nrowspadded+idx] = a;
-                }
+                vec[v*nrowspadded+idx] = a;
             }
         }
     }

{
    int idx = blockIdx.x*blockDim.x+threadIdx.x;

    for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
        if (rowmask[idx]) {
            ghost_idx_t v;
            for (v=0; v<ncols; v++) {
                if (colmask[v]) {
                    v1[v*nrowspadded+idx] = axpby<T>(v2[v*nrowspadded+idx],v1[v*nrowspadded+idx],a[v],b[v]);
                }
            }
        }
    }
}

template<typename T>  

        {
            if (vec->traits.datatype & GHOST_DT_DOUBLE)
            {
                CUBLAS_CALL_GOTO(cublasZdotc(ghost_cublas_handle,vec->traits.nrows,
                            (const cuDoubleComplex *)v1,1,(const cuDoubleComplex *)v2,1,&((cuDoubleComplex *)res)[v]),err,ret);
            } 
            else 
            {
                CUBLAS_CALL_GOTO(cublasCdotc(ghost_cublas_handle,vec->traits.nrows,
                            (const cuFloatComplex *)v1,1,(const cuFloatComplex *)v2,1,&((cuFloatComplex *)res)[v]),err,ret);
            }
        }
        else
        {
            if (vec->traits.datatype & GHOST_DT_DOUBLE)
            {
                CUBLAS_CALL_GOTO(cublasDdot(ghost_cublas_handle,vec->traits.nrows,
                            (const double *)v1,1,(const double *)v2,1,&((double *)res)[v]),err,ret);
            } 
            else 
            {
                CUBLAS_CALL_GOTO(cublasSdot(ghost_cublas_handle,vec->traits.nrows,
                            (const float *)v1,1,(const float *)v2,1,&((float *)res)[v]),err,ret);
            }
        }
    }

    goto out;
err:
out:
    if (!ghost_bitmap_iscompact(vec->ldmask) || 
            !ghost_bitmap_iscompact(vec->trmask)) {
        vecclone->destroy(vecclone);
    }
    
    if (!ghost_bitmap_iscompact(vec2->ldmask) || 
            !ghost_bitmap_iscompact(vec2->trmask)) {
        vec2clone->destroy(vec2clone);
    }
    GHOST_INSTR_STOP(dot);

    return ret;
}

extern "C" ghost_error_t ghost_densemat_cm_cu_axpy(ghost_densemat_t *vec, ghost_densemat_t *vec2, void *a)
