From a72016b6e7e6115472e91b25661d162560d1303d Mon Sep 17 00:00:00 2001
From: Mathieu Faverge <mathieu.faverge@inria.fr>
Date: Mon, 23 Apr 2018 17:22:53 +0200
Subject: [PATCH] Integrate changes made in pastix to print the norm2 of the
residual as done in iterative refinements processes
---
src/z_spm_genrhs.c | 21 ++++++++++++++++-----
1 file changed, 16 insertions(+), 5 deletions(-)
diff --git a/src/z_spm_genrhs.c b/src/z_spm_genrhs.c
index 8e1227a4..7016a0c0 100644
--- a/src/z_spm_genrhs.c
+++ b/src/z_spm_genrhs.c
@@ -131,7 +131,7 @@ z_spmRndVect( double scale, int m, int n, spm_complex64_t *A, int lda,
for (i = 0; i < m; ++i) {
*tmp = (0.5f - ran * RndF_Mul) * scale;
ran = Rnd64_A * ran + Rnd64_C;
-#ifdef COMPLEX
+#if defined(PRECISION_z) || defined(PRECISION_c)
*tmp += (I*(0.5f - ran * RndF_Mul)) * scale;
ran = Rnd64_A * ran + Rnd64_C;
#endif
@@ -368,7 +368,8 @@ z_spmCheckAxb( spm_fixdbl_t eps, int nrhs,
const spm_complex64_t *zx = (const spm_complex64_t *)x;
spm_complex64_t *zx0 = (spm_complex64_t *)x0;
spm_complex64_t *zb = (spm_complex64_t *)b;
- double normA, normB, normX, normX0, normR;
+ double *nb2 = malloc( nrhs * sizeof(double) );
+ double normA, normB, normX, normX0, normR, normR2;
double backward, forward;
int failure = 0;
int i;
@@ -394,6 +395,8 @@ z_spmCheckAxb( spm_fixdbl_t eps, int nrhs,
normB = (norm > normB ) ? norm : normB;
norm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'I', spm->n, 1, zx + i * ldx, ldx );
normX = (norm > normX ) ? norm : normX;
+
+ nb2[i] = cblas_dznrm2( spm->n, zb + i * ldb, 1 );
}
printf( " || A ||_1 %e\n"
" max(|| b_i ||_oo) %e\n"
@@ -406,22 +409,27 @@ z_spmCheckAxb( spm_fixdbl_t eps, int nrhs,
spmMatMat( SpmNoTrans, nrhs, &mzone, spm, x, ldx, &zone, b, ldb );
normR = 0.;
+ normR2 = 0.;
backward = 0.;
failure = 0;
for( i=0; i<nrhs; i++ ) {
double nx = cblas_dzasum( spm->n, zx + i * ldx, 1 );
double nr = cblas_dzasum( spm->n, zb + i * ldb, 1 );
+ double nr2 = cblas_dznrm2( spm->n, zb + i * ldb, 1 ) / nb2[i];
double back = ((nr / normA) / nx) / eps;
int fail = 0;
normR = (nr > normR ) ? nr : normR;
+ normR2 = (nr2 > normR2 ) ? nr2 : normR2;
backward = (back > backward) ? back : backward;
fail = isnan(nr) || isinf(nr) || isnan(back) || isinf(back) || (back > 1.e2);
if ( fail ) {
- printf( " || b_%d - A x_%d ||_1 %e\n"
+ printf( " || b_%d - A x_%d ||_2 / || b_%d ||_2 %e\n"
+ " || b_%d - A x_%d ||_1 %e\n"
" || b_%d - A x_%d ||_1 / (||A||_1 * ||x_%d||_oo * eps) %e (%s)\n",
+ i, i, i, nr2,
i, i, nr,
i, i, i, back,
fail ? "FAILED" : "SUCCESS" );
@@ -430,11 +438,14 @@ z_spmCheckAxb( spm_fixdbl_t eps, int nrhs,
failure = failure || fail;
}
- printf( " max(|| b_i - A x_i ||_1) %e\n"
+ printf( " max(|| b_i - A x_i ||_2 / || b_i ||_2) %e\n"
+ " 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,
+ normR2, normR, backward,
failure ? "FAILED" : "SUCCESS" );
+ free(nb2);
+
/**
* Compute r = x0 - x
*/
--
GitLab