diff --git a/src/z_spm_dof_extend.c b/src/z_spm_dof_extend.c index 4ab207636bdaa71d227f39239dbefb6aab9358b7..0a77ab8d150d6f696db70da29cf8f9f34333cb6a 100644 --- a/src/z_spm_dof_extend.c +++ b/src/z_spm_dof_extend.c @@ -79,7 +79,7 @@ z_spm_dof_extend_update_values( spm_complex64_t *newval, * * @ingroup spm_dev_dof * - * @brief Extend a multi-dof CSX sparse matrix to a single dof sparse matrix. + * @brief Extend a single dof CSX sparse matrix to a multi-dof sparse matrix. * ******************************************************************************* * @@ -128,7 +128,7 @@ z_spm_dof_extend_csx( spmatrix_t *spm ) * * @ingroup spm_dev_dof * - * @brief Extend a multi-dof CSX sparse matrix to a single dof sparse matrix. + * @brief Extend a single dof IJV sparse matrix to a multi-dof sparse matrix. * ******************************************************************************* * @@ -172,7 +172,7 @@ z_spm_dof_extend_ijv( spmatrix_t *spm ) * * @ingroup spm_dev_dof * - * @brief Extend a multi-dof sparse matrix to a single dof sparse matrix. + * @brief Extend a single dof sparse matrix to a multi-dof sparse matrix. * ******************************************************************************* * diff --git a/src/z_spm_genrhs.c b/src/z_spm_genrhs.c index c09fc2f4224ea8c8615848c8bb696e15e5acb468..db4c626d1179cd1111dc637a4900f4f4955e699c 100644 --- a/src/z_spm_genrhs.c +++ b/src/z_spm_genrhs.c @@ -455,10 +455,6 @@ z_spmGenRHS( spm_rhstype_t type, int nrhs, return SPM_ERR_BADPARAMETER; } - /* if( spm->dof != 1 ) { */ - /* return SPM_ERR_BADPARAMETER; */ - /* } */ - if (nrhs == 1) { ldb = spm->nexp; ldx = spm->nexp; diff --git a/src/z_spm_norm.c b/src/z_spm_norm.c index c94daf14ad2cf1f6a55063d4691fd342982387c2..432f77009076e0007568394aca266256d7ed573b 100644 --- a/src/z_spm_norm.c +++ b/src/z_spm_norm.c @@ -621,6 +621,103 @@ z_spm_oneinf_elt( const spm_mtxtype_t mtxtype, } } +/** + * @brief Compute the one/inf norm of an spm CSX structure. + */ +static inline void +z_spmOneInfNorm_csx( spm_normtype_t ntype, + const spmatrix_t *spm, + double *sumtab, + spm_int_t baseval ) +{ + spm_int_t i, j, ig, jg, col, row; + spm_int_t dofi, dofj, dof; + spm_int_t *colptr, *rowptr, *loc2glob, *dofs; + spm_complex64_t *valptr; + + colptr = (spm->fmttype == SpmCSC) ? spm->colptr : spm->rowptr; + rowptr = (spm->fmttype == SpmCSC) ? spm->rowptr : spm->colptr; + valptr = (spm_complex64_t*)(spm->values); + loc2glob = spm->loc2glob; + dofs = spm->dofs; + dof = spm->dof; + for(j=0; j<spm->n; j++, colptr++, loc2glob++) + { + jg = (spm->loc2glob == NULL) ? j : (*loc2glob) - baseval; + if ( dof > 0 ) { + dofj = dof; + col = dof * jg; + } + else { + dofj = dofs[jg+1] - dofs[jg]; + col = dofs[jg] - baseval; + } + + for(i=colptr[0]; i<colptr[1]; i++, rowptr++) + { + ig = (*rowptr - baseval); + if ( dof > 0 ) { + dofi = dof; + row = dof * ig; + } + else { + dofi = dofs[ig+1] - dofs[ig]; + row = dofs[ig] - baseval; + } + + z_spm_oneinf_elt( spm->mtxtype, spm->layout, + row, dofi, col, dofj, valptr, + ntype, sumtab ); + valptr += dofi * dofj; + } + } +} + +/** + * @brief Compute the one/inf norm of an spm IJV structure. + */ +static inline void +z_spmOneInfNorm_ijv( spm_normtype_t ntype, + const spmatrix_t *spm, + double *sumtab, + spm_int_t baseval ) +{ + spm_int_t k, ig, jg, col, row; + spm_int_t dofi, dofj, dof; + spm_int_t *colptr, *rowptr, *dofs; + spm_complex64_t *valptr; + + colptr = spm->colptr; + rowptr = spm->rowptr; + valptr = (spm_complex64_t*)(spm->values); + dofs = spm->dofs; + dof = spm->dof; + + for(k=0; k<spm->nnz; k++, rowptr++, colptr++) + { + ig = *rowptr - baseval; + jg = *colptr - baseval; + + if ( dof > 0 ) { + dofi = dof; + row = dof * ig; + dofj = dof; + col = dof * jg; + } + else { + dofi = dofs[ig+1] - dofs[ig]; + row = dofs[ig] - baseval; + dofj = dofs[jg+1] - dofs[jg]; + col = dofs[jg] - baseval; + } + + z_spm_oneinf_elt( spm->mtxtype, spm->layout, + row, dofi, col, dofj, valptr, + ntype, sumtab ); + valptr += dofi * dofj; + } +} + /** ******************************************************************************* * @@ -647,120 +744,17 @@ static inline double z_spmOneInfNorm( spm_normtype_t ntype, const spmatrix_t *spm ) { - spm_int_t i, j, k, baseval; - spm_int_t ig, dofi, row; - spm_int_t jg, dofj, col; - const spm_int_t *colptr; - const spm_int_t *rowptr; - const spm_int_t *dofs; - const spm_int_t *loc2glob; - const spm_complex64_t *valptr; + spm_int_t k, baseval; double *sumtab = calloc( spm->gNexp, sizeof(double) ); double norm = 0.; baseval = spmFindBase( spm ); - colptr = spm->colptr; - rowptr = spm->rowptr; - valptr = (spm_complex64_t*)(spm->values); - dofs = spm->dofs; - loc2glob = spm->loc2glob; - - switch( spm->fmttype ){ - case SpmCSC: - for(j=0; j<spm->n; j++, colptr++, loc2glob++) - { - jg = (spm->loc2glob == NULL) ? j : (*loc2glob) - baseval; - if ( spm->dof > 0 ) { - dofj = spm->dof; - col = spm->dof * jg; - } - else { - dofj = dofs[jg+1] - dofs[jg]; - col = dofs[jg] - baseval; - } - - for(k=colptr[0]; k<colptr[1]; k++, rowptr++) - { - ig = (*rowptr - baseval); - if ( spm->dof > 0 ) { - dofi = spm->dof; - row = spm->dof * ig; - } - else { - dofi = dofs[ig+1] - dofs[ig]; - row = dofs[ig] - baseval; - } - - z_spm_oneinf_elt( spm->mtxtype, spm->layout, - row, dofi, col, dofj, valptr, - ntype, sumtab ); - valptr += dofi * dofj; - } - } - break; - - case SpmCSR: - for(i=0; i<spm->n; i++, rowptr++, loc2glob++) - { - ig = (spm->loc2glob == NULL) ? i : (*loc2glob) - baseval; - if ( spm->dof > 0 ) { - dofi = spm->dof; - row = spm->dof * ig; - } - else { - dofi = dofs[ig+1] - dofs[ig]; - row = dofs[ig] - baseval; - } - - for(k=rowptr[0]; k<rowptr[1]; k++, colptr++) - { - jg = (*colptr - baseval); - if ( spm->dof > 0 ) { - dofj = spm->dof; - col = spm->dof * jg; - } - else { - dofj = dofs[jg+1] - dofs[jg]; - col = dofs[jg] - baseval; - } - - z_spm_oneinf_elt( spm->mtxtype, spm->layout, - row, dofi, col, dofj, valptr, - ntype, sumtab ); - valptr += dofi * dofj; - } - } - break; - - case SpmIJV: - for(k=0; k<spm->nnz; k++, rowptr++, colptr++) - { - ig = *rowptr - baseval; - jg = *colptr - baseval; - - if ( spm->dof > 0 ) { - dofi = spm->dof; - row = spm->dof * ig; - dofj = spm->dof; - col = spm->dof * jg; - } - else { - dofi = dofs[ig+1] - dofs[ig]; - row = dofs[ig] - baseval; - dofj = dofs[jg+1] - dofs[jg]; - col = dofs[jg] - baseval; - } - - z_spm_oneinf_elt( spm->mtxtype, spm->layout, - row, dofi, col, dofj, valptr, - ntype, sumtab ); - valptr += dofi * dofj; - } - break; - - default: - return SPM_ERR_BADPARAMETER; + if( spm->fmttype != SpmIJV ){ + z_spmOneInfNorm_csx( ntype, spm, sumtab, baseval ); + } + else{ + z_spmOneInfNorm_ijv( ntype, spm, sumtab, baseval ); } #if defined(SPM_WITH_MPI) @@ -772,7 +766,7 @@ z_spmOneInfNorm( spm_normtype_t ntype, /* Look for the maximum */ { const double *sumtmp = sumtab; - for( i=0; i<spm->gNexp; i++, sumtmp++ ) + for( k=0; k<spm->gNexp; k++, sumtmp++ ) { if( norm < *sumtmp ) { norm = *sumtmp;