Commit 279613b2 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Use runtime refinement settings in CLUSTER and INCLUSION

parent e744c4dd

src/cluster/cluster.cpp

+9 −28
Original line number Diff line number Diff line
@@ -244,8 +244,10 @@ void cluster(const string& config_file, const string& data_file, const string& o 
	message = "INFO: using RBT for inversion.\n";

	logger->log(message);

      } else if(gconf->invert_mode == RuntimeSettings::INV_MODE_SVD) {

	message = "INFO: using SVD for inversion.\n";

	//message = "INFO: using SVD for inversion.\n";

	message = "ERROR: SVD inversion mode not yet implemented!\n";

	logger->log(message);

	exit(1);

      }

      // Overlapping spheres test

      double tolerance = gconf->tolerance;
@@ -389,7 +391,6 @@ void cluster(const string& config_file, const string& data_file, const string& o 
      }



      int jxi488;

      int initialmaxrefiters = cid->maxrefiters;

      int jer = 0;



      //==================================================
@@ -840,7 +841,6 @@ int cluster_jxi488_cycle( 
	cid->update_orders(sconf->_rcf, new_li, new_le);

	is_first_scale = true;

	jaw = 1;

	cid->refinemode = 2;

      }

    }

  }
@@ -943,11 +943,7 @@ int cluster_jxi488_cycle( 
#ifdef USE_NVTX

  nvtxRangePush("Invert the matrix");

#endif

  // we put the accuracygoal in, get the actual accuracy back out

  double actualaccuracy = cid->accuracygoal;

  invert_matrix(cid->am, ndit, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, output_path, jxi488, mxndm, cid->proc_device, rs);

  // in principle, we should check whether the returned actualaccuracy is indeed lower than the accuracygoal, and do something about it if not

  cid->refinemode = 0;

  invert_matrix(cid->am, ndit, jer, output_path, jxi488, mxndm, cid->proc_device, rs);

#ifdef DEBUG_AM

  VirtualAsciiFile *outam2 = new VirtualAsciiFile();

  string outam2_name = output_path + "/c_AM2_JXI" + to_string(jxi488) + ".txt";
@@ -1746,11 +1742,6 @@ ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, Scatter 


  // the first time the object is used, it will be a "first iteration", to ensure data structures are properly initialised.

  is_first_scale = 1;

  // In the first iteration, if refinement is enabled, determine the number of refinement iterations required to arrive at the target accuracy (if achievable in a reasonable number of iterations)

  refinemode = 2;

  // maxrefiters and accuracygoal should be configurable and preferably set somewhere else

  maxrefiters = 20;

  accuracygoal = 1e-6;

}



ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
@@ -1895,9 +1886,6 @@ ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) { 


  proc_device = rhs.proc_device;

  is_first_scale = rhs.is_first_scale;

  refinemode = rhs.refinemode;

  maxrefiters = rhs.maxrefiters;

  accuracygoal = rhs.accuracygoal;

}



#ifdef MPI_VERSION
@@ -2031,10 +2019,7 @@ ClusterIterationData::ClusterIterationData(const mixMPI *mpidata, const int devi 
#else

  proc_device = 0;

#endif

  MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&is_first_scale, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);

  MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

}



void ClusterIterationData::mpibcast(const mixMPI *mpidata) {
@@ -2103,10 +2088,7 @@ void ClusterIterationData::mpibcast(const mixMPI *mpidata) { 
  MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&is_first_scale, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);

  MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

}

#endif
@@ -2192,14 +2174,13 @@ long ClusterIterationData::get_size(GeometryConfiguration *gconf, ScattererConfi 
  37 root pointers



  double scan; double cfmp; double sfmp; double cfsp; double sfsp; double sqsfi; double vk; double wn; double xip;

  double accuracygoal;

  10 double values

  9 double values

  

  dcomplex arg;

  1 dcomplex value



  int number_of_scales; int xiblock; int firstxi; int lastxi; int proc_device; int refinemode; int maxrefiters;

  7 int values

  int number_of_scales; int xiblock; int firstxi; int lastxi; int proc_device;

  5 int values

  

  bool is_first_scale;

  1 boolean value
@@ -2233,9 +2214,9 @@ long ClusterIterationData::get_size(GeometryConfiguration *gconf, ScattererConfi 
  const np_int ndit = 2 * ndi;

  const int lm = (gconf->li > gconf->le) ? gconf->li : gconf->le;

  long result = sizeof(long) * 37;

  result += sizeof(double) * 10;

  result += sizeof(double) * 9;

  result += sizeof(dcomplex);

  result += sizeof(int) * 7;

  result += sizeof(int) * 5;

  result += sizeof(bool);

  result += sizeof(long) * (44 + 6 * lm + ndit);

  result += sizeof(double) * (132 + 5 * nsph + 12 * lm);

src/inclusion/inclusion.cpp

+22 −38
Original line number Diff line number Diff line
@@ -234,6 +234,21 @@ void inclusion(const string& config_file, const string& data_file, const string& 
      sprintf(virtual_line, "%.5lg m.\n", sconf->get_particle_radius(gconf));

      message = "INFO: particle radius is " + (string)virtual_line;

      logger->log(message);

      if (gconf->invert_mode == RuntimeSettings::INV_MODE_LU) {

	message = "INFO: using LU factorization for inversion.\n";

	logger->log(message);

      } else if(gconf->invert_mode == RuntimeSettings::INV_MODE_GESV) {

	message = "INFO: using GESV system solver for inversion.\n";

	logger->log(message);

      } else if(gconf->invert_mode == RuntimeSettings::INV_MODE_RBT) {

	message = "INFO: using RBT for inversion.\n";

	logger->log(message);

      } else if(gconf->invert_mode == RuntimeSettings::INV_MODE_SVD) {

	//message = "INFO: using SVD for inversion.\n";

	message = "ERROR: SVD inversion mode not yet implemented!\n";

	logger->log(message);

	exit(1);

      }

      // Overlapping spheres test

      double tolerance = gconf->tolerance;

      if (tolerance < 0.0) {
@@ -371,7 +386,6 @@ void inclusion(const string& config_file, const string& data_file, const string& 
      }

      

      int jxi488;

      int initialmaxrefiters = cid->maxrefiters;

      int jer = 0;



      //==================================================
@@ -730,6 +744,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
  string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";

  Logger *logger = new Logger(LOG_DEBG);

  logger->log(message);

  RuntimeSettings rs(gconf, logger);

  chrono::duration<double> elapsed;

  chrono::time_point<chrono::high_resolution_clock> interval_start, interval_end;

  int jer = 0;
@@ -816,7 +831,6 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
	cid->update_orders(sconf->_rcf, new_li, new_le);

	is_first_scale = true;

	jaw = 1;

	cid->refinemode = 2;

      }

    }

  }
@@ -923,22 +937,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
#ifdef USE_NVTX

  nvtxRangePush("Invert the matrix");

#endif

  // we the accuracygoal in, get the actual accuracy back out

  double actualaccuracy = cid->accuracygoal;

  invert_matrix(cid->am, cid->c1->ndm, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, output_path, jxi488, mxndm, cid->proc_device);

  // in principle, we should check whether the returned actualaccuracy is indeed lower than the accuracygoal, and do something about it if not

  if (gconf->refine_flag > 0) {

    if (cid->refinemode==2) {

      message = "DEBUG: iterative refinement enabled at run-time.\n";

      logger->log(message, LOG_DEBG);

      message = "INFO: calibration obtained accuracy " + to_string(actualaccuracy) + " (" + to_string(cid->accuracygoal) + " requested) in " + to_string(cid->maxrefiters) + " refinement iterations\n";

      logger->log(message);

      if (actualaccuracy > 1e-2) {

	printf("Accuracy worse than 0.01, stopping");

	exit(1);

      }

    }

  }

  invert_matrix(cid->am, cid->c1->ndm, jer, output_path, jxi488, mxndm, cid->proc_device, rs);

#ifdef USE_NVTX

  nvtxRangePop();

#endif
@@ -1596,11 +1595,6 @@ InclusionIterationData::InclusionIterationData(GeometryConfiguration *gconf, Sca 


  // the first time the object is used, it will be a "first iteration", to ensure data structures are properly initialised.

  is_first_scale = 1;

  // In the first iteration, if refinement is enabled, determine the number of refinement iterations required to arrive at the target accuracy (if achievable in a reasonable number of iterations)

  refinemode = 2;

  // maxrefiters and accuracygoal should be configurable and preferably set somewhere else

  maxrefiters = 20;

  accuracygoal = 1e-6;

}



InclusionIterationData::InclusionIterationData(const InclusionIterationData& rhs) {
@@ -1748,9 +1742,6 @@ InclusionIterationData::InclusionIterationData(const InclusionIterationData& rhs 
  

  proc_device = rhs.proc_device;

  is_first_scale = rhs.is_first_scale;

  refinemode = rhs.refinemode;

  maxrefiters = rhs.maxrefiters;

  accuracygoal = rhs.accuracygoal;

}



#ifdef MPI_VERSION
@@ -1885,10 +1876,7 @@ InclusionIterationData::InclusionIterationData(const mixMPI *mpidata, const int 
#else

  proc_device = 0;

#endif

  MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&is_first_scale, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);

  MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

}



void InclusionIterationData::mpibcast(const mixMPI *mpidata) {
@@ -1951,10 +1939,7 @@ void InclusionIterationData::mpibcast(const mixMPI *mpidata) { 
  MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&nimd, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&extr, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&is_first_scale, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);

  MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

}

#endif
@@ -2040,15 +2025,14 @@ long InclusionIterationData::get_size(GeometryConfiguration *gconf, ScattererCon 
  37 root pointers

 

  double extr; double vk; double wn; double xip; double scan; double cfmp; double sfmp; double cfsp;

  double sfsp; double sqsfi; double accuracygoal;

  11 double values

  double sfsp; double sqsfi;

  10 double values



  dcomplex arg;

  1 dcomplex value

  

  int nimd; int number_of_scales; int xiblock; int firstxi; int lastxi; int proc_device;

  int refinemode; int maxrefiters;

  8 int values

  6 int values



  bool is_first_scale;

  1 boolean value
@@ -2083,9 +2067,9 @@ long InclusionIterationData::get_size(GeometryConfiguration *gconf, ScattererCon 
  const int ndm = 2 * (nsph * nlim + nlem);

  const int lm = (gconf->li > gconf->le) ? gconf->li : gconf->le;

  long result = sizeof(long) * 37;

  result += sizeof(double) * 11;

  result += sizeof(double) * 10;

  result += sizeof(dcomplex);

  result += sizeof(int) * 8;

  result += sizeof(int) * 6;

  result += sizeof(bool);

  result += sizeof(long) * (44 + 6 * lm + ndm);

  result += sizeof(double) * (132 + 5 * nsph + 12 * lm);