     onevec->fromScalar(onevec,one);
 
     one[1] = 0.;
-
-    ghost_lidx_t v;
-    for (v=0; v<vec->traits.ncols; v++)
-    {
-        if (vec->traits.datatype & GHOST_DT_COMPLEX)
-        {
-            if (vec->traits.datatype & GHOST_DT_DOUBLE)
-            {
-                CURAND_CALL_GOTO(curandGenerateUniformDouble(gen,
-                            &((double *)(vec->cu_val))[v*vec->traits.nrowspadded],
-                            vec->traits.nrows*2),err,ret);
-            } 
-            else 
-            {
-                CURAND_CALL_GOTO(curandGenerateUniform(gen,
-                            &((float *)(vec->cu_val))[v*vec->traits.nrowspadded],
-                            vec->traits.nrows*2),err,ret);
-            }
-        }
-        else
-        {
-            if (vec->traits.datatype & GHOST_DT_DOUBLE)
-            {
-                CURAND_CALL_GOTO(curandGenerateUniformDouble(gen,
-                            &((double *)(vec->cu_val))[v*vec->traits.nrowspadded],
-                            vec->traits.nrows),err,ret);
-            } 
-            else 
-            {
-                CURAND_CALL_GOTO(curandGenerateUniform(gen,
-                            &((float *)(vec->cu_val))[v*vec->traits.nrowspadded],
-                            vec->traits.nrows),err,ret);
-            }
-        }
-    }
-    if (vec->traits.datatype & GHOST_DT_DOUBLE) {
-        vec->axpby(vec,onevec,minusahalf,one);
-    } else {
-        vec->axpby(vec,onevec,fminusahalf,fone);
+    void *valptr;
+    ghost_densemat_t *compactvec;
+
+    if ((vec->traits.ncolsorig != vec->traits.ncols) || (vec->traits.flags & GHOST_DENSEMAT_SCATTERED)) {
+        INFO_LOG("Cloning (and compressing) vec before operation");
+        vec->clone(vec,&compactvec,vec->traits.nrows,0,vec->traits.ncols,0);
+    } else {
+        compactvec = vec;
+    }
+    ghost_densemat_cu_valptr(compactvec,&valptr);
+
+
+    if (vec->traits.datatype & GHOST_DT_COMPLEX)
+    {
+        if (vec->traits.datatype & GHOST_DT_DOUBLE)
+        {
+            CURAND_CALL_GOTO(curandGenerateUniformDouble(gen,
+                        (double *)valptr,
+                        compactvec->traits.nrowspadded*compactvec->traits.ncols*2),err,ret);
+        } 
+        else 
+        {
+            CURAND_CALL_GOTO(curandGenerateUniform(gen,
+                        (float *)valptr,
+                        compactvec->traits.nrowspadded*compactvec->traits.ncols*2),err,ret);
+        }
+    }
+    else
+    {
+        if (vec->traits.datatype & GHOST_DT_DOUBLE)
+        {
+            CURAND_CALL_GOTO(curandGenerateUniformDouble(gen,
+                        (double *)valptr,
+                        compactvec->traits.nrowspadded*compactvec->traits.ncols),err,ret);
+        } 
+        else 
+        {
+            CURAND_CALL_GOTO(curandGenerateUniform(gen,
+                        (float *)valptr,
+                        compactvec->traits.nrowspadded*compactvec->traits.ncols),err,ret);
+        }
+    }
+    if (compactvec->traits.datatype & GHOST_DT_DOUBLE) {
+        compactvec->axpby(compactvec,onevec,minusahalf,one);
+    } else {
+        compactvec->axpby(compactvec,onevec,fminusahalf,fone);
+    }
+    if (compactvec != vec) {
+        vec->fromVec(vec,compactvec,0,0);
+        compactvec->destroy(compactvec);
     }
     goto out;
 err:

extern "C" ghost_error_t ghost_densemat_cm_cu_vaxpy(ghost_densemat_t *v1, ghost_densemat_t *v2, void *a)
{
    GHOST_FUNC_ENTER(GHOST_FUNCTYPE_MATH);
    ghost_error_t ret = GHOST_SUCCESS;
    void *d_a;
    size_t sizeofdt;
    char colfield[v1->traits.ncolsorig];
    char rowfield[v1->traits.nrowsorig];
    char *cucolfield = NULL, *curowfield = NULL;
    ghost_idx_t c,v=0;
    ghost_datatype_size(&sizeofdt,v1->traits.datatype);
    
    GHOST_CALL_GOTO(ghost_cu_malloc(&d_a,v1->traits.ncols*sizeofdt),err,ret);
    GHOST_CALL_GOTO(ghost_cu_memset(d_a,0,v1->traits.ncolsorig*sizeofdt),err,ret);
    
    for (c=0; c<v1->traits.ncolsorig; c++) {
        if (ghost_bitmap_isset(v1->ldmask,c)) {
            GHOST_CALL_GOTO(ghost_cu_upload(&((char *)d_a)[c*sizeofdt],&((char *)a)[v*sizeofdt],sizeofdt),err,ret);
            v++;
        }
    }
    
    if (v1->traits.datatype != v2->traits.datatype)
    {
        ERROR_LOG("Cannot VAXPY vectors with different data types");
        return GHOST_ERR_NOT_IMPLEMENTED;
    }

    
    if ((v1->traits.flags & GHOST_DENSEMAT_VIEW) || (v2->traits.flags & GHOST_DENSEMAT_VIEW)) {
        WARNING_LOG("Potentially slow VAXPY operation because some rows or columns may be masked out!");
        
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&cucolfield,v1->traits.ncolsorig),err,ret);
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&curowfield,v1->traits.nrowsorig),err,ret);

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

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

    if (v1->traits.datatype & GHOST_DT_COMPLEX)
    {
        if (v1->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_vaxpy_kernel<cuDoubleComplex><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>((cuDoubleComplex *)v1->cu_val, (cuDoubleComplex *)v2->cu_val,(cuDoubleComplex *)d_a,v1->traits.nrowsorig,curowfield,v1->traits.ncolsorig,cucolfield,v1->traits.nrowspadded);
        } 
        else 
        {
            cu_vaxpy_kernel<cuFloatComplex><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>((cuFloatComplex *)v1->cu_val, (cuFloatComplex *)v2->cu_val,(cuFloatComplex *)d_a,v1->traits.nrowsorig,curowfield,v1->traits.ncolsorig,cucolfield,v1->traits.nrowspadded);
        }
    }
    else
    {
        if (v1->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_vaxpy_kernel<double><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>((double *)v1->cu_val, (double *)v2->cu_val,(double *)d_a,v1->traits.nrowsorig,curowfield,v1->traits.ncolsorig,cucolfield,v1->traits.nrowspadded);
        } 
        else 
        {
            cu_vaxpy_kernel<float><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>((float *)v1->cu_val, (float *)v2->cu_val,(float *)d_a,v1->traits.nrowsorig,curowfield,v1->traits.ncolsorig,cucolfield,v1->traits.nrowspadded);
        }
    }
    
    goto out;
err:
out:
    GHOST_CALL_RETURN(ghost_cu_free(cucolfield));
    GHOST_CALL_RETURN(ghost_cu_free(curowfield));
    GHOST_CALL_RETURN(ghost_cu_free(d_a));
    cudaDeviceSynchronize();
    GHOST_FUNC_EXIT(GHOST_FUNCTYPE_MATH);

    return ret;
}
    
extern "C" ghost_error_t ghost_densemat_cm_cu_vaxpby(ghost_densemat_t *v1, ghost_densemat_t *v2, void *a, void *b)
{
    GHOST_FUNC_ENTER(GHOST_FUNCTYPE_MATH);
    ghost_error_t ret = GHOST_SUCCESS;

    void *d_a;
    void *d_b;
    size_t sizeofdt;
    char colfield[v1->traits.ncolsorig];
    char rowfield[v1->traits.nrowsorig];
    ghost_idx_t c,v=0;

    char *cucolfield = NULL, *curowfield = NULL;
    
    ghost_datatype_size(&sizeofdt,v1->traits.datatype);
    
    GHOST_CALL_GOTO(ghost_cu_malloc(&d_a,v1->traits.ncols*sizeofdt),err,ret);
    GHOST_CALL_GOTO(ghost_cu_malloc(&d_b,v1->traits.ncols*sizeofdt),err,ret);
    GHOST_CALL_GOTO(ghost_cu_memset(d_a,0,v1->traits.ncolsorig*sizeofdt),err,ret);
    GHOST_CALL_GOTO(ghost_cu_memset(d_b,0,v1->traits.ncolsorig*sizeofdt),err,ret);
    
    for (c=0; c<v1->traits.ncolsorig; c++) {
        if (ghost_bitmap_isset(v1->ldmask,c)) {
            GHOST_CALL_GOTO(ghost_cu_upload(&((char *)d_a)[c*sizeofdt],&((char *)a)[v*sizeofdt],sizeofdt),err,ret);
            GHOST_CALL_GOTO(ghost_cu_upload(&((char *)d_b)[c*sizeofdt],&((char *)b)[v*sizeofdt],sizeofdt),err,ret);
            v++;
        }
    }
    
    
    if (v1->traits.datatype != v2->traits.datatype)
    {
        ERROR_LOG("Cannot VAXPBY vectors with different data types");
        return GHOST_ERR_NOT_IMPLEMENTED;
    }
    
    if ((v1->traits.flags & GHOST_DENSEMAT_VIEW) || (v2->traits.flags & GHOST_DENSEMAT_VIEW)) {
        
        WARNING_LOG("Potentially slow VAXPBY operation because some rows or columns may be masked out!");
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&cucolfield,v1->traits.ncolsorig),err,ret);
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&curowfield,v1->traits.nrowsorig),err,ret);

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

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

    if (v1->traits.datatype & GHOST_DT_COMPLEX)
    {
        if (v1->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_vaxpby_kernel<cuDoubleComplex><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                (cuDoubleComplex *)v1->cu_val, (cuDoubleComplex *)v2->cu_val,(cuDoubleComplex *)d_a,(cuDoubleComplex *)d_b,
                 v1->traits.nrowsorig,curowfield,v1->traits.ncolsorig,cucolfield,v1->traits.nrowspadded);
        } 
        else 
        {
            cu_vaxpby_kernel<cuFloatComplex><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                (cuFloatComplex *)v1->cu_val, (cuFloatComplex *)v2->cu_val,(cuFloatComplex *)d_a,(cuFloatComplex *)d_b,
                 v1->traits.nrowsorig,curowfield,v1->traits.ncolsorig,cucolfield,v1->traits.nrowspadded);
        }
    }
    else
    {
        if (v1->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_vaxpby_kernel<double><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                 (double *)v1->cu_val, (double *)v2->cu_val,(double *)d_a,(double *)d_b,
                 v1->traits.nrowsorig,curowfield,v1->traits.ncolsorig,cucolfield,v1->traits.nrowspadded);
        } 
        else 
        {
            cu_vaxpby_kernel<float><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                (float *)v1->cu_val, (float *)v2->cu_val,(float *)d_a,(float *)d_b,
                 v1->traits.nrowsorig,curowfield,v1->traits.ncolsorig,cucolfield,v1->traits.nrowspadded);
        }
    }
    goto out;
err:
out:
    GHOST_CALL_RETURN(ghost_cu_free(cucolfield));
    GHOST_CALL_RETURN(ghost_cu_free(curowfield));
    GHOST_CALL_RETURN(ghost_cu_free(d_a));
    GHOST_CALL_RETURN(ghost_cu_free(d_a));
    cudaDeviceSynchronize();
    GHOST_FUNC_EXIT(GHOST_FUNCTYPE_MATH);




template<typename T>  
__global__ static void cu_vaxpby_kernel(T *v1, T *v2, T *a, T *b, ghost_lidx_t nrows, ghost_lidx_t ncols, ghost_lidx_t ld1, ghost_lidx_t ld2)
{
    int idx = blockIdx.x*blockDim.x+threadIdx.x;

    for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
        ghost_lidx_t v;
        for (v=0; v<ncols; v++) {
            v1[v*ld1+idx] = axpby<T>(v2[v*ld2+idx],v1[v*ld1+idx],a[v],b[v]);
        }
    }
}

template<typename T>  
__global__ static void cu_axpby_kernel(T *v1, T *v2, T a, T b, ghost_lidx_t nrows, ghost_lidx_t ncols, ghost_lidx_t ld1, ghost_lidx_t ld2) 
{
    int idx = blockIdx.x*blockDim.x+threadIdx.x;
    for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
        ghost_lidx_t v;
        for (v=0; v<ncols; v++) {
            v1[v*ld1+idx] = axpby<T>(v2[v*ld2+idx],v1[v*ld1+idx],a,b);
        }
    }
}

template<typename T>  
__global__ static void cu_scale_kernel(T *vec, T a, ghost_lidx_t nrows, ghost_lidx_t ncols, ghost_lidx_t ld)
{
    int idx = blockIdx.x*blockDim.x+threadIdx.x;

    for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
        ghost_lidx_t v;
        for (v=0; v<ncols; v++) {
            vec[v*ld+idx] = scale<T>(a,vec[v*ld+idx]);
        }
    }

}

template<typename T>  
__global__ static void cu_vscale_kernel(T *vec, T *a, ghost_lidx_t nrows, ghost_lidx_t ncols, ghost_lidx_t ld)
{
    int idx = blockIdx.x*blockDim.x+threadIdx.x;

    for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
        ghost_lidx_t v;
        for (v=0; v<ncols; v++) {
            vec[v*ld+idx] = scale<T>(a[v],vec[v*ld+idx]);
        }
    }
}

template<typename T>  
__global__ static void cu_fromscalar_kernel(T *vec, T a, ghost_lidx_t nrows, ghost_lidx_t ncols, ghost_lidx_t ld)
{
    int idx = blockIdx.x*blockDim.x+threadIdx.x;

    for (;idx < nrows; idx+=gridDim.x*blockDim.x) {
        ghost_lidx_t v;
        for (v=0; v<ncols; v++) {
            vec[v*ld+idx] = a;
        }
    }
}
