diff --git a/z_spm_genrhs.c b/z_spm_genrhs.c index 0876c435c6633856c3a682d77f4f63713679190e..a9b9b83328f2088972771862a88cd634ae31a319 100644 --- a/z_spm_genrhs.c +++ b/z_spm_genrhs.c @@ -209,23 +209,29 @@ z_spmGenRHS( pastix_rhstype_t type, int nrhs, } /* Other format not supported for now */ - if( spm->fmttype != PastixCSC ) + if( spm->fmttype != PastixCSC ) { return PASTIX_ERR_BADPARAMETER; + } - if( spm->gN <= 0 ) + if( spm->gN <= 0 ) { return PASTIX_ERR_BADPARAMETER; + } - if( nrhs <= 0 ) + if( nrhs <= 0 ) { return PASTIX_ERR_BADPARAMETER; + } - if( (nrhs > 1) && (ldx < spm->n) ) + if( (nrhs > 1) && (ldx < spm->n) ) { return PASTIX_ERR_BADPARAMETER; + } - if( (nrhs > 1) && (ldb < spm->n) ) + if( (nrhs > 1) && (ldb < spm->n) ) { return PASTIX_ERR_BADPARAMETER; + } - if( spm->dof != 1 ) + if( spm->dof != 1 ) { return PASTIX_ERR_BADPARAMETER; + } if (nrhs == 1) { ldb = spm->n; @@ -292,7 +298,8 @@ z_spmGenRHS( pastix_rhstype_t type, int nrhs, spm->gN, 0, 0, 24356 ); } - rc = z_spmCSCMatVec( PastixNoTrans, &zone, spm, xptr, &zzero, bptr ); + /* Compute B */ + rc = z_spmCSCMatMat( PastixNoTrans, nrhs, &zone, spm, xptr, ldx, &zzero, bptr, ldb ); if ( x == NULL ) { memFree_null(xptr); @@ -360,6 +367,7 @@ z_spmCheckAxb( int nrhs, double normA, normB, normX, normX0, normR; double backward, forward, eps; int failure = 0; + int i; assert( spm->nexp == spm->n ); assert( spm->dof == 1 ); @@ -369,26 +377,43 @@ z_spmCheckAxb( int nrhs, /** * Compute the starting norms */ - normA = spmNorm( PastixFrobeniusNorm, spm ); - normX = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, nrhs, x, ldx ); - normB = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, nrhs, b, ldb ); + normA = spmNorm( PastixOneNorm, spm ); - printf( " || A ||_oo %e\n" - " || b ||_oo %e\n" - " || x ||_oo %e\n", + normB = 0.; + normX = 0.; + for( i=0; i<nrhs; i++ ) { + double norm; + + norm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, 1, b + i * ldb, ldb ); + normB = (norm > normB ) ? norm : normB; + norm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, 1, b + i * ldb, ldb ); + normX = (norm > normB ) ? norm : normB; + } + printf( " || A ||_1 %e\n" + " max(|| b_i ||_oo) %e\n" + " max(|| x_i ||_oo) %e\n", normA, normB, normX ); /** - * Compute r = A * x - b + * Compute r = b - A * x */ - spmMatVec( PastixNoTrans, &mzone, spm, x, &zone, b ); - normR = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, nrhs, b, ldb ); + spmMatMat( PastixNoTrans, nrhs, &mzone, spm, x, ldx, &zone, b, ldb ); + + normR = 0.; + backward = 0.; + for( i=0; i<nrhs; i++ ) { + double nx = cblas_dzasum( spm->n, x + i * ldx, 1 ); + double nb = cblas_dzasum( spm->n, b + i * ldb, 1 ); + double r = ((nb / normA) / nx) / eps; + + normR = (nb > normR) ? nb : normR; + backward = (r > backward) ? r : backward; + } - backward = normR / ( normA * normX + normB ); - failure = isnan(normX) || isinf(normX) || isnan(backward) || isinf(backward) || ((backward / eps) > 1.e3); + failure = isnan(normR) || isinf(normR) || isnan(backward) || isinf(backward) || (backward > 1.e3); - printf( " ||b-Ax||_oo %e\n" - " ||b-Ax||_oo / (||A||_oo ||x||_oo + ||b||_oo) %e (%s)\n", + printf( " max(|| b_i - A x_i ||_1) %e\n" + " max(|| b_i - A x_i ||_1) / (||A||_1 * ||x_i||_oo * eps) %e (%s)\n", normR, backward, failure ? "FAILED" : "SUCCESS" ); @@ -397,8 +422,7 @@ z_spmCheckAxb( int nrhs, */ if ( x0 != NULL ) { const pastix_complex64_t *zx = (const pastix_complex64_t *)x; - pastix_complex64_t *zx0 = (pastix_complex64_t *)x0;; - pastix_int_t i; + pastix_complex64_t *zx0 = (pastix_complex64_t *)x0; normX0 = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, nrhs, x0, ldx0 );