Remove HMAT options for computing residual norm and iterative refinement

include/main.h

@@ -45,17 +45,6 @@ extern ScalarType stype;
 /*! \brief  Wavelength (for oscillatory kernels). */
 extern double lambda;
 
-/*! \brief  Flag to compute the residual norm after the Hmat solve (default is 0).
-
-  It can be changed with --hmat_residual (read in readFlagsTestHMAT() in hmat.c)
-*/
-extern int hmat_residual;
-
-/*! \brief  Flag to use Hmat iterative refinement (default is 0).
-
-  It can be changed with --hmat_refine (read in readFlagsTestHMAT() in hmat.c)
-*/
-extern int hmat_refine;
 /*! \brief  Flag to use Hmat shuffle clustering (default is 0).
 
   It can be changed with --use_hmat_shuffle (read in readFlagsTestHMAT() in hmat.c)

src/help.c

@@ -42,8 +42,6 @@ int printHelp() {
          "     --hmat-ldlt: Force LDLt decomposition (for symmetric matrices).\n"
          "     --hmat-coarsening: Coarsening of the HMatrix.\n"
          "     --hmat-recompress: Recompress the HMatrix.\n"
-         "     --hmat_residual: Compute the residual norm.\n"
-         "     --hmat_refine: Use iterative refinement.\n"
          "     --hmat-val-null: Validate the detection of null rows and columns.\n"
          "     --hmat-val: Validate the rk-matrices after compression.\n"
          "     --hmat-val-threshold: Error threshold for the compression validation.\n"

src/hmat.c

@@ -20,19 +20,6 @@ int init_hmat_interface() {
 }
 
 int readFlagsTestHMAT(int *argc, char ***argv) {
-  /* Flag to compute the residual norm after the Hmat solve */
-  if (MpfArgHasName(argc, *argv, 1, "--hmat_residual")) {
-    hmat_residual=1;
-    printf("Computing the residual norm after the Hmat solve\n") ;
-  }
-
-  /* Flag to use Hmat iterative refinement */
-  if (MpfArgHasName(argc, *argv, 1, "--hmat_refine")) {
-    hmat_refine=1;
-    printf("Using Hmat iterative refinement\n") ;
-    hmat_residual=1; // Iterative refinement needs the residual computation
-  }
-
   /* Flag to use Hmat shuffle clustering */
   if (MpfArgHasName(argc, *argv, 1, "--use_hmat_shuffle")) {
     use_hmat_shuffle=1;

src/main.c

@@ -19,8 +19,6 @@ int simplePrec = 0;
 int complexALGO = 1;
 ScalarType stype = DOUBLE_COMPLEX;
 double lambda;
-int hmat_residual = 0;
-int hmat_refine = 0;
 int use_hmat_shuffle = 0;
 int divider_min = 2;
 int divider_max = 4;

src/testHMAT.c

@@ -60,23 +60,6 @@ int testHMAT(double * relative_error) {
   /* Perform a HMAT solve */
   printf("\n") ;
 
-  hmat_matrix_t* hmatrix_orig = NULL;
-  void *solCLA_orig = NULL;
-  void *solCLA_sum = NULL;
-  double normRes, normOrig, temps_refine=-getTime();
-  int nbIter = 0;
-  if (hmat_residual) {
-    hmatrix_orig = hi->copy(hmatrix); // ideally, we do a copy with lower accuracy for factorisation
-    /* Vector to duplicate solCLA */
-    solCLA_orig = MpfCalloc(vectorSize, sizeof(D_type)) ; CHKPTRQ(solCLA_orig) ;
-    /* Norm of the originaal RHS */
-    ierr = computeVecNorm(solCLA, &normOrig); CHKERRQ(ierr);
-  }
-  if (hmat_refine) {
-    /* Vector to accumulate solCLA */
-    solCLA_sum = MpfCalloc(vectorSize, sizeof(D_type)) ; CHKPTRQ(solCLA_sum) ;
-  }
-
   /* Factorize the H-matrix */
 
   printf("\n**** Factorizing HMAT...\n") ;
@@ -93,8 +76,6 @@ int testHMAT(double * relative_error) {
   temps_cpu = (temps_final - temps_initial) ;
   printf("<PERFTESTS> TpsCpuFacto%s = %f \n", postfix_async, temps_cpu) ;
 
-  temps_refine -= temps_cpu; // We do not count facto in the iterative refinement timer
-
   /* display some informations (synchronous) */
   if (hmat_get_sync_exec()) {
     hmat_info_t info;
@@ -107,14 +88,7 @@ int testHMAT(double * relative_error) {
 
   /* Solve the system : A-1.solCLA -> solCLA */
 
-  printf("\n**** Solve HMAT%s...\n", hmat_refine ? " with Iterative Refinement" : "") ;
-
-_refineLoop: // Iterative Refinement loop
-  nbIter++ ;
-  if (hmat_residual) {
-    /* Duplicate the RHS solCLA into solCLA_orig */
-    memcpy(solCLA_orig, solCLA, vectorSize*sizeof(D_type));
-  }
+  printf("\n**** Solve HMAT...\n") ;
 
   temps_initial = getTime ();
 
@@ -130,58 +104,6 @@ _refineLoop: // Iterative Refinement loop
   temps_cpu = (temps_final - temps_initial) ;
   printf("<PERFTESTS> TpsCpuSolve%s = %f \n", postfix_async, temps_cpu) ;
 
-  /* Compute the residual solCLA_orig-A_orig.solCLA */
-
-  if (hmat_residual) {
-    temps_initial = getTime ();
-    if (simplePrec) {
-      doubleToSimple(solCLA, vectorSize);
-      doubleToSimple(solCLA_orig, vectorSize);
-    }
-    ierr = hi->gemv('N', Mpf_mone[stype], hmatrix_orig, solCLA, Mpf_pone[stype], solCLA_orig, nbRHS);
-    if (simplePrec) {
-      simpleToDouble(solCLA, vectorSize);
-      simpleToDouble(solCLA_orig, vectorSize);
-    }
-    temps_final = getTime ();
-    temps_cpu = (temps_final - temps_initial) ;
-    printf("<PERFTESTS> TpsCpuGEMV%s = %f \n", postfix_async, temps_cpu) ;
-    ierr = computeVecNorm(solCLA_orig, &normRes); CHKERRQ(ierr);
-    printf("<PERFTESTS> Residual_%d = %.4e \n", nbIter, normRes/normOrig);
-    // Test if convergence is achieved or if iteration number is too large
-    if (normRes/normOrig < mpf_hmat_settings.acaEpsilon/10 || nbIter > 20) {
-      // TODO stop the computation if the residual norm increases
-      printf("\n") ;
-      printf("<PERFTESTS> NbIterRefinement = %d \n", nbIter);
-      hmat_refine=0;
-    }
-  }
-  if (solCLA_sum) {
-    /* Compute the updated solution solCLA_sum += solCLA */
-    size_t k;
-    for (k=0 ; k<vectorSize ; k++)
-      ((double*)solCLA_sum)[k] += ((double*)solCLA)[k];
-    ierr=computeRelativeError(solCLA_sum, rhs, relative_error) ; CHKERRQ(ierr) ;
-    printf("<PERFTESTS> Error_%d = %.4e \n", nbIter, *relative_error);
-  }
-  if (hmat_refine) {
-    /* Copy the updated RHS of the system solCLA_orig into solCLA */
-    memcpy(solCLA, solCLA_orig, vectorSize*sizeof(D_type));
-    goto _refineLoop;
-  }
-
-  if (hmatrix_orig) hi->destroy(hmatrix_orig);
-  hmatrix_orig=NULL;
-  if (solCLA_orig) MpfFree(solCLA_orig);
-  solCLA_orig=NULL ;
-  if (solCLA_sum) { // If solCLA_sum has been allocated, I use it as a replacement for solCLA
-    MpfFree(solCLA) ;
-    solCLA=solCLA_sum ;
-    temps_refine += getTime();
-    printf("<PERFTESTS> TpsCpuRefinement%s = %f \n", postfix_async, temps_refine) ;
-
-  }
-
   /* Compare the two vectors solCLA and rhs */
 
   printf("\n**** Comparing results...\n") ;