diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 6ccfa80dd59ee3f895090331576768c4d5c8724d..84eeb8c03d13369bf1eb25c4a27c99718cd7a600 100644
--- a/src/cluster/cluster.cpp
+++ b/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, 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, 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, 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);
/*! \brief C++ implementation of CLU
*
@@ -63,6 +63,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
*/
void cluster(const string& config_file, const string& data_file, const string& output_path) {
chrono::time_point<chrono::high_resolution_clock> t_start = chrono::high_resolution_clock::now();
+ chrono::duration<double> elapsed;
+ string message;
string timing_name = output_path + "/c_timing.log";
FILE *timing_file = fopen(timing_name.c_str(), "w");
Logger *time_logger = new Logger(LOG_DEBG, timing_file);
@@ -103,18 +105,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
int isam = (int)gconf->get_param("meridional_type");
int nxi = (int)sconf->get_param("number_of_scales");
int idfc = (int)sconf->get_param("idfc");
- double th = gconf->get_param("in_theta_start");
- double thstp = gconf->get_param("in_theta_step");
- double thlst = gconf->get_param("in_theta_end");
- double ths = gconf->get_param("sc_theta_start");
- double thsstp = gconf->get_param("sc_theta_step");
- double thslst = gconf->get_param("sc_theta_end");
- double ph = gconf->get_param("in_phi_start");
- double phstp = gconf->get_param("in_phi_step");
- double phlst = gconf->get_param("in_phi_end");
- double phs = gconf->get_param("sc_phi_start");
- double phsstp = gconf->get_param("sc_phi_step");
- double phslst = gconf->get_param("sc_phi_end");
+ 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");
@@ -216,12 +207,16 @@ void cluster(const string& config_file, const string& data_file, const string& o
fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
fprintf(
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
- th, thstp, thlst, ths, thsstp, thslst
+ p_scattering_angles->th, p_scattering_angles->thstp,
+ p_scattering_angles->thlst, p_scattering_angles->ths,
+ p_scattering_angles->thsstp, p_scattering_angles->thslst
);
fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
fprintf(
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
- ph, phstp, phlst, phs, phsstp, phslst
+ p_scattering_angles->ph, p_scattering_angles->phstp,
+ p_scattering_angles->phlst, p_scattering_angles->phs,
+ p_scattering_angles->phsstp, p_scattering_angles->phslst
);
fprintf(output, " READ(IR,*)JWTM\n");
fprintf(output, " %5d\n", jwtm);
@@ -232,33 +227,6 @@ void cluster(const string& config_file, const string& data_file, const string& o
fprintf(output, " READ(ITIN)NSHL(I),ROS(I)\n");
fprintf(output, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
fprintf(output, " \n");
- double small = 1.0e-3;
- int nth = 0, nph = 0;
- if (thstp != 0.0) nth = int((thlst - th) / thstp + small);
- nth++;
- if (phstp != 0.0) nph = int((phlst - ph) / phstp + small);
- nph++;
- int nths = 0, nphs = 0;
- double thsca = 0.0;
- if (isam > 1) {
- nths = 1;
- thsca = ths - th;
- } else { // ISAM <= 1
- if (thsstp == 0.0) nths = 0;
- else nths = int ((thslst - ths) / thsstp + small);
- nths++;
- }
- if (isam >= 1) {
- nphs = 1;
- } else {
- if (phsstp == 0.0) nphs = 0;
- else nphs = int((phslst - phs) / phsstp + small);
- nphs++;
- }
- int nk = nth * nph;
- int nks = nths * nphs;
- int nkks = nk * nks;
- double th1 = th, ph1 = ph, ths1 = ths, phs1 = phs;
str(sconf, c1, c1ao, c3, c4, c6);
double ****zpv = new double***[c4->lm]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
for (int zi = 0; zi < c4->lm; zi++) {
@@ -285,6 +253,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 nth = p_scattering_angles->nth;
+ int nths = p_scattering_angles->nths;
+ int nph = p_scattering_angles->nph;
+ int nphs = p_scattering_angles->nphs;
tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&inpol), sizeof(int));
@@ -303,7 +275,14 @@ 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;
- 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, isq, ibf, logger);
+ 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);
+ 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
int ompnumthreads = 1;
@@ -362,23 +341,6 @@ void cluster(const string& config_file, const string& data_file, const string& o
int npntts_2 = npntts;
int isam_2 = isam;
int lm_2 = lm;
- double th_2 = th;
- double thstp_2 = thstp;
- double thlst_2 = thlst;
- double ths_2 = ths;
- double thsstp_2 = thsstp;
- double thslst_2 = thslst;
- double ph_2 = ph;
- double phstp_2 = phstp;
- double phlst_2 = phlst;
- double phs_2 = phs;
- double phsstp_2 = phsstp;
- double phslst_2 = phslst;
- double th1_2 = th1;
- double ph1_2 = ph1;
- double ths1_2 = ths1;
- double phs1_2 = phs1;
- double thsca_2 = thsca;
double *u_2 = NULL;
double *us_2 = NULL;
double *un_2 = NULL;
@@ -403,13 +365,6 @@ void cluster(const string& config_file, const string& data_file, const string& o
np_int ndit_2 = ndit;
int isq_2 = isq;
int ibf_2 = ibf;
- int nth_2 = nth;
- int nths_2 = nths;
- int nph_2 = nph;
- int nphs_2 = nph;
- int nk_2 = nk;
- int nks_2 = nks;
- int nkks_2 = nkks;
// 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;
@@ -590,7 +545,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, 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, isq_2, ibf_2, logger);
+ 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);
}
#pragma omp barrier
@@ -732,6 +687,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
delete[] zpv[zi];
}
delete[] zpv;
+ delete p_scattering_angles;
delete c1;
delete c1ao;
delete c2;
@@ -802,8 +758,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
delete sconf;
delete gconf;
chrono::time_point<chrono::high_resolution_clock> t_end = chrono::high_resolution_clock::now();
- const chrono::duration<double> elapsed = t_end - t_start;
- string message = "Calculation lasted " + to_string(elapsed.count()) + ".\n";
+ elapsed = t_end - t_start;
+ message = "INFO: Calculation lasted " + to_string(elapsed.count()) + "s.\n";
logger->log(message);
logger->log("Finished: output written to " + output_path + "/c_OCLU\n");
time_logger->log(message);
@@ -813,7 +769,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
}
-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, 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, 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)
{
logger->log("INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
int jer = 0;
@@ -867,35 +823,11 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
);
if (jer != 0) {
fprintf(output, " STOP IN DME\n");
- // delete[] u;
- // delete[] us;
- // delete[] un;
- // delete[] uns;
- // delete[] up;
- // delete[] ups;
- // delete[] unmp;
- // delete[] unsmp;
- // delete[] upmp;
- // delete[] upsmp;
- // delete[] am_vector;
- // delete[] am;
return jer;
//break;
}
}
if (jer != 0) {
- // delete[] u;
- // delete[] us;
- // delete[] un;
- // delete[] uns;
- // delete[] up;
- // delete[] ups;
- // delete[] unmp;
- // delete[] unsmp;
- // delete[] upmp;
- // delete[] upsmp;
- // delete[] am_vector;
- // delete[] am;
return jer;
//break;
}
@@ -908,16 +840,6 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
cms(am, c1, c1ao, c4, c6);
invert_matrix(am, ndit, jer, mxndm);
if (jer != 0) {
- // delete[] u;
- // delete[] us;
- // delete[] un;
- // delete[] uns;
- // delete[] up;
- // delete[] ups;
- // delete[] unmp;
- // delete[] unsmp;
- // delete[] upmp;
- // delete[] upsmp;
delete[] am_vector;
delete[] am;
return jer;
@@ -992,13 +914,13 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
pcrsm0(vk, exri, inpol, c1, c1ao, c4);
apcra(zpv, c4->le, c1ao->am0m, inpol, sqk, gapm, gappm);
- th = th1;
- for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP portable?
- ph = ph1;
+ double th = sa->th;
+ for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable?
+ double ph = sa->ph;
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
- for (int jph484 = 1; jph484 <= nph; jph484++) {
+ for (int jph484 = 1; jph484 <= sa->nph; jph484++) {
int jw = 0;
- if (nk != 1 || jxi488 <= 1) {
+ if (sa->nk != 1 || jxi488 <= 1) {
upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
if (isam >= 0) {
wmamp(
@@ -1019,12 +941,12 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
}
}
// label 184
- double thsl = ths1;
+ double thsl = sa->ths;
double phsph = 0.0;
- for (int jths = 1; jths <= nths; jths++) {
- ths = thsl;
+ for (int jths = 1; jths <= sa->nths; jths++) {
+ double ths = thsl;
int icspnv = 0;
- if (isam > 1) ths += thsca;
+ if (isam > 1) ths += sa->thsca;
if (isam >= 1) {
phsph = 0.0;
if (ths < 0.0 || ths > 180.0) phsph = 180.0;
@@ -1033,15 +955,15 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
if (phsph != 0.0) icspnv = 1;
}
// label 186
- phs = phs1;
- for (int jphs = 1; jphs <= nphs; jphs++) {
+ double phs = sa->phs;
+ for (int jphs = 1; jphs <= sa->nphs; jphs++) {
double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0;
if (isam >= 1) {
- phs = ph + phsph;
+ phs = sa->ph + phsph;
if (phs > 360.0) phs -= 360.0;
}
// label 188
- bool goto190 = (nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
+ bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
if (!goto190) {
upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
if (isam >= 0)
@@ -1051,7 +973,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
);
}
// label 190
- if (nkks != 1 || jxi488 <= 1) {
+ if (sa->nkks != 1 || jxi488 <= 1) {
upvsp(
u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns,
duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp
@@ -1468,13 +1390,13 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
}
}
// label 420, continues jphs loop
- if (isam < 1) phs += phsstp;
+ if (isam < 1) phs += sa->phsstp;
} // jphs loop, labeled 480
- if (isam <= 1) thsl += thsstp;
+ if (isam <= 1) thsl += sa->thsstp;
} // jths loop, labeled 482
- ph += phstp;
+ ph += sa->phstp;
} // jph484 loop
- th += thstp;
+ th += sa->thstp;
} // jth486 loop
// delete[] u;