 {
     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] = scale<T>(a,vec[v*nrowspadded+idx]);
+                    }
+                }
+            }
+        }
+    } 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] = scale<T>(a,vec[v*nrowspadded+idx]);
-                }
-            }
-        }
-    }
+                vec[v*nrowspadded+idx] = scale<T>(a,vec[v*nrowspadded+idx]);
+            }
+        }
+    }
+
 }
 
 template<typename T>  
 __global__ static void cu_vscale_kernel(T *vec, T *a, ghost_idx_t nrows, char *rowmask, ghost_idx_t ncols, char *colmask, ghost_idx_t nrowspadded)
 {
     int idx = blockIdx.x*blockDim.x+threadIdx.x;
-    int c;
-
-    for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
-        if (rowmask[idx]) {
+
+    if (rowmask || colmask) {
+        int c;
+        for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
+            if (rowmask[idx]) {
+                ghost_idx_t v;
+                for (c=0,v=0; v<ncols; v++) {
+                    if (colmask[v]) {
+                        vec[v*nrowspadded+idx] = scale<T>(a[c],vec[v*nrowspadded+idx]);
+                        c++;
+                    }
+                }
+            }
+        }
+    } else {
+        for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
             ghost_idx_t v;
-            for (c=0,v=0; v<ncols; v++) {
-                if (colmask[v]) {
-                    vec[v*nrowspadded+idx] = scale<T>(a[c],vec[v*nrowspadded+idx]);
-                    c++;
-                }
+            for (v=0; v<ncols; v++) {
+                vec[v*nrowspadded+idx] = scale<T>(a[v],vec[v*nrowspadded+idx]);
             }
         }
     }

                        (float *)vec->cu_val,1));
        }
    }*/
    
    GHOST_CALL_RETURN(ghost_cu_free(cucolfield));
    GHOST_CALL_RETURN(ghost_cu_free(curowfield));

    return GHOST_SUCCESS;
}

extern "C" ghost_error_t ghost_densemat_cm_cu_axpby(ghost_densemat_t *v1, ghost_densemat_t *v2, void *a, void *b)
{
    if (v1->traits.datatype != v2->traits.datatype)
    {
        ERROR_LOG("Cannot AXPY vectors with different data types");
        return GHOST_ERR_NOT_IMPLEMENTED;
    }
    char colfield[v1->traits.ncolsorig];
    char rowfield[v1->traits.nrowsorig];

    char *cucolfield, *curowfield;
    GHOST_CALL_RETURN(ghost_cu_malloc((void **)&cucolfield,v1->traits.ncolsorig));
    GHOST_CALL_RETURN(ghost_cu_malloc((void **)&curowfield,v1->traits.nrowsorig));

    ghost_densemat_mask2charfield(v1->trmask,v1->traits.ncolsorig,colfield);
    ghost_densemat_mask2charfield(v1->ldmask,v1->traits.nrowsorig,rowfield);

    GHOST_CALL_RETURN(ghost_cu_upload(cucolfield,colfield,v1->traits.ncolsorig));
    GHOST_CALL_RETURN(ghost_cu_upload(curowfield,rowfield,v1->traits.nrowsorig));
    if (v1->traits.datatype & GHOST_DT_COMPLEX)
    {
        if (v1->traits.datatype & GHOST_DT_DOUBLE)

                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        }
    }
    goto out;

err:

out:
    GHOST_CALL_RETURN(ghost_cu_free(cucolfield));
    GHOST_CALL_RETURN(ghost_cu_free(curowfield));
    GHOST_INSTR_STOP(scale);

    
    return ret;
}

extern "C" ghost_error_t ghost_densemat_cm_cu_vscale(ghost_densemat_t *vec, void *a)
{
    GHOST_INSTR_START(vscale);
    ghost_error_t ret = GHOST_SUCCESS;

    void *d_a;
    size_t sizeofdt;
    char colfield[vec->traits.ncolsorig];
    char rowfield[vec->traits.nrowsorig];

    char *cucolfield = NULL, *curowfield = NULL;
    
    ghost_datatype_size(&sizeofdt,vec->traits.datatype);
    GHOST_CALL_GOTO(ghost_cu_malloc(&d_a,vec->traits.ncols*sizeofdt),err,ret);
    ghost_cu_upload(d_a,a,vec->traits.ncols*sizeofdt);
    
    if (!ghost_bitmap_iscompact(vec->ldmask) || 
            !ghost_bitmap_iscompact(vec->trmask)) {
        
        WARNING_LOG("Potentially slow VSCALE operation because some rows or columns are masked out!");
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&cucolfield,vec->traits.ncolsorig),err,ret);
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&curowfield,vec->traits.nrowsorig),err,ret);

        ghost_densemat_mask2charfield(vec->trmask,vec->traits.ncolsorig,colfield);
        ghost_densemat_mask2charfield(vec->ldmask,vec->traits.nrowsorig,rowfield);

        GHOST_CALL_GOTO(ghost_cu_upload(cucolfield,colfield,vec->traits.ncolsorig),err,ret);
        GHOST_CALL_GOTO(ghost_cu_upload(curowfield,rowfield,vec->traits.nrowsorig),err,ret);
    }

    if (vec->traits.datatype & GHOST_DT_COMPLEX)
    {
