Skip to content
Commits on Source (2)
Hide whitespace changes
Inline Side-by-side
src/cluster/cluster.cpp
View file @ ffabe9cd
...@@ -53,7 +53,7 @@ using namespace std; ...@@ -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... // 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 /*! \brief C++ implementation of CLU
* *
...@@ -63,6 +63,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf ...@@ -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) { 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::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"; string timing_name = output_path + "/c_timing.log";
FILE *timing_file = fopen(timing_name.c_str(), "w"); FILE *timing_file = fopen(timing_name.c_str(), "w");
Logger *time_logger = new Logger(LOG_DEBG, timing_file); 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 ...@@ -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 isam = (int)gconf->get_param("meridional_type");
int nxi = (int)sconf->get_param("number_of_scales"); int nxi = (int)sconf->get_param("number_of_scales");
int idfc = (int)sconf->get_param("idfc"); int idfc = (int)sconf->get_param("idfc");
double th = gconf->get_param("in_theta_start"); ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);
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");
double wp = sconf->get_param("wp"); double wp = sconf->get_param("wp");
// Global variables for CLU // Global variables for CLU
int lm = (int)gconf->get_param("l_max"); 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 ...@@ -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, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
fprintf( fprintf(
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n", 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, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
fprintf( fprintf(
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n", 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, " READ(IR,*)JWTM\n");
fprintf(output, " %5d\n", jwtm); fprintf(output, " %5d\n", jwtm);
...@@ -232,33 +227,6 @@ void cluster(const string& config_file, const string& data_file, const string& o ...@@ -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)NSHL(I),ROS(I)\n");
fprintf(output, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n"); fprintf(output, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
fprintf(output, " \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); 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] 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++) { 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 ...@@ -285,6 +253,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
#else #else
logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO); logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);
#endif #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 *>(&iavm), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int)); tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&inpol), 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 ...@@ -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 // do the first iteration on jxi488 separately, since it seems to be different from the others
int jxi488 = 1; 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 // 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; int ompnumthreads = 1;
...@@ -362,23 +341,6 @@ void cluster(const string& config_file, const string& data_file, const string& o ...@@ -362,23 +341,6 @@ void cluster(const string& config_file, const string& data_file, const string& o
int npntts_2 = npntts; int npntts_2 = npntts;
int isam_2 = isam; int isam_2 = isam;
int lm_2 = lm; 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 *u_2 = NULL;
double *us_2 = NULL; double *us_2 = NULL;
double *un_2 = NULL; double *un_2 = NULL;
...@@ -403,13 +365,6 @@ void cluster(const string& config_file, const string& data_file, const string& o ...@@ -403,13 +365,6 @@ void cluster(const string& config_file, const string& data_file, const string& o
np_int ndit_2 = ndit; np_int ndit_2 = ndit;
int isq_2 = isq; int isq_2 = isq;
int ibf_2 = ibf; 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 // 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) { if (myompthread == 0) {
sconf_2 = sconf; sconf_2 = sconf;
...@@ -590,7 +545,7 @@ void cluster(const string& config_file, const string& data_file, const string& o ...@@ -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 // ok, now I can actually start the parallel calculations
#pragma omp for #pragma omp for
for (jxi488 = 2; jxi488 <= nxi; jxi488++) { 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 #pragma omp barrier
...@@ -732,6 +687,7 @@ void cluster(const string& config_file, const string& data_file, const string& o ...@@ -732,6 +687,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
delete[] zpv[zi]; delete[] zpv[zi];
} }
delete[] zpv; delete[] zpv;
delete p_scattering_angles;
delete c1; delete c1;
delete c1ao; delete c1ao;
delete c2; delete c2;
...@@ -802,8 +758,8 @@ void cluster(const string& config_file, const string& data_file, const string& o ...@@ -802,8 +758,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
delete sconf; delete sconf;
delete gconf; delete gconf;
chrono::time_point<chrono::high_resolution_clock> t_end = chrono::high_resolution_clock::now(); chrono::time_point<chrono::high_resolution_clock> t_end = chrono::high_resolution_clock::now();
const chrono::duration<double> elapsed = t_end - t_start; elapsed = t_end - t_start;
string message = "Calculation lasted " + to_string(elapsed.count()) + ".\n"; message = "INFO: Calculation lasted " + to_string(elapsed.count()) + "s.\n";
logger->log(message); logger->log(message);
logger->log("Finished: output written to " + output_path + "/c_OCLU\n"); logger->log("Finished: output written to " + output_path + "/c_OCLU\n");
time_logger->log(message); time_logger->log(message);
...@@ -813,7 +769,7 @@ void cluster(const string& config_file, const string& data_file, const string& o ...@@ -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"); logger->log("INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
int jer = 0; int jer = 0;
...@@ -867,35 +823,11 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf ...@@ -867,35 +823,11 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
); );
if (jer != 0) { if (jer != 0) {
fprintf(output, " STOP IN DME\n"); 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; return jer;
//break; //break;
} }
} }
if (jer != 0) { 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; return jer;
//break; //break;
} }
...@@ -908,16 +840,6 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf ...@@ -908,16 +840,6 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
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) { 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_vector;
delete[] am; delete[] am;
return jer; return jer;
...@@ -992,13 +914,13 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf ...@@ -992,13 +914,13 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double)); tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
pcrsm0(vk, exri, inpol, c1, c1ao, c4); pcrsm0(vk, exri, inpol, c1, c1ao, c4);
apcra(zpv, c4->le, c1ao->am0m, inpol, sqk, gapm, gappm); apcra(zpv, c4->le, c1ao->am0m, inpol, sqk, gapm, gappm);
th = th1; double th = sa->th;
for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP portable? for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable?
ph = ph1; double ph = sa->ph;
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0; 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; 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); upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
if (isam >= 0) { if (isam >= 0) {
wmamp( wmamp(
...@@ -1019,12 +941,12 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf ...@@ -1019,12 +941,12 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
} }
} }
// label 184 // label 184
double thsl = ths1; double thsl = sa->ths;
double phsph = 0.0; double phsph = 0.0;
for (int jths = 1; jths <= nths; jths++) { for (int jths = 1; jths <= sa->nths; jths++) {
ths = thsl; double ths = thsl;
int icspnv = 0; int icspnv = 0;
if (isam > 1) ths += thsca; if (isam > 1) ths += sa->thsca;
if (isam >= 1) { if (isam >= 1) {
phsph = 0.0; phsph = 0.0;
if (ths < 0.0 || ths > 180.0) phsph = 180.0; if (ths < 0.0 || ths > 180.0) phsph = 180.0;
...@@ -1033,15 +955,15 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf ...@@ -1033,15 +955,15 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
if (phsph != 0.0) icspnv = 1; if (phsph != 0.0) icspnv = 1;
} }
// label 186 // label 186
phs = phs1; double phs = sa->phs;
for (int jphs = 1; jphs <= nphs; jphs++) { for (int jphs = 1; jphs <= sa->nphs; jphs++) {
double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0; double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0;
if (isam >= 1) { if (isam >= 1) {
phs = ph + phsph; phs = sa->ph + phsph;
if (phs > 360.0) phs -= 360.0; if (phs > 360.0) phs -= 360.0;
} }
// label 188 // 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) { if (!goto190) {
upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp); upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
if (isam >= 0) if (isam >= 0)
...@@ -1051,7 +973,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf ...@@ -1051,7 +973,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
); );
} }
// label 190 // label 190
if (nkks != 1 || jxi488 <= 1) { if (sa->nkks != 1 || jxi488 <= 1) {
upvsp( upvsp(
u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns, u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns,
duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp
...@@ -1468,13 +1390,13 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf ...@@ -1468,13 +1390,13 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
} }
} }
// label 420, continues jphs loop // label 420, continues jphs loop
if (isam < 1) phs += phsstp; if (isam < 1) phs += sa->phsstp;
} // jphs loop, labeled 480 } // jphs loop, labeled 480
if (isam <= 1) thsl += thsstp; if (isam <= 1) thsl += sa->thsstp;
} // jths loop, labeled 482 } // jths loop, labeled 482
ph += phstp; ph += sa->phstp;
} // jph484 loop } // jph484 loop
th += thstp; th += sa->thstp;
} // jth486 loop } // jth486 loop
// delete[] u; // delete[] u;
......
src/include/Commons.h
View file @ ffabe9cd
...@@ -385,4 +385,106 @@ public: ...@@ -385,4 +385,106 @@ public:
*/ */
~C9(); ~C9();
}; };
/*! \brief A data structure representing the angles to be evaluated in the problem.
*
*/
class ScatteringAngles {
protected:
//! \brief Number of incident field azimuth angles.
int _nth;
//! \brief Number of scattered field azimuth angles.
int _nths;
//! \brief Number of incident field elevation angles.
int _nph;
//! \brief Number of scattered field elevation angles.
int _nphs;
//! \brief Number of incident field propagation angles.
int _nk;
//! \brief Number of scattered field propagation angles.
int _nks;
//! \brief Total number of field propagation angles.
int _nkks;
//! \brief First incident field azimuth angle.
double _th;
//! \brief Incident field azimuth angle increment.
double _thstp;
//! \brief Last incident field azimuth angle.
double _thlst;
//! \brief First scattered field azimuth angle.
double _ths;
//! \brief Scattered field azimuth angle increment.
double _thsstp;
//! \brief Last scattered field azimuth angle.
double _thslst;
//! \brief First incident field elevation angle.
double _ph;
//! \brief Incident field elevation angle increment.
double _phstp;
//! \brief Last incident field elevation angle.
double _phlst;
//! \brief First scattered field elevation angle.
double _phs;
//! \brief Scattered field elevation angle increment.
double _phsstp;
//! \brief Last scattered field elevation angle.
double _phslst;
//! \brief Azimuth scattering deflection.
double _thsca;
public:
//! \brief Read only view of `_nth`.
const int& nth = _nth;
//! \brief Read only view of `_nths`.
const int& nths = _nths;
//! \brief Read only view of `_nph`.
const int& nph = _nph;
//! \brief Read only view of `_nphs`.
const int& nphs = _nphs;
//! \brief Read only view of `_nk`.
const int& nk = _nk;
//! \brief Read only view of `_nks`.
const int& nks = _nks;
//! \brief Read only view of `_nkks`.
const int& nkks = _nkks;
//! \brief Read only view of `_th`.
const double& th = _th;
//! \brief Read only view of `_thstp`.
const double& thstp = _thstp;
//! \brief Read only view of `_thlst`.
const double& thlst = _thlst;
//! \brief Read only view of `_ths`.
const double& ths = _ths;
//! \brief Read only view of `_thsstp`.
const double& thsstp = _thsstp;
//! \brief Read only view of `_thslst`.
const double& thslst = _thslst;
//! \brief Read only view of `_ph`.
const double& ph = _ph;
//! \brief Read only view of `_phstp`.
const double& phstp = _phstp;
//! \brief Read only view of `_phlst`.
const double& phlst = _phlst;
//! \brief Read only view of `_phs`.
const double& phs = _phs;
//! \brief Read only view of `_phsstp`.
const double& phsstp = _phsstp;
//! \brief Read only view of `_phslst`.
const double& phslst = _phslst;
//! \brief Read only view of `_thsca`.
const double& thsca = _thsca;
/*! \brief ScatteringAngles instance constructor.
*
* \param gconf: `GeometryConfiguration*` Pointer to a GeometryConfiguration object.
*/
ScatteringAngles(GeometryConfiguration *gconf);
/*! \brief ScatteringAngles copy constructor.
*
* \param rhs: `ScatteringAngles&` Reference to the ScatteringAngles object to be copied.
*/
ScatteringAngles(const ScatteringAngles &rhs);
};
#endif #endif
src/libnptm/Commons.cpp
View file @ ffabe9cd
...@@ -2,15 +2,7 @@ ...@@ -2,15 +2,7 @@
/*! \file Commons.cpp /*! \file Commons.cpp
* *
* DEVELOPMENT NOTE: * \brief Implementation of the common data structures.
* The construction of common blocks requires some information
* that is stored in configuration objects and is needed to compute
* the allocation size of vectors and matrices. Currently, this
* information is passed as arguments to the constructors of the
* common blocks. A simpler and more logical way to operate would
* be to design the constructors to take as arguments only pointers
* to the configuration objects. These, on their turn, need to
* expose methods to access the relevant data in read-only mode.
*/ */
#ifndef INCLUDE_TYPES_H_ #ifndef INCLUDE_TYPES_H_
#include "../include/types.h" #include "../include/types.h"
...@@ -536,3 +528,70 @@ C9::~C9() { ...@@ -536,3 +528,70 @@ C9::~C9() {
for (int si = sam_size_0 - 1; si > -1; si--) delete[] sam[si]; for (int si = sam_size_0 - 1; si > -1; si--) delete[] sam[si];
delete[] sam; delete[] sam;
} }
ScatteringAngles::ScatteringAngles(GeometryConfiguration *gconf) {
int isam = (int)gconf->get_param("meridional_type");
_th = gconf->get_param("in_theta_start");
_thstp = gconf->get_param("in_theta_step");
_thlst = gconf->get_param("in_theta_end");
_ths = gconf->get_param("sc_theta_start");
_thsstp = gconf->get_param("sc_theta_step");
_thslst = gconf->get_param("sc_theta_end");
_ph = gconf->get_param("in_phi_start");
_phstp = gconf->get_param("in_phi_step");
_phlst = gconf->get_param("in_phi_end");
_phs = gconf->get_param("sc_phi_start");
_phsstp = gconf->get_param("sc_phi_step");
_phslst = gconf->get_param("sc_phi_end");
double small = 1.0e-3;
_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++;
_nths = 0;
_nphs = 0;
_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++;
}
_nk = nth * nph;
_nks = nths * nphs;
_nkks = nk * nks;
}
ScatteringAngles::ScatteringAngles(const ScatteringAngles &rhs) {
_th = rhs._th;
_thstp = rhs._thstp;
_thlst = rhs._thlst;
_ths = rhs._ths;
_thsstp = rhs._thsstp;
_thslst = rhs._thslst;
_ph = rhs._ph;
_phstp = rhs._phstp;
_phlst = rhs._phlst;
_phs = rhs._phs;
_phsstp = rhs._phsstp;
_phslst = rhs._phslst;
_nth = rhs._nth;
_nph = rhs._nph;
_nths = rhs._nths;
_nphs = rhs._nphs;
_thsca = rhs._thsca;
_nk = rhs._nk;
_nks = rhs._nks;
_nkks = rhs._nkks;
}