 
 extern "C" ghost_error_t ghost_densemat_cm_cu_axpy(ghost_densemat_t *v1, ghost_densemat_t *v2, void *a)
 {
+    GHOST_FUNC_ENTER(GHOST_FUNCTYPE_MATH);
+    ghost_error_t ret = GHOST_SUCCESS; 
+    
+    if (v1->traits.datatype & GHOST_DT_COMPLEX)
+    {
+        if (v1->traits.datatype & GHOST_DT_DOUBLE)
+        {
+            const cuDoubleComplex one = make_cuDoubleComplex(1.,0);
+            ret =  ghost_densemat_cm_cu_axpby(v1,v2,a,(void *)&one);
+        } 
+        else 
+        {
+            const cuFloatComplex one = make_cuFloatComplex(1.,0.);
+            ret = ghost_densemat_cm_cu_axpby(v1,v2,a,(void *)&one);
+        }
+    }
+    else
+    {
+        if (v1->traits.datatype & GHOST_DT_DOUBLE)
+        {
+            const double one = 1.;
+            ret = ghost_densemat_cm_cu_axpby(v1,v2,a,(void *)&one);
+        } 
+        else 
+        {
+            const float one = 1.f;
+            ret = ghost_densemat_cm_cu_axpby(v1,v2,a,(void *)&one);
+        }
+    }
+    
+    GHOST_FUNC_EXIT(GHOST_FUNCTYPE_MATH);
+
+    return ret;
+}
+
+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");
+        ERROR_LOG("Cannot AXPBY vectors with different data types");
         return GHOST_ERR_NOT_IMPLEMENTED;
     }
     GHOST_FUNC_ENTER(GHOST_FUNCTYPE_MATH);
     ghost_error_t ret = GHOST_SUCCESS;
-    
-    char colfield1[v1->traits.ncolsorig];
-    char rowfield1[v1->traits.nrowsorig];
-    char *cucolfield1 = NULL, *curowfield1 = NULL;
-
-    char colfield2[v2->traits.ncolsorig];
-    char rowfield2[v2->traits.nrowsorig];
-    char *cucolfield2 = NULL, *curowfield2 = NULL;
-    
-    if ((v1->traits.flags & GHOST_DENSEMAT_VIEW) || (v2->traits.flags & GHOST_DENSEMAT_VIEW)) {
-        
-/*        if (!ghost_bitmap_isequal(vec->ldmask,vec2->ldmask) || !ghost_bitmap_isequal(vec->trmask,vec2->trmask)) {
-            ERROR_LOG("The masks have to be equal!");
-            ret = GHOST_ERR_INVALID_ARG;
-            goto err;
-        }*/
-       
-        WARNING_LOG("Potentially slow AXPY operation because some rows or columns may be masked out!");
-        
-        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&cucolfield1,v1->traits.ncolsorig),err,ret);
-        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&curowfield1,v1->traits.nrowsorig),err,ret);
-
-        ghost_densemat_mask2charfield(v1->trmask,v1->traits.ncolsorig,colfield1);
-        ghost_densemat_mask2charfield(v1->ldmask,v1->traits.nrowsorig,rowfield1);
-
-        GHOST_CALL_GOTO(ghost_cu_upload(cucolfield1,colfield1,v1->traits.ncolsorig),err,ret);
-        GHOST_CALL_GOTO(ghost_cu_upload(curowfield1,rowfield1,v1->traits.nrowsorig),err,ret);
-        
-        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&cucolfield2,v2->traits.ncolsorig),err,ret);
-        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&curowfield2,v2->traits.nrowsorig),err,ret);
-
-        ghost_densemat_mask2charfield(v2->trmask,v2->traits.ncolsorig,colfield2);
-        ghost_densemat_mask2charfield(v2->ldmask,v2->traits.nrowsorig,rowfield2);
-
-        GHOST_CALL_GOTO(ghost_cu_upload(cucolfield2,colfield2,v2->traits.ncolsorig),err,ret);
-        GHOST_CALL_GOTO(ghost_cu_upload(curowfield2,rowfield2,v2->traits.nrowsorig),err,ret);
-    }
-
-    
+    void *v1val, *v2val;
+    ghost_densemat_t *v1compact, *v2compact;
+    
+    if (v1->traits.flags & GHOST_DENSEMAT_SCATTERED) {
+        INFO_LOG("Cloning (and compressing) v1 before operation");
+        GHOST_CALL_GOTO(v1->clone(v1,&v1compact,v1->traits.nrows,0,v1->traits.ncols,0),err,ret);
+    } else {
+        v1compact = v1;
+    }
+    if (v2->traits.flags & GHOST_DENSEMAT_SCATTERED) {
+        INFO_LOG("Cloning (and compressing) v2 before operation");
+        GHOST_CALL_GOTO(v2->clone(v2,&v2compact,v2->traits.nrows,0,v2->traits.ncols,0),err,ret);
+    } else {
+        v2compact = v2;
+    }
+    GHOST_CALL_GOTO(ghost_densemat_cu_valptr(v1compact,&v1val),err,ret);
+    GHOST_CALL_GOTO(ghost_densemat_cu_valptr(v2compact,&v2val),err,ret);
+
     if (v1->traits.datatype & GHOST_DT_COMPLEX)
     {
         if (v1->traits.datatype & GHOST_DT_DOUBLE)
         {
-            const cuDoubleComplex one = make_cuDoubleComplex(1.,1.);
             cu_axpby_kernel<cuDoubleComplex><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>
-                ((cuDoubleComplex *)v1->cu_val, (cuDoubleComplex *)v2->cu_val,*((cuDoubleComplex *)a),one,
-                 v1->traits.nrowsorig,curowfield1,v1->traits.ncolsorig,cucolfield1,v1->traits.nrowspadded,
-                 v2->traits.nrowsorig,curowfield2,v2->traits.ncolsorig,cucolfield2,v2->traits.nrowspadded);
-        } 
-        else 
-        {
-            const cuFloatComplex one = make_cuFloatComplex(1.,1.);
+                ((cuDoubleComplex *)v1val, (cuDoubleComplex *)v2val,*((cuDoubleComplex *)a),*((cuDoubleComplex *)b),v1->traits.nrows,v1->traits.ncols,v1->traits.nrowspadded,v2->traits.nrowspadded);
+        } 
+        else 
+        {
             cu_axpby_kernel<cuFloatComplex><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>
-                ((cuFloatComplex *)v1->cu_val, (cuFloatComplex *)v2->cu_val,*((cuFloatComplex *)a),one,
-                 v1->traits.nrowsorig,curowfield1,v1->traits.ncolsorig,cucolfield1,v1->traits.nrowspadded,
-                 v2->traits.nrowsorig,curowfield2,v2->traits.ncolsorig,cucolfield2,v2->traits.nrowspadded);
+                ((cuFloatComplex *)v1val, (cuFloatComplex *)v2val,*((cuFloatComplex *)a),*((cuFloatComplex *)b),v1->traits.nrows,v1->traits.ncols,v1->traits.nrowspadded,v2->traits.nrowspadded);
+            
         }
     }
     else

        if (v1->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_axpby_kernel<double><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>
                ((double *)v1->cu_val, (double *)v2->cu_val,*((double *)a),(double)1.,
                 v1->traits.nrowsorig,curowfield1,v1->traits.ncolsorig,cucolfield1,v1->traits.nrowspadded,
                 v2->traits.nrowsorig,curowfield2,v2->traits.ncolsorig,cucolfield2,v2->traits.nrowspadded);
        } 
        else 
        {
            cu_axpby_kernel<float><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>
                ((float *)v1->cu_val, (float *)v2->cu_val,*((float *)a),(float)1.,
                 v1->traits.nrowsorig,curowfield1,v1->traits.ncolsorig,cucolfield1,v1->traits.nrowspadded,
                 v2->traits.nrowsorig,curowfield2,v2->traits.ncolsorig,cucolfield2,v2->traits.nrowspadded);
        }
    }

    
    /*cublasHandle_t ghost_cublas_handle;
    GHOST_CALL_RETURN(ghost_cu_cublas_handle(&ghost_cublas_handle)); 
    if (vec->traits.datatype & GHOST_DT_COMPLEX)
    {
        if (vec->traits.datatype & GHOST_DT_DOUBLE)
        {
            CUBLAS_CALL_RETURN(cublasZaxpy(ghost_cublas_handle,vec->traits.nrows,
                        (const cuDoubleComplex *)a,
                        (const cuDoubleComplex *)vec2->cu_val,1,
                        (cuDoubleComplex *)vec->cu_val,1));
        } 
        else 
        {
            CUBLAS_CALL_RETURN(cublasCaxpy(ghost_cublas_handle,vec->traits.nrows,
                        (const cuFloatComplex *)a,
                        (const cuFloatComplex *)vec2->cu_val,1,
                        (cuFloatComplex *)vec->cu_val,1));
        }
    }
    else
    {
        if (vec->traits.datatype & GHOST_DT_DOUBLE)
        {
            CUBLAS_CALL_RETURN(cublasDaxpy(ghost_cublas_handle,vec->traits.nrows,
                        (const double *)a,
                        (const double *)vec2->cu_val,1,
                        (double *)vec->cu_val,1));
        } 
        else 
        {
            CUBLAS_CALL_RETURN(cublasSaxpy(ghost_cublas_handle,vec->traits.nrows,
                        (const float *)a,
                        (const float *)vec2->cu_val,1,
                        (float *)vec->cu_val,1));
        }
    }*/
    goto out;
err:
out:
    
    GHOST_CALL_RETURN(ghost_cu_free(cucolfield1));
    GHOST_CALL_RETURN(ghost_cu_free(curowfield1));
    GHOST_CALL_RETURN(ghost_cu_free(cucolfield2));
    GHOST_CALL_RETURN(ghost_cu_free(curowfield2));
    cudaDeviceSynchronize();
    GHOST_FUNC_EXIT(GHOST_FUNCTYPE_MATH);

    return ret;
}

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 AXPBY vectors with different data types");
        return GHOST_ERR_NOT_IMPLEMENTED;
    }
    GHOST_FUNC_ENTER(GHOST_FUNCTYPE_MATH);
    ghost_error_t ret = GHOST_SUCCESS;

    char colfield1[v1->traits.ncolsorig];
    char rowfield1[v1->traits.nrowsorig];
    char *cucolfield1 = NULL, *curowfield1 = NULL;
    
    char colfield2[v2->traits.ncolsorig];
    char rowfield2[v2->traits.nrowsorig];
    char *cucolfield2 = NULL, *curowfield2 = NULL;
    
    if ((v1->traits.flags & GHOST_DENSEMAT_VIEW) || (v2->traits.flags & GHOST_DENSEMAT_VIEW)) {
        
/*        if (!ghost_bitmap_isequal(v1->ldmask,v2->ldmask) || !ghost_bitmap_isequal(v1->trmask,v2->trmask)) {
            ERROR_LOG("The masks have to be equal!");
            ret = GHOST_ERR_INVALID_ARG;
            goto err;
        }*/
       
        WARNING_LOG("Potentially slow AXPBY operation because some rows or columns may be masked out!");
        
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&cucolfield1,v1->traits.ncolsorig),err,ret);
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&curowfield1,v1->traits.nrowsorig),err,ret);

        ghost_densemat_mask2charfield(v1->trmask,v1->traits.ncolsorig,colfield1);
        ghost_densemat_mask2charfield(v1->ldmask,v1->traits.nrowsorig,rowfield1);

        GHOST_CALL_GOTO(ghost_cu_upload(cucolfield1,colfield1,v1->traits.ncolsorig),err,ret);
        GHOST_CALL_GOTO(ghost_cu_upload(curowfield1,rowfield1,v1->traits.nrowsorig),err,ret);
        
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&cucolfield2,v2->traits.ncolsorig),err,ret);
        GHOST_CALL_GOTO(ghost_cu_malloc((void **)&curowfield2,v2->traits.nrowsorig),err,ret);

        ghost_densemat_mask2charfield(v1->trmask,v2->traits.ncolsorig,colfield2);
        ghost_densemat_mask2charfield(v1->ldmask,v2->traits.nrowsorig,rowfield2);

        GHOST_CALL_GOTO(ghost_cu_upload(cucolfield2,colfield2,v2->traits.ncolsorig),err,ret);
        GHOST_CALL_GOTO(ghost_cu_upload(curowfield2,rowfield2,v2->traits.nrowsorig),err,ret);
    }

    if (v1->traits.datatype & GHOST_DT_COMPLEX)
    {
        if (v1->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_axpby_kernel<cuDoubleComplex><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>
                ((cuDoubleComplex *)v1->cu_val, (cuDoubleComplex *)v2->cu_val,*((cuDoubleComplex *)a),*((cuDoubleComplex *)b),
                 v1->traits.nrowsorig,curowfield1,v1->traits.ncolsorig,cucolfield1,v1->traits.nrowspadded,
                 v2->traits.nrowsorig,curowfield2,v2->traits.ncolsorig,cucolfield2,v2->traits.nrowspadded);
        } 
        else 
        {
            cu_axpby_kernel<cuFloatComplex><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>
                ((cuFloatComplex *)v1->cu_val, (cuFloatComplex *)v2->cu_val,*((cuFloatComplex *)a),*((cuFloatComplex *)b),
                 v1->traits.nrowsorig,curowfield1,v1->traits.ncolsorig,cucolfield1,v1->traits.nrowspadded,
                 v2->traits.nrowsorig,curowfield2,v2->traits.ncolsorig,cucolfield2,v2->traits.nrowspadded);
        }
    }
    else
    {
        if (v1->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_axpby_kernel<double><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>
                ((double *)v1->cu_val, (double *)v2->cu_val,*((double *)a),*((double *)b),
                 v1->traits.nrowsorig,curowfield1,v1->traits.ncolsorig,cucolfield1,v1->traits.nrowspadded,
                 v2->traits.nrowsorig,curowfield2,v2->traits.ncolsorig,cucolfield2,v2->traits.nrowspadded);
        } 
        else 
        {
            cu_axpby_kernel<float><<< (int)ceil((double)v1->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>
                ((float *)v1->cu_val, (float *)v2->cu_val,*((float *)a),*((float *)b),
                 v1->traits.nrowsorig,curowfield1,v1->traits.ncolsorig,cucolfield1,v1->traits.nrowspadded,
                 v2->traits.nrowsorig,curowfield2,v2->traits.ncolsorig,cucolfield2,v2->traits.nrowspadded);
        }
    }
    
    goto out;
err:
out:
    GHOST_CALL_RETURN(ghost_cu_free(cucolfield1));
    GHOST_CALL_RETURN(ghost_cu_free(cucolfield2));
    GHOST_CALL_RETURN(ghost_cu_free(curowfield1));
    GHOST_CALL_RETURN(ghost_cu_free(curowfield2));
    cudaDeviceSynchronize();
    GHOST_FUNC_EXIT(GHOST_FUNCTYPE_MATH);

    return ret;
}

extern "C" ghost_error_t ghost_densemat_cm_cu_scale(ghost_densemat_t *vec, void *a)
{
    GHOST_FUNC_ENTER(GHOST_FUNCTYPE_MATH);
    ghost_error_t ret = GHOST_SUCCESS;
    
    char colfield[vec->traits.ncolsorig];
    char rowfield[vec->traits.nrowsorig];

    char *cucolfield = NULL, *curowfield = NULL;
    
    if (vec->traits.flags & GHOST_DENSEMAT_VIEW) {
        WARNING_LOG("Potentially slow SCAL operation because some rows or columns may be 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)
    {
        if (vec->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_scale_kernel<cuDoubleComplex><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (cuDoubleComplex *)vec->cu_val, *(cuDoubleComplex *)a,
                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        } 
        else 
        {
            cu_scale_kernel<cuFloatComplex><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (cuFloatComplex *)vec->cu_val, *(cuFloatComplex *)a,
                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        }
    }
    else
    {
        if (vec->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_scale_kernel<double><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (double *)vec->cu_val, *(double *)a,
                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        } 
        else 
        {
            cu_scale_kernel<float><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (float *)vec->cu_val, *(float *)a,
                    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));
    cudaDeviceSynchronize();
    GHOST_FUNC_EXIT(GHOST_FUNCTYPE_MATH);

    
    return ret;
}

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

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

    char *cucolfield = NULL, *curowfield = NULL;
    
    ghost_datatype_size(&sizeofdt,vec->traits.datatype);
    GHOST_CALL_GOTO(ghost_cu_malloc(&d_a,vec->traits.ncolsorig*sizeofdt),err,ret);
    GHOST_CALL_GOTO(ghost_cu_memset(d_a,0,vec->traits.ncolsorig*sizeofdt),err,ret);
    
    for (c=0; c<vec->traits.ncolsorig; c++) {
        if (ghost_bitmap_isset(vec->ldmask,c)) {
            GHOST_CALL_GOTO(ghost_cu_upload(&((char *)d_a)[c*sizeofdt],&((char *)a)[v*sizeofdt],sizeofdt),err,ret);
            v++;
        }
    }
    
    if (vec->traits.flags & GHOST_DENSEMAT_VIEW) {
        
        WARNING_LOG("Potentially slow VSCALE operation because some rows or columns may be 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)
    {
        if (vec->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_vscale_kernel<cuDoubleComplex><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (cuDoubleComplex *)vec->cu_val, (cuDoubleComplex *)d_a,
                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        } 
        else 
        {
            cu_vscale_kernel<cuFloatComplex><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (cuFloatComplex *)vec->cu_val, (cuFloatComplex *)d_a,
                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        }
    }
    else
    {
        if (vec->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_vscale_kernel<double><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (double *)vec->cu_val, (double *)d_a,
                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        } 
        else 
        {
            cu_vscale_kernel<float><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (float *)vec->cu_val, (float *)d_a,
                    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_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_fromScalar(ghost_densemat_t *vec, void *a)
{
    ghost_error_t ret = GHOST_SUCCESS;
    
    char colfield[vec->traits.ncolsorig];
    char rowfield[vec->traits.nrowsorig];

    char *cucolfield = NULL, *curowfield = NULL;

    if (vec->traits.flags & GHOST_DENSEMAT_VIEW) {
        
        WARNING_LOG("Potentially slow fromScalar operation because some rows or columns may be 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);
    }
    
    ghost_densemat_cm_malloc(vec);
    if (vec->traits.datatype & GHOST_DT_COMPLEX)
    {
        if (vec->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_fromscalar_kernel<cuDoubleComplex><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (cuDoubleComplex *)vec->cu_val, *(cuDoubleComplex *)a,
                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        } 
        else 
        {
            cu_fromscalar_kernel<cuFloatComplex><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (cuFloatComplex *)vec->cu_val, *(cuFloatComplex *)a,
                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        }
    }
    else
    {
        if (vec->traits.datatype & GHOST_DT_DOUBLE)
        {
            cu_fromscalar_kernel<double><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (double *)vec->cu_val, *(double *)a,
                    vec->traits.nrowsorig,curowfield,vec->traits.ncolsorig,cucolfield,vec->traits.nrowspadded);
        } 
        else 
        {
            cu_fromscalar_kernel<float><<< (int)ceil((double)vec->traits.nrows/THREADSPERBLOCK),THREADSPERBLOCK >>>(
                    (float *)vec->cu_val, *(float *)a,
                    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));

    return ret;
}

extern "C" ghost_error_t ghost_densemat_cm_cu_fromRand(ghost_densemat_t *vec)
{
    ghost_error_t ret = GHOST_SUCCESS;
    if (vec->traits.flags & GHOST_DENSEMAT_VIEW) {
        ERROR_LOG("fromRand does currently not consider vector views!");
        return GHOST_ERR_NOT_IMPLEMENTED;
    }

    ghost_densemat_t *onevec;
    long pid = getpid();



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;
        }
    }
}
