diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8faf7143e8abbf9d543019c3231114ed375098e1..e0293612cdf9027ace07b05a24a6bef62874ac98 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -23,6 +23,7 @@ set(SOURCES
z_spm_convert_to_ijv.c
z_spm_dofs2flat.c
z_spm_expand.c
+ z_spm_dof_extend.c
z_spm_genrhs.c
z_spm_integer.c
z_spm_laplacian.c
@@ -47,7 +48,7 @@ set(spm_sources
spm.c
spm_io.c
spm_integers.c
- spm_expand.c
+ spm_dofs.c
spm_read_driver.c
drivers/skitf.f
drivers/iohb.c
diff --git a/spm.h b/spm.h
index d46a446dcb800e5c6c169b6c6b05a0b63623501a..c91de942ecb26fcdb04b5f172215bd223e01a682 100644
--- a/spm.h
+++ b/spm.h
@@ -79,12 +79,13 @@ struct pastix_spm_s {
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) */
+ pastix_int_t *dofs; /**< Array of the first column of each element in the
+ expanded matrix [+baseval] */
pastix_order_t layout; /**< PastixColMajor, or PastixRowMajor */
- 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 */
+ pastix_int_t *colptr; /**< List of indirections to rows for each vertex [+baseval] */
+ pastix_int_t *rowptr; /**< List of edges for each vertex [+baseval] */
+ pastix_int_t *loc2glob;/**< Corresponding numbering from local to global [+baseval] */
void *values; /**< Values stored in the matrix */
};
diff --git a/spm_dofs.c b/spm_dofs.c
new file mode 100644
index 0000000000000000000000000000000000000000..1732564327a6717b2c7784883dcf10d5dca71483
--- /dev/null
+++ b/spm_dofs.c
@@ -0,0 +1,147 @@
+/**
+ *
+ * @file spm_dofs.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"
+#include "c_spm.h"
+#include "d_spm.h"
+#include "s_spm.h"
+#include "p_spm.h"
+
+
+pastix_spm_t *
+spmDofExtend( const int type,
+ const int dof,
+ const pastix_spm_t *spm )
+{
+ pastix_spm_t *newspm;
+
+ /* Quick return */
+ if ( dof == 1 )
+ return (pastix_spm_t *)spm;
+
+ if ( spm->dof != 1 ) {
+ pastix_error_print( "Cannot extend spm including dofs already\n" );
+ return (pastix_spm_t *)spm;
+ }
+
+ newspm = spmCopy( spm );
+
+ /**
+ * Generate constant dof
+ */
+ if (type == 0) {
+ newspm->dof = dof;
+ newspm->nexp = spm->n * dof;
+ newspm->gNexp = spm->gN * dof;
+ newspm->nnzexp = spm->nnz * dof * dof;
+ newspm->gnnzexp = spm->gnnz * dof * dof;
+ }
+ else {
+ pastix_int_t i, j, k, dofi, dofj, baseval;
+ pastix_int_t *dofptr, *colptr, *rowptr;
+
+ baseval = spmFindBase( spm );
+
+ newspm->dof = -1;
+ newspm->dofs = malloc( (spm->n+1) * sizeof(pastix_int_t) );
+ dofptr = newspm->dofs;
+
+ /**
+ * Initialize the dofs array where the degree of freedom of vertex i is
+ * dof[i+1] - dof[i]
+ */
+ *dofptr = 0;
+ for(i=0; i<spm->n; i++, dofptr++) {
+ dofi = 1 + ( rand() % dof );
+ dofptr[1] = dofptr[0] + dofi;
+ }
+
+ newspm->nexp = *dofptr;
+ newspm->gNexp = *dofptr;
+
+ /**
+ * Count the number of non zeroes
+ */
+ newspm->nnzexp = 0;
+ colptr = newspm->colptr;
+ rowptr = newspm->rowptr;
+ dofptr = newspm->dofs;
+
+ switch(spm->fmttype)
+ {
+ case PastixCSC:
+ for(j=0; j<newspm->n; j++, colptr++) {
+ dofj = dofptr[j+1] - dofptr[j];
+
+ for(k=colptr[0]; k<colptr[1]; k++, rowptr++) {
+ i = *rowptr - baseval;
+ dofi = dofptr[i+1] - dofptr[i];
+
+ newspm->nnzexp += dofi * dofj;
+ }
+ }
+ break;
+ case PastixCSR:
+ for(i=0; i<newspm->n; i++, rowptr++) {
+ dofi = dofptr[i+1] - dofptr[i];
+
+ for(k=rowptr[0]; k<rowptr[1]; k++, rowptr++) {
+ j = *colptr - baseval;
+ dofj = dofptr[j+1] - dofptr[j];
+
+ newspm->nnzexp += dofi * dofj;
+ }
+ }
+ break;
+ case PastixIJV:
+ for(k=0; k<newspm->nnz; k++, rowptr++, colptr++)
+ {
+ i = *rowptr - baseval;
+ j = *colptr - baseval;
+ dofi = dofptr[i+1] - dofptr[i];
+ dofj = dofptr[j+1] - dofptr[j];
+
+ newspm->nnzexp += dofi * dofj;
+ }
+ }
+ newspm->gnnzexp = newspm->nnzexp;
+ }
+
+ switch (spm->flttype) {
+ case PastixFloat:
+ s_spmDofExtend( newspm );
+ break;
+
+ case PastixDouble:
+ d_spmDofExtend( newspm );
+ break;
+
+ case PastixComplex32:
+ c_spmDofExtend( newspm );
+ break;
+
+ case PastixComplex64:
+ z_spmDofExtend( newspm );
+ break;
+
+ case PastixPattern:
+ ;
+ }
+
+ return newspm;
+}
diff --git a/spm_expand.c b/spm_expand.c
deleted file mode 100644
index 2957feb995540c162e4958092eff76b6fcc495e3..0000000000000000000000000000000000000000
--- a/spm_expand.c
+++ /dev/null
@@ -1,191 +0,0 @@
-/**
- *
- * @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"
-#include "c_spm.h"
-#include "d_spm.h"
-#include "s_spm.h"
-#include "p_spm.h"
-
-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;
-
- spmExpand(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;
- spmExpand(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);
- */
-}
-
-
diff --git a/z_spm.h b/z_spm.h
index c846eeca1380e82be614cdc8ef9a5eeecf8f324d..fc8321ddc6e3cd1e7ddb6adbf2521bf242a7c5f0 100644
--- a/z_spm.h
+++ b/z_spm.h
@@ -66,6 +66,6 @@ void z_spmPrint( const pastix_spm_t *spm );
pastix_spm_t *z_spmExpand(const pastix_spm_t *spm);
-void z_spmDofs2Flat(pastix_spm_t *spm);
+void z_spmDofExtend(pastix_spm_t *spm);
#endif /* _z_spm_H_ */
diff --git a/z_spm_dof_extend.c b/z_spm_dof_extend.c
new file mode 100644
index 0000000000000000000000000000000000000000..c9e2e0c317bb5ba45d2c48483e3b10f55fc0fdd7
--- /dev/null
+++ b/z_spm_dof_extend.c
@@ -0,0 +1,111 @@
+/**
+ *
+ * @file z_spm_dof_extend.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 Alban Bellot
+ * @date 2015-01-01
+ *
+ * @precisions normal z -> c d s
+ **/
+#include "common.h"
+#include "spm.h"
+#include "z_spm.h"
+
+void
+z_spmDofExtend(pastix_spm_t *spm)
+{
+ pastix_int_t i, j, k, ii, jj, dofi, dofj, baseval;
+ pastix_int_t *colptr, *rowptr, *dofs;
+ pastix_complex64_t *newval, *oldval;
+
+ oldval = (pastix_complex64_t*)(spm->values);
+ newval = spm->values = malloc( spm->nnzexp * sizeof(pastix_complex64_t) );
+
+ baseval = spmFindBase( spm );
+ colptr = spm->colptr;
+ rowptr = spm->rowptr;
+ dofs = spm->dofs;
+
+ switch(spm->fmttype)
+ {
+ case PastixCSC:
+ /**
+ * Loop on col
+ */
+ for(j=0; j<spm->n; j++, colptr++)
+ {
+ dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+
+ /**
+ * Loop on rows
+ */
+ for(k=colptr[0]; k<colptr[1]; k++, rowptr++, oldval++)
+ {
+ i = *rowptr - baseval;
+ dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+
+ for(jj=0; jj<dofj; jj++)
+ {
+ for(ii=0; ii<dofi; ii++, newval++)
+ {
+ *newval = *oldval;
+ }
+ }
+ }
+ }
+ break;
+ case PastixCSR:
+ /**
+ * Loop on row
+ */
+ for(i=0; i<spm->n; i++, rowptr++)
+ {
+ dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+
+ /**
+ * Loop on cols
+ */
+ for(k=rowptr[0]; k<rowptr[1]; k++, colptr++, oldval++)
+ {
+ j = *colptr - baseval;
+ dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+
+ for(jj=0; jj<dofj; jj++)
+ {
+ for(ii=0; ii<dofi; ii++, newval++)
+ {
+ *newval = *oldval;
+ }
+ }
+ }
+ }
+ break;
+ case PastixIJV:
+ /**
+ * Loop on coordinates
+ */
+ for(k=0; k<spm->nnz; k++, rowptr++, colptr++, oldval++)
+ {
+ i = *rowptr - baseval;
+ j = *colptr - 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++, newval++)
+ {
+ *newval = *oldval;
+ }
+ }
+ }
+ break;
+ }
+ return;
+}