Integrate changes made in pastix to print the norm2 of the residual as done in...

Integrate changes made in pastix to print the norm2 of the residual as done in iterative refinements processes

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
      */