diff --git a/include/spm.h b/include/spm.h
index 9d4ac2f90958503c5f63129762f8891ba4e3de3c..4c7f9ff5bd2e724a5acafa1db5cc0564d62140c2 100644
--- a/include/spm.h
+++ b/include/spm.h
@@ -172,7 +172,8 @@ int spmParseLaplacianInfo( const char * filename,
spm_int_t * dim2,
spm_int_t * dim3,
double * alpha,
- double * beta );
+ double * beta,
+ spm_int_t * dof );
/**
* @}
diff --git a/src/drivers/laplacian.c b/src/drivers/laplacian.c
index e2845506069b1fd6e1beb27d4baca063723a4b05..e280d2bd10320f642c4e7367aaf7b1f7e9644957 100644
--- a/src/drivers/laplacian.c
+++ b/src/drivers/laplacian.c
@@ -85,6 +85,9 @@ laplacian_usage(void)
* @param[out] beta
* The beta coefficient for the adjacency matrix
*
+ * @param[out] dof
+ * The degree of freedom of each unknown
+ *
*******************************************************************************
*
* @retval SPM_SUCCESS if the matrix has been generated successfully
@@ -98,13 +101,15 @@ spmParseLaplacianInfo( const char *filename,
spm_int_t *dim2,
spm_int_t *dim3,
double *alpha,
- double *beta )
+ double *beta,
+ spm_int_t *dof )
{
double val1, val2;
- long tmp1, tmp2, tmp3;
+ long tmp1, tmp2, tmp3, tmp4;
*alpha = 1.;
*beta = 1.;
+ *dof = 1;
/* Look for the datatype */
{
@@ -173,7 +178,16 @@ spmParseLaplacianInfo( const char *filename,
/* Scan the dimensions */
*dim1 = *dim2 = *dim3 = 1;
- if ( sscanf( filename, "%ld:%ld:%ld:%lf:%lf", &tmp1, &tmp2, &tmp3, &val1, &val2 ) == 5 ) {
+ if ( sscanf( filename, "%ld:%ld:%ld:%lf:%lf:%ld",
+ &tmp1, &tmp2, &tmp3, &val1, &val2, &tmp4 ) == 6 ) {
+ *dim1 = (spm_int_t)tmp1;
+ *dim2 = (spm_int_t)tmp2;
+ *dim3 = (spm_int_t)tmp3;
+ *alpha = val1;
+ *beta = val2;
+ *dof = (spm_int_t)tmp4;
+ }
+ else if ( sscanf( filename, "%ld:%ld:%ld:%lf:%lf", &tmp1, &tmp2, &tmp3, &val1, &val2 ) == 5 ) {
*dim1 = (spm_int_t)tmp1;
*dim2 = (spm_int_t)tmp2;
*dim3 = (spm_int_t)tmp3;
@@ -203,6 +217,8 @@ spmParseLaplacianInfo( const char *filename,
return SPM_ERR_BADPARAMETER;
}
+ assert( *dof != 0 );
+
/* One of the dimension was set to 0 */
if ( (*dim1 == 0) || (*dim2 == 0) || (*dim3 == 0) ) {
laplacian_usage();
@@ -269,20 +285,34 @@ genLaplacian( const char *filename,
spmatrix_t *spm )
{
spm_coeftype_t flttype;
- spm_int_t dim1, dim2, dim3;
+ spm_int_t dim1, dim2, dim3, dof;
double alpha = 1.;
double beta = 1.;
int rc;
- rc = spmParseLaplacianInfo(filename, &flttype, &dim1, &dim2, &dim3, &alpha, &beta );
+ rc = spmParseLaplacianInfo(filename, &flttype, &dim1, &dim2, &dim3, &alpha, &beta, &dof );
if (rc != SPM_SUCCESS)
return rc;
spm->flttype = flttype;
spm->n = dim1 * dim2 * dim3;
+ spm->dof = dof;
laplacian_7points[spm->flttype](spm, dim1, dim2, dim3, alpha, beta);
+ if ( dof != 1 ) {
+ spmatrix_t *spmtmp;
+ spmUpdateComputedFields( spm );
+ if ( dof < 1 ) {
+ spmtmp = spmDofExtend( spm, 1, -dof );
+ }
+ else {
+ spmtmp = spmDofExtend( spm, 0, dof );
+ }
+ spmExit( spm );
+ memcpy( spm, spmtmp, sizeof(spmatrix_t) );
+ }
+
return SPM_SUCCESS;
}
@@ -323,19 +353,33 @@ genExtendedLaplacian( const char *filename,
spmatrix_t *spm )
{
spm_coeftype_t flttype;
- spm_int_t dim1, dim2, dim3;
+ spm_int_t dim1, dim2, dim3, dof;
double alpha = 1.;
double beta = 1.;
int rc;
- rc = spmParseLaplacianInfo(filename, &flttype, &dim1, &dim2, &dim3, &alpha, &beta);
+ rc = spmParseLaplacianInfo(filename, &flttype, &dim1, &dim2, &dim3, &alpha, &beta, &dof );
if (rc != SPM_SUCCESS)
return rc;
spm->flttype = flttype;
spm->n = dim1 * dim2 * dim3;
+ spm->dof = dof;
laplacian_27points[spm->flttype](spm, dim1, dim2, dim3, alpha, beta);
+ if ( dof != 1 ) {
+ spmatrix_t *spmtmp;
+ spmUpdateComputedFields( spm );
+ if ( dof < 1 ) {
+ spmtmp = spmDofExtend( spm, 1, -dof );
+ }
+ else {
+ spmtmp = spmDofExtend( spm, 0, dof );
+ }
+ spmExit( spm );
+ memcpy( spm, spmtmp, sizeof(spmatrix_t) );
+ }
+
return SPM_SUCCESS;
}
diff --git a/src/drivers/readijv.c b/src/drivers/readijv.c
index e7ae19170ef04f6eb76de05113ad26debe7fdcce..a587788fa32b4e0fe08ba96df7f794b5e701fd63 100644
--- a/src/drivers/readijv.c
+++ b/src/drivers/readijv.c
@@ -103,7 +103,7 @@ readIJV( const char *dirname,
double *tempval;
spm_int_t i, Nrow, Ncol, Nnzero;
- filename = malloc(strlen(dirname)+10);
+ filename = malloc(strlen(dirname)+20);
spm->flttype = SpmDouble;
spm->mtxtype = SpmGeneral;
@@ -113,7 +113,7 @@ readIJV( const char *dirname,
/* Read the header information */
{
- sprintf(filename,"%s/header",dirname);
+ sprintf( filename, "%s/header", dirname );
hdrfile = fopen (filename,"r");
if (hdrfile == NULL)
{
@@ -134,7 +134,7 @@ readIJV( const char *dirname,
spm->values = (double *) malloc(Nnzero*sizeof(double));
/* Open the 3 files */
- sprintf(filename,"%s/ia_threeFiles",dirname);
+ sprintf( filename, "%s/ia_threeFiles", dirname );
iafile = fopen(filename,"r");
if (iafile == NULL)
{
@@ -143,7 +143,7 @@ readIJV( const char *dirname,
return SPM_ERR_BADPARAMETER;
}
- sprintf(filename,"%s/ja_threeFiles",dirname);
+ sprintf( filename, "%s/ja_threeFiles", dirname );
jafile = fopen(filename,"r");
if (jafile == NULL)
{
@@ -153,7 +153,7 @@ readIJV( const char *dirname,
return SPM_ERR_BADPARAMETER;
}
- sprintf(filename,"%s/ra_threeFiles",dirname);
+ sprintf( filename, "%s/ra_threeFiles", dirname );
rafile = fopen(filename,"r");
if (rafile == NULL)
{
diff --git a/src/spm.c b/src/spm.c
index fe4800ede8a9a9f6e08b9babc3c728f587850b28..25f5ab3bbb60d963f021c4d46e0ae7a15cbef5cb 100644
--- a/src/spm.c
+++ b/src/spm.c
@@ -128,7 +128,6 @@ spmInit( spmatrix_t *spm )
void
spmUpdateComputedFields( spmatrix_t *spm )
{
-
/*
* Compute the local expended field for multi-dofs
*/
@@ -399,6 +398,8 @@ spmFindBase( const spmatrix_t *spm )
if ( ( baseval != 0 ) &&
( baseval != 1 ) )
{
+ assert( spm->fmttype == SpmIJV );
+
baseval = spm->n;
tmp = spm->colptr;
for(i=0; i<spm->nnz; i++, tmp++){
@@ -594,8 +595,8 @@ int
spmSort( spmatrix_t *spm )
{
if ( (spm->dof != 1) && (spm->flttype != SpmPattern) ) {
- fprintf(stderr, "WARNING: spm expanded due to non implemented sort for non-expanded spm with values\n");
- spm = spmExpand( spm );
+ assert( 0 );
+ fprintf(stderr, "ERROR: spmSort should not be called with non expanded matrices including values\n");
}
switch (spm->flttype) {
case SpmPattern:
@@ -644,9 +645,9 @@ spmSort( spmatrix_t *spm )
spm_int_t
spmMergeDuplicate( spmatrix_t *spm )
{
- if ( (spm->dof != 1) && (spm->flttype != SpmPattern) ) {
- fprintf(stderr, "WARNING: spm expanded due to non implemented merge for non-expanded spm with values\n");
- spm = spmExpand( spm );
+ if ( (spm->dof < 1) && (spm->flttype != SpmPattern) ) {
+ assert( 0 );
+ fprintf(stderr, "Error: spmMergeDuplicate should not be called with non expanded matrices with variadic degrees of freedom and values\n" );
}
switch (spm->flttype) {
case SpmPattern:
@@ -695,8 +696,8 @@ spm_int_t
spmSymmetrize( spmatrix_t *spm )
{
if ( (spm->dof != 1) && (spm->flttype != SpmPattern) ) {
- fprintf(stderr, "WARNING: spm expanded due to non implemented symmetrize for non-expanded spm with values\n");
- spm = spmExpand( spm );
+ assert( 0 );
+ fprintf(stderr, "ERROR: spmSymmetrize should not be called with non expanded matrices including values\n");
}
switch (spm->flttype) {
case SpmPattern:
@@ -755,14 +756,14 @@ spmCheckAndCorrect( const spmatrix_t *spm_in,
/* Let's work on a copy */
newspm = spmCopy( spm_in );
- /* PaStiX works on CSC matrices */
- spmConvert( SpmCSC, newspm );
-
if ( (newspm->dof != 1) && (newspm->flttype != SpmPattern) ) {
- fprintf(stderr, "WARNING: newspm expanded due to missing check functions implementations\n");
+ fprintf(stderr, "spmCheckAndCorrect: spm is expanded due to multiple degrees of freedom\n");
newspm = spmExpand( newspm );
}
+ /* PaStiX works on CSC matrices */
+ spmConvert( SpmCSC, newspm );
+
/* Sort the rowptr for each column */
spmSort( newspm );
@@ -858,19 +859,19 @@ spmCopy( const spmatrix_t *spm )
}
if(spm->colptr != NULL) {
- newspm->colptr = (spm_int_t*)malloc( colsize * sizeof(spm_int_t));
- memcpy( newspm->colptr, spm->colptr, colsize * sizeof(spm_int_t));
+ newspm->colptr = (spm_int_t*)malloc( colsize * sizeof(spm_int_t) );
+ memcpy( newspm->colptr, spm->colptr, colsize * sizeof(spm_int_t) );
}
if(spm->rowptr != NULL) {
- newspm->rowptr = (spm_int_t*)malloc(rowsize * sizeof(spm_int_t));
- memcpy( newspm->rowptr, spm->rowptr, rowsize * sizeof(spm_int_t));
+ newspm->rowptr = (spm_int_t*)malloc( rowsize * sizeof(spm_int_t) );
+ memcpy( newspm->rowptr, spm->rowptr, rowsize * sizeof(spm_int_t) );
}
if(spm->loc2glob != NULL) {
- newspm->loc2glob = (spm_int_t*)malloc(dofsize * sizeof(spm_int_t));
- memcpy( newspm->loc2glob, spm->loc2glob, dofsize * sizeof(spm_int_t));
+ newspm->loc2glob = (spm_int_t*)malloc( dofsize * sizeof(spm_int_t) );
+ memcpy( newspm->loc2glob, spm->loc2glob, dofsize * sizeof(spm_int_t) );
}
if(spm->dofs != NULL) {
- newspm->dofs = (spm_int_t*)malloc(dofsize * sizeof(spm_int_t));
+ newspm->dofs = (spm_int_t*)malloc( dofsize * sizeof(spm_int_t) );
memcpy( newspm->dofs, spm->dofs, dofsize * sizeof(spm_int_t) );
}
if(spm->values != NULL) {
diff --git a/src/z_spm.c b/src/z_spm.c
index 9ce44ee2fa4fc58f9e7d269629cd950d8f6208cc..28cabf49e028f69da046360ebee318bbc492bc48 100644
--- a/src/z_spm.c
+++ b/src/z_spm.c
@@ -138,17 +138,23 @@ z_spmMergeDuplicate( spmatrix_t *spm )
spm_int_t n = spm->n;
spm_int_t baseval = spm->colptr[0];
spm_int_t dof2 = spm->dof * spm->dof;
- spm_int_t idx, i, j, size;
+ spm_int_t idx, i, j, size, savedcolptr;
spm_int_t merge = 0;
#if !defined(PRECISION_p)
spm_int_t d;
#endif
+ assert( spm->dof >= 1 );
+ assert( spm->fmttype == SpmCSC );
+
if ( spm->fmttype == SpmCSC ) {
- idx = 0;
+ idx = baseval;
+ savedcolptr = colptr[0];
for (i=0; i<n; i++, colptr++)
{
- size = colptr[1] - colptr[0];
+ size = colptr[1] - savedcolptr;
+ savedcolptr = colptr[1];
+
for (j = 0; j < size;
j++, oldrow++, oldval+=dof2, newrow++, newval+=dof2, idx++)
{
@@ -172,13 +178,13 @@ z_spmMergeDuplicate( spmatrix_t *spm )
merge++;
}
}
- assert( ((merge == 0) && ( colptr[1] == idx+baseval)) ||
- ((merge != 0) && ( colptr[1] > idx+baseval)) );
+ assert( ( (merge == 0) && (colptr[1] == idx) ) ||
+ ( (merge != 0) && (colptr[1] > idx) ) );
- colptr[1] = idx + baseval;
+ colptr[1] = idx;
}
- assert( ((merge == 0) && (spm->nnz == idx)) ||
- ((merge != 0) && (spm->nnz - merge == idx)) );
+ assert( ((merge == 0) && (spm->nnz == (idx-baseval))) ||
+ ((merge != 0) && (spm->nnz - merge == (idx-baseval))) );
if (merge > 0) {
spm->nnz = spm->nnz - merge;
diff --git a/src/z_spm_dof_extend.c b/src/z_spm_dof_extend.c
index 2a397a5bcb652f82f80a7a8b7d96d6ffea1d7fef..b51c32ba6e3f59ee1b1bc9c5cbbb2517a035994b 100644
--- a/src/z_spm_dof_extend.c
+++ b/src/z_spm_dof_extend.c
@@ -31,7 +31,7 @@
*
*******************************************************************************/
void
-z_spmDofExtend(spmatrix_t *spm)
+z_spmDofExtend( spmatrix_t *spm )
{
spm_int_t i, j, k, ii, jj, dofi, dofj, baseval;
spm_int_t *colptr, *rowptr, *dofs;
@@ -74,38 +74,17 @@ z_spmDofExtend(spmatrix_t *spm)
{
for(ii=0; ii<dofi; ii++, newval++)
{
- *newval = *oldval;
+ if ( i == j ) {
+ *newval = *oldval / (abs(ii - jj) + 1.);
+ }
+ else {
+ *newval = *oldval;
+ }
}
}
}
}
break;
- /* case SpmCSR: */
- /* /\** */
- /* * 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 SpmIJV:
/**
* Loop on coordinates
@@ -121,7 +100,12 @@ z_spmDofExtend(spmatrix_t *spm)
{
for(ii=0; ii<dofi; ii++, newval++)
{
- *newval = *oldval;
+ if ( i == j ) {
+ *newval = *oldval / (abs(ii - jj) + 1.);
+ }
+ else {
+ *newval = *oldval;
+ }
}
}
}
diff --git a/src/z_spm_norm.c b/src/z_spm_norm.c
index 14203318ecf5acf464ed649da05f06de29bc5db4..861a4f536b1042f44c139ad9699f2c577387db01 100644
--- a/src/z_spm_norm.c
+++ b/src/z_spm_norm.c
@@ -340,8 +340,9 @@ z_spmOneNorm( const spmatrix_t *spm )
{
for(i=0; i < spm->nnz; i++)
{
- if(spm->rowptr[i] != spm->colptr[i])
+ if( spm->rowptr[i] != spm->colptr[i] ) {
sumcol[spm->rowptr[i]-baseval] += cabs( valptr[i] );
+ }
}
}
break;
diff --git a/tests/spm_tests.h b/tests/spm_tests.h
index 7d79984f5407899dbc2dee1fa8343c398aed8630..624841317a0badff845e083fd4f3111b6bb52057 100644
--- a/tests/spm_tests.h
+++ b/tests/spm_tests.h
@@ -95,7 +95,7 @@ spm_norm_print_result( double norms, double normd, double result )
printf("SUCCESS !\n");
} else {
printf("FAILED !\n");
- rc=1;
+ rc = 1;
}
printf(" Nsparse = %e, Ndense = %e\n", norms, normd );
diff --git a/tests/z_spm_tests.c b/tests/z_spm_tests.c
index 17f0bb8c3d32cc552a617026bd60310d21fe8448..511c0d5c3bf75022187762c3e10500a9d8a9efee 100644
--- a/tests/z_spm_tests.c
+++ b/tests/z_spm_tests.c
@@ -207,18 +207,23 @@ z_spm_norm_check( const spmatrix_t *spm )
{
spm_complex64_t *A;
double norms, normd;
- double eps, result;
+ double eps, result, nnz;
int ret = 0;
eps = LAPACKE_dlamch_work('e');
+ nnz = spm->gnnzexp;
+ if ( spm->mtxtype != SpmGeneral ) {
+ nnz *= 2.;
+ }
+
/* Create a dense backup of spm */
A = z_spm2dense( spm );
/**
* Test Norm Max
*/
- printf(" -- Test norm Max :");
+ printf(" -- Test norm Max : ");
norms = spmNorm( SpmMaxNorm, spm );
normd = LAPACKE_zlange( LAPACK_COL_MAJOR, 'M', spm->gNexp, spm->gNexp, A, spm->gNexp );
result = fabs(norms - normd) / (normd * eps);
@@ -227,31 +232,31 @@ z_spm_norm_check( const spmatrix_t *spm )
/**
* Test Norm Inf
*/
- printf(" -- Test norm Inf :");
+ printf(" -- Test norm Inf : ");
norms = spmNorm( SpmInfNorm, spm );
normd = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->gNexp, spm->gNexp, A, spm->gNexp );
result = fabs(norms - normd) / (normd * eps);
- result = result * ((double)(spm->gNexp)) / ((double)(spm->gnnzexp));
+ result = result * ((double)(spm->gNexp)) / nnz;
ret += spm_norm_print_result( norms, normd, result );
/**
* Test Norm One
*/
- printf(" -- Test norm One :");
+ printf(" -- Test norm One : ");
norms = spmNorm( SpmOneNorm, spm );
normd = LAPACKE_zlange( LAPACK_COL_MAJOR, 'O', spm->gNexp, spm->gNexp, A, spm->gNexp );
result = fabs(norms - normd) / (normd * eps);
- result = result * ((double)(spm->gNexp)) / ((double)(spm->gnnzexp));
+ result = result * ((double)(spm->gNexp)) / nnz;
ret += spm_norm_print_result( norms, normd, result );
/**
* Test Norm Frobenius
*/
- printf(" -- Test norm Frb :");
+ printf(" -- Test norm Frb : ");
norms = spmNorm( SpmFrobeniusNorm, spm );
normd = LAPACKE_zlange( LAPACK_COL_MAJOR, 'F', spm->gNexp, spm->gNexp, A, spm->gNexp );
result = fabs(norms - normd) / (normd * eps);
- result = result / ((double)spm->gnnzexp);
+ result = result / nnz;
ret += spm_norm_print_result( norms, normd, result );
free(A);
diff --git a/wrappers/fortran90/src/spmf.f90 b/wrappers/fortran90/src/spmf.f90
index ac0ce6d69b98cb3e65e68debaccacbe0810f8c62..e5d4cb2923ad094f694dd8e454f3ceb50dfd9c1f 100644
--- a/wrappers/fortran90/src/spmf.f90
+++ b/wrappers/fortran90/src/spmf.f90
@@ -344,7 +344,7 @@ module spmf
interface
function spmParseLaplacianInfo_c(filename, flttype, dim1, dim2, dim3, &
- alpha, beta) &
+ alpha, beta, dof) &
bind(c, name='spmParseLaplacianInfo')
use iso_c_binding
import spm_int_t
@@ -357,6 +357,7 @@ module spmf
type(c_ptr), value :: dim3
type(c_ptr), value :: alpha
type(c_ptr), value :: beta
+ type(c_ptr), value :: dof
end function spmParseLaplacianInfo_c
end interface
@@ -667,7 +668,7 @@ contains
end subroutine spmReadDriver
subroutine spmParseLaplacianInfo(filename, flttype, dim1, dim2, dim3, alpha, &
- beta, info)
+ beta, dof, info)
use iso_c_binding
implicit none
character(kind=c_char), intent(in), target :: filename
@@ -677,10 +678,12 @@ contains
integer(kind=spm_int_t), intent(inout), target :: dim3
real(kind=c_double), intent(inout), target :: alpha
real(kind=c_double), intent(inout), target :: beta
+ integer(kind=spm_int_t), intent(inout), target :: dof
integer(kind=c_int), intent(out) :: info
info = spmParseLaplacianInfo_c(c_loc(filename), c_loc(flttype), &
- c_loc(dim1), c_loc(dim2), c_loc(dim3), c_loc(alpha), c_loc(beta))
+ c_loc(dim1), c_loc(dim2), c_loc(dim3), c_loc(alpha), c_loc(beta), &
+ c_loc(dof))
end subroutine spmParseLaplacianInfo
subroutine spm2Dense(spm)
diff --git a/wrappers/python/spm/__spm__.py b/wrappers/python/spm/__spm__.py
index 97148fc180fda384e3528f2f8a1384f8ed4c93a1..d72fef8676cb197e9826047b2a61915911622227 100644
--- a/wrappers/python/spm/__spm__.py
+++ b/wrappers/python/spm/__spm__.py
@@ -165,16 +165,17 @@ def pyspm_spmReadDriver( driver, filename, spm ):
return libspm.spmReadDriver( driver, filename, spm )
def pyspm_spmParseLaplacianInfo( filename, flttype, dim1, dim2, dim3, alpha,
- beta ):
+ beta, dof ):
libspm.spmParseLaplacianInfo.argtypes = [ c_char_p, c_int_p,
POINTER(__spm_int__),
POINTER(__spm_int__),
POINTER(__spm_int__),
POINTER(c_double),
- POINTER(c_double) ]
+ POINTER(c_double),
+ POINTER(__spm_int__) ]
libspm.spmParseLaplacianInfo.restype = c_int
return libspm.spmParseLaplacianInfo( filename, flttype, dim1, dim2, dim3,
- alpha, beta )
+ alpha, beta, dof )
def pyspm_spm2Dense( spm ):
libspm.spm2Dense.argtypes = [ POINTER(pyspm_spmatrix_t) ]