From 097d268223ae6f4a21b9605cfffdd9b289981123 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Sat, 23 Mar 2024 21:00:37 +0100
Subject: [PATCH] Remove redundant C1_AddOns vectors and set up logging for
 parallel execution

---
 src/cluster/cluster.cpp     | 66 ++++++-------------------------------
 src/include/Commons.h       |  8 ++---
 src/include/Configuration.h |  5 +--
 src/libnptm/Commons.cpp     |  6 ++--
 4 files changed, 17 insertions(+), 68 deletions(-)

diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index aae0e743..9609ebe2 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -53,53 +53,8 @@ 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...
 //
-// struct clu_jxi488_cycle_args {
-//   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 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, Logger &logger);
 
 /*! \brief C++ implementation of CLU
  *
@@ -351,7 +306,7 @@ void cluster(string config_file, string data_file, string output_path) {
       }
       // do the first iteration on jxi488 separately, since it seems to be different from the others
       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);
+      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, logger);
 
       // 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;
@@ -549,7 +504,7 @@ void cluster(string config_file, string data_file, string output_path) {
 #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);
+	  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, logger);
 	}
 
 #pragma omp barrier
@@ -762,12 +717,11 @@ void cluster(string config_file, string data_file, string output_path) {
   }
   delete sconf;
   delete gconf;
-  logger.flush(LOG_INFO);
   logger.log("Finished: output written to " + output_path + "/c_OCLU\n");
 }
 
 
-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, Logger &logger)
 {
 
   // int nxi = sconf->number_of_scales;
@@ -1225,10 +1179,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 		      real(c1->sas[i226 - 1][0][1]), imag(c1->sas[i226 - 1][0][1]),
 		      real(c1->sas[i226 - 1][1][1]), imag(c1->sas[i226 - 1][1][1])
 		      );
-	      for (int j225 = 0; j225 < 16; j225++) { // QUESTION: check that 16 is a fixed dimension
-		c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];
+	      for (int j225 = 0; j225 < 16; j225++) {
+		c1->vint[j225] = c1->vints[i226 - 1][j225];
 	      } // j225 loop
-	      mmulc(c1ao->vint, cmullr, cmul);
+	      mmulc(c1->vint, cmullr, cmul);
 	      fprintf(output, "  MULS\n");
 	      for (int i1 = 0; i1 < 4; i1++) {
 		fprintf(
@@ -1258,7 +1212,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 	  fprintf(output, "     CLUSTER\n");
 	  pcros(vk, exri, c1, c1ao, c4);
 	  mextc(vk, exri, c1ao->fsac, cextlr, cext);
-	  mmulc(c1ao->vint, cmullr, cmul);
+	  mmulc(c1->vint, cmullr, cmul);
 	  if (jw != 0) {
 	    jw = 0;
 	    // Some implicit loops writing to binary.
@@ -1323,9 +1277,9 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 	  }
 	  // label 254
 	  for (int i = 0; i < 16; i++) {
-	    double value = real(c1ao->vint[i]);
+	    double value = real(c1->vint[i]);
 	    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	    value = imag(c1ao->vint[i]);
+	    value = imag(c1->vint[i]);
 	    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 	  }
 	  for (int i = 0; i < 4; i++) {
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 381f9ec0..90a911eb 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -63,10 +63,10 @@ public:
   //! \brief QUESTION: definition?
   dcomplex *fsas;
   //! \brief QUESTION: definition?
-  dcomplex **vints;
-  //! \brief QUESTION: definition?
   dcomplex *vint;
   //! \brief QUESTION: definition?
+  dcomplex **vints;
+  //! \brief QUESTION: definition?
   double *sscs;
   //! \brief QUESTION: definition?
   double *sexs;
@@ -270,10 +270,6 @@ public:
   //! \brief QUESTION: definition?
   dcomplex *vintm;
   //! \brief QUESTION: definition?
-  dcomplex *vint;
-  //! \brief QUESTION: definition?
-  dcomplex *vintm;
-  //! \brief QUESTION: definition?
   dcomplex *vintt;
   //! \brief QUESTION: definition?
   dcomplex **fsac;
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index 591fc0ba..efa4986e 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -32,6 +32,7 @@
 
 class ScattererConfiguration;
 class GeometryConfiguration;
+class Logger;
 class C1;
 class C1_AddOns;
 class C2;
@@ -53,7 +54,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, Logger &);
   friend void sphere(std::string, std::string, std::string);
 protected:
   //! \brief Number of spherical components.
@@ -185,7 +186,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, Logger &);
   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].
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index b76e3cc1..d6b0c558 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -102,6 +102,8 @@ C1::C1(const C1& rhs) {
     for (int wj=0; wj<4; wj++) w[wi][wj] = rhs.w[wi][wj];
   }
   configurations = rhs.configurations;
+  vint = new dcomplex[16]();
+  for (int vi = 0; vi < 16; vi++) vint[vi] = rhs.vint[vi];
   vec_vints = new dcomplex[nsph * 16]();
   vints = new dcomplex*[nsph];
   rc = new double*[nsph];
@@ -211,7 +213,6 @@ C1_AddOns::C1_AddOns(C4 *c4) {
   for (int ai = 0; ai < nlemt; ai++) {
     am0m[ai] = new dcomplex[nlemt]();
   }
-  vint = new dcomplex[16]();
   vintm = new dcomplex[16]();
   vintt = new dcomplex[16]();
   // This calculates the size of v3j0
@@ -270,11 +271,9 @@ C1_AddOns::C1_AddOns(const C1_AddOns& rhs) {
     am0m[ai] = new dcomplex[nlemt]();
     for (int aj = 0; aj < nlemt; aj++) am0m[ai][aj] = rhs.am0m[ai][aj];
   }
-  vint = new dcomplex[16](); // QUESTION: is dimension 16 generally fixed?
   vintm = new dcomplex[16]();
   vintt = new dcomplex[16]();
   for (int hi = 0; hi < 16; hi++) {
-    vint[hi] = rhs.vint[hi];
     vintm[hi] = rhs.vintm[hi];
     vintt[hi] = rhs.vintt[hi];
   }
@@ -333,7 +332,6 @@ C1_AddOns::~C1_AddOns() {
     delete[] am0m[ai];
   }
   delete[] am0m;
-  delete[] vint;
   delete[] vintm;
   delete[] vintt;
   for (int fi = 1; fi > -1; fi--) {
-- 
GitLab