Multi-dof fin de stage

spm.c

@@ -813,3 +813,5 @@ spmCheckAxb( int nrhs,
         return ptrfunc[id](nrhs, spm, x0, ldx0, b, ldb, x, ldx );
     }
 }
+
+

spm.h

@@ -41,7 +41,7 @@ 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
+    int           colmajor;  /*< If > 0, column major with dofs
                                  otherwise, row major                                        */
 
     pastix_int_t *colptr;    /*< List of indirections to rows for each vertex                */

spm_expand.c

+/**
+ *
+ * @file spm_expand.c
+ *
+ *  PaStiX spm routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @author Alban Bellot
+ * @date 2015-01-01
+ *
+ **/
+#include "common.h"
+#include "spm.h"
+#include "z_spm.h"
+
+void
+spmExpandCSC(pastix_spm_t* spm)
+{
+    switch(spm->flttype)
+    {
+    case PastixFloat:
+        s_extandCSC(spm);
+        break;
+    case PastixComplex32:
+        c_extandCSC(spm);
+        break;
+    case PastixComplex64:
+        z_extandCSC(spm);
+        break;
+    case PastixDouble:
+    default:
+        d_extandCSC(spm);
+        break;
+    }
+}
+
+void print_tab_int(pastix_int_t* tab, pastix_int_t size)
+{
+    int i;
+    for(i=0;i<size;i++)
+        printf("%d ",tab[i]);
+    printf("\n");
+}
+
+void print_tab_cpx(double* tab, pastix_int_t size)
+{
+    int i;
+    for(i=0;i<size;i++)
+        printf("%f ",tab[i]);
+    printf("\n");
+}
+
+
+void init_rand()
+{
+  srand(time(NULL));
+}
+
+int randminmax(pastix_int_t min, pastix_int_t max)
+{
+  return min + rand()%(max-min+1);
+}
+
+pastix_int_t* computeCorrespondIndex(const pastix_int_t *dofs, pastix_int_t n)
+{
+    pastix_int_t* dofs_coef=malloc(sizeof(pastix_int_t)*(n+1));
+    dofs_coef[0]=0;
+    int i;
+    for(i=1;i<n+1;i++)
+    {
+        dofs_coef[i]=dofs_coef[i-1]+dofs[i-1];
+    }
+    return dofs_coef;
+}
+
+void genDof(pastix_int_t step,pastix_int_t* dofs,pastix_int_t size)
+{
+  int i;
+  for(i=0;i<size;i++)
+      dofs[i]=i%step+1;
+}
+
+void genRandDof(pastix_int_t min,pastix_int_t max,pastix_int_t* dofs,pastix_int_t size)
+{
+  int i;
+  for(i=0;i<size;i++)
+    dofs[i]=randminmax(min,max);
+}
+
+pastix_int_t* computeEltPerCol(const pastix_spm_t *spm)
+{
+  pastix_int_t col,row,baseval;
+
+  pastix_int_t* coef=malloc(sizeof(pastix_int_t)*spm->n);
+  int i;
+  for(i=0;i<spm->n;i++)
+      coef[i]=0;
+  baseval=spmFindBase( spm );
+  row=0;
+  for( col=0; col < spm->n; col++ )
+    {
+      for( row=spm->colptr[col]-baseval; row<spm->colptr[col+1]-baseval; row++)
+	{
+            coef[col] += spm->dofs[col] * spm->dofs[spm->rowptr[row]-baseval];
+	}
+    }
+  return coef;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup spm_internal
+ *
+ * dofCst - Expand a CSC matrix without dofs into a CSC matrix with constant dofs
+ *
+ *
+ *******************************************************************************
+ *
+ * @param[in] spm
+ *           The spm to expand
+ *
+ *******************************************************************************/
+void
+dofCst(pastix_spm_t* spm,pastix_int_t dof)
+{
+    int i;
+    pastix_int_t* dofs = malloc((spm->n+1)*sizeof(pastix_int_t));
+
+    spm->dof      = dof;
+    spm->gNexp    = spm->n*dof;
+    spm->nexp     = spm->n*dof;
+    spm->nnzexp   = spm->n*dof*dof;
+    spm->gnnzexp  = spm->n*dof*dof;
+    spm->dofs     = dofs;
+
+    for(i=0;i<spm->n+1;i++)
+        dofs[i]=dof*i;
+
+    spmExpandCSC(spm);
+}
+
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup spm_internal
+ *
+ * dofVar - Expand a CSC matrix without dofs into a CSC matrix with variable dofs
+ *
+ *
+ *******************************************************************************
+ *
+ * @param[in] spm
+ *           The spm to expand
+ *
+ *******************************************************************************/
+void
+dofVar(pastix_spm_t* spm)
+{
+    pastix_int_t *dofs=malloc(sizeof(pastix_int_t)*spm->n);
+    pastix_int_t i,nnzexp,nexp;
+
+    spm->gNexp    = spm->gN;
+    spm->nexp     = spm->n;
+    spm->gnnzexp  = spm->gnnz;
+    spm->nnzexp   = spm->nnz;
+    spm->dof      = -1;
+
+    genDof(3,dofs,spm->n); //cyclique, pas de 3
+    //genRandDof(/*min*/1,/*max*/5 ,spm->dofs,spm->n); //random dofs, entre 1 et 5
+
+    spm->dofs = dofs;
+    pastix_int_t *coef = computeEltPerCol(spm); // free coef
+
+    nnzexp = 0;
+    nexp   = 0;
+
+    for(i=0;i<spm->n;i++)
+    {
+        nnzexp += coef[i];
+        nexp   += spm->dofs[i];
+    }
+
+    spm->nnzexp  = nnzexp;
+    spm->gnnzexp = nnzexp;
+    spm->gNexp   = nexp;
+    spm->nexp    = nexp;
+
+    pastix_int_t* tab = computeCorrespondIndex(dofs,spm->n);
+
+    spm->dofs   = tab;
+    spmExpandCSC(spm);
+
+    /*
+    print_tab_int(spm->colptr,spm->n);
+    print_tab_int(spm->rowptr,spm->nnz);
+    print_tab_cpx(spm->values,spm->nnzexp);
+    print_tab_int(spm->dofs,spm->n+1);
+     */
+}
+
+

z_spm.h

@@ -41,4 +41,11 @@ pastix_complex64_t *z_spm2dense( const pastix_spm_t *spm );
 void z_spmDensePrint( pastix_int_t m, pastix_int_t n, pastix_complex64_t *A, pastix_int_t lda );
 void z_spmPrint( const pastix_spm_t *spm );
 
+void z_spmConvertColMaj2RowMaj(pastix_spm_t *spm);
+void z_spmConvertRowMaj2ColMaj(pastix_spm_t *spm);
+
+void z_extandCSC(pastix_spm_t *spm);
+
+void z_spmDofs2Flat(pastix_spm_t *spm);
+
 #endif /* _z_spm_H_ */

z_spm_convert_col_row_major.c

+/**
+ *
+ * @file z_spm_convert_col_row_major.c
+ *
+ *  PaStiX spm routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @author Alban Bellot
+ * @date 2015-01-01
+ *
+ * @precisions normal z -> c d s p
+ **/
+#include "common.h"
+#include "spm.h"
+#include "z_spm.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_spm
+ *
+ * z_spmConvertColMaj2RowMaj - convert a matrix in Column Major format to a
+ * matrix in Row Major format.
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          The colmaj matrix at enter,
+ *          the rowmaj matrix at exit.
+ *
+ *******************************************************************************/
+void
+z_spmConvertColMaj2RowMaj(pastix_spm_t *spm)
+{
+    assert(spm->colmajor > 0);
+    assert(spm->dof != 1);
+    assert(spm->dofs);
+
+    spm->colmajor = -1;
+
+    pastix_int_t dofi, dofj, ii, jj, i, j, k, baseval;
+    pastix_int_t cpt=0;
+    pastix_int_t *dofs = spm->dofs;
+    pastix_int_t *tmp;
+    pastix_complex64_t *oavals, *navals;
+
+    baseval = spmFindBase( spm );
+
+    oavals = (pastix_complex64_t*)spm->values;
+    navals = malloc(sizeof(pastix_complex64_t)*spm->nnzexp);
+
+    switch(spm->fmttype)
+    {
+    case PastixCSC:
+        for(i=0; i<spm->n; i++)
+        {
+            //dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            dofi = dofs[i+1] - dofs[i];
+            for(k=spm->colptr[i]; k<spm->colptr[i+1]; k++)
+            {
+                j = spm->rowptr[k-baseval] - baseval;
+                //dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                dofj = 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;
+        break;
+
+    case PastixCSR:
+        tmp           = spm->rowptr;
+        spm->rowptr   = spm->colptr;
+        spm->colptr   = tmp;
+        spm->fmttype  = PastixCSC;
+
+        z_spmConvertRowMaj2ColMaj(spm);
+
+        spm->colmajor = -1;
+        tmp          = spm->rowptr;
+        spm->rowptr  = spm->colptr;
+        spm->colptr  = tmp;
+        spm->fmttype = PastixCSR;
+        break;
+
+    case PastixIJV:
+        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];
+            dofi = dofs[i+1] - dofs[i];
+            dofj = 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;
+    default:
+        break;
+    }
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_spm
+ *
+ * z_spmConvertColMaj2RowMaj - convert a matrix in Row Major format to a matrix
+ * in Column Major format.
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          The rowmaj matrix at enter,
+ *          the colmaj matrix at exit.
+ *
+ *******************************************************************************/
+void
+z_spmConvertRowMaj2ColMaj(pastix_spm_t *spm)
+{
+    assert(spm->colmajor < 0);
+    assert(spm->dof != 1);
+    assert(spm->dofs);
+
+    spm->colmajor = 1;
+    pastix_int_t dofi, dofj, ii, jj, i, j, k, baseval;
+    pastix_int_t cpt=0;
+    pastix_int_t *dofs = spm->dofs;
+    pastix_int_t *tmp;
+    pastix_complex64_t *oavals, *navals;
+
+    baseval = spmFindBase( spm );
+
+    oavals = (pastix_complex64_t*)spm->values;
+    navals = malloc(sizeof(pastix_complex64_t)*spm->nnzexp);
+
+    switch(spm->fmttype)
+    {
+    case PastixCSC :
+        for(i=0; i<spm->n; i++)
+        {
+            //dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            dofi = dofs[i+1] - dofs[i];
+            for(k=spm->colptr[i]; k<spm->colptr[i+1]; k++)
+            {
+                j = spm->rowptr[k-baseval]-baseval;
+                //dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                dofj = 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;
+        break;
+
+    case PastixCSR :
+        tmp          = spm->rowptr;
+        spm->rowptr  = spm->colptr;
+        spm->colptr  = tmp;
+        spm->fmttype = PastixCSC;
+
+        z_spmConvertColMaj2RowMaj(spm);
+
+        spm->colmajor = 1;
+        tmp          = spm->rowptr;
+        spm->rowptr  = spm->colptr;
+        spm->colptr  = tmp;
+        spm->fmttype = PastixCSR;
+        break;
+
+    case PastixIJV:
+        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];
+            dofi = dofs[i+1] - dofs[i];
+            dofj = 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;
+    default:
+        break;
+    }
+}

z_spm_convert_to_csc.c

@@ -162,7 +162,7 @@ z_spmConvertIJV2CSC( pastix_spm_t *spm )
     spmExit( &oldspm );
 
     spm->fmttype = PastixCSC;
-
+    spm->colmajor = 1;
     return PASTIX_SUCCESS;
 }
 
@@ -354,6 +354,6 @@ z_spmConvertCSR2CSC( pastix_spm_t *spm )
 
     if(spm-> dof != 1)
         spm->mtxtype = type;
-
+    spm->colmajor = 1;
     return PASTIX_SUCCESS;
 }

z_spm_convert_to_csr.c

@@ -97,6 +97,7 @@ z_spmConvertCSC2CSR( pastix_spm_t *spm )
 
     if(spm-> dof != 1)
         spm->mtxtype = type;
+    spm->colmajor = -1;
 
     return result;
 }
@@ -144,10 +145,10 @@ z_spmConvertIJV2CSR( pastix_spm_t *spm )
     baseval = spmFindBase( spm );
 
 #if !defined(PRECISION_p)
+    pastix_int_t ii, jj, k;
     /* Transpose values in row major format */
     if( spm->colmajor )
     {
-        pastix_int_t ii, jj, k;
         pastix_int_t cpt=0;
         pastix_int_t* dofs = spm->dofs;
 
@@ -278,6 +279,7 @@ z_spmConvertIJV2CSR( pastix_spm_t *spm )
 
     spmExit( &oldspm );
 
+    spm->colmajor = -1;
     spm->fmttype = PastixCSR;
 
     return PASTIX_SUCCESS;

z_spm_convert_to_ijv.c

@@ -64,7 +64,8 @@ z_spmConvertCSC2IJV( pastix_spm_t *spm )
     }
 
     /* Transpose values in row major format */
-    if( !spm->colmajor ) //A test
+    /*
+    if( spm->colmajor < 1 ) //A test
     {
         int k, ii, jj, dofi, dofj;
         int cpt = 0;
@@ -86,11 +87,13 @@ z_spmConvertCSC2IJV( pastix_spm_t *spm )
                     oavals++;
                 }
             }
-            cpt + = dofi * dofj;
+            cpt += dofi * dofj;
         }
         spm->values = navals;
     }
 
+    */
+
     memFree_null(spm->colptr);
     spm->colptr = col_ijv;
 
@@ -143,7 +146,8 @@ z_spmConvertCSR2IJV( pastix_spm_t *spm )
     }
 
     /* Transpose values in column major format */
-    if( spm->colmajor )
+    /*
+    if( spm->colmajor > 1 )
     {
         pastix_int_t k, ii, jj, dofi, dofj;
         pastix_int_t cpt=0;
@@ -168,6 +172,7 @@ z_spmConvertCSR2IJV( pastix_spm_t *spm )
         }
         spm->values = navals;
     }
+     */
 
     memFree_null(spm->rowptr);
     spm->rowptr = row_ijv;

z_spm_dofs2flat.c

@@ -44,83 +44,190 @@
 void
 z_spmDofs2Flat(pastix_spm_t *spm)
 {
-    pastix_int_t i, j, k, ii, jj, dofi, dofj, col, row, baseval;
+    pastix_int_t i, j, k, ii, jj, dofi, dofj, col, row, baseval, cpt;
     baseval = spmFindBase( spm );
 
     pastix_int_t       *new_col  = calloc(spm->nexp+1,sizeof(pastix_int_t));
-    pastix_int_t       *new_row  = malloc(sizeof(pastix_int_t)*spm->nnzexp);
+    pastix_int_t       *new_row;
     pastix_int_t       *dofs     = spm->dofs;
 
-    pastix_complex64_t *new_vals = malloc(sizeof(pastix_int_t)*spm->nnzexp);
+    pastix_complex64_t *new_vals;
     pastix_complex64_t *vals     = (pastix_complex64_t*)spm->values;
 
-    for(i=0;i<spm->n ; i++)
+    switch(spm->mtxtype)
     {
-        col = ( spm->dof > 0 ) ? i : dofs[i];
-        for(k=spm->colptr[i]-baseval; k<spm->colptr[i+1]-baseval; k++)
-	{
-            j = spm->rowptr[k]-baseval;
-            row  = ( spm->dof > 0 ) ? j        : dofs[j];
+    case PastixGeneral:
+        new_row  = malloc(sizeof(pastix_int_t)*spm->nnzexp);
+        new_vals = malloc(sizeof(pastix_complex64_t)*spm->nnzexp);
+        for(i=0;i<spm->n ; i++)
+        {
+            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];
-            for(ii=0; ii<dofi; ii++)
-	    {
-                new_col[col+ii+1] +=  dofj;
-	    }
-	}
-    }
+            for(k=spm->colptr[i]-baseval; k<spm->colptr[i+1]-baseval; k++)
+            {
+                j = spm->rowptr[k]-baseval;
+                row  = ( spm->dof > 0 ) ? j        : dofs[j];
+                dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                for(ii=0; ii<dofi; ii++)
+                {
+                    new_col[col+ii+1] +=  dofj;
+                }
+            }
+        }
 
-    for(i=0; i<spm->nexp; i++)
-    {
-        new_col[i+1]+=new_col[i];
-    }
+        for(i=0; i<spm->nexp; i++)
+        {
+            new_col[i+1]+=new_col[i];
+        }
 
-    int cpt = 0;
-    for(i=0; i < spm->n;i++)
-    {
-        col  = ( spm->dof > 0 ) ? i        : dofs[i];
-        dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
-        for(k=spm->colptr[i]-baseval ; k<spm->colptr[i+1]-baseval ;k++)
-	{
-            j = spm->rowptr[k] - baseval;
-            row  = ( spm->dof > 0 ) ? j        : dofs[j];
-            dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
-            for(ii=0;ii < dofi; ii++)
-	    {
-                for(jj=0;jj < dofj ; jj++)
-		{
-                    new_vals[new_col[col+ii]] = vals[cpt];
-                    new_row[new_col[col+ii]]  = row + jj + baseval;
-                    cpt++;
-                    new_col[col+ii]++;
-		}
-	    }
-	}
-    }
+        cpt = 0;
+        for(i=0; i < spm->n;i++)
+        {
+            col  = ( spm->dof > 0 ) ? i        : dofs[i];
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            for(k=spm->colptr[i]-baseval ; k<spm->colptr[i+1]-baseval ;k++)
+            {
+                j = spm->rowptr[k] - baseval;
+                row  = ( spm->dof > 0 ) ? j        : dofs[j];
+                dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                for(ii=0;ii < dofi; ii++)
+                {
+                    for(jj=0;jj < dofj ; jj++)
+                    {
+                        new_vals[new_col[col+ii]] = vals[cpt];
+                        new_row[new_col[col+ii]]  = row + jj + baseval;
+                        new_col[col+ii]++;
+                        cpt++;
+                    }
+                }
+            }
+        }
 
-    {
-        int tmp;
-        int tmp1 = 0;
+        {
+            int tmp;
+            int tmp1 = 0;
+            for(i=0; i<spm->nexp; i++)
+            {
+                tmp = new_col[i];
+                new_col[i] = tmp1+baseval;
+                tmp1 = tmp;
+            }
+            new_col[i] += baseval;
+        }
+        spm->gN   = spm->gNexp;
+        spm->n    = spm->nexp;
+        spm->gnnz = spm->gnnzexp;
+        spm->nnz  = spm->nnzexp;
+
+        spm->dof       = 1;
+        spm->dofs      = NULL;
+        spm->colmajor  = 1;
+
+        spm->colptr   = new_col;
+        spm->rowptr   = new_row;
+        //spm->loc2glob = NULL; // ?
+        spm->values   = new_vals;
+        break;
+
+    case PastixSymmetric:
+        for(i=0;i<spm->n ; i++)
+        {
+            col  = ( spm->dof > 0 ) ? i        : dofs[i];
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            for(k=spm->colptr[i]-baseval; k<spm->colptr[i+1]-baseval; k++)
+            {
+                j = spm->rowptr[k]-baseval;
+                row  = ( spm->dof > 0 ) ? j        : dofs[j];
+                dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                for(ii=0; ii<dofi; ii++)
+                {
+                    for(jj=0; jj<dofj; jj++)
+                    {
+                        if( i != j )
+                            new_col[col+ii+1] +=  1;
+                        else
+                            if(ii <= jj )
+                                new_col[col+ii+1] += 1;
+                    }
+                }
+            }
+        }
         for(i=0; i<spm->nexp; i++)
         {
-            tmp = new_col[i];
-            new_col[i] = tmp1+baseval;
-            tmp1 = tmp;
+            new_col[i+1] += new_col[i];
         }
-        new_col[i] += baseval;
-    }
+        pastix_int_t nnz = new_col[spm->nexp];
+        new_row  = malloc(sizeof(pastix_int_t)*nnz);
+        new_vals = malloc(sizeof(pastix_complex64_t)*nnz);
 
-    spm->gN   = spm->gNexp;
-    spm->n    = spm->nexp;
-    spm->gnnz = spm->gnnzexp;
-    spm->nnz  = spm->nnzexp;
+        cpt = 0;
+        for(i=0; i < spm->n;i++)
+        {
+            col  = ( spm->dof > 0 ) ? i        : dofs[i];
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            for(k=spm->colptr[i]-baseval ; k<spm->colptr[i+1]-baseval ;k++)
+            {
+                j = spm->rowptr[k] - baseval;
+                row  = ( spm->dof > 0 ) ? j        : dofs[j];
+                dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                for(ii=0;ii < dofi; ii++)
+                {
+                    for(jj=0;jj < dofj ; jj++)
+                    {
+                        if( i == j )
+                        {
+                            if ( ii <= jj )
+                            {
+                                /* diagonal dominant for spd matrix
+                                if( ii == jj)
+                                    new_vals[new_col[col+ii]] = 2*vals[cpt];
+                                 else
+                                */
+                                new_vals[new_col[col+ii]] = vals[cpt];
+                                new_row[new_col[col+ii]]  = row + jj + baseval;
+                                new_col[col+ii]++;
+                            }
+                        }
+                        else
+                        {
+                            new_vals[new_col[col+ii]] = vals[cpt];
+                            new_row[new_col[col+ii]]  = row + jj + baseval;
+                            new_col[col+ii]++;
 
-    spm->dof       = 1;
-    spm->dofs      = NULL;
-    spm->colmajor  = 1;
+                        }
+                        cpt++;
+                    }
+                }
+            }
+        }
+        {
+            int tmp;
+            int tmp1 = 0;
+            for(i=0; i<spm->nexp; i++)
+            {
+                tmp = new_col[i];
+                new_col[i] = tmp1+baseval;
+                tmp1 = tmp;
+            }
+            new_col[i] += baseval;
+        }
+        spm->gN   = spm->gNexp;
+        spm->n    = spm->nexp;
+        spm->gnnz = nnz;
+        spm->nnz  = nnz;
+        spm->gnnzexp = nnz;
+        spm->nnzexp  = nnz;
 
-    spm->colptr   = new_col;
-    spm->rowptr   = new_row;
-    spm->loc2glob = NULL; // ?
-    spm->values   = new_vals;
+        spm->dof       = 1;
+        spm->dofs      = NULL;
+        spm->colmajor  = 1;
+
+        spm->colptr   = new_col;
+        spm->rowptr   = new_row;
+        //spm->loc2glob = NULL; //
+        spm->values   = new_vals;
+
+        break;
+    }
+    assert(spm->loc2glob == NULL);//to do
 }

z_spm_expand.c

+/**
+ *
+ * @file z_spm_expand.c
+ *
+ *  PaStiX spm routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @author Alban Bellot
+ * @date 2015-01-01
+ *
+ * @precisions normal z -> c d s p
+ **/
+#include "common.h"
+#include "spm.h"
+#include "z_spm.h"
+
+
+void z_extandCSC(pastix_spm_t* spm)
+{
+  pastix_int_t col, row, i, cpt, dofj, dofi, baseval;
+  cpt=0;
+
+  pastix_complex64_t *valptr = (pastix_complex64_t*)spm->values;
+  pastix_complex64_t *newval = calloc(spm->nnzexp,sizeof(pastix_complex64_t));
+
+  baseval=spmFindBase( spm );
+
+  for( col=0; col<spm->n; col++)
+  {
+      for( row=spm->colptr[col]-baseval; row<spm->colptr[col+1]-baseval; row++)
+      {
+          dofi = ( spm->dof > 0 ) ? spm->dof : spm->dofs[col+1] - spm->dofs[col];
+          dofj = ( spm->dof > 0 ) ? spm->dof : spm->dofs[spm->rowptr[row]-baseval+1] - spm->dofs[spm->rowptr[row]-baseval];
+	  for( i=0; i<dofi*dofj; i++)
+          {
+              newval[cpt] = valptr[row] / ((i/dofj) + (i%dofj) + 2); // Col major
+	      cpt++;
+          }
+      }
+  }
+  spm->values=newval;
+}