Make INCLUSION a minimal executable to test memory structures

 * \param mpidata: `mixMPI *` Pointer to an instance of MPI data settings.

 */

void inclusion(const string& config_file, const string& data_file, const string& output_path, const mixMPI *mpidata) {

  chrono::time_point<chrono::high_resolution_clock> t_start = chrono::high_resolution_clock::now();

  chrono::duration<double> elapsed;

  string message;

  string timing_name;

  FILE *timing_file;

  Logger *time_logger;

  if (mpidata->rank == 0) {

    timing_name = output_path + "/c_timing_mpi"+ to_string(mpidata->rank) +".log";

    timing_file = fopen(timing_name.c_str(), "w");

    time_logger = new Logger(LOG_DEBG, timing_file);

  ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file);

  GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file);

  ParticleDescriptorSphere *pds = NULL;

  ParticleDescriptorCluster *pdc = NULL;

  ParticleDescriptorInclusion *pdi = NULL;

  if (gconf->number_of_spheres == 1) pds = new ParticleDescriptorSphere(gconf, sconf);

  else {

    if (sconf->use_external_sphere) pdi = new ParticleDescriptorInclusion(gconf, sconf);

    else pdc = new ParticleDescriptorCluster(gconf, sconf);

  }

  Logger *logger = new Logger(LOG_DEBG);

  int device_count = 0;

  //===========

  // Initialise MAGMA

  //===========

#ifdef USE_MAGMA

  cudaGetDeviceCount(&device_count);

  logger->log("DEBUG: Proc-" + to_string(mpidata->rank) + " found " + to_string(device_count) + " CUDA devices.\n", LOG_DEBG);

  logger->log("INFO: Process " + to_string(mpidata->rank) + " initializes MAGMA.\n");

  magma_int_t magma_result = magma_init();

  if (magma_result != MAGMA_SUCCESS) {

    logger->err("ERROR: Process " + to_string(mpidata->rank) + " failed to initilize MAGMA.\n");

    logger->err("PROC-" + to_string(mpidata->rank) + ": MAGMA error code " + to_string(magma_result) + "\n");

    if (mpidata->rank == 0) {

      fclose(timing_file);

      delete time_logger;

    }

    delete logger;

    return;

  }

#endif

  // end MAGMA initialisation

  //===========================

  // the following only happens on MPI process 0

  //===========================

  if (mpidata->rank == 0) {

#ifdef USE_NVTX

    nvtxRangePush("Set up");

#endif

    //=======================

    // Initialise sconf from configuration file

    //=======================

    logger->log("INFO: making legacy configuration...", LOG_INFO);

    ScattererConfiguration *sconf = NULL;

    try {

      sconf = ScattererConfiguration::from_dedfb(config_file);

    } catch(const OpenConfigurationFileException &ex) {

      logger->err("\nERROR: failed to open scatterer configuration file.\n");

      string message = "FILE: " + string(ex.what()) + "\n";

      logger->err(message);

      fclose(timing_file);

      delete time_logger;

      delete logger;

      return;

    }

    sconf->write_formatted(output_path + "/c_OEDFB");

    sconf->write_binary(output_path + "/c_TEDF");

    sconf->write_binary(output_path + "/c_TEDF.hd5", "HDF5");

    // end logger initialisation

    //========================

    // Initialise gconf from configuration files

    //========================

    GeometryConfiguration *gconf = NULL;

    try {

      gconf = GeometryConfiguration::from_legacy(data_file);

    } catch (const OpenConfigurationFileException &ex) {

      logger->err("\nERROR: failed to open geometry configuration file.\n");

      string message = "FILE: " + string(ex.what()) + "\n";

      logger->err(message);

      if (sconf) delete sconf;

      fclose(timing_file);

      delete time_logger;

      delete logger;

      return;

    }

    logger->log(" done.\n", LOG_INFO);

    //end gconf initialisation

#ifdef USE_NVTX

    nvtxRangePop();

#endif

    int s_nsph = sconf->number_of_spheres;

    int nsph = gconf->number_of_spheres;

    // Sanity check on number of sphere consistency, should always be verified

    if (s_nsph == nsph) {

      // Shortcuts to variables stored in configuration objects

      ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);

      double wp = sconf->wp;

      // Open empty virtual ascii file for output

      VirtualAsciiFile *p_output = new VirtualAsciiFile();

      // for the time being, this is ok. When we can, add some logic in the sprintf calls that checks if a longer buffer would be needed, and in case expand it

      // in any case, replace all sprintf() with snprintf(), to avoid in any case writing more than the available buffer size

      char virtual_line[256];

      // Create and initialise pristine cid for MPI proc 0 and thread 0

      ClusterIterationData *cid = new ClusterIterationData(gconf, sconf, mpidata, device_count);

      const int ndi = cid->c4->nsph * cid->c4->nlim;

      np_int ndit = 2 * ndi;

      logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");

      time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");

      //==========================

      // Write a block of info to the ascii output file

      //==========================

      sprintf(virtual_line, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");

      p_output->append_line(virtual_line);

#ifdef USE_ILP64

      sprintf(virtual_line, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",

	      nsph, cid->c4->li, cid->c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,

	      gconf->npntts, gconf->iavm, gconf->iavm

	      );

#else

      sprintf(virtual_line, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",

	      nsph, cid->c4->li, cid->c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,

	      gconf->npntts, gconf->iavm, gconf->iavm

	      );

#endif      

      p_output->append_line(virtual_line);

      sprintf(virtual_line, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");

      p_output->append_line(virtual_line);

      for (int ri = 0; ri < nsph; ri++) {

	sprintf(virtual_line, "%17.8lE%17.8lE%17.8lE\n",

		gconf->get_sph_x(ri), gconf->get_sph_y(ri), gconf->get_sph_z(ri)

		);

	p_output->append_line(virtual_line);

      }

      sprintf(virtual_line, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");

      p_output->append_line(virtual_line);

      sprintf(

	      virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",

	      p_scattering_angles->th, p_scattering_angles->thstp,

	      p_scattering_angles->thlst, p_scattering_angles->ths,

	      p_scattering_angles->thsstp, p_scattering_angles->thslst

	      );

      p_output->append_line(virtual_line);

      sprintf(virtual_line, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");

      p_output->append_line(virtual_line);

      sprintf(

	      virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",

	      p_scattering_angles->ph, p_scattering_angles->phstp,

	      p_scattering_angles->phlst, p_scattering_angles->phs,

	      p_scattering_angles->phsstp, p_scattering_angles->phslst

	      );

      p_output->append_line(virtual_line);

      sprintf(virtual_line, " READ(IR,*)JWTM\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, " %5d\n", gconf->jwtm);

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)NSPHT\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)(IOG(I),I=1,NSPH)\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)(XIV(I),I=1,NXI)\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)NSHL(I),ROS(I)\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, " \n");

      p_output->append_line(virtual_line);

      str(sconf, cid->c1, cid->c1ao, cid->c3, cid->c4, cid->c6);

      thdps(cid->c4->lm, cid->zpv);

      double exdc = sconf->exdc;

      double exri = sqrt(exdc);

      sprintf(virtual_line, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);

      p_output->append_line(virtual_line);

      // Create empty virtual binary file

      VirtualBinaryFile *vtppoanp = new VirtualBinaryFile();

      string tppoan_name = output_path + "/c_TPPOAN";

#ifdef USE_MAGMA

      logger->log("INFO: using MAGMA calls.\n", LOG_INFO);

#elif defined USE_LAPACK

      logger->log("INFO: using LAPACK calls.\n", LOG_INFO);

#else

      logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);

#endif

      int iavm = gconf->iavm;

      int isam = gconf->isam;

      int inpol = gconf->in_pol;

      int nxi = sconf->number_of_scales;

      int nth = p_scattering_angles->nth;

      int nths = p_scattering_angles->nths;

      int nph = p_scattering_angles->nph;

      int nphs = p_scattering_angles->nphs;

      //========================

      // write a block of info to virtual binary file

      //========================

      vtppoanp->append_line(VirtualBinaryLine(iavm));

      vtppoanp->append_line(VirtualBinaryLine(isam));

      vtppoanp->append_line(VirtualBinaryLine(inpol));

      vtppoanp->append_line(VirtualBinaryLine(nxi));

      vtppoanp->append_line(VirtualBinaryLine(nth));

      vtppoanp->append_line(VirtualBinaryLine(nph));

      vtppoanp->append_line(VirtualBinaryLine(nths));

      vtppoanp->append_line(VirtualBinaryLine(nphs));

      if (sconf->idfc < 0) {

	cid->vk = cid->xip * cid->wn;

	sprintf(virtual_line, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", cid->vk);

	p_output->append_line(virtual_line);

	sprintf(virtual_line, " \n");

	p_output->append_line(virtual_line);

      }

      //==================================================

      // do the first outputs here, so that I open here the new files, afterwards I only append

      //==================================================

      p_output->write_to_disk(output_path + "/c_OINCLU");

      delete p_output;

      vtppoanp->write_to_disk(output_path + "/c_TPPOAN");

      delete vtppoanp;

      delete cid;

      delete p_scattering_angles;

    }

    delete sconf;

  ParticleDescriptor *pd;

  if (pds != NULL) pd = pds;

  else if (pdc != NULL) pd = pdc;

  else pd = pdi;

  printf("INFO: %s\n", (pd->get_descriptor_type()).c_str());

  delete pd;

  delete gconf;

#ifdef USE_MAGMA

    logger->log("INFO: Process " + to_string(mpidata->rank) + " finalizes MAGMA.\n");

    magma_finalize();

#endif

    chrono::time_point<chrono::high_resolution_clock> t_end = chrono::high_resolution_clock::now();

    elapsed = t_end - t_start;

    string message = "INFO: Calculation lasted " + to_string(elapsed.count()) + "s.\n";

    logger->log(message);

    logger->log("Finished: output written to " + output_path + "/c_OINCLU\n");

    time_logger->log(message);

    fclose(timing_file);

    delete time_logger;

    delete logger;

  } // end instructions block of MPI process 0

  delete sconf;

}