Remove duplicate geometry configuration arguments from OpneMP loop function calls

src/cluster/cluster.cpp

@@ -53,7 +53,7 @@ using namespace std;
 
 // I would like to put it all in a struct, but then I'd have to write a constructor for it, due to members defined as references, creating a worse nightmare than the one I'd like to simplify...
 
-int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, 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, 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 *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, int isq, int ibf, Logger *logger);
+int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, 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, 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, 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, dcomplex arg, double wn, double vk, Logger *logger);
 
 /*! \brief C++ implementation of CLU
  *
@@ -97,22 +97,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
   int nsph = gconf->get_param("nsph");
   if (s_nsph == nsph) {
     // Shortcuts to variables stored in configuration objects
-    np_int mxndm = (np_int)gconf->get_param("mxndm");
-    int inpol = (int)gconf->get_param("in_pol");
-    int npnt = (int)gconf->get_param("npnt");
-    int npntts = (int)gconf->get_param("npntts");
-    int iavm = (int)gconf->get_param("iavm");
-    int isam = (int)gconf->get_param("meridional_type");
-    int nxi = (int)sconf->get_param("number_of_scales");
-    int idfc = (int)sconf->get_param("idfc");
     ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);
     double wp = sconf->get_param("wp");
     // Global variables for CLU
-    int lm = (int)gconf->get_param("l_max");
-    int li = (int)gconf->get_param("li");
-    int le = (int)gconf->get_param("le");
-    if (li > lm) lm = li;
-    if (le > lm) lm = le;
     C1 *c1 = new C1(gconf, sconf);
     C3 *c3 = new C3();
     C4 *c4 = new C4(gconf);
@@ -120,18 +107,15 @@ void cluster(const string& config_file, const string& data_file, const string& o
     // End of add-ons initialization
     C6 *c6 = new C6(c4->lmtpo);
     FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
-    int jer = 0;
-    int lcalc = 0;
     dcomplex arg = 0.0 + 0.0 * I;
     dcomplex ccsam = 0.0 + 0.0 * I;
-    int configurations = (int)sconf->get_param("configurations");
     C2 *c2 = new C2(gconf, sconf);
-    np_int ndit = 2 * nsph * c4->nlim;
+    const int ndi = c4->nsph * 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");
-    const int ndi = c4->nsph * c4->nlim;
     C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
-    double *gaps = new double[nsph]();
+    double *gaps = new double[(int)gconf->get_param("nsph")]();
     double *tqev = new double[3]();
     double *tqsv = new double[3]();
     double **tqse, **tqss, **tqce, **tqcs;
@@ -145,10 +129,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
     tqcs = new double*[2];
     tqcps = new dcomplex*[2];
     for (int ti = 0; ti < 2; ti++) {
-      tqse[ti] = new double[nsph]();
-      tqspe[ti] = new dcomplex[nsph]();
-      tqss[ti] = new double[nsph]();
-      tqsps[ti] = new dcomplex[nsph]();
+      tqse[ti] = new double[(int)gconf->get_param("nsph")]();
+      tqspe[ti] = new dcomplex[(int)gconf->get_param("nsph")]();
+      tqss[ti] = new double[(int)gconf->get_param("nsph")]();
+      tqsps[ti] = new dcomplex[(int)gconf->get_param("nsph")]();
       tqce[ti] = new double[3]();
       tqcpe[ti] = new dcomplex[3]();
       tqcs[ti] = new double[3]();
@@ -192,13 +176,14 @@ void cluster(const string& config_file, const string& data_file, const string& o
       cmullr[ci] = new double[4]();
       cmul[ci] = new double[4]();
     }
-    int isq, ibf;
-    int jwtm = (int)gconf->get_param("jwtm");
     double scan, cfmp, sfmp, cfsp, sfsp;
     // End of global variables for CLU
     fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
     fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
-	    nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam
+	    nsph, c4->li, c4->le, (np_int)gconf->get_param("mxndm"),
+	    (int)gconf->get_param("in_pol"), (int)gconf->get_param("npnt"),
+	    (int)gconf->get_param("npntts"), (int)gconf->get_param("iavm"),
+	    (int)gconf->get_param("meridional_type")
 	    );
     fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
     for (int ri = 0; ri < nsph; ri++) fprintf(output, "%17.8lE%17.8lE%17.8lE\n",
@@ -219,7 +204,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	    p_scattering_angles->phsstp, p_scattering_angles->phslst
 	    );
     fprintf(output, " READ(IR,*)JWTM\n");
-    fprintf(output, " %5d\n", jwtm);
+    fprintf(output, " %5d\n", (int)gconf->get_param("jwtm"));
     fprintf(output, "  READ(ITIN)NSPHT\n");
     fprintf(output, "  READ(ITIN)(IOG(I),I=1,NSPH)\n");
     fprintf(output, "  READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
@@ -241,7 +226,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
     thdps(c4->lm, zpv);
     double exdc = sconf->get_param("exdc");
     double exri = sqrt(exdc);
-    double vk = 0.0; // NOTE: Not really sure it should be initialized at 0
+    double vk = 0.0;
     fprintf(output, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
     fstream *tppoanp = new fstream;
     fstream &tppoan = *tppoanp;
@@ -253,6 +238,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
 #else
       logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);
 #endif
+      int iavm = (int)gconf->get_param("iavm");
+      int isam = (int)gconf->get_param("meridional_type");
+      int inpol = (int)gconf->get_param("in_pol");
+      int nxi = (int)sconf->get_param("nxi");
       int nth = p_scattering_angles->nth;
       int nths = p_scattering_angles->nths;
       int nph = p_scattering_angles->nph;
@@ -268,7 +257,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
       double wn = wp / 3.0e8;
       double xip = sconf->get_param("xip");
       double sqsfi = 1.0;
-      if (idfc < 0) {
+      if (sconf->get_param("idfc") < 0.0) {
 	vk = xip * wn;
 	fprintf(output, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
 	fprintf(output, " \n");
@@ -276,12 +265,11 @@ void cluster(const string& config_file, const string& data_file, const string& o
       // do the first iteration on jxi488 separately, since it seems to be different from the others
       int jxi488 = 1;
       chrono::time_point<chrono::high_resolution_clock> start_iter_1 = chrono::high_resolution_clock::now();
-      jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, 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, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, isq, ibf, logger);
+      int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, argi, args, gap, gapp, tqce, tqcpe, tqcs, tqcps, duk, tppoan, cextlr, cext, cmullr, cmul, gapv, tqev, tqsv, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, arg, wn, vk, logger);
       chrono::time_point<chrono::high_resolution_clock> end_iter_1 = chrono::high_resolution_clock::now();
       elapsed = end_iter_1 - start_iter_1;
       message = "INFO: First iteration took " + to_string(elapsed.count()) + "s.\n";
       logger->log(message);
-      logger->log("Finished: output written to " + output_path + "/c_OCLU\n");
       time_logger->log(message);
 
       // Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
@@ -332,15 +320,6 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	double *gapv_2 = NULL;
 	double *tqev_2 = NULL;
 	double *tqsv_2 = NULL;
-	int nxi_2 = nxi;
-	int nsph_2 = nsph;
-	np_int mxndm_2 = mxndm;
-	int inpol_2 = inpol;
-	int iavm_2 = iavm;
-	int npnt_2 = npnt;
-	int npntts_2 = npntts;
-	int isam_2 = isam;
-	int lm_2 = lm;
 	double *u_2 = NULL;
 	double *us_2 = NULL;
 	double *un_2 = NULL;
@@ -358,13 +337,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	double sfsp_2 = sfsp;
 	double sqsfi_2 = sqsfi;
 	double exri_2 = exri;
-	int lcalc_2 = lcalc;
 	dcomplex arg_2 = arg;
 	double wn_2 = wn;
 	double vk_2 = vk;
-	np_int ndit_2 = ndit;
-	int isq_2 = isq;
-	int ibf_2 = ibf;
 	// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
 	if (myompthread == 0) {
 	  sconf_2 = sconf;
@@ -545,7 +520,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	// ok, now I can actually start the parallel calculations
 #pragma omp for
 	for (jxi488 = 2; jxi488 <= nxi; jxi488++) {
-	  jer = cluster_jxi488_cycle(jxi488, sconf_2, gconf_2, p_scattering_angles, 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, 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, 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, isq_2, ibf_2, logger);
+	  int jer = cluster_jxi488_cycle(jxi488, sconf_2, gconf_2, p_scattering_angles, 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, 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, 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, arg_2, wn_2, vk_2, logger);
 	}
 
 #pragma omp barrier
@@ -631,7 +606,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	}
 #pragma omp barrier
 	{
-	  string message = "Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
+	  message = "Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
 	  logger->log(message);
 	}
       } // closes pragma omp parallel
@@ -642,7 +617,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	for (int ri = 1; 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);
-	  string message = "Copying ASCII output of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
+	  message = "Copying ASCII output of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
 	  logger->log(message);
 	  FILE *partial_output = fopen(partial_file_name.c_str(), "r");
 	  char c = fgetc(partial_output);
@@ -769,14 +744,29 @@ void cluster(const string& config_file, const string& data_file, const string& o
 }
 
 
-int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, 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, 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 *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, int isq, int ibf, Logger *logger)
+int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, 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, 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, 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, dcomplex arg, double wn, double vk, Logger *logger)
 {
+  int nxi = (int)sconf->get_param("nxi");
   logger->log("INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
   int jer = 0;
+  int lcalc = 0;
   int jaw = 1;
-  int jwtm = (int)gconf->get_param("jwtm");
   int li = (int)gconf->get_param("li");
   int le = (int)gconf->get_param("le");
+  int lm = 0;
+  if (le > lm) lm = le;
+  if (li > lm) lm = li;
+  int nsph = gconf->get_param("nsph");
+  int mxndm = gconf->get_param("mxndm");
+  int iavm = gconf->get_param("iavm");
+  int inpol = gconf->get_param("in_pol");
+  int npnt = gconf->get_param("npnt");
+  int npntts = gconf->get_param("npntts");
+  int isam = gconf->get_param("meridional_type");
+  int jwtm = (int)gconf->get_param("jwtm");
+  np_int ndit = 2 * nsph * c4->nlim;
+  int isq, ibf;
+  
   fprintf(output, "========== JXI =%3d ====================\n", jxi488);
   double xi = sconf->get_scale(jxi488 - 1);
   double exdc = sconf->get_param("exdc");