Skip to content
Snippets Groups Projects
Commit ef512fd2 authored by Mulas, Giacomo's avatar Mulas, Giacomo
Browse files

Move the content of the loop over scales in a separate function, so it can be...

Move the content of the loop over scales in a separate function, so it can be called with fresh arguments every time to check parallelism and dependencies.
parent cb809b82
Branches
Tags
No related merge requests found
/* Distributed under the terms of GPLv3 or later. See COPYING for details. */ /* Distributed under the terms of GPLv3 or later. See COPYING for details. */
/*! \file cluster.cpp /*! \file cluster.cp
* *
* \brief Implementation of the calculation for a cluster of spheres. * \brief Implementation of the calculation for a cluster of spheres.
*/ */
...@@ -43,6 +43,8 @@ ...@@ -43,6 +43,8 @@
using namespace std; using namespace std;
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);
/*! \brief C++ implementation of CLU /*! \brief C++ implementation of CLU
* *
* \param config_file: `string` Name of the configuration file. * \param config_file: `string` Name of the configuration file.
...@@ -112,20 +114,20 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -112,20 +114,20 @@ void cluster(string config_file, string data_file, string output_path) {
// End of add-ons initialization // End of add-ons initialization
C6 *c6 = new C6(c4->lmtpo); C6 *c6 = new C6(c4->lmtpo);
FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w"); FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
int jer = 0, lcalc = 0; int jer = 0;
dcomplex arg = 0.0 + 0.0 * I; // int lcalc = 0;
dcomplex ccsam = 0.0 + 0.0 * I; dcomplex ccsam = 0.0 + 0.0 * I;
int configurations = 0; int configurations = 0;
for (int ci = 1; ci <= nsph; ci++) { for (int ci = 1; ci <= nsph; ci++) {
if (sconf->iog_vec[ci -1] >= ci) configurations++; if (sconf->iog_vec[ci -1] >= ci) configurations++;
} }
C2 *c2 = new C2(nsph, configurations, npnt, npntts); C2 *c2 = new C2(nsph, configurations, npnt, npntts);
np_int ndit = 2 * nsph * c4->nlim; // np_int ndit = 2 * nsph * c4->nlim;
dcomplex *am_vector = new dcomplex[ndit * ndit](); // dcomplex *am_vector = new dcomplex[ndit * ndit]();
dcomplex **am = new dcomplex*[ndit]; // dcomplex **am = new dcomplex*[ndit];
for (int ai = 0; ai < ndit; ai++) { // for (int ai = 0; ai < ndit; ai++) {
am[ai] = (am_vector + ai * ndit); // am[ai] = (am_vector + ai * ndit);
} // }
const int ndi = c4->nsph * c4->nlim; const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem); C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
double *gaps = new double[nsph](); double *gaps = new double[nsph]();
...@@ -164,16 +166,6 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -164,16 +166,6 @@ void cluster(string config_file, string data_file, string output_path) {
gap[gi] = new double[2](); gap[gi] = new double[2]();
gapm[gi] = new double[2](); gapm[gi] = new double[2]();
} }
double *u = new double[3]();
double *us = new double[3]();
double *un = new double[3]();
double *uns = new double[3]();
double *up = new double[3]();
double *ups = new double[3]();
double *unmp = new double[3]();
double *unsmp = new double[3]();
double *upmp = new double[3]();
double *upsmp = new double[3]();
double *argi = new double[1](); double *argi = new double[1]();
double *args = new double[1](); double *args = new double[1]();
double *duk = new double[3](); double *duk = new double[3]();
...@@ -189,8 +181,8 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -189,8 +181,8 @@ void cluster(string config_file, string data_file, string output_path) {
cmullr[ci] = new double[4](); cmullr[ci] = new double[4]();
cmul[ci] = new double[4](); cmul[ci] = new double[4]();
} }
int isq, ibf; //int isq, ibf;
double scan, cfmp, sfmp, cfsp, sfsp; //double scan, cfmp, sfmp, cfsp, sfsp;
// End of global variables for CLU // End of global variables for CLU
fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n"); 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", fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
...@@ -279,13 +271,160 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -279,13 +271,160 @@ void cluster(string config_file, string data_file, string output_path) {
tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int)); tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int)); tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
double wn = wp / 3.0e8; double wn = wp / 3.0e8;
double sqsfi = 1.0;
if (sconf->idfc < 0) { if (sconf->idfc < 0) {
vk = sconf->xip * wn; vk = sconf->xip * wn;
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk); fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n"); fprintf(output, " \n");
} }
for (int jxi488 = 1; jxi488 <= nxi; jxi488++) { for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
// 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
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);
} // jxi488 loop
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
printf("\nERROR: failed to open TPPOAN file.\n");
}
fclose(output);
// Clean memory
delete c1;
delete c1ao;
delete c3;
delete c4;
delete c6;
delete c9;
for (int zi = c4->lm - 1; zi > -1; zi--) {
for (int zj = 2; zj > -1; zj--) {
delete[] zpv[zi][zj][1];
delete[] zpv[zi][zj][0];
delete[] zpv[zi][zj];
}
delete[] zpv[zi];
}
delete[] zpv;
//delete[] am_vector;
//delete[] am;
//delete[] tam;
delete[] gaps;
for (int ti = 1; ti > -1; ti--) {
delete[] tqse[ti];
delete[] tqss[ti];
delete[] tqspe[ti];
delete[] tqsps[ti];
delete[] tqce[ti];
delete[] tqcpe[ti];
delete[] tqcs[ti];
delete[] tqcps[ti];
}
delete[] tqse;
delete[] tqss;
delete[] tqspe;
delete[] tqsps;
delete[] tqce;
delete[] tqcpe;
delete[] tqcs;
delete[] tqcps;
delete[] tqev;
delete[] tqsv;
for (int gi = 2; gi > -1; gi--) {
delete[] gapp[gi];
delete[] gappm[gi];
delete[] gap[gi];
delete[] gapm[gi];
}
delete[] gapp;
delete[] gappm;
delete[] gap;
delete[] gapm;
delete[] gapv;
delete[] argi;
delete[] args;
delete[] duk;
for (int ci = 3; ci > -1; ci--) {
delete[] cextlr[ci];
delete[] cext[ci];
delete[] cmullr[ci];
delete[] cmul[ci];
}
delete[] cextlr;
delete[] cext;
delete[] cmullr;
delete[] cmul;
} else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations."
);
}
delete sconf;
delete gconf;
printf("Finished: output written to %s.\n", (output_path + "/c_OCLU").c_str());
}
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 = sconf->number_of_scales;
int nsph = gconf->number_of_spheres;
np_int mxndm = (np_int)gconf->mxndm;
int inpol = gconf->in_pol;
int iavm = gconf->iavm;
int npnt = gconf->npnt;
int npntts = gconf->npntts;
int isam = gconf->meridional_type;
int lm = gconf->l_max;
double th = gconf->in_theta_start;
double thstp = gconf->in_theta_step;
double thlst = gconf->in_theta_end;
double ths = gconf->sc_theta_start;
double thsstp = gconf->sc_theta_step;
double thslst = gconf->sc_theta_end;
double ph = gconf->in_phi_start;
double phstp = gconf->in_phi_step;
double phlst = gconf->in_phi_end;
double phs = gconf->sc_phi_start;
double phsstp = gconf->sc_phi_step;
double phslst = gconf->sc_phi_end;
double th1 = th;
double ph1 = ph;
double ths1 = ths;
double phs1 = phs;
double thsca = 0.0;
if (isam > 1) thsca = ths - th;
double *u = new double[3]();
double *us = new double[3]();
double *un = new double[3]();
double *uns = new double[3]();
double *up = new double[3]();
double *ups = new double[3]();
double *unmp = new double[3]();
double *unsmp = new double[3]();
double *upmp = new double[3]();
double *upsmp = new double[3]();
double scan;
double cfmp;
double sfmp;
double cfsp;
double sfsp;
double sqsfi = 1.0;
double exri = sqrt(sconf->exdc);
int jer = 0;
int lcalc = 0;
dcomplex arg = 0.0 + 0.0 * I;
double wp = sconf->wp;
double wn = wp / 3.0e8;
double vk = 0.0;
if (sconf->idfc < 0) vk = sconf->xip * wn;
np_int ndit = 2 * nsph * c4->nlim;
dcomplex *am_vector = new dcomplex[ndit * ndit]();
dcomplex **am = new dcomplex*[ndit];
for (int ai = 0; ai < ndit; ai++) {
am[ai] = (am_vector + ai * ndit);
}
int isq;
int ibf;
printf("INFO: running scale iteration %d of %d...", jxi488, nxi); printf("INFO: running scale iteration %d of %d...", jxi488, nxi);
int jaw = 1; int jaw = 1;
fprintf(output, "========== JXI =%3d ====================\n", jxi488); fprintf(output, "========== JXI =%3d ====================\n", jxi488);
...@@ -303,7 +442,20 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -303,7 +442,20 @@ void cluster(string config_file, string data_file, string output_path) {
hjv(exri, vkarg, jer, lcalc, arg, c1, c1ao, c4); hjv(exri, vkarg, jer, lcalc, arg, c1, c1ao, c4);
if (jer != 0) { if (jer != 0) {
fprintf(output, " STOP IN HJV\n"); fprintf(output, " STOP IN HJV\n");
break; // jxi488 loop: goes to memory cleaning and return 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; // rewrite this to go to the end of the function, to free locally allocated variables and return jer
} }
for (int i132 = 1; i132 <= nsph; i132++) { for (int i132 = 1; i132 <= nsph; i132++) {
int iogi = c1->iog[i132 - 1]; int iogi = c1->iog[i132 - 1];
...@@ -331,14 +483,57 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -331,14 +483,57 @@ void cluster(string config_file, string data_file, string output_path) {
); );
if (jer != 0) { if (jer != 0) {
fprintf(output, " STOP IN DME\n"); fprintf(output, " STOP IN DME\n");
break; 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) 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;
}
} // i132 loop } // i132 loop
cms(am, c1, c1ao, c4, c6); cms(am, c1, c1ao, c4, c6);
invert_matrix(am, ndit, jer, mxndm); invert_matrix(am, ndit, jer, mxndm);
if (jer != 0) break; // jxi488 loop: goes to memory clean 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; // jxi488 loop: goes to memory clean
}
ztm(am, c1, c1ao, c4, c6, c9); ztm(am, c1, c1ao, c4, c6, c9);
if (sconf->idfc >= 0) { if (sconf->idfc >= 0) {
if (jxi488 == gconf->jwtm) { if (jxi488 == gconf->jwtm) {
...@@ -891,64 +1086,7 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -891,64 +1086,7 @@ void cluster(string config_file, string data_file, string output_path) {
} // jph484 loop } // jph484 loop
th += thstp; th += thstp;
} // jth486 loop } // jth486 loop
printf(" done.\n");
} // jxi488 loop
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
printf("\nERROR: failed to open TPPOAN file.\n");
}
fclose(output);
// Clean memory
delete c1;
delete c1ao;
delete c3;
delete c4;
delete c6;
delete c9;
for (int zi = c4->lm - 1; zi > -1; zi--) {
for (int zj = 2; zj > -1; zj--) {
delete[] zpv[zi][zj][1];
delete[] zpv[zi][zj][0];
delete[] zpv[zi][zj];
}
delete[] zpv[zi];
}
delete[] zpv;
delete[] am_vector;
delete[] am;
//delete[] tam;
delete[] gaps;
for (int ti = 1; ti > -1; ti--) {
delete[] tqse[ti];
delete[] tqss[ti];
delete[] tqspe[ti];
delete[] tqsps[ti];
delete[] tqce[ti];
delete[] tqcpe[ti];
delete[] tqcs[ti];
delete[] tqcps[ti];
}
delete[] tqse;
delete[] tqss;
delete[] tqspe;
delete[] tqsps;
delete[] tqce;
delete[] tqcpe;
delete[] tqcs;
delete[] tqcps;
delete[] tqev;
delete[] tqsv;
for (int gi = 2; gi > -1; gi--) {
delete[] gapp[gi];
delete[] gappm[gi];
delete[] gap[gi];
delete[] gapm[gi];
}
delete[] gapp;
delete[] gappm;
delete[] gap;
delete[] gapm;
delete[] gapv;
delete[] u; delete[] u;
delete[] us; delete[] us;
delete[] un; delete[] un;
...@@ -959,25 +1097,11 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -959,25 +1097,11 @@ void cluster(string config_file, string data_file, string output_path) {
delete[] unsmp; delete[] unsmp;
delete[] upmp; delete[] upmp;
delete[] upsmp; delete[] upsmp;
delete[] argi; delete[] am_vector;
delete[] args; delete[] am;
delete[] duk;
for (int ci = 3; ci > -1; ci--) { printf(" done.\n");
delete[] cextlr[ci];
delete[] cext[ci]; return jer;
delete[] cmullr[ci];
delete[] cmul[ci];
}
delete[] cextlr;
delete[] cext;
delete[] cmullr;
delete[] cmul;
} else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations."
);
}
delete sconf;
delete gconf;
printf("Finished: output written to %s.\n", (output_path + "/c_OCLU").c_str());
} }
...@@ -30,6 +30,18 @@ ...@@ -30,6 +30,18 @@
#ifndef INCLUDE_CONFIGURATION_H_ #ifndef INCLUDE_CONFIGURATION_H_
#define INCLUDE_CONFIGURATION_H_ #define INCLUDE_CONFIGURATION_H_
class ScattererConfiguration;
class GeometryConfiguration;
class C1;
class C1_AddOns;
class C2;
class C3;
class C4;
class C6;
class C9;
#include <fstream>
/** /**
* \brief A class to represent the configuration of the scattering geometry. * \brief A class to represent the configuration of the scattering geometry.
* *
...@@ -41,6 +53,7 @@ ...@@ -41,6 +53,7 @@
class GeometryConfiguration { class GeometryConfiguration {
//! Temporary work-around to allow cluster() and sphere() peeking in. //! Temporary work-around to allow cluster() and sphere() peeking in.
friend void cluster(std::string, std::string, std::string); 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 *);
friend void sphere(std::string, std::string, std::string); friend void sphere(std::string, std::string, std::string);
protected: protected:
//! \brief Number of spherical components. //! \brief Number of spherical components.
...@@ -167,6 +180,7 @@ public: ...@@ -167,6 +180,7 @@ public:
class ScattererConfiguration { class ScattererConfiguration {
//! Temporary work-around to allow cluster() and sphere() peeking in. //! Temporary work-around to allow cluster() and sphere() peeking in.
friend void cluster(std::string, std::string, std::string); 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 *);
friend void sphere(std::string, std::string, std::string); friend void sphere(std::string, std::string, std::string);
protected: protected:
//! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES]. //! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment