 
 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

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


#include <curand.h>
#include <sys/types.h>
#include <unistd.h>
#include <complex.h>

#include "ghost/cu_complex.h"

