diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..21caac4ee5dcc3c861e917b2b25660fd20cd08aa
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,55 @@
+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}
+ )
diff --git a/spm.c b/spm.c
index c08db2fd553d1468c401f00eb0fceeef2733ba29..7e6b362bc61c2402bdf64e2142f3827b58bca9bc 100644
--- a/spm.c
+++ b/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.
*
diff --git a/spm.h b/spm.h
index 18fd91d2bccc60f8233e298cf976a90ad6b8c0b1..71fd46984fa014bdee1eb995297d11adb207c9b2 100644
--- a/spm.h
+++ b/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_ */
diff --git a/spm_expand.c b/spm_expand.c
new file mode 100644
index 0000000000000000000000000000000000000000..6dead986a5be667ed3d5cbd0327960fd32b7e4e5
--- /dev/null
+++ b/spm_expand.c
@@ -0,0 +1,207 @@
+/**
+ *
+ * @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);
+ */
+}
+
+
diff --git a/z_spm.h b/z_spm.h
index 877da184bbc1a4e9a638d9e3f3d20889513373bb..886f39e90f5bc44a54728afc376c50d684e13d27 100644
--- a/z_spm.h
+++ b/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_ */
diff --git a/z_spm_2dense.c b/z_spm_2dense.c
new file mode 100644
index 0000000000000000000000000000000000000000..ff3dddfa9e7cc8ffc8ca0ec5cf708a44aa6643fd
--- /dev/null
+++ b/z_spm_2dense.c
@@ -0,0 +1,191 @@
+/**
+ *
+ * @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;
+}
+
diff --git a/z_spm_convert_col_row_major.c b/z_spm_convert_col_row_major.c
new file mode 100644
index 0000000000000000000000000000000000000000..e13b28ff81b1ca64faf352a9e0d164a445e44270
--- /dev/null
+++ b/z_spm_convert_col_row_major.c
@@ -0,0 +1,220 @@
+/**
+ *
+ * @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;
+ }
+}
diff --git a/z_spm_convert_to_csc.c b/z_spm_convert_to_csc.c
index bd76019cc490bee69ae093d8693344a48c6ab963..5ec2702ca00c13cb679e5f63ec68d1d6940e60d5 100644
--- a/z_spm_convert_to_csc.c
+++ b/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;
}
diff --git a/z_spm_convert_to_csr.c b/z_spm_convert_to_csr.c
index c1d0d4582f0d2edd50c74a231cb6a1419dbfcd50..73822a86722ff523d42c41e67fcbef29b10ddf28 100644
--- a/z_spm_convert_to_csr.c
+++ b/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;
}
+
diff --git a/z_spm_convert_to_ijv.c b/z_spm_convert_to_ijv.c
index 6ab25d095da59c6dc59bcf576e9afcc1c5544472..457f353169a6c87eb36c90e28883bc378c3185a5 100644
--- a/z_spm_convert_to_ijv.c
+++ b/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;
diff --git a/z_spm_dofs2flat.c b/z_spm_dofs2flat.c
new file mode 100644
index 0000000000000000000000000000000000000000..b4db59f69be193c6f9eaadbef6fcfe78e1bf8ce5
--- /dev/null
+++ b/z_spm_dofs2flat.c
@@ -0,0 +1,233 @@
+/**
+ *
+ * @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
+}
diff --git a/z_spm_expand.c b/z_spm_expand.c
new file mode 100644
index 0000000000000000000000000000000000000000..8167b0c537ae154e49c418a2dcb9c14253769042
--- /dev/null
+++ b/z_spm_expand.c
@@ -0,0 +1,46 @@
+/**
+ *
+ * @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;
+}
diff --git a/z_spm_matrixvector.c b/z_spm_matrixvector.c
index e5346ebf9c67c735626e2a43394b9f99133dcc5d..93734bbdf80399195ff8f1eabfc1c4c2fe893b6f 100644
--- a/z_spm_matrixvector.c
+++ b/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
diff --git a/z_spm_norm.c b/z_spm_norm.c
index 422ace0b4cb31e7eeb46a987f731ebdb3f340622..60880281fb44b686916caeab192301c3da756ae4 100644
--- a/z_spm_norm.c
+++ b/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])
{
diff --git a/z_spm_print.c b/z_spm_print.c
new file mode 100644
index 0000000000000000000000000000000000000000..6e4f99298b45d6ba5d4188550fc4c76932e18ce5
--- /dev/null
+++ b/z_spm_print.c
@@ -0,0 +1,129 @@
+/**
+ *
+ * @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;
+}
+