Commit 71a70bde authored by Mulas, Giacomo's avatar Mulas, Giacomo
Browse files

- For safety, since some loops are murky, overallocate some arrays of size... 

- For safety, since some loops are murky, overallocate some arrays of size nspheres instead of configurations
- correct a race in the definition of myompthread
- add some omp barriers to force synchronisations
parent 037cc8c1

src/cluster/cluster.cpp

+52 −36
Original line number Diff line number Diff line
@@ -95,7 +95,7 @@ using namespace std; 
//   double *tqsv;

// }



int cluster_jxi488_cycle(int &jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, int &nth, int &nths, int &nph, int &nphs, int &nk, int &nks, int &nkks, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, int &nxi, int &nsph, np_int &mxndm, int &inpol, int &iavm, int &npnt, int &npntts, int &isam, int &lm, double &th, double &thstp, double &thlst, double &ths, double &thsstp, double &thslst, double &ph, double &phstp, double &phlst, double &phs, double &phsstp, double &phslst, double &th1, double &ph1, double &ths1, double &phs1, double &thsca, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double &sqsfi, double &exri, int &lcalc, dcomplex &arg, double &wn, double &vk, np_int &ndit, dcomplex **am, int &isq, int &ibf);

int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, int nth, int nths, int nph, int nphs, int nk, int nks, int nkks, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, int nxi, int nsph, np_int mxndm, int inpol, int iavm, int npnt, int npntts, int isam, int lm, double th, double thstp, double thlst, double ths, double thsstp, double thslst, double ph, double phstp, double phlst, double phs, double phsstp, double phslst, double th1, double ph1, double ths1, double phs1, double thsca, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double sqsfi, double exri, int lcalc, dcomplex arg, double wn, double vk, np_int ndit, dcomplex **am, int isq, int ibf);



/*! \brief C++ implementation of CLU

 *
@@ -159,6 +159,8 @@ void cluster(string config_file, string data_file, string output_path) { 
      c1->rxx[c1i] = gconf->sph_x[c1i];

      c1->ryy[c1i] = gconf->sph_y[c1i];

      c1->rzz[c1i] = gconf->sph_z[c1i];

    }

    for (int c1i = 0; c1i < c1->configurations; c1i++) {

      c1->ros[c1i] = sconf->radii_of_spheres[c1i];

    }

    C4 *c4 = new C4(gconf->li, gconf->le, nsph);
@@ -344,8 +346,9 @@ void cluster(string config_file, string data_file, string output_path) { 
      int jxi488 = 1;

      jer = cluster_jxi488_cycle(jxi488, sconf, gconf, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, nth, nths, nph, nphs, nk, nks, nkks, argi, args, gap, gapp, tqce, tqcpe, tqcs, tqcps, duk, tppoan, cextlr, cext, cmullr, cmul, gapv, tqev, tqsv, nxi, nsph, mxndm, inpol, iavm, npnt, npntts, isam, lm, th, thstp, thlst, ths, thsstp, thslst, ph, phstp, phlst, phs, phsstp, phslst, th1, ph1, ths1, phs1, thsca, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, am, isq, ibf);



      int myompthread = 0;

      // Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled

      int ompnumthreads = 1;



#pragma omp parallel

      {

	// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one
@@ -358,13 +361,17 @@ void cluster(string config_file, string data_file, string output_path) { 
	C4 *c4_2 = new C4(*c4);

	C6 *c6_2 = new C6(*c6);

	C9 *c9_2 = new C9(*c9);

	// Create and initialise this variable here, so that if OpenMP is enabled it is local to the thread, and if OpenMP is not enabled it has a well-defiled value anyway

	int myompthread = 0;

#ifdef _OPENMP

	// If OpenMP is enabled, give actual values to myompthread and ompnumthreads, and open thread-local output files

	myompthread = omp_get_thread_num();

	if (myompthread == 0) ompnumthreads = omp_get_num_threads();

	FILE *output_2 = fopen((output_path + "/c_OCLU_" + to_string(myompthread)).c_str(), "w");

	fstream tppoan_2;

	tppoan_2.open((output_path + "/c_TPPOAN_" + to_string(myompthread)).c_str(), ios::out | ios::binary);

#else

	// If OpenMP is not enabled, just keep doing the output to the global files

	FILE *output_2 = output;

	fstream &tppoan_2 = tppoan;  

#endif
@@ -529,11 +536,14 @@ void cluster(string config_file, string data_file, string output_path) { 
	int nks_2 = nks;

	int nkks_2 = nkks;



	

#pragma omp barrier

#pragma omp for

	for (jxi488 = 2; jxi488 <= nxi; jxi488++) {

	jer = cluster_jxi488_cycle(jxi488, sconf_2, gconf_2, c1_2, c1ao_2, c2_2, c3_2, c4_2, c6_2, c9_2, output_2, output_path, gaps_2, tqse_2, tqspe_2, tqss_2, tqsps_2, zpv_2, gapm_2, gappm_2, nth_2, nths_2, nph_2, nphs_2, nk_2, nks_2, nkks_2, argi_2, args_2, gap_2, gapp_2, tqce_2, tqcpe_2, tqcs_2, tqcps_2, duk_2, tppoan_2, cextlr_2, cext_2, cmullr_2, cmul_2, gapv_2, tqev_2, tqsv_2, nxi_2, nsph_2, mxndm_2, inpol_2, iavm_2, npnt_2, npntts_2, isam_2, lm_2, th_2, thstp_2, thlst_2, ths_2, thsstp_2, thslst_2, ph_2, phstp_2, phlst_2, phs_2, phsstp_2, phslst_2, th1_2, ph1_2, ths1_2, phs1_2, thsca_2, u_2, us_2, un_2, uns_2, up_2, ups_2, unmp_2, unsmp_2, upmp_2, upsmp_2, scan_2, cfmp_2, sfmp_2, cfsp_2, sfsp_2, sqsfi_2, exri_2, lcalc_2, arg_2, wn_2, vk_2, ndit_2, am_2, isq_2, ibf_2);

	}



#pragma omp barrier

	delete sconf_2;

	delete gconf_2;

	delete c1_2;
@@ -546,6 +556,10 @@ void cluster(string config_file, string data_file, string output_path) { 
#ifdef _OPENMP

	fclose(output_2);

	tppoan_2.close();

#pragma omp barrier

	{

	  printf("Closing thread-local output files of thread %d and syncing threads\n", myompthread);

	}

#endif

	delete[] gaps_2;

	for (int ti = 0; ti <2 -1; ti++) {
@@ -614,9 +628,10 @@ void cluster(string config_file, string data_file, string output_path) { 
	delete[] upsmp_2;

	delete[] am_vector_2;

	delete[] am_2;

	

      } // jxi488 loop and omp parallel

      } // closes pragma omp parallel

#ifdef _OPENMP

#pragma omp barrier

      {

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

	  // Giovanni, please add here in this loop code to reopen the temporary files, reread them and append them respectively to the global output and tppoan, before closing them

	  string partial_file_name = output_path + "/c_OCLU_" + to_string(ri);
@@ -647,6 +662,7 @@ void cluster(string config_file, string data_file, string output_path) { 
	  remove(partial_file_name.c_str());

	  printf("done.\n");

	}

      }

#endif

      tppoan.close();

    } else { // In case TPPOAN could not be opened. Should never happen.
@@ -738,7 +754,7 @@ void cluster(string config_file, string data_file, string output_path) { 
}





int cluster_jxi488_cycle(int &jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, int &nth, int &nths, int &nph, int &nphs, int &nk, int &nks, int &nkks, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, int &nxi, int &nsph, np_int &mxndm, int &inpol, int &iavm, int &npnt, int &npntts, int &isam, int &lm, double &th, double &thstp, double &thlst, double &ths, double &thsstp, double &thslst, double &ph, double &phstp, double &phlst, double &phs, double &phsstp, double &phslst, double &th1, double &ph1, double &ths1, double &phs1, double &thsca, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double &sqsfi, double &exri, int &lcalc, dcomplex &arg, double &wn, double &vk, np_int &ndit, dcomplex **am, int &isq, int &ibf)

int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, int nth, int nths, int nph, int nphs, int nk, int nks, int nkks, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, int nxi, int nsph, np_int mxndm, int inpol, int iavm, int npnt, int npntts, int isam, int lm, double th, double thstp, double thlst, double ths, double thsstp, double thslst, double ph, double phstp, double phlst, double phs, double phsstp, double phslst, double th1, double ph1, double ths1, double phs1, double thsca, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double sqsfi, double exri, int lcalc, dcomplex arg, double wn, double vk, np_int ndit, dcomplex **am, int isq, int ibf)

{



  // int nxi = sconf->number_of_scales;

src/include/Commons.h

+1 −1
Original line number Diff line number Diff line
@@ -42,9 +42,9 @@ protected: 
  int lm;

  //! \brief NLMMT. QUESTION: definition?

  int nlmmt;

public:

  //! \brief Number of configurations

  int configurations;

public:

  //! \brief QUESTION: definition?

  dcomplex **rmi;

  //! \brief QUESTION: definition?

src/include/Configuration.h

+2 −2
Original line number Diff line number Diff line
@@ -53,7 +53,7 @@ class C9; 
class GeometryConfiguration {

  //! Temporary work-around to allow cluster() and sphere() peeking in.

  friend void cluster(std::string, std::string, std::string);

  friend int cluster_jxi488_cycle(int &, ScattererConfiguration *, GeometryConfiguration *, C1 *, C1_AddOns *, C2 *,  C3 *,  C4 *,  C6 *,  C9 *, std::FILE *, std::string, double *, double **, dcomplex **, double **, dcomplex **, double ****, double **, dcomplex **, int &, int &, int &, int &, int &, int &, int &, double *, double *, double **, dcomplex **, double **, dcomplex **, double **, dcomplex **, double *, std::fstream &, double **, double **, double **, double **, double *, double *, double *, int &, int &, np_int &, int &, int &, int &, int &, int &, int &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double &, double &, double &, double &, double &, double &, double &, int &, dcomplex &, double &, double &, np_int &, dcomplex **, int &, int &);

  friend int cluster_jxi488_cycle(int, ScattererConfiguration *, GeometryConfiguration *, C1 *, C1_AddOns *, C2 *,  C3 *,  C4 *,  C6 *,  C9 *, std::FILE *, std::string, double *, double **, dcomplex **, double **, dcomplex **, double ****, double **, dcomplex **, int, int, int, int, int, int, int, double *, double *, double **, dcomplex **, double **, dcomplex **, double **, dcomplex **, double *, std::fstream &, double **, double **, double **, double **, double *, double *, double *, int, int, np_int, int, int, int, int, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double &, double &, double &, double &, double &, double, double, int, dcomplex, double, double, np_int, dcomplex **, int, int);

  friend void sphere(std::string, std::string, std::string);

protected:

  //! \brief Number of spherical components.
@@ -185,7 +185,7 @@ public: 
class ScattererConfiguration {

  //! Temporary work-around to allow cluster() and sphere() peeking in.

  friend void cluster(std::string, std::string, std::string);

  friend int cluster_jxi488_cycle(int &, ScattererConfiguration *, GeometryConfiguration *, C1 *, C1_AddOns *, C2 *,  C3 *,  C4 *,  C6 *,  C9 *, std::FILE *, std::string, double *, double **, dcomplex **, double **, dcomplex **, double ****, double **, dcomplex **, int &, int &, int &, int &, int &, int &, int &, double *, double *, double **, dcomplex **, double **, dcomplex **, double **, dcomplex **, double *, std::fstream &, double **, double **, double **, double **, double *, double *, double *, int &, int &, np_int &, int &, int &, int &, int &, int &, int &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double &, double &, double &, double &, double &, double &, double &, int &, dcomplex &, double &, double &, np_int &, dcomplex **, int &, int &);

  friend int cluster_jxi488_cycle(int, ScattererConfiguration *, GeometryConfiguration *, C1 *, C1_AddOns *, C2 *,  C3 *,  C4 *,  C6 *,  C9 *, std::FILE *, std::string, double *, double **, dcomplex **, double **, dcomplex **, double ****, double **, dcomplex **, int, int, int, int, int, int, int, double *, double *, double **, dcomplex **, double **, dcomplex **, double **, dcomplex **, double *, std::fstream &, double **, double **, double **, double **, double *, double *, double *, int, int, np_int, int, int, int, int, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double &, double &, double &, double &, double &, double, double, int, dcomplex, double, double, np_int, dcomplex **, int, int);

  friend void sphere(std::string, std::string, std::string);

protected:

  //! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].

src/libnptm/Commons.cpp

+7 −5
Original line number Diff line number Diff line
@@ -38,8 +38,8 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) { 
    if (_iog[ci - 1] >= ci) configurations++;

  }

  vints = new dcomplex*[nsph];

  rc = new double*[configurations];

  nshl = new int[configurations]();

  rc = new double*[nsph];

  nshl = new int[nsph]();

  iog = new int[nsph]();

  int conf_index = 0;

  for (int vi = 0; vi < nsph; vi++) {
@@ -95,9 +95,9 @@ C1::C1(const C1& rhs) { 
  }

  configurations = rhs.configurations;

  vints = new dcomplex*[nsph];

  rc = new double*[configurations];

  nshl = new int[configurations]();

  for (int ni=0; ni<configurations; ni++) nshl[ni] = rhs.nshl[ni];

  rc = new double*[nsph];

  nshl = new int[nsph]();

  for (int ni=0; ni<nsph; ni++) nshl[ni] = rhs.nshl[ni];

  iog = new int[nsph]();

  for (int ni=0; ni<nsph; ni++) iog[ni] = rhs.iog[ni];

  int conf_index = 0;
@@ -140,6 +140,8 @@ C1::C1(const C1& rhs) { 
    rxx[si] = rhs.rxx[si];

    ryy[si] = rhs.ryy[si];

    rzz[si] = rhs.rzz[si];

  }

  for (int si = 0; si < configurations; si++) {

    ros[si] = rhs.ros[si];    

  }

}

src/libnptm/Configuration.cpp

+6 −6
Original line number Diff line number Diff line
@@ -523,9 +523,9 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam 
    }

  }

  int index = 0;

  double *ros_vector = new double[configurations]();

  int *nshl_vector = new int[configurations]();

  double **rcf_vector = new double*[configurations];

  double *ros_vector = new double[nsph]();

  int *nshl_vector = new int[nsph]();

  double **rcf_vector = new double*[nsph];

  for (int i113 = 1; i113 <= nsph; i113++) {

    if (iog_vector[i113 - 1] < i113) continue;

    str_target = file_lines[++last_read_line];
@@ -736,9 +736,9 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) { 
  for (int i = 1; i <= _nsph; i++) {

    if (_iog_vec[i - 1] >= i) configurations++;

  }

  _nshl_vec = new int[configurations]();

  _ros_vec = new double[configurations]();

  _rcf_vec = new double*[configurations];

  _nshl_vec = new int[_nsph]();

  _ros_vec = new double[_nsph]();

  _rcf_vec = new double*[_nsph];

  for (int i115 = 1; i115 <= _nsph; i115++) {

    if (_iog_vec[i115 - 1] < i115) continue;

    input.read(reinterpret_cast<char *>(&(_nshl_vec[i115 - 1])), sizeof(int));