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

Add first OpenMP parallelisation in cluster.cpp. YAY!!

parent 825b43d0
Branches
Tags
No related merge requests found
...@@ -8,6 +8,9 @@ ...@@ -8,6 +8,9 @@
#include <exception> #include <exception>
#include <fstream> #include <fstream>
#include <string> #include <string>
#ifdef _OPENMP
#include <omp.h>
#endif
#ifndef INCLUDE_TYPES_H_ #ifndef INCLUDE_TYPES_H_
#include "../include/types.h" #include "../include/types.h"
...@@ -92,7 +95,7 @@ using namespace std; ...@@ -92,7 +95,7 @@ using namespace std;
// double *tqsv; // 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);
/*! \brief C++ implementation of CLU /*! \brief C++ implementation of CLU
* *
...@@ -340,7 +343,11 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -340,7 +343,11 @@ 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 // 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, 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);
for (jxi488 = 2; jxi488 <= nxi; jxi488++) {
int myompthread = 0;
int ompnumthreads = 1;
#pragma omp parallel
{
// 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 // 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
ScattererConfiguration *sconf_2 = new ScattererConfiguration(*sconf); ScattererConfiguration *sconf_2 = new ScattererConfiguration(*sconf);
GeometryConfiguration *gconf_2 = new GeometryConfiguration(*gconf); GeometryConfiguration *gconf_2 = new GeometryConfiguration(*gconf);
...@@ -351,9 +358,13 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -351,9 +358,13 @@ void cluster(string config_file, string data_file, string output_path) {
C4 *c4_2 = new C4(*c4); C4 *c4_2 = new C4(*c4);
C6 *c6_2 = new C6(*c6); C6 *c6_2 = new C6(*c6);
C9 *c9_2 = new C9(*c9); C9 *c9_2 = new C9(*c9);
// FILE *output_2 = fopen((output_path + "/c_OCLU_" + to_string(jxi488)).c_str(), "w"); #ifdef _OPENMP
// fstream tppoan_2; myompthread = omp_get_thread_num();
// tppoan_2.open((output_path + "/c_TPPOAN_" + to_string(jxi488)).c_str(), ios::out | ios::binary); if (myompthread == 0) ompnumthreads = omp_get_num_threads();
#endif
FILE *output_2 = fopen((output_path + "/c_OCLU_" + to_string(myompthread)).c_str(), "w");
fstream tppoan_2;
tppoan_2.open((output_path + "/c_TPPOAN_" + to_string(myompthread)).c_str(), ios::out | ios::binary);
double *gaps_2 = new double[nsph](); double *gaps_2 = new double[nsph]();
double **tqse_2 = new double*[2]; double **tqse_2 = new double*[2];
double **tqce_2 = new double*[2]; double **tqce_2 = new double*[2];
...@@ -440,32 +451,32 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -440,32 +451,32 @@ void cluster(string config_file, string data_file, string output_path) {
tqev_2[ti] = tqev[ti]; tqev_2[ti] = tqev[ti];
tqsv_2[ti] = tqsv[ti]; tqsv_2[ti] = tqsv[ti];
} }
int nxi_2 = nxi; // int nxi_2 = nxi;
int nsph_2 = nsph; // int nsph_2 = nsph;
np_int mxndm_2 = mxndm; // np_int mxndm_2 = mxndm;
int inpol_2 = inpol; // int inpol_2 = inpol;
int iavm_2 = iavm; // int iavm_2 = iavm;
int npnt_2 = npnt; // int npnt_2 = npnt;
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 th_2 = th;
double thstp_2 = thstp; // double thstp_2 = thstp;
double thlst_2 = thlst; // double thlst_2 = thlst;
double ths_2 = ths; // double ths_2 = ths;
double thsstp_2 = thsstp; // double thsstp_2 = thsstp;
double thslst_2 = thslst; // double thslst_2 = thslst;
double ph_2 = ph; // double ph_2 = ph;
double phstp_2 = phstp; // double phstp_2 = phstp;
double phlst_2 = phlst; // double phlst_2 = phlst;
double phs_2 = phs; // double phs_2 = phs;
double phsstp_2 = phsstp; // double phsstp_2 = phsstp;
double phslst_2 = phslst; // double phslst_2 = phslst;
double th1_2 = th1; // double th1_2 = th1;
double ph1_2 = ph1; // double ph1_2 = ph1;
double ths1_2 = ths1; // double ths1_2 = ths1;
double phs1_2 = phs1; // double phs1_2 = phs1;
double thsca_2 = thsca; // double thsca_2 = thsca;
double *u_2 = new double[3](); double *u_2 = new double[3]();
double *us_2 = new double[3](); double *us_2 = new double[3]();
double *un_2 = new double[3](); double *un_2 = new double[3]();
...@@ -488,27 +499,30 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -488,27 +499,30 @@ void cluster(string config_file, string data_file, string output_path) {
upmp_2[ti] = upmp[ti]; upmp_2[ti] = upmp[ti];
upsmp_2[ti] = upsmp[ti]; upsmp_2[ti] = upsmp[ti];
} }
double scan_2 = scan; // double scan_2 = scan;
double cfmp_2 = cfmp; // double cfmp_2 = cfmp;
double sfmp_2 = sfmp; // double sfmp_2 = sfmp;
double cfsp_2 = cfsp; // double cfsp_2 = cfsp;
double sfsp_2 = sfsp; // double sfsp_2 = sfsp;
double sqsfi_2 = sqsfi; // double sqsfi_2 = sqsfi;
double exri_2 = exri; // double exri_2 = exri;
int lcalc_2 = lcalc; // int lcalc_2 = lcalc;
dcomplex arg_2 = arg; // dcomplex arg_2 = arg;
double wn_2 = wn; // double wn_2 = wn;
double vk_2 = vk; // double vk_2 = vk;
np_int ndit_2 = ndit; // np_int ndit_2 = ndit;
dcomplex **am_2 = new dcomplex*[ndit]; dcomplex **am_2 = new dcomplex*[ndit];
dcomplex *am_vector_2 = new dcomplex[ndit * ndit](); dcomplex *am_vector_2 = new dcomplex[ndit * ndit]();
for (int ai = 0; ai < ndit; ai++) { for (int ai = 0; ai < ndit; ai++) {
am_2[ai] = (am_vector_2 + ai * ndit); am_2[ai] = (am_vector_2 + ai * ndit);
} }
int isq_2 = isq; // int isq_2 = isq;
int ibf_2 = ibf; // int ibf_2 = ibf;
jer = cluster_jxi488_cycle(jxi488, sconf_2, gconf_2, c1_2, c1ao_2, c2_2, c3_2, c4_2, c6_2, c9_2, output, output_path, gaps_2, tqse_2, tqspe_2, tqss_2, tqsps_2, zpv_2, gapm_2, gappm_2, nth, nths, nph, nphs, nk, nks, nkks, argi_2, args_2, gap_2, gapp_2, tqce_2, tqcpe_2, tqcs_2, tqcps_2, duk_2, tppoan, 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); #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, nths, nph, nphs, nk, nks, nkks, 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, 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_2, us_2, un_2, uns_2, up_2, ups_2, unmp_2, unsmp_2, upmp_2, upsmp_2, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, am_2, isq, ibf);
}
delete sconf_2; delete sconf_2;
delete gconf_2; delete gconf_2;
...@@ -519,8 +533,8 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -519,8 +533,8 @@ void cluster(string config_file, string data_file, string output_path) {
delete c4_2; delete c4_2;
delete c6_2; delete c6_2;
delete c9_2; delete c9_2;
// fclose(output_2); fclose(output_2);
// tppoan_2.close(); tppoan_2.close();
delete[] gaps_2; delete[] gaps_2;
for (int ti = 0; ti <2 -1; ti++) { for (int ti = 0; ti <2 -1; ti++) {
delete[] tqse_2[ti]; delete[] tqse_2[ti];
...@@ -590,6 +604,9 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -590,6 +604,9 @@ void cluster(string config_file, string data_file, string output_path) {
delete[] am_2; delete[] am_2;
} // jxi488 loop } // jxi488 loop
for (int ri=0; 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
}
tppoan.close(); tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen. } else { // In case TPPOAN could not be opened. Should never happen.
printf("\nERROR: failed to open TPPOAN file.\n"); printf("\nERROR: failed to open TPPOAN file.\n");
...@@ -680,7 +697,7 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -680,7 +697,7 @@ void cluster(string config_file, string data_file, string output_path) {
} }
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)
{ {
// int nxi = sconf->number_of_scales; // int nxi = sconf->number_of_scales;
......
...@@ -53,7 +53,7 @@ class C9; ...@@ -53,7 +53,7 @@ class C9;
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 *, 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);
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.
...@@ -185,7 +185,7 @@ public: ...@@ -185,7 +185,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 *, 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);
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