Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 53289a3e authored by Alban Bellot's avatar Alban Bellot
Browse files

matvec test with dofs

parent 99beb4d7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 2
− 0
View file @ 53289a3e
......@@ -21,6 +21,8 @@ set(SOURCES
z_spm_genrhs.c
z_spm_matrixvector.c
z_spm_norm.c
z_spm_2dense.c
z_spm_print.c
)
set(spm_headers
......
z_spm.h
+ 4
− 0
View file @ 53289a3e
......@@ -37,4 +37,8 @@ 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 );
#endif /* _z_spm_H_ */
z_spm_convert_to_csc.c
+ 1
− 0
View file @ 53289a3e
......@@ -204,6 +204,7 @@ z_spmConvertCSR2CSC( pastix_spm_t *spm )
assert( val_csc );
#endif
/* Count the number of elements per column */
for (j=0; j<nnz; j++) {
col = spm->colptr[j] - baseval;
......
z_spm_matrixvector.c
+ 91
− 25
View file @ 53289a3e
......@@ -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=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=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=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=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 )
{
......@@ -194,7 +217,10 @@ 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 +236,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 = 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 = 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 +319,10 @@ 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,17 +336,43 @@ 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( i=0; i < csc->gN; i++ )
{
dofi = ( csc->dof > 0 ) ? csc->dof : dofs[i+1] - dofs[i];
col = dofs[i];
for( k=csc->colptr[i]; k<csc->colptr[i+1]; k++ )
{
j = csc->rowptr[k-baseval]-baseval;
dofj = ( csc->dof > 0 ) ? csc->dof : dofs[j+1] - dofs[j];
row = 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];
}
}
}
}
}
}
/*
if( alpha != 0. ) {
for( col=0; col < csc->gN; col++ )
{
......@@ -320,6 +385,7 @@ z_spmHeCSCv( pastix_complex64_t alpha,
}
}
}
*/
return PASTIX_SUCCESS;
}
......
z_spm_norm.c
+ 47
− 26
View file @ 53289a3e
......@@ -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,31 @@ z_spmFrobeniusNorm( const pastix_spm_t *spm )
else {
pastix_int_t *colptr = spm->colptr;
pastix_int_t *rowptr = spm->rowptr;
int nb;
int nb, dofi, dofj, ii, jj;
pastix_int_t *dofs=spm->dofs;
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;
......@@ -211,7 +222,7 @@ z_spmInfNorm( const pastix_spm_t *spm )
*/
/* Dofs */
for( i=0; i < spm->n; i++ )
for( i=0; i < spm->n; i++)
{
dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
for(k=spm->colptr[i]; k < spm->colptr[i+1]; k++)
......@@ -230,7 +241,6 @@ z_spmInfNorm( const pastix_spm_t *spm )
}
}
}
/* Add the symmetric/hermitian part */
if ( (spm->mtxtype == PastixHermitian) ||
(spm->mtxtype == PastixSymmetric) )
......@@ -244,16 +254,20 @@ z_spmInfNorm( const pastix_spm_t *spm )
{
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++)
if(i != j)
{
for(ii=0; ii < dofi; ii++)
{
for(jj=0; jj < dofj; jj++, valptr++)
{
if(i != j)
{
sumcol[col + ii] += cabs( *valptr );
}
sumcol[col + ii] += cabs( *valptr );
}
}
}
else
{
valptr += dofi * dofj;
}
}
}
}
......@@ -368,18 +382,20 @@ z_spmOneNorm( const pastix_spm_t *spm )
for(i=0; i < spm->n; i++)
{
col = dofs[i];
dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - spm->dofs[i];
dofi = ( spm->dof > 0 ) ? spm->dof : 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] - spm->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++, valptr++)
{
{
sumrow[col + ii] += cabs( *valptr );
}
}
}
}
}
/* Add the symmetric/hermitian part */
......@@ -389,22 +405,27 @@ z_spmOneNorm( const pastix_spm_t *spm )
valptr = (pastix_complex64_t*)spm->values;
for(i=0; i < spm->n; i++)
{
dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - spm->dofs[i];
dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
for(k=spm->colptr[i]; k < spm->colptr[i+1]; k++)
{
j = spm->rowptr[k - baseval] - baseval;
row = dofs[j];
dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - spm->dofs[j];
for(ii=0; ii < dofi; ii++)
for(jj=0; jj < dofj; jj++, valptr++)
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++)
{
if(i != j)
{
sumrow[row + jj] += cabs( *valptr );
}
sumrow[row + jj] += cabs( *valptr );
}
}
}
else
{
valptr += dofi * dofj;
}
}
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment