diff --git a/src/spm.c b/src/spm.c index 82a7fff45721e5b5dddfdf1e778f209874bd7337..5028953c13b5423e55aa857b42a861bbeceb3ea5 100644 --- a/src/spm.c +++ b/src/spm.c @@ -165,36 +165,19 @@ spmInit( spmatrix_t *spm ) void spmAlloc( spmatrix_t *spm ) { - spm_int_t colsize, rowsize, valsize, dofsize; - - switch(spm->fmttype){ - case SpmCSC: - colsize = spm->n + 1; - rowsize = spm->nnz; - valsize = spm->nnzexp; - dofsize = spm->gN + 1; - break; - case SpmCSR: - colsize = spm->nnz; - rowsize = spm->n + 1; - valsize = spm->nnzexp; - dofsize = spm->gN + 1; - break; - case SpmIJV: - default: - colsize = spm->nnz; - rowsize = spm->nnz; - valsize = spm->nnzexp; - dofsize = spm->gN + 1; - } + spm_int_t colsize = (spm->fmttype == SpmCSC) ? spm->n + 1 : spm->nnz; + spm_int_t rowsize = (spm->fmttype == SpmCSR) ? spm->n + 1 : spm->nnz; spm->colptr = (spm_int_t*)malloc( colsize * sizeof(spm_int_t) ); spm->rowptr = (spm_int_t*)malloc( rowsize * sizeof(spm_int_t) ); + if ( spm->dof < 1 ) { + spm_int_t dofsize = spm->gN + 1; spm->dofs = (spm_int_t*)malloc( dofsize * sizeof(spm_int_t) ); } + if(spm->flttype != SpmPattern) { - valsize = valsize * spm_size_of( spm->flttype ); + spm_int_t valsize = spm->nnzexp * spm_size_of( spm->flttype ); spm->values = malloc(valsize); } } @@ -260,7 +243,7 @@ spmBase( spmatrix_t *spm, int baseval ) { spm_int_t baseadj; - spm_int_t i, n, nnz; + spm_int_t i, n, nnz, colsize, rowsize; /* Parameter checks */ if ( spm == NULL ) { @@ -281,39 +264,20 @@ spmBase( spmatrix_t *spm, } baseadj = baseval - spmFindBase( spm ); - if (baseadj == 0) + if (baseadj == 0) { return; + } - n = spm->n; - nnz = spm->nnz; - - switch(spm->fmttype) - { - case SpmCSC: - assert( nnz == (spm->colptr[n] - spm->colptr[0]) ); - - for (i = 0; i <= n; i++) { - spm->colptr[i] += baseadj; - } - for (i = 0; i < nnz; i++) { - spm->rowptr[i] += baseadj; - } - break; + n = spm->n; + nnz = spm->nnz; + colsize = (spm->fmttype == SpmCSC) ? n + 1 : nnz; + rowsize = (spm->fmttype == SpmCSR) ? n + 1 : nnz; - case SpmCSR: - assert( nnz == (spm->rowptr[n] - spm->rowptr[0]) ); - for (i = 0; i <= n; i++) { - spm->rowptr[i] += baseadj; - } - for (i = 0; i < nnz; i++) { - spm->colptr[i] += baseadj; - } - break; - case SpmIJV: - for (i = 0; i < nnz; i++) { - spm->rowptr[i] += baseadj; - spm->colptr[i] += baseadj; - } + for (i = 0; i < colsize; i++) { + spm->colptr[i] += baseadj; + } + for (i = 0; i < rowsize; i++) { + spm->rowptr[i] += baseadj; } if (spm->loc2glob != NULL) { @@ -508,10 +472,6 @@ spmNorm( spm_normtype_t ntype, { double norm = -1.; - if ( spm->flttype == SpmPattern ) { - return norm; - } - switch (spm->flttype) { case SpmFloat: norm = (double)s_spmNorm( ntype, spm ); @@ -531,7 +491,7 @@ spmNorm( spm_normtype_t ntype, case SpmPattern: default: - ; + return norm; } return norm; @@ -807,26 +767,10 @@ spmCopy( const spmatrix_t *spm ) memcpy( newspm, spm, sizeof(spmatrix_t)); - switch(spm->fmttype){ - case SpmCSC: - colsize = spm->n + 1; - rowsize = spm->nnz; - valsize = spm->nnzexp; - dofsize = spm->gN + 1; - break; - case SpmCSR: - colsize = spm->nnz; - rowsize = spm->n + 1; - valsize = spm->nnzexp; - dofsize = spm->gN + 1; - break; - case SpmIJV: - default: - colsize = spm->nnz; - rowsize = spm->nnz; - valsize = spm->nnzexp; - dofsize = spm->gN + 1; - } + colsize = (spm->fmttype == SpmCSC) ? spm->n + 1 : spm->nnz; + rowsize = (spm->fmttype == SpmCSR) ? spm->n + 1 : spm->nnz; + valsize = spm->nnzexp; + dofsize = spm->gN + 1; if(spm->colptr != NULL) { newspm->colptr = (spm_int_t*)malloc( colsize * sizeof(spm_int_t) ); @@ -1038,6 +982,7 @@ spmPrintRHS( const spmatrix_t *spm, z_spmPrintRHS( stream, spm, nrhs, x, ldx ); break; case SpmDouble: + default: d_spmPrintRHS( stream, spm, nrhs, x, ldx ); } @@ -1145,10 +1090,6 @@ spmMatVec( spm_trans_t trans, return SPM_ERR_BADPARAMETER; } - if ( spm->dof != 1 ) { - espm = malloc( sizeof(spmatrix_t) ); - spmExpand( spm, espm ); - } switch (spm->flttype) { case SpmFloat: rc = spm_sspmv( trans, alpha, espm, x, 1, beta, y, 1 ); @@ -1631,50 +1572,50 @@ spm_get_distribution( const spmatrix_t *spm ) { int distribution = 0; + /* The matrix is not distributed */ if( (spm->loc2glob == NULL) || (spm->n == spm->gN) ) { distribution = ( SpmDistByColumn | SpmDistByRow ); + return distribution; + } + if( spm->fmttype == SpmCSC ){ + distribution = SpmDistByColumn; + } + else if ( spm->fmttype == SpmCSR ) { + distribution = SpmDistByRow; } else { - if( spm->fmttype == SpmCSC ){ - distribution = SpmDistByColumn; - } - else if ( spm->fmttype == SpmCSR ) { - distribution = SpmDistByRow; - } - else { - spm_int_t i, baseval; - spm_int_t *colptr = spm->colptr; - spm_int_t *glob2loc = spm->glob2loc; - - baseval = spmFindBase( spm ); - distribution = 1; - assert( glob2loc != NULL ); - for ( i = 0; i < spm->nnz; i++, colptr++ ) - { - /* - * If the global index is not in the local colptr - * -> row distribution - */ - if( glob2loc[ *colptr - baseval ] < 0 ) { - distribution = SpmDistByRow; - break; - } + spm_int_t i, baseval; + spm_int_t *colptr = spm->colptr; + spm_int_t *glob2loc = spm->glob2loc; + + baseval = spmFindBase( spm ); + distribution = 1; + assert( glob2loc != NULL ); + for ( i = 0; i < spm->nnz; i++, colptr++ ) + { + /* + * If the global index is not in the local colptr + * -> row distribution + */ + if( glob2loc[ *colptr - baseval ] < 0 ) { + distribution = SpmDistByRow; + break; } + } - #if defined(SPM_WITH_MPI) - { - int check = 0; - MPI_Allreduce( &distribution, &check, 1, MPI_INT, - MPI_BOR, spm->comm ); - /* - * If a matrix is distributed - * it cannot be distributed by row AND column - */ - assert( check != ( SpmDistByColumn | SpmDistByRow ) ); - assert( distribution == check ); - } - #endif +#if defined(SPM_WITH_MPI) + { + int check = 0; + MPI_Allreduce( &distribution, &check, 1, MPI_INT, + MPI_BOR, spm->comm ); + /* + * If a matrix is distributed + * it cannot be distributed by row AND column + */ + assert( check != ( SpmDistByColumn | SpmDistByRow ) ); + assert( distribution == check ); } +#endif } assert(distribution > 0); return distribution; diff --git a/src/spm_scatter.c b/src/spm_scatter.c index 1d59b0162baa9540778e902734d9708cc9d27de7..7077e46db05f597424845fbb5f7d70c43770994c 100644 --- a/src/spm_scatter.c +++ b/src/spm_scatter.c @@ -231,7 +231,7 @@ spm_scatter_csx_get_locals( const spmatrix_t *oldspm, root, newspm->comm ); } - /* Syore the local information */ + /* Store the local information */ newspm->nnz = allcounts[ newspm->clustnum * 3 + 1 ]; newspm->nnzexp = allcounts[ newspm->clustnum * 3 + 2 ]; @@ -339,7 +339,7 @@ spm_scatter_ijv_get_locals( const spmatrix_t *oldspm, root, newspm->comm ); } - /* Syore the local information */ + /* Store the local information */ newspm->nnz = allcounts[ newspm->clustnum * 3 + 1 ]; newspm->nnzexp = allcounts[ newspm->clustnum * 3 + 2 ]; diff --git a/src/z_spm_dof_extend.c b/src/z_spm_dof_extend.c index 0a77ab8d150d6f696db70da29cf8f9f34333cb6a..23d75bcc007aedabe725cfcfc93ba65b7abc1656 100644 --- a/src/z_spm_dof_extend.c +++ b/src/z_spm_dof_extend.c @@ -62,7 +62,6 @@ z_spm_dof_extend_update_values( spm_complex64_t *newval, } } } - else { for(jj=0; jj<dofj; jj++) { @@ -189,6 +188,4 @@ z_spmDofExtend( spmatrix_t *spm ) else { z_spm_dof_extend_ijv( spm ); } - - return; } diff --git a/src/z_spm_norm.c b/src/z_spm_norm.c index 432f77009076e0007568394aca266256d7ed573b..72b1997ae857f55d88812f0085690c43e89563f3 100644 --- a/src/z_spm_norm.c +++ b/src/z_spm_norm.c @@ -355,6 +355,7 @@ z_spmFrobeniusNorm( const spmatrix_t *spm ) break; case SpmIJV: + default: z_spmFrobeniusNorm_ijv( spm, data ); } } @@ -622,10 +623,10 @@ z_spm_oneinf_elt( const spm_mtxtype_t mtxtype, } /** - * @brief Compute the one/inf norm of an spm CSX structure. + * @brief Compute the one/inf norm of an spm CSC structure. */ static inline void -z_spmOneInfNorm_csx( spm_normtype_t ntype, +z_spmOneInfNorm_csc( spm_normtype_t ntype, const spmatrix_t *spm, double *sumtab, spm_int_t baseval ) @@ -635,8 +636,8 @@ z_spmOneInfNorm_csx( spm_normtype_t ntype, 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; + colptr = spm->colptr; + rowptr = spm->rowptr; valptr = (spm_complex64_t*)(spm->values); loc2glob = spm->loc2glob; dofs = spm->dofs; @@ -673,6 +674,58 @@ z_spmOneInfNorm_csx( spm_normtype_t ntype, } } +/** + * @brief Compute the one/inf norm of an spm CSR structure. + */ +static inline void +z_spmOneInfNorm_csr( 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->colptr; + rowptr = spm->rowptr; + valptr = (spm_complex64_t*)(spm->values); + loc2glob = spm->loc2glob; + dofs = spm->dofs; + dof = spm->dof; + for(i=0; i<spm->n; i++, rowptr++, loc2glob++) + { + ig = (spm->loc2glob == NULL) ? i : (*loc2glob) - baseval; + if ( dof > 0 ) { + dofi = dof; + row = dof * ig; + } + else { + dofi = dofs[ig+1] - dofs[ig]; + row = dofs[ig] - baseval; + } + + for(j=rowptr[0]; j<rowptr[1]; j++, colptr++) + { + jg = (*colptr - baseval); + if ( dof > 0 ) { + dofj = dof; + col = 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; + } + } +} + /** * @brief Compute the one/inf norm of an spm IJV structure. */ @@ -750,11 +803,19 @@ z_spmOneInfNorm( spm_normtype_t ntype, baseval = spmFindBase( spm ); - if( spm->fmttype != SpmIJV ){ - z_spmOneInfNorm_csx( ntype, spm, sumtab, baseval ); - } - else{ - z_spmOneInfNorm_ijv( ntype, spm, sumtab, baseval ); + switch (spm->fmttype) + { + case SpmCSC: + z_spmOneInfNorm_csc( ntype, spm, sumtab, baseval ); + break; + + case SpmCSR: + z_spmOneInfNorm_csr( ntype, spm, sumtab, baseval ); + break; + + case SpmIJV: + default: + z_spmOneInfNorm_ijv( ntype, spm, sumtab, baseval ); } #if defined(SPM_WITH_MPI)