From ce4b1162c8a6deda9f538f2b212a84adf4f9b48d Mon Sep 17 00:00:00 2001
From: Alban Bellot <alban.bellot@inria.fr>
Date: Thu, 21 Jul 2016 16:09:56 +0200
Subject: [PATCH] test convert (some bugs with symmetric)

---
 spm.c                  |   5 +-
 spm.h                  |   2 +
 z_spm_convert_to_csc.c | 117 +++++++++++++++++++++++++++++++++++++----
 z_spm_convert_to_csr.c |  91 ++++++++++++++++++++++++++++++--
 z_spm_convert_to_ijv.c |  54 +++++++++++++++++++
 5 files changed, 251 insertions(+), 18 deletions(-)

diff --git a/spm.c b/spm.c
index 9a4132e4..1ab5af08 100644
--- a/spm.c
+++ b/spm.c
@@ -98,8 +98,9 @@ spmInit( pastix_spm_t *spm )
     spm->gnnzexp = 0;
     spm->nnzexp  = 0;
 
-    spm->dof     = 1;
-    spm->dofs    = NULL;
+    spm->dof       = 1;
+    spm->dofs      = NULL;
+    spm->colmajor  = 1;
 
     spm->colptr   = NULL;
     spm->rowptr   = NULL;
diff --git a/spm.h b/spm.h
index 54144fa1..cb18f6d3 100644
--- a/spm.h
+++ b/spm.h
@@ -41,6 +41,8 @@ struct pastix_spm_s {
                                  if > 0, constant degree of freedom
                                  otherwise, irregular degree of freedom (refer to dofs)      */
     pastix_int_t *dofs;      /*< Number of degrees of freedom per unknown (NULL, if dof > 0) */
+    int           colmajor;  /*< If > 0, column major when PastixIJV with dofs
+                                 otherwise, row major                                        */
 
     pastix_int_t *colptr;    /*< List of indirections to rows for each vertex                */
     pastix_int_t *rowptr;    /*< List of edges for each vertex                               */
diff --git a/z_spm_convert_to_csc.c b/z_spm_convert_to_csc.c
index d07996f4..43f8367f 100644
--- a/z_spm_convert_to_csc.c
+++ b/z_spm_convert_to_csc.c
@@ -80,11 +80,46 @@ z_spmConvertIJV2CSC( pastix_spm_t *spm )
 
     /* Sort the rows and avals arrays by column */
     spm->rowptr = malloc(spm->nnz * sizeof(pastix_int_t));
+    pastix_int_t* node = calloc(spm->nnz+1,sizeof(pastix_int_t));
+    pastix_int_t* old_node = calloc(spm->nnz+1,sizeof(pastix_int_t));
+    pastix_int_t row, col, dofi, dofj;
+    pastix_int_t *dofs=spm->dofs;
+
+
+    for(i=0; i<spm->nnz; i++)
+    {
+        row = oldspm.rowptr[i]-baseval;
+        col = oldspm.colptr[i]-baseval;
+        dofi = ( spm->dof > 0 ) ? spm->dof : dofs[col+1] - dofs[col];
+        dofj = ( spm->dof > 0 ) ? spm->dof : dofs[row+1] - dofs[row];
+        old_node[i+1] = dofi * dofj;
+        node[spm->colptr[col]+1] += dofi * dofj;
+        spm->colptr[col]++;
+    }
+
+    for(i=0;i<spm->nnz;i++)
+    {
+        node[i+1]+=node[i];
+        old_node[i+1]+=old_node[i];
+    }
+
+    /* Restore the colptr indexes */
+    {
+        pastix_int_t tmp, tmp2;
+        tmp = spm->colptr[0];
+        spm->colptr[0] = 0;
+        for (j=0; j<spm->n; j++) {
+            tmp2 = spm->colptr[j+1];
+            spm->colptr[j+1] = tmp;
+            tmp = tmp2;
+        }
+    }
 
 #if defined(PRECISION_p)
     spm->values = NULL;
 #else
-    spm->values = malloc(spm->nnz * sizeof(pastix_complex64_t));
+    pastix_int_t ii, jj;
+    spm->values = malloc(spm->nnzexp * sizeof(pastix_complex64_t));
     navals = (pastix_complex64_t*)(spm->values);
     oavals = (pastix_complex64_t*)(oldspm.values);
 #endif
@@ -96,7 +131,18 @@ z_spmConvertIJV2CSC( pastix_spm_t *spm )
         spm->rowptr[ spm->colptr[i] ] = oldspm.rowptr[j];
 
 #if !defined(PRECISION_p)
-        navals[ spm->colptr[i] ] = oavals[j];
+        dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+        dofj = ( spm->dof > 0 ) ? spm->dof : dofs[oldspm.rowptr[j]-baseval+1] - dofs[oldspm.rowptr[j]-baseval];
+
+        for(ii=0; ii<dofi; ii++)
+        {
+            for(jj=0; jj<dofj; jj++)
+            {
+                navals[node[spm->colptr[i]]] = oavals[ old_node[j]];
+                old_node[j]++;
+                node[spm->colptr[i]]++;
+            }
+        }
 #endif
         (spm->colptr[i])++;
 
@@ -184,27 +230,39 @@ z_spmConvertCSR2CSC( pastix_spm_t *spm )
     {
         pastix_int_t       *row_csc;
         pastix_int_t       *col_csc;
+        pastix_int_t       *node;
+        pastix_int_t       *dofs;
+
 #if !defined(PRECISION_p)
         pastix_complex64_t *val_csc;
         pastix_complex64_t *valptr = (pastix_complex64_t*)(spm->values);
+        pastix_int_t ii, jj;
+        int    cpt = 0;
 #endif
-        pastix_int_t j, k, col, row, nnz, baseval;
+        pastix_int_t  i, j, k, col, row, nnz, baseval;
+        pastix_int_t dofi, dofj;
+
 
         baseval = spmFindBase( spm );
         nnz = spm->nnz;
 
         row_csc = malloc(nnz * sizeof(pastix_int_t));
         col_csc = calloc(spm->n+1,sizeof(pastix_int_t));
+        node    = calloc(spm->nnz+1,sizeof(pastix_int_t));
+        dofs    = spm->dofs;
 
         assert( row_csc );
         assert( col_csc );
+        assert( node );
+        assert( ( (dofs) && !(spm->dof > 0) ) ||
+                ( !(dofs) && (spm->dof > 0) ) ); // (dofs) xor (spm->dof > 0)
+
 
 #if !defined(PRECISION_p)
-        val_csc = malloc(nnz*sizeof(pastix_complex64_t));
+        val_csc = malloc(spm->nnzexp*sizeof(pastix_complex64_t));
         assert( val_csc );
 #endif
 
-
         /* Count the number of elements per column */
         for (j=0; j<nnz; j++) {
             col = spm->colptr[j] - baseval;
@@ -218,6 +276,35 @@ z_spmConvertCSR2CSC( pastix_spm_t *spm )
             col_csc[j+1] += col_csc[j];
         }
 
+        for(i=0; i<spm->n; i++)
+        {
+            for(k=spm->rowptr[i]; k<spm->rowptr[i+1]; k++)
+            {
+                j = spm->colptr[k-baseval] - baseval;
+                row =  ( spm->dof > 0 ) ? j        : dofs[j];
+                col =  ( spm->dof > 0)  ? i        : dofs[i];
+                dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+                dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                node[col_csc[j]+1] += dofi * dofj;
+                col_csc[j]++;
+            }
+        }
+
+        for(i=0;i <spm->nnz; i++)
+            node[i+1] += node[i];
+
+        /* Restore the colptr indexes */
+        {
+            pastix_int_t tmp, tmp2;
+            tmp = col_csc[0];
+            col_csc[0] = 0;
+            for (j=0; j<spm->n; j++) {
+                tmp2 = col_csc[j+1];
+                col_csc[j+1] = tmp;
+                tmp = tmp2;
+            }
+        }
+
         assert( (col_csc[spm->gN]) == nnz );
 
         for (row=0; row<spm->n; row++) {
@@ -228,27 +315,35 @@ z_spmConvertCSR2CSC( pastix_spm_t *spm )
                 col = spm->colptr[k] - baseval;
                 j = col_csc[col];
                 row_csc[j] = row + baseval;
-
 #if !defined(PRECISION_p)
-                val_csc[j] = valptr[k];
+                dofi = ( spm->dof > 0 ) ? spm->dof : dofs[col+1] - dofs[col];
+                dofj = ( spm->dof > 0 ) ? spm->dof : dofs[row+1] - dofs[row];
+                //printf("dof dof %d %d\n",dofi,dofj);
+                for(jj=0; jj < dofj ; jj++)
+                {
+                    for(ii=0; ii < dofi ; ii++)
+                    {
+                        val_csc[node[j]+ii*dofj+jj] = valptr[cpt];
+                        cpt++;
+                    }
+                }
 #endif
-                col_csc[col] ++;
+                col_csc[col]++;
             }
         }
 
         /* Restore the colptr indexes */
         {
             pastix_int_t tmp, tmp2;
-
             tmp = col_csc[0];
             col_csc[0] = baseval;
-            for (j=0; j<spm->n; j++) {
+            for (j=0; j<spm->n; j++)
+            {
                 tmp2 = col_csc[j+1];
                 col_csc[j+1] = tmp + baseval;
                 tmp = tmp2;
             }
         }
-
         spmExit( spm );
         spm->colptr = col_csc;
         spm->rowptr = row_csc;
diff --git a/z_spm_convert_to_csr.c b/z_spm_convert_to_csr.c
index ef7b8f35..792a6ae5 100644
--- a/z_spm_convert_to_csr.c
+++ b/z_spm_convert_to_csr.c
@@ -123,16 +123,57 @@ z_spmConvertIJV2CSR( pastix_spm_t *spm )
 #endif
     pastix_int_t       *spmptx, *otmp;
     pastix_int_t i, j, tmp, baseval, total;
-    pastix_spm_t oldspm;
+    pastix_int_t row, col, dofi, dofj;
 
-    /* Backup the input */
-    memcpy( &oldspm, spm, sizeof(pastix_spm_t) );
+    pastix_int_t *node     = calloc(spm->nnz+1,sizeof(pastix_int_t));
+    pastix_int_t *old_node = calloc(spm->nnz+1,sizeof(pastix_int_t));
+    pastix_int_t *dofs     =spm->dofs;
+
+    pastix_spm_t oldspm;
 
     /*
      * Check the baseval, we consider that arrays are sorted by columns or rows
      */
     baseval = spmFindBase( spm );
 
+
+#if !defined(PRECISION_p)
+    pastix_int_t ii, jj;
+    /* Transpose values in row major format */
+    if( spm->colmajor ) //A test
+    {
+        int k;
+        int cpt=0;
+        oavals = (pastix_complex64_t*)spm->values;
+        navals = malloc(sizeof(pastix_complex64_t)*spm->nnzexp);
+        pastix_int_t* dofs=spm->dofs;
+
+        for(k=0; k<spm->nnz; k++)
+        {
+            j = spm->rowptr[k]-baseval;
+            i = spm->colptr[k]-baseval;
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+            for(ii=0; ii<dofi; ii++)
+            {
+                for(jj=0; jj<dofj; jj++)
+                {
+                    navals[cpt+jj*dofi+ii]=*oavals;
+                    oavals++;
+                }
+            }
+            cpt += dofi*dofj;
+        }
+        spm->values=navals;
+    }
+#endif
+
+
+
+    /* Backup the input */
+    memcpy( &oldspm, spm, sizeof(pastix_spm_t) );
+
+
     /* Compute the new rowptr */
     spm->rowptr = (pastix_int_t *) calloc(spm->n+1,sizeof(pastix_int_t));
 
@@ -158,10 +199,39 @@ z_spmConvertIJV2CSR( pastix_spm_t *spm )
     /* Sort the colptr and avals arrays by rows */
     spm->colptr  = malloc(spm->nnz * sizeof(pastix_int_t));
 
+    for(i=0; i<spm->nnz; i++)
+    {
+        row = oldspm.rowptr[i]-baseval;
+        col = oldspm.colptr[i]-baseval;
+        dofi = ( spm->dof > 0 ) ? spm->dof : dofs[col+1] - dofs[col];
+        dofj = ( spm->dof > 0 ) ? spm->dof : dofs[row+1] - dofs[row];
+        old_node[i+1] = dofi * dofj;
+        node[spm->rowptr[row]+1] += dofi * dofj;
+        spm->rowptr[row]++;
+    }
+
+    for(i=0;i<spm->nnz;i++)
+    {
+        node[i+1]+=node[i];
+        old_node[i+1]+=old_node[i];
+    }
+
+    /* Restore the rowptr indexes */
+    {
+        pastix_int_t tmp, tmp2;
+        tmp = spm->rowptr[0];
+        spm->rowptr[0] = 0;
+        for (j=0; j<spm->n; j++) {
+            tmp2 = spm->rowptr[j+1];
+            spm->rowptr[j+1] = tmp;
+            tmp = tmp2;
+        }
+    }
+
 #if defined(PRECISION_p)
     spm->values = NULL;
 #else
-    spm->values = malloc(spm->nnz * sizeof(pastix_complex64_t));
+    spm->values = malloc(spm->nnzexp * sizeof(pastix_complex64_t));
     navals = (pastix_complex64_t*)(spm->values);
     oavals = (pastix_complex64_t*)(oldspm.values);
 #endif
@@ -173,7 +243,17 @@ z_spmConvertIJV2CSR( pastix_spm_t *spm )
         spm->colptr[ spm->rowptr[i] ] = oldspm.colptr[j];
 
 #if !defined(PRECISION_p)
-        navals[ spm->rowptr[i] ] = oavals[j];
+        dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+        dofj = ( spm->dof > 0 ) ? spm->dof : dofs[oldspm.colptr[j]-baseval+1] - dofs[oldspm.colptr[j]-baseval];
+        for(ii=0; ii<dofi; ii++)
+        {
+            for(jj=0; jj<dofj; jj++)
+            {
+                navals[node[spm->rowptr[i]]] = oavals[ old_node[j]];
+                old_node[j]++;
+                node[spm->rowptr[i]]++;
+            }
+        }
 #endif
         (spm->rowptr[i])++;
 
@@ -197,5 +277,6 @@ z_spmConvertIJV2CSR( pastix_spm_t *spm )
 
     spm->fmttype = PastixCSR;
 
+
     return PASTIX_SUCCESS;
 }
diff --git a/z_spm_convert_to_ijv.c b/z_spm_convert_to_ijv.c
index 68fa025e..d81ab2d4 100644
--- a/z_spm_convert_to_ijv.c
+++ b/z_spm_convert_to_ijv.c
@@ -62,6 +62,33 @@ z_spmConvertCSC2IJV( pastix_spm_t *spm )
         }
     }
 
+    /* Transpose values in row major format */
+    if( !spm->colmajor ) //A test
+    {
+        int k,ii,jj,dofi,dofj;
+        int cpt=0;
+        pastix_complex64_t* oavals = (pastix_complex64_t*)spm->values;
+        pastix_complex64_t* navals = malloc(sizeof(pastix_complex64_t)*spm->nnzexp);
+        pastix_int_t* dofs=spm->dofs;
+        for(k=0; k<spm->nnz; k++)
+        {
+            j = spm->rowptr[k]-baseval;
+            i = col_ijv[k]-baseval;
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+            for(ii=0; ii<dofi; ii++)
+            {
+                for(jj=0; jj<dofj; jj++)
+                {
+                    navals[cpt+jj*dofi+ii]=*oavals;
+                    oavals++;
+                }
+            }
+            cpt+=dofi*dofj;
+        }
+        spm->values=navals;
+    }
+
     memFree_null(spm->colptr);
     spm->colptr = col_ijv;
 
@@ -113,6 +140,33 @@ z_spmConvertCSR2IJV( pastix_spm_t *spm )
         }
     }
 
+    /* Transpose values in column major format */
+    if( spm->colmajor )
+    {
+        int k,ii,jj,dofi,dofj;
+        int cpt=0;
+        pastix_complex64_t* oavals = (pastix_complex64_t*)spm->values;
+        pastix_complex64_t* navals = malloc(sizeof(pastix_complex64_t)*spm->nnzexp);
+        pastix_int_t* dofs=spm->dofs;
+        for(k=0; k<spm->nnz; k++)
+        {
+            i = spm->colptr[k]-baseval;
+            j = row_ijv[k]-baseval;
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+            for(jj=0; jj<dofj; jj++)
+            {
+                for(ii=0; ii<dofi; ii++)
+                {
+                    navals[cpt+ii*dofj+jj]=*oavals;
+                    oavals++;
+                }
+            }
+            cpt+=dofi*dofj;
+        }
+        spm->values=navals;
+    }
+
     memFree_null(spm->rowptr);
     spm->rowptr = row_ijv;
 
-- 
GitLab