Merge default

CMakeLists.txt

+include(RulesPrecisions)
+
+### reset variables
+set(generated_sources "")
+set(generated_headers "")
+
+set(HEADERS
+  z_spm.h
+)
+
+### Generate the headers in all precisions
+precisions_rules_py(generated_headers
+                    "${HEADERS}"
+                    PRECISIONS "p;s;d;c;z")
+
+set(SOURCES
+  z_spm.c
+  z_spm_convert_to_csc.c
+  z_spm_convert_to_csr.c
+  z_spm_convert_to_ijv.c
+  z_spm_genrhs.c
+  z_spm_matrixvector.c
+  z_spm_norm.c
+  z_spm_2dense.c
+  z_spm_print.c
+  )
+
+set(spm_headers
+  ${generated_headers}
+  spm.h
+  )
+
+### Generate the sources in all precisions
+precisions_rules_py(generated_sources
+                    "${SOURCES}"
+                    PRECISIONS "p;s;d;c;z")
+
+set(spm_sources
+  ${generated_sources}
+  spm.c
+  spm_io.c
+  )
+
+add_custom_target(spm_headers_tgt
+  DEPENDS ${spm_headers} )
+
+set_source_files_properties(
+  spm.c
+  PROPERTIES DEPENDS spm_headers_tgt
+  )
+
+set(PASTIX_LIB_SRCS
+  ${PASTIX_LIB_SRCS}
+  ${spm_sources}
+  )

spm.c

@@ -68,6 +68,46 @@ static int (*conversionTable[3][3][6])(pastix_spm_t*) = {
 };
 
 
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_spm
+ *
+ * spmInit - Init the spm structure given as parameter
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          The sparse matrix to init.
+ *
+ *******************************************************************************/
+void
+spmInit( pastix_spm_t *spm )
+{
+    spm->mtxtype = PastixGeneral;
+    spm->flttype = PastixComplex64;
+    spm->fmttype = PastixCSC;
+
+    spm->gN   = 0;
+    spm->n    = 0;
+    spm->gnnz = 0;
+    spm->nnz  = 0;
+
+    spm->gNexp   = 0;
+    spm->nexp    = 0;
+    spm->gnnzexp = 0;
+    spm->nnzexp  = 0;
+
+    spm->dof       = 1;
+    spm->dofs      = NULL;
+    spm->colmajor  = 1;
+
+    spm->colptr   = NULL;
+    spm->rowptr   = NULL;
+    spm->loc2glob = NULL;
+    spm->values   = NULL;
+}
+
 /**
  *******************************************************************************
  *
@@ -106,6 +146,10 @@ spmExit( pastix_spm_t *spm )
  *
  *******************************************************************************
  *
+ * @param[in] ofmttype
+ *          The output format of the sparse matrix. It might be PastixCSC,
+ *          PastixCSR, or PastixIJV.
+ *
  * @param[in,out] spm
  *          The sparse matrix to rebase.
  *

spm.h

@@ -36,14 +36,24 @@ struct pastix_spm_s {
     int           flttype;   /**< avals datatype: PastixPattern, PastixFloat, PastixDouble,
                                   PastixComplex32 or PastixComplex64                          */
     int           fmttype;   /**< Matrix storage format: PastixCSC, PastixCSR, PastixIJV      */
+
     pastix_int_t  gN;        /**< Global number of vertices in the compressed graph           */
     pastix_int_t  n;         /**< Local number of vertices in the compressed graph            */
     pastix_int_t  gnnz;      /**< Global number of non zeroes in the compressed graph         */
     pastix_int_t  nnz;       /**< Local number of non zeroes in the compressed graph          */
+
+    pastix_int_t  gNexp;     /**< Global number of vertices in the compressed graph           */
+    pastix_int_t  nexp;      /**< Local number of vertices in the compressed graph            */
+    pastix_int_t  gnnzexp;   /**< Global number of non zeroes in the compressed graph         */
+    pastix_int_t  nnzexp;    /**< Local number of non zeroes in the compressed graph          */
+
     pastix_int_t  dof;       /**< Number of degrees of freedom per unknown,
                                   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 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                               */
     pastix_int_t *loc2glob;  /**< Corresponding numbering from local to global                */
@@ -89,4 +99,6 @@ pastix_int_t  spmSymmetrize( pastix_spm_t *spm );
 
 pastix_spm_t *spmCheckAndCorrect( pastix_spm_t *spm );
 
+void dofVar(pastix_spm_t* spm);//tmp
+
 #endif /* _SPM_H_ */

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

@@ -37,4 +37,15 @@ pastix_int_t z_spmSymmetrize( pastix_spm_t *csc );
 int z_spmGenRHS(int type, int nrhs, const pastix_spm_t *spm, void *x, int ldx, void *b, int ldb );
 int z_spmCheckAxb( int nrhs, const pastix_spm_t *spm, void *x0, int ldx0, void *b, int ldb, const void *x, int ldx );
 
+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_2dense.c

+/**
+ *
+ * @file z_spm_2dense.c
+ *
+ * Convert a sparse matrix into a dense matrix.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @author Alban Bellot
+ * @date 2015-01-01
+ *
+ * @precisions normal z -> c d s
+ *
+ **/
+#include <stdint.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include "pastix.h"
+#include "common.h"
+#include "spm.h"
+#include "z_spm.h"
+
+pastix_complex64_t *
+z_spm2dense( const pastix_spm_t *spm )
+{
+    pastix_int_t i, j, lda, baseval;
+    pastix_complex64_t *A, *valptr;
+    pastix_int_t *colptr, *rowptr;
+
+    assert( spm->fmttype == PastixCSC );
+    assert( spm->flttype == PastixComplex64 );
+
+    lda = spm->gNexp;
+    A = (pastix_complex64_t*)malloc(lda * lda * sizeof(pastix_complex64_t));
+    memset( A, 0, lda * lda * sizeof(pastix_complex64_t));
+
+    baseval = spmFindBase( spm );
+    i = 0; j = 0;
+
+    colptr = spm->colptr;
+    rowptr = spm->rowptr;
+    valptr = (pastix_complex64_t*)(spm->values);
+
+    /**
+     * Constant degree of fredom of 1
+     */
+    if ( spm->dof == 1 ) {
+        switch( spm->mtxtype ){
+#if defined(PRECISION_z) || defined(PRECISION_c)
+        case PastixHermitian:
+            for(i=0; i<spm->n; i++, colptr++)
+            {
+                for(j=colptr[0]; j<colptr[1]; j++, rowptr++, valptr++)
+                {
+                    A[ i * lda + (*rowptr - baseval) ] = *valptr;
+                    A[ (*rowptr - baseval) * lda + i ] = conj(*valptr);
+                }
+            }
+            break;
+#endif
+        case PastixSymmetric:
+            for(i=0; i<spm->n; i++, colptr++)
+            {
+                for(j=colptr[0]; j<colptr[1]; j++, rowptr++, valptr++)
+                {
+                    A[ i * lda + (*rowptr - baseval) ] = *valptr;
+                    A[ (*rowptr - baseval) * lda + i ] = *valptr;
+                }
+            }
+            break;
+        case PastixGeneral:
+        default:
+            for(i=0; i<spm->n; i++, colptr++)
+            {
+                for(j=colptr[0]; j<colptr[1]; j++, rowptr++, valptr++)
+                {
+                    A[ i * lda + (*rowptr - baseval) ] = *valptr;
+                }
+            }
+        }
+    }
+    /**
+     * General degree of freedom (constant or variable)
+     */
+    else {
+        pastix_int_t  k, ii, jj, dofi, dofj, col, row;
+        pastix_int_t *dofs = spm->dofs;
+
+        switch( spm->mtxtype ){
+#if defined(PRECISION_z) || defined(PRECISION_c)
+        case PastixHermitian:
+            for(i=0; i<spm->n; i++, colptr++)
+            {
+                dofi = ( spm->dof > 1 ) ? spm->dof : dofs[i+1] - dofs[i];
+                col = dofs[i];
+
+                for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
+                {
+                    j = (*rowptr - baseval);
+                    dofj = ( spm->dof > 1 ) ? spm->dof : dofs[j+1] - dofs[j];
+                    row = dofs[j];
+
+                    for(ii=0; ii<dofi; ii++)
+                    {
+                        for(jj=0; jj<dofj; jj++, valptr++)
+                        {
+                            A[ (col + ii) * lda + (row + jj) ] = *valptr;
+                            A[ (row + jj) * lda + (col + ii) ] = conj(*valptr);
+                        }
+                    }
+                }
+            }
+            break;
+#endif
+        case PastixSymmetric:
+            for(i=0; i<spm->n; i++, colptr++)
+            {
+                dofi = ( spm->dof > 1 ) ? spm->dof : dofs[i+1] - dofs[i];
+                col = dofs[i];
+
+                for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
+                {
+                    j = (*rowptr - baseval);
+                    dofj = ( spm->dof > 1 ) ? spm->dof : dofs[j+1] - dofs[j];
+                    row = dofs[j];
+
+                    for(ii=0; ii<dofi; ii++)
+                    {
+                        for(jj=0; jj<dofj; jj++, valptr++)
+                        {
+                            A[ (col + ii) * lda + (row + jj) ] = *valptr;
+                            A[ (row + jj) * lda + (col + ii) ] = *valptr;
+                        }
+                    }
+                }
+            }
+            break;
+        case PastixGeneral:
+        default:
+            for(i=0; i<spm->n; i++, colptr++)
+            {
+                dofi = ( spm->dof > 1 ) ? spm->dof : dofs[i+1] - dofs[i];
+                col = dofs[i];
+
+                for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
+                {
+                    j = (*rowptr - baseval);
+                    dofj = ( spm->dof > 1 ) ? spm->dof : dofs[j+1] - dofs[j];
+                    row = dofs[j];
+
+                    for(ii=0; ii<dofi; ii++)
+                    {
+                        for(jj=0; jj<dofj; jj++, valptr++)
+                        {
+                            A[ (col + ii) * lda + (row + jj) ] = *valptr;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    return A;
+}
+
+void
+z_spmDensePrint( pastix_int_t m, pastix_int_t n, pastix_complex64_t *A, pastix_int_t lda )
+{
+    pastix_int_t i, j;
+
+    for(i=0; i<m; i++)
+    {
+        for(j=0; j<n; j++)
+        {
+            if ( cabs( A[ lda *j + i] ) != 0. ) {
+#if defined(PRECISION_z) || defined(PRECISION_c)
+                fprintf( stderr, "%ld %ld (%e, %e)\n",
+                         i, j, creal(A[lda * j + i]), cimag(A[lda * j + i]) );
+#else
+                fprintf( stderr, "%ld %ld %e\n",
+                         i, j, A[lda * j + i] );
+#endif
+            }
+        }
+    }
+    return;
+}
+

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

@@ -9,6 +9,7 @@
  * @version 5.1.0
  * @author Mathieu Faverge
  * @author Theophile Terraz
+ * @author Alban Bellot
  * @date 2015-01-01
  *
  * @precisions normal z -> c d s p
@@ -80,11 +81,45 @@ 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 *dofs     = spm->dofs;
+    pastix_int_t row, col, dofi, dofj;
+
+    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
@@ -92,11 +127,20 @@ z_spmConvertIJV2CSC( pastix_spm_t *spm )
     for (j=0; j<spm->nnz; j++)
     {
         i = oldspm.colptr[j] - baseval;
-
         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])++;
 
@@ -106,7 +150,6 @@ z_spmConvertIJV2CSC( pastix_spm_t *spm )
     /* Rebuild the colptr with the correct baseval */
     tmp = spm->colptr[0];
     spm->colptr[0] = baseval;
-
     spmptx = spm->colptr + 1;
     for (i=1; i<(spm->n+1); i++, spmptx++)
     {
@@ -119,7 +162,7 @@ z_spmConvertIJV2CSC( pastix_spm_t *spm )
     spmExit( &oldspm );
 
     spm->fmttype = PastixCSC;
-
+    spm->colmajor = 1;
     return PASTIX_SUCCESS;
 }
 
@@ -152,6 +195,10 @@ z_spmConvertCSR2CSC( pastix_spm_t *spm )
 
     spm->fmttype = PastixCSC;
 
+    pastix_int_t type = spm->mtxtype;
+    if(spm->dof != 1)
+        spm->mtxtype = PastixGeneral;
+
     switch( spm->mtxtype ) {
 #if defined(PRECISION_z) || defined(PRECISION_c)
     case PastixHermitian:
@@ -184,26 +231,36 @@ 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;
@@ -217,6 +274,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++) {
@@ -227,27 +313,34 @@ 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];
+                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;
@@ -259,5 +352,8 @@ 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

@@ -9,6 +9,7 @@
  * @version 5.1.0
  * @author Mathieu Faverge
  * @author Theophile Terraz
+ * @authot Alban Bellot
  * @date 2015-01-01
  *
  * @precisions normal z -> c d s p
@@ -43,6 +44,10 @@ z_spmConvertCSC2CSR( pastix_spm_t *spm )
 {
     pastix_int_t *tmp;
     pastix_int_t  result;
+    pastix_int_t  type = spm->mtxtype;
+
+    if(spm->dof != 1)
+        spm->mtxtype = PastixGeneral;
 
     switch( spm->mtxtype ) {
 #if defined(PRECISION_z) || defined(PRECISION_c)
@@ -66,7 +71,6 @@ z_spmConvertCSC2CSR( pastix_spm_t *spm )
         spm->rowptr  = spm->colptr;
         spm->colptr  = tmp;
         spm->fmttype = PastixCSR;
-
         return PASTIX_SUCCESS;
     }
     break;
@@ -91,6 +95,10 @@ z_spmConvertCSC2CSR( pastix_spm_t *spm )
     }
     }
 
+    if(spm-> dof != 1)
+        spm->mtxtype = type;
+    spm->colmajor = -1;
+
     return result;
 }
 
@@ -123,16 +131,53 @@ 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, k;
+    /* Transpose values in row major format */
+    if( spm->colmajor )
+    {
+        pastix_int_t cpt=0;
+        pastix_int_t* dofs = spm->dofs;
+
+        oavals = (pastix_complex64_t*)spm->values;
+        navals = malloc(sizeof(pastix_complex64_t)*spm->nnzexp);
+
+        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 +203,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 +247,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])++;
 
@@ -195,7 +279,9 @@ z_spmConvertIJV2CSR( pastix_spm_t *spm )
 
     spmExit( &oldspm );
 
+    spm->colmajor = -1;
     spm->fmttype = PastixCSR;
 
     return PASTIX_SUCCESS;
 }
+

z_spm_convert_to_ijv.c

@@ -9,6 +9,7 @@
  * @version 5.1.0
  * @author Mathieu Faverge
  * @author Theophile Terraz
+ * @author Alban Bellot
  * @date 2015-01-01
  *
  * @precisions normal z -> c d s p
@@ -62,6 +63,37 @@ z_spmConvertCSC2IJV( pastix_spm_t *spm )
         }
     }
 
+    /* Transpose values in row major format */
+    /*
+    if( spm->colmajor < 1 ) //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 +145,35 @@ z_spmConvertCSR2IJV( pastix_spm_t *spm )
         }
     }
 
+    /* Transpose values in column major format */
+    /*
+    if( spm->colmajor > 1 )
+    {
+        pastix_int_t k, ii, jj, dofi, dofj;
+        pastix_int_t 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;
 

z_spm_dofs2flat.c

+/**
+ *
+ * @file z_spm_dofs2flat.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
+ *
+ **/
+#include <stdint.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include "pastix.h"
+#include "common.h"
+#include "spm.h"
+#include "z_spm.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup spm_internal
+ *
+ * z_spmDofs2Flat - Convert a sparse matrix with dofs into a sparse matrix without dofs.
+ *
+ *
+ *******************************************************************************
+ *
+ * @param[in] spm
+ *           The sparse matrix which needs to be converted
+ *
+ *******************************************************************************
+ **/
+void
+z_spmDofs2Flat(pastix_spm_t *spm)
+{
+    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;
+    pastix_int_t       *dofs     = spm->dofs;
+
+    pastix_complex64_t *new_vals;
+    pastix_complex64_t *vals     = (pastix_complex64_t*)spm->values;
+
+    switch(spm->mtxtype)
+    {
+    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];
+            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];
+        }
+
+        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;
+            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++)
+        {
+            new_col[i+1] += new_col[i];
+        }
+        pastix_int_t nnz = new_col[spm->nexp];
+        new_row  = malloc(sizeof(pastix_int_t)*nnz);
+        new_vals = malloc(sizeof(pastix_complex64_t)*nnz);
+
+        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]++;
+
+                        }
+                        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->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;
+}

z_spm_matrixvector.c

@@ -73,7 +73,9 @@ z_spmGeCSCv(      int                 trans,
     const pastix_complex64_t *valptr = (pastix_complex64_t*)csc->values;
     const pastix_complex64_t *xptr   = x;
     pastix_complex64_t *yptr = y;
-    pastix_int_t col, row, i, baseval;
+    pastix_int_t col, row, i, j, baseval;
+    pastix_int_t ii, jj, k, dofi, dofj;
+    pastix_int_t *dofs = csc->dofs;
 
     if ( (csc == NULL) || (x == NULL) || (y == NULL ) )
     {
@@ -89,10 +91,10 @@ z_spmGeCSCv(      int                 trans,
 
     /* first, y = beta*y */
     if( beta == 0. ) {
-        memset( yptr, 0, csc->gN * sizeof(pastix_complex64_t) );
+        memset( yptr, 0, csc->gNexp * sizeof(pastix_complex64_t) );
     }
     else {
-        for( i=0; i<csc->gN; i++, yptr++ ) {
+        for( i=0; i<csc->gNexp; i++, yptr++ ) {
             (*yptr) *= beta;
         }
         yptr = y;
@@ -104,12 +106,22 @@ z_spmGeCSCv(      int                 trans,
          */
         if( trans == PastixNoTrans )
         {
-            for( col=0; col < csc->gN; col++ )
+            for( i=0; i < csc->gN; i++ )
             {
-                for( i=csc->colptr[col]; i<csc->colptr[col+1]; i++ )
+                dofi = ( csc->dof > 0 ) ? csc->dof : dofs[i+1] - dofs[i];
+                col  = ( csc->dof > 0 ) ? i        : dofs[i];
+                for( k=csc->colptr[i]; k<csc->colptr[i+1]; k++ )
                 {
-                    row = csc->rowptr[i-baseval]-baseval;
-                    yptr[row] += alpha * valptr[i-baseval] * xptr[col];
+                    j = csc->rowptr[k-baseval]-baseval;
+                    dofj = ( csc->dof > 0 ) ? csc->dof : dofs[j+1] - dofs[j];
+                    row  = ( csc->dof > 0 ) ? j        : dofs[j];
+                    for(ii=0; ii<dofi; ii++)
+                    {
+                        for(jj=0; jj<dofj; jj++, valptr++)
+                        {
+                            yptr[row+jj] += alpha * (*valptr) * xptr[col+ii];
+                        }
+                    }
                 }
             }
         }
@@ -118,15 +130,26 @@ z_spmGeCSCv(      int                 trans,
          */
         else if( trans == PastixTrans )
         {
-            for( col=0; col < csc->gN; col++ )
+            for( i=0; i < csc->gN; i++ )
             {
-                for( i=csc->colptr[col]; i<csc->colptr[col+1]; i++ )
+                dofi = ( csc->dof > 0 ) ? csc->dof : dofs[i+1] - dofs[i];
+                col  = ( csc->dof > 0 ) ? i        : dofs[i];
+                for( k=csc->colptr[i]; k<csc->colptr[i+1]; k++ )
                 {
-                    row = csc->rowptr[i-baseval]-baseval;
-                    yptr[col] += alpha * valptr[i-baseval] * xptr[row];
+                    j = csc->rowptr[k-baseval]-baseval;
+                    dofj = ( csc->dof > 0 ) ? csc->dof : dofs[j+1] - dofs[j];
+                    row  = ( csc->dof > 0 ) ? j        : dofs[j];
+                    for(ii=0; ii<dofi; ii++)
+                    {
+                        for(jj=0; jj<dofj; jj++, valptr++)
+                        {
+                            yptr[col+ii] += alpha * (*valptr) * xptr[row+jj];
+                        }
+                    }
                 }
             }
         }
+
 #if defined(PRECISION_c) || defined(PRECISION_z)
         else if( trans == PastixConjTrans )
         {
@@ -145,7 +168,6 @@ z_spmGeCSCv(      int                 trans,
             return PASTIX_ERR_BADPARAMETER;
         }
     }
-
     return PASTIX_SUCCESS;
 }
 
@@ -194,7 +216,9 @@ z_spmSyCSCv(      pastix_complex64_t  alpha,
     const pastix_complex64_t *valptr = (pastix_complex64_t*)csc->values;
     const pastix_complex64_t *xptr   = x;
     pastix_complex64_t *yptr = y;
-    pastix_int_t col, row, i, baseval;
+    pastix_int_t col, row, i, j, baseval;
+    pastix_int_t ii, jj, k, dofi, dofj;
+    pastix_int_t *dofs = csc->dofs;
 
     if ( (csc == NULL) || (x == NULL) || (y == NULL ) )
     {
@@ -210,30 +234,40 @@ z_spmSyCSCv(      pastix_complex64_t  alpha,
 
     /* First, y = beta*y */
     if( beta == 0. ) {
-        memset( yptr, 0, csc->gN * sizeof(pastix_complex64_t) );
+        memset( yptr, 0, csc->gNexp * sizeof(pastix_complex64_t) );
     }
     else {
-        for( i=0; i<csc->gN; i++, yptr++ ) {
+        for( i=0; i<csc->gNexp; i++, yptr++ ) {
             (*yptr) *= beta;
         }
         yptr = y;
     }
 
-    if( alpha != 0. ) {
-        for( col=0; col < csc->gN; col++ )
+    if(alpha != 0.)
+    {
+        for( i=0; i < csc->gN; i++ )
         {
-            for( i=csc->colptr[col]; i < csc->colptr[col+1]; i++ )
+            dofi = ( csc->dof > 0 ) ? csc->dof : dofs[i+1] - dofs[i];
+            col  = ( csc->dof > 0 )  ? i        : dofs[i];
+            for( k=csc->colptr[i]; k<csc->colptr[i+1]; k++ )
             {
-                row = csc->rowptr[i-baseval]-baseval;
-                yptr[row] += alpha * valptr[i-baseval] * xptr[col];
-                if( col != row )
+                j = csc->rowptr[k-baseval]-baseval;
+                dofj = ( csc->dof > 0 ) ? csc->dof : dofs[j+1] - dofs[j];
+                row  = ( csc->dof > 0 ) ? j        : dofs[j];
+                for(ii=0; ii<dofi; ii++)
                 {
-                    yptr[col] += alpha * valptr[i-baseval] * xptr[row];
+                    for(jj=0; jj<dofj; jj++, valptr++)
+                    {
+                        yptr[row+jj] += alpha * (*valptr) * xptr[col+ii];
+                        if( i != j )
+                        {
+                            yptr[col+ii] += alpha * (*valptr) * xptr[row+jj];
+                        }
+                    }
                 }
             }
         }
     }
-
     return PASTIX_SUCCESS;
 }
 
@@ -283,7 +317,9 @@ z_spmHeCSCv(      pastix_complex64_t  alpha,
     const pastix_complex64_t *valptr = (pastix_complex64_t*)csc->values;
     const pastix_complex64_t *xptr   = x;
     pastix_complex64_t *yptr = y;
-    pastix_int_t col, row, i, baseval;
+    pastix_int_t col, row, i, j, baseval;
+    pastix_int_t ii, jj, k, dofi, dofj;
+    pastix_int_t *dofs = csc->dofs;
 
     if ( (csc == NULL) || (x == NULL) || (y == NULL ) )
     {
@@ -297,30 +333,41 @@ z_spmHeCSCv(      pastix_complex64_t  alpha,
 
     /* First, y = beta*y */
     if( beta == 0. ) {
-        memset( yptr, 0, csc->gN * sizeof(pastix_complex64_t) );
+        memset( yptr, 0, csc->gNexp * sizeof(pastix_complex64_t) );
     }
     else {
-        for( i=0; i<csc->gN; i++, yptr++ ) {
+        for( i=0; i<csc->gNexp; i++, yptr++ ) {
             (*yptr) *= beta;
         }
         yptr = y;
     }
 
     baseval = spmFindBase( csc );
-
-    if( alpha != 0. ) {
-        for( col=0; col < csc->gN; col++ )
+    if( alpha != 0.)
+    {
+        for( i=0; i < csc->gN; i++ )
         {
-            for( i=csc->colptr[col]; i < csc->colptr[col+1]; i++ )
+            dofi = ( csc->dof > 0 ) ? csc->dof : dofs[i+1] - dofs[i];
+            col  = ( csc->dof > 0 ) ? i        : dofs[i];
+            for( k=csc->colptr[i]; k<csc->colptr[i+1]; k++ )
             {
-                row=csc->rowptr[i-baseval]-baseval;
-                yptr[row] += alpha * valptr[i-baseval] * xptr[col];
-                if( col != row )
-                    yptr[col] += alpha * conj( valptr[i-baseval] ) * xptr[row];
+                j = csc->rowptr[k-baseval]-baseval;
+                dofj = ( csc->dof > 0 ) ? csc->dof : dofs[j+1] - dofs[j];
+                row  = ( csc->dof > 0 ) ? j        : dofs[j];
+                for(ii=0; ii<dofi; ii++)
+                {
+                    for(jj=0; jj<dofj; jj++, valptr++)
+                    {
+                        yptr[row+jj] += alpha * (*valptr) * xptr[col+ii];
+                        if( i != j )
+                        {
+                            yptr[col+ii] += alpha * conj( *valptr ) * xptr[row+jj];
+                        }
+                    }
+                }
             }
         }
     }
-
     return PASTIX_SUCCESS;
 }
 #endif

z_spm_norm.c

@@ -50,7 +50,7 @@ z_spmFrobeniusNorm( const pastix_spm_t *spm )
     double sumsq = 0.;
 
     if (spm->mtxtype == PastixGeneral) {
-        for(i=0; i <spm->nnz; i++, valptr++) {
+        for(i=0; i <spm->nnzexp; i++, valptr++) {
             frobenius_update( 1, &scale, &sumsq, valptr );
 
 #if defined(PRECISION_z) || defined(PRECISION_c)
@@ -62,20 +62,32 @@ z_spmFrobeniusNorm( const pastix_spm_t *spm )
     else {
         pastix_int_t *colptr = spm->colptr;
         pastix_int_t *rowptr = spm->rowptr;
-        int nb;
+        pastix_int_t *dofs   = spm->dofs;
+
+        int nb, dofi, dofj, ii, jj;
         baseval = spmFindBase( spm );
 
         switch( spm->fmttype ){
         case PastixCSC:
-            for(i=0; i<spm->n; i++, colptr++) {
-                for(j=colptr[0]; j<colptr[1]; j++, rowptr++, valptr++) {
-                    nb = ( i == (*rowptr-baseval) ) ? 1 : 2;
-                    frobenius_update( nb, &scale, &sumsq, valptr );
+            for(i=0; i<spm->n; i++, colptr++)
+            {
+                dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+                for(j=colptr[0]; j<colptr[1]; j++, rowptr++)
+                {
+                    dofj = ( spm->dof > 0 ) ? spm->dof : dofs[(*rowptr-baseval)+1] - dofs[(*rowptr-baseval)];
+                    for(ii=0 ; ii < dofi; ii++)
+                    {
+                        for(jj=0; jj < dofj ;jj++,valptr++)
+                        {
+                            nb = ( i == (*rowptr-baseval) ) ? 1 : 2;
+                            frobenius_update( nb, &scale, &sumsq, valptr );
 
 #if defined(PRECISION_z) || defined(PRECISION_c)
-                    valptr++;
-                    frobenius_update( nb, &scale, &sumsq, valptr );
+                            valptr++;
+                            frobenius_update( nb, &scale, &sumsq, valptr );
 #endif
+                        }
+                    }
                 }
             }
             break;
@@ -139,11 +151,10 @@ z_spmMaxNorm( const pastix_spm_t *spm )
     pastix_complex64_t *valptr = (pastix_complex64_t*)spm->values;
     double tmp, norm = 0.;
 
-    for(i=0; i <spm->nnz; i++, valptr++) {
+    for(i=0; i < spm->nnzexp; i++, valptr++) {
         tmp = cabs( *valptr );
         norm = norm > tmp ? norm : tmp;
     }
-
     return norm;
 }
 
@@ -171,22 +182,37 @@ z_spmMaxNorm( const pastix_spm_t *spm )
 double
 z_spmInfNorm( const pastix_spm_t *spm )
 {
-    pastix_int_t col, row, i, baseval;
+    pastix_int_t col, row, i, j, k, ii, jj, baseval;
+    pastix_int_t dofi, dofj;
+    pastix_int_t *dofs = spm->dofs;
+
     pastix_complex64_t *valptr = (pastix_complex64_t*)spm->values;
     double norm = 0.;
     double *sumcol;
 
-    MALLOC_INTERN( sumcol, spm->gN, double );
-    memset( sumcol, 0, spm->gN * sizeof(double) );
+    MALLOC_INTERN( sumcol, spm->gNexp, double );
+    memset( sumcol, 0, spm->gNexp * sizeof(double) );
     baseval = spmFindBase( spm );
 
     switch( spm->fmttype ){
     case PastixCSC:
-        for( col=0; col < spm->gN; col++ )
+        for( i=0; i < spm->n; i++)
         {
-            for( i=spm->colptr[col]-baseval; i<spm->colptr[col+1]-baseval; i++ )
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            for(k=spm->colptr[i]; k < spm->colptr[i+1]; k++)
             {
-                sumcol[spm->rowptr[i]-baseval] += cabs( valptr[i] );
+                j = spm->rowptr[k - baseval] - baseval;
+                dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                row  = ( spm->dof > 0 ) ? j        : dofs[j];
+                for(ii=0; ii < dofi; ii++)
+                {
+                    for(jj=0; jj < dofj; jj++, valptr++)
+                    {
+                        {
+                            sumcol[row + jj] += cabs( *valptr );
+                        }
+                    }
+                }
             }
         }
 
@@ -194,11 +220,29 @@ z_spmInfNorm( const pastix_spm_t *spm )
         if ( (spm->mtxtype == PastixHermitian) ||
              (spm->mtxtype == PastixSymmetric) )
         {
-            for( col=0; col < spm->gN; col++ )
+            valptr = (pastix_complex64_t*)spm->values;
+            for(i=0; i < spm->n; i++)
             {
-                for( i=spm->colptr[col]-baseval+1; i<spm->colptr[col+1]-baseval; i++ )
+                col  = ( spm->dof > 0 ) ? i        : dofs[i];
+                dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+                for(k=spm->colptr[i]; k < spm->colptr[i+1]; k++)
                 {
-                    sumcol[col] += cabs( valptr[i] );
+                    j = spm->rowptr[k - baseval] - baseval;
+                    dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                    if(i != j)
+                    {
+                        for(ii=0; ii < dofi; ii++)
+                        {
+                            for(jj=0; jj < dofj; jj++, valptr++)
+                            {
+                                sumcol[col + ii] += cabs( *valptr );
+                            }
+                        }
+                    }
+                    else
+                    {
+                        valptr += dofi * dofj;
+                    }
                 }
             }
         }
@@ -250,7 +294,7 @@ z_spmInfNorm( const pastix_spm_t *spm )
         return PASTIX_ERR_BADPARAMETER;
     }
 
-    for( i=0; i<spm->gN; i++)
+    for( i=0; i<spm->gNexp; i++)
     {
         if(norm < sumcol[i])
         {
@@ -286,22 +330,37 @@ z_spmInfNorm( const pastix_spm_t *spm )
 double
 z_spmOneNorm( const pastix_spm_t *spm )
 {
-    pastix_int_t col, row, i, baseval;
+    pastix_int_t col, row, i, j, k, ii, jj, baseval;
+    pastix_int_t dofi, dofj;
+    pastix_int_t* dofs = spm->dofs;
+
     pastix_complex64_t *valptr = (pastix_complex64_t*)spm->values;
     double norm = 0.;
     double *sumrow;
 
-    MALLOC_INTERN( sumrow, spm->gN, double );
-    memset( sumrow, 0, spm->gN * sizeof(double) );
+    MALLOC_INTERN( sumrow, spm->gNexp, double );
+    memset( sumrow, 0, spm->gNexp * sizeof(double) );
     baseval = spmFindBase( spm );
 
     switch( spm->fmttype ){
     case PastixCSC:
-        for( col=0; col < spm->gN; col++ )
+        for(i=0; i < spm->n; i++)
         {
-            for( i=spm->colptr[col]-baseval; i<spm->colptr[col+1]-baseval; i++ )
+            col  = ( spm->dof > 0 ) ? i        : dofs[i];
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            for(k=spm->colptr[i]; k < spm->colptr[i+1]; k++)
             {
-                sumrow[col] += cabs( valptr[i] );
+                j = spm->rowptr[k - baseval] - baseval;
+                dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                for(ii=0; ii < dofi; ii++)
+                {
+                    for(jj=0; jj < dofj; jj++, valptr++)
+                    {
+                        {
+                            sumrow[col + ii] += cabs( *valptr );
+                        }
+                    }
+                }
             }
         }
 
@@ -309,11 +368,29 @@ z_spmOneNorm( const pastix_spm_t *spm )
         if ( (spm->mtxtype == PastixHermitian) ||
              (spm->mtxtype == PastixSymmetric) )
         {
-            for( col=0; col < spm->gN; col++ )
+            valptr = (pastix_complex64_t*)spm->values;
+            for(i=0; i < spm->n; i++)
             {
-                for( i=spm->colptr[col]-baseval+1; i<spm->colptr[col+1]-baseval; i++ )
+                dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+                for(k=spm->colptr[i]; k < spm->colptr[i+1]; k++)
                 {
-                    sumrow[spm->rowptr[i]-baseval] += cabs( valptr[i] );
+                    j = spm->rowptr[k - baseval] - baseval;
+                    row  = ( spm->dof > 0 ) ? j        : dofs[j];
+                    dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+                    if(i != j)
+                    {
+                        for(ii=0; ii < dofi; ii++)
+                        {
+                            for(jj=0; jj < dofj; jj++, valptr++)
+                            {
+                                sumrow[row + jj] += cabs( *valptr );
+                            }
+                        }
+                    }
+                    else
+                    {
+                        valptr += dofi * dofj;
+                    }
                 }
             }
         }
@@ -365,7 +442,7 @@ z_spmOneNorm( const pastix_spm_t *spm )
         return PASTIX_ERR_BADPARAMETER;
     }
 
-    for( i=0; i<spm->gN; i++)
+    for( i=0; i<spm->gNexp; i++)
     {
         if(norm < sumrow[i])
         {

z_spm_print.c

+/**
+ *
+ * @file z_spm_print.c
+ *
+ * Convert a sparse matrix into a dense matrix.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @author Alban Bellot
+ * @date 2015-01-01
+ *
+ * @precisions normal z -> c d s
+ *
+ **/
+#include <stdint.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include "pastix.h"
+#include "common.h"
+#include "spm.h"
+#include "z_spm.h"
+
+void
+z_spmPrint( const pastix_spm_t *spm )
+{
+    pastix_int_t i, j, baseval;
+    pastix_int_t k, ii, jj, dofi, dofj, col, row;
+    pastix_complex64_t *valptr;
+    pastix_int_t *colptr, *rowptr, *dofs;
+
+    assert( spm->fmttype == PastixCSC );
+    assert( spm->flttype == PastixComplex64 );
+
+    baseval = spmFindBase( spm );
+    i = 0; j = 0;
+
+    colptr = spm->colptr;
+    rowptr = spm->rowptr;
+    valptr = (pastix_complex64_t*)(spm->values);
+    dofs   = spm->dofs;
+
+    switch( spm->mtxtype ){
+#if defined(PRECISION_z) || defined(PRECISION_c)
+    case PastixHermitian:
+        for(i=0; i<spm->n; i++, colptr++)
+        {
+            dofi = ( spm->dof > 1 ) ? spm->dof : dofs[i+1] - dofs[i];
+            col = dofs[i];
+
+            for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
+            {
+                j = (*rowptr - baseval);
+                dofj = ( spm->dof > 1 ) ? spm->dof : dofs[j+1] - dofs[j];
+                row = dofs[j];
+
+                for(ii=0; ii<dofi; ii++)
+                {
+                    for(jj=0; jj<dofj; jj++, valptr++)
+                    {
+                        fprintf( stderr, "%ld %ld %e\n",
+                                 row + jj, col + ii, *valptr );
+                        if (i != j) {
+                            fprintf( stderr, "%ld %ld %e\n",
+                                     col + ii, row + jj, conj(*valptr) );
+                        }
+                    }
+                }
+            }
+        }
+        break;
+#endif
+    case PastixSymmetric:
+        for(i=0; i<spm->n; i++, colptr++)
+        {
+            dofi = ( spm->dof > 1 ) ? spm->dof : dofs[i+1] - dofs[i];
+            col = dofs[i];
+
+            for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
+            {
+                j = (*rowptr - baseval);
+                dofj = ( spm->dof > 1 ) ? spm->dof : dofs[j+1] - dofs[j];
+                row = dofs[j];
+
+                for(ii=0; ii<dofi; ii++)
+                {
+                    for(jj=0; jj<dofj; jj++, valptr++)
+                    {
+                        fprintf( stderr, "%ld %ld %e\n",
+                                 row + jj, col + ii, *valptr );
+                        if (i != j) {
+                            fprintf( stderr, "%ld %ld %e\n",
+                                     col + ii, row + jj, *valptr );
+                        }
+                    }
+                }
+            }
+        }
+        break;
+    case PastixGeneral:
+    default:
+        for(i=0; i<spm->n; i++, colptr++)
+        {
+            dofi = ( spm->dof > 1 ) ? spm->dof : dofs[i+1] - dofs[i];
+            col = dofs[i];
+
+            for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
+            {
+                j = (*rowptr - baseval);
+                dofj = ( spm->dof > 1 ) ? spm->dof : dofs[j+1] - dofs[j];
+                row = dofs[j];
+
+                for(ii=0; ii<dofi; ii++)
+                {
+                    for(jj=0; jj<dofj; jj++, valptr++)
+                    {
+                        fprintf( stderr, "%ld %ld %e\n",
+                                 row + jj, col + ii, *valptr );
+                    }
+                }
+            }
+        }
+    }
+    return;
+}
+