From 6007657a5e0a245c1b162a6a4075a9eab7e60c4a Mon Sep 17 00:00:00 2001
From: Alban Bellot <alban.bellot@inria.fr>
Date: Fri, 29 Jul 2016 17:14:17 +0200
Subject: [PATCH] Multi-dof fin de stage
---
spm.c | 2 +
spm.h | 2 +-
spm_expand.c | 207 ++++++++++++++++++++++++++++++
z_spm.h | 7 ++
z_spm_convert_col_row_major.c | 220 ++++++++++++++++++++++++++++++++
z_spm_convert_to_csc.c | 4 +-
z_spm_convert_to_csr.c | 4 +-
z_spm_convert_to_ijv.c | 11 +-
z_spm_dofs2flat.c | 229 +++++++++++++++++++++++++---------
z_spm_expand.c | 46 +++++++
10 files changed, 664 insertions(+), 68 deletions(-)
create mode 100644 spm_expand.c
create mode 100644 z_spm_convert_col_row_major.c
create mode 100644 z_spm_expand.c
diff --git a/spm.c b/spm.c
index 1ab5af08..f5d7bfb3 100644
--- a/spm.c
+++ b/spm.c
@@ -813,3 +813,5 @@ spmCheckAxb( int nrhs,
return ptrfunc[id](nrhs, spm, x0, ldx0, b, ldb, x, ldx );
}
}
+
+
diff --git a/spm.h b/spm.h
index 46b80c3a..7d0d63db 100644
--- a/spm.h
+++ b/spm.h
@@ -41,7 +41,7 @@ struct pastix_spm_s {
if > 0, constant degree of freedom
otherwise, irregular degree of freedom (refer to dofs) */
pastix_int_t *dofs; /*< Number of degrees of freedom per unknown (NULL, if dof > 0) */
- int colmajor; /*< If > 0, column major when PastixIJV with dofs
+ int colmajor; /*< If > 0, column major with dofs
otherwise, row major */
pastix_int_t *colptr; /*< List of indirections to rows for each vertex */
diff --git a/spm_expand.c b/spm_expand.c
new file mode 100644
index 00000000..6dead986
--- /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 73b31693..886f39e9 100644
--- a/z_spm.h
+++ b/z_spm.h
@@ -41,4 +41,11 @@ pastix_complex64_t *z_spm2dense( const pastix_spm_t *spm );
void z_spmDensePrint( pastix_int_t m, pastix_int_t n, pastix_complex64_t *A, pastix_int_t lda );
void z_spmPrint( const pastix_spm_t *spm );
+void z_spmConvertColMaj2RowMaj(pastix_spm_t *spm);
+void z_spmConvertRowMaj2ColMaj(pastix_spm_t *spm);
+
+void z_extandCSC(pastix_spm_t *spm);
+
+void z_spmDofs2Flat(pastix_spm_t *spm);
+
#endif /* _z_spm_H_ */
diff --git a/z_spm_convert_col_row_major.c b/z_spm_convert_col_row_major.c
new file mode 100644
index 00000000..e13b28ff
--- /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 435ecc50..5c535651 100644
--- a/z_spm_convert_to_csc.c
+++ b/z_spm_convert_to_csc.c
@@ -162,7 +162,7 @@ z_spmConvertIJV2CSC( pastix_spm_t *spm )
spmExit( &oldspm );
spm->fmttype = PastixCSC;
-
+ spm->colmajor = 1;
return PASTIX_SUCCESS;
}
@@ -354,6 +354,6 @@ z_spmConvertCSR2CSC( pastix_spm_t *spm )
if(spm-> dof != 1)
spm->mtxtype = type;
-
+ spm->colmajor = 1;
return PASTIX_SUCCESS;
}
diff --git a/z_spm_convert_to_csr.c b/z_spm_convert_to_csr.c
index a3eca787..71d3cf2c 100644
--- a/z_spm_convert_to_csr.c
+++ b/z_spm_convert_to_csr.c
@@ -97,6 +97,7 @@ z_spmConvertCSC2CSR( pastix_spm_t *spm )
if(spm-> dof != 1)
spm->mtxtype = type;
+ spm->colmajor = -1;
return result;
}
@@ -144,10 +145,10 @@ z_spmConvertIJV2CSR( pastix_spm_t *spm )
baseval = spmFindBase( spm );
#if !defined(PRECISION_p)
+ pastix_int_t ii, jj, k;
/* Transpose values in row major format */
if( spm->colmajor )
{
- pastix_int_t ii, jj, k;
pastix_int_t cpt=0;
pastix_int_t* dofs = spm->dofs;
@@ -278,6 +279,7 @@ z_spmConvertIJV2CSR( pastix_spm_t *spm )
spmExit( &oldspm );
+ spm->colmajor = -1;
spm->fmttype = PastixCSR;
return PASTIX_SUCCESS;
diff --git a/z_spm_convert_to_ijv.c b/z_spm_convert_to_ijv.c
index 4c518b98..4f2a6017 100644
--- a/z_spm_convert_to_ijv.c
+++ b/z_spm_convert_to_ijv.c
@@ -64,7 +64,8 @@ z_spmConvertCSC2IJV( pastix_spm_t *spm )
}
/* Transpose values in row major format */
- if( !spm->colmajor ) //A test
+ /*
+ if( spm->colmajor < 1 ) //A test
{
int k, ii, jj, dofi, dofj;
int cpt = 0;
@@ -86,11 +87,13 @@ z_spmConvertCSC2IJV( pastix_spm_t *spm )
oavals++;
}
}
- cpt + = dofi * dofj;
+ cpt += dofi * dofj;
}
spm->values = navals;
}
+ */
+
memFree_null(spm->colptr);
spm->colptr = col_ijv;
@@ -143,7 +146,8 @@ z_spmConvertCSR2IJV( pastix_spm_t *spm )
}
/* Transpose values in column major format */
- if( spm->colmajor )
+ /*
+ if( spm->colmajor > 1 )
{
pastix_int_t k, ii, jj, dofi, dofj;
pastix_int_t cpt=0;
@@ -168,6 +172,7 @@ z_spmConvertCSR2IJV( pastix_spm_t *spm )
}
spm->values = navals;
}
+ */
memFree_null(spm->rowptr);
spm->rowptr = row_ijv;
diff --git a/z_spm_dofs2flat.c b/z_spm_dofs2flat.c
index 42b9805a..b4db59f6 100644
--- a/z_spm_dofs2flat.c
+++ b/z_spm_dofs2flat.c
@@ -44,83 +44,190 @@
void
z_spmDofs2Flat(pastix_spm_t *spm)
{
- pastix_int_t i, j, k, ii, jj, dofi, dofj, col, row, baseval;
+ pastix_int_t i, j, k, ii, jj, dofi, dofj, col, row, baseval, cpt;
baseval = spmFindBase( spm );
pastix_int_t *new_col = calloc(spm->nexp+1,sizeof(pastix_int_t));
- pastix_int_t *new_row = malloc(sizeof(pastix_int_t)*spm->nnzexp);
+ pastix_int_t *new_row;
pastix_int_t *dofs = spm->dofs;
- pastix_complex64_t *new_vals = malloc(sizeof(pastix_int_t)*spm->nnzexp);
+ pastix_complex64_t *new_vals;
pastix_complex64_t *vals = (pastix_complex64_t*)spm->values;
- for(i=0;i<spm->n ; i++)
+ switch(spm->mtxtype)
{
- col = ( spm->dof > 0 ) ? i : dofs[i];
- for(k=spm->colptr[i]-baseval; k<spm->colptr[i+1]-baseval; k++)
- {
- j = spm->rowptr[k]-baseval;
- row = ( spm->dof > 0 ) ? j : dofs[j];
+ case PastixGeneral:
+ new_row = malloc(sizeof(pastix_int_t)*spm->nnzexp);
+ new_vals = malloc(sizeof(pastix_complex64_t)*spm->nnzexp);
+ for(i=0;i<spm->n ; i++)
+ {
+ col = ( spm->dof > 0 ) ? i : dofs[i];
dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
- dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
- for(ii=0; ii<dofi; ii++)
- {
- new_col[col+ii+1] += dofj;
- }
- }
- }
+ for(k=spm->colptr[i]-baseval; k<spm->colptr[i+1]-baseval; k++)
+ {
+ j = spm->rowptr[k]-baseval;
+ row = ( spm->dof > 0 ) ? j : dofs[j];
+ dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+ for(ii=0; ii<dofi; ii++)
+ {
+ new_col[col+ii+1] += dofj;
+ }
+ }
+ }
- for(i=0; i<spm->nexp; i++)
- {
- new_col[i+1]+=new_col[i];
- }
+ for(i=0; i<spm->nexp; i++)
+ {
+ new_col[i+1]+=new_col[i];
+ }
- int cpt = 0;
- for(i=0; i < spm->n;i++)
- {
- col = ( spm->dof > 0 ) ? i : dofs[i];
- dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
- for(k=spm->colptr[i]-baseval ; k<spm->colptr[i+1]-baseval ;k++)
- {
- j = spm->rowptr[k] - baseval;
- row = ( spm->dof > 0 ) ? j : dofs[j];
- dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
- for(ii=0;ii < dofi; ii++)
- {
- for(jj=0;jj < dofj ; jj++)
- {
- new_vals[new_col[col+ii]] = vals[cpt];
- new_row[new_col[col+ii]] = row + jj + baseval;
- cpt++;
- new_col[col+ii]++;
- }
- }
- }
- }
+ cpt = 0;
+ for(i=0; i < spm->n;i++)
+ {
+ col = ( spm->dof > 0 ) ? i : dofs[i];
+ dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+ for(k=spm->colptr[i]-baseval ; k<spm->colptr[i+1]-baseval ;k++)
+ {
+ j = spm->rowptr[k] - baseval;
+ row = ( spm->dof > 0 ) ? j : dofs[j];
+ dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+ for(ii=0;ii < dofi; ii++)
+ {
+ for(jj=0;jj < dofj ; jj++)
+ {
+ new_vals[new_col[col+ii]] = vals[cpt];
+ new_row[new_col[col+ii]] = row + jj + baseval;
+ new_col[col+ii]++;
+ cpt++;
+ }
+ }
+ }
+ }
- {
- int tmp;
- int tmp1 = 0;
+ {
+ int tmp;
+ int tmp1 = 0;
+ for(i=0; i<spm->nexp; i++)
+ {
+ tmp = new_col[i];
+ new_col[i] = tmp1+baseval;
+ tmp1 = tmp;
+ }
+ new_col[i] += baseval;
+ }
+ spm->gN = spm->gNexp;
+ spm->n = spm->nexp;
+ spm->gnnz = spm->gnnzexp;
+ spm->nnz = spm->nnzexp;
+
+ spm->dof = 1;
+ spm->dofs = NULL;
+ spm->colmajor = 1;
+
+ spm->colptr = new_col;
+ spm->rowptr = new_row;
+ //spm->loc2glob = NULL; // ?
+ spm->values = new_vals;
+ break;
+
+ case PastixSymmetric:
+ for(i=0;i<spm->n ; i++)
+ {
+ col = ( spm->dof > 0 ) ? i : dofs[i];
+ dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+ for(k=spm->colptr[i]-baseval; k<spm->colptr[i+1]-baseval; k++)
+ {
+ j = spm->rowptr[k]-baseval;
+ row = ( spm->dof > 0 ) ? j : dofs[j];
+ dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+ for(ii=0; ii<dofi; ii++)
+ {
+ for(jj=0; jj<dofj; jj++)
+ {
+ if( i != j )
+ new_col[col+ii+1] += 1;
+ else
+ if(ii <= jj )
+ new_col[col+ii+1] += 1;
+ }
+ }
+ }
+ }
for(i=0; i<spm->nexp; i++)
{
- tmp = new_col[i];
- new_col[i] = tmp1+baseval;
- tmp1 = tmp;
+ new_col[i+1] += new_col[i];
}
- new_col[i] += baseval;
- }
+ pastix_int_t nnz = new_col[spm->nexp];
+ new_row = malloc(sizeof(pastix_int_t)*nnz);
+ new_vals = malloc(sizeof(pastix_complex64_t)*nnz);
- spm->gN = spm->gNexp;
- spm->n = spm->nexp;
- spm->gnnz = spm->gnnzexp;
- spm->nnz = spm->nnzexp;
+ cpt = 0;
+ for(i=0; i < spm->n;i++)
+ {
+ col = ( spm->dof > 0 ) ? i : dofs[i];
+ dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+ for(k=spm->colptr[i]-baseval ; k<spm->colptr[i+1]-baseval ;k++)
+ {
+ j = spm->rowptr[k] - baseval;
+ row = ( spm->dof > 0 ) ? j : dofs[j];
+ dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+ for(ii=0;ii < dofi; ii++)
+ {
+ for(jj=0;jj < dofj ; jj++)
+ {
+ if( i == j )
+ {
+ if ( ii <= jj )
+ {
+ /* diagonal dominant for spd matrix
+ if( ii == jj)
+ new_vals[new_col[col+ii]] = 2*vals[cpt];
+ else
+ */
+ new_vals[new_col[col+ii]] = vals[cpt];
+ new_row[new_col[col+ii]] = row + jj + baseval;
+ new_col[col+ii]++;
+ }
+ }
+ else
+ {
+ new_vals[new_col[col+ii]] = vals[cpt];
+ new_row[new_col[col+ii]] = row + jj + baseval;
+ new_col[col+ii]++;
- spm->dof = 1;
- spm->dofs = NULL;
- spm->colmajor = 1;
+ }
+ cpt++;
+ }
+ }
+ }
+ }
+ {
+ int tmp;
+ int tmp1 = 0;
+ for(i=0; i<spm->nexp; i++)
+ {
+ tmp = new_col[i];
+ new_col[i] = tmp1+baseval;
+ tmp1 = tmp;
+ }
+ new_col[i] += baseval;
+ }
+ spm->gN = spm->gNexp;
+ spm->n = spm->nexp;
+ spm->gnnz = nnz;
+ spm->nnz = nnz;
+ spm->gnnzexp = nnz;
+ spm->nnzexp = nnz;
- spm->colptr = new_col;
- spm->rowptr = new_row;
- spm->loc2glob = NULL; // ?
- spm->values = new_vals;
+ spm->dof = 1;
+ spm->dofs = NULL;
+ spm->colmajor = 1;
+
+ spm->colptr = new_col;
+ spm->rowptr = new_row;
+ //spm->loc2glob = NULL; //
+ spm->values = new_vals;
+
+ break;
+ }
+ assert(spm->loc2glob == NULL);//to do
}
diff --git a/z_spm_expand.c b/z_spm_expand.c
new file mode 100644
index 00000000..8167b0c5
--- /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;
+}
--
GitLab