diff --git a/src/Makefile b/src/Makefile
index dd65769a321c70e4883df53b2e981f74936ee568..e6032377ec1b3c6a829d90d64967579767de7bf2 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -8,9 +8,13 @@ ifndef BUILDDIR_NPTM
override BUILDDIR_NPTM=$(BUILDDIR)/libnptm
endif
ifndef LIBNPTM
-# choose one of the two following lines, depending on whether a static or dynamic libnptm is wanted
-#override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.a
+ifdef STATIC_NPTM
+override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.a
+override STATICFLAG="-lsz -lz -laec -static"
+else
override LIBNPTM=$(BUILDDIR_NPTM)/libnptm.so
+override STATICFLAG=""
+endif
endif
DOCSDIR=$(SRCDIR)/../doc
@@ -28,7 +32,7 @@ $(LIBNPTM):
BUILDDIR=$(BUILDDIR) BUILDDIR_NPTM=$(BUILDDIR_NPTM) LIBNPTM=$(LIBNPTM) $(MAKE) -C libnptm $@
$(SUBDIRS): $(BUILDDIR) $(LIBNPTM)
- BUILDDIR=$(BUILDDIR) BUILDDIR_NPTM=$(BUILDDIR_NPTM) LIBNPTM=$(LIBNPTM) $(MAKE) -C $@
+ BUILDDIR=$(BUILDDIR) BUILDDIR_NPTM=$(BUILDDIR_NPTM) LIBNPTM=$(LIBNPTM) STATICFLAG=$(STATICFLAG) $(MAKE) -C $@
clean: $(BUILDDIR)
BUILDDIR=$(BUILDDIR) $(MAKE) -C cluster clean
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 93737f22beedbeb90e70dbe8f0ce92094abb2cbf..8195ce0932a685ad0f41b5e01aec0058a632add7 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -20,6 +20,10 @@
#include "../include/errors.h"
#endif
+#ifndef INCLUDE_LOGGING_H_
+#include "../include/logging.h"
+#endif
+
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
@@ -48,7 +52,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, 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
*
@@ -57,13 +61,15 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
* \param output_path: `string` Directory to write the output files in.
*/
void cluster(string config_file, string data_file, string output_path) {
- printf("INFO: making legacy configuration...");
+ Logger *logger = new Logger(LOG_INFO);
+ logger->log("INFO: making legacy configuration...", LOG_INFO);
ScattererConfiguration *sconf = NULL;
try {
sconf = ScattererConfiguration::from_dedfb(config_file);
} catch(const OpenConfigurationFileException &ex) {
- printf("\nERROR: failed to open scatterer configuration file.\n");
- printf("FILE: %s\n", ex.what());
+ logger->err("\nERROR: failed to open scatterer configuration file.\n");
+ string message = "FILE: " + string(ex.what()) + "\n";
+ logger->err(message);
exit(1);
}
sconf->write_formatted(output_path + "/c_OEDFB");
@@ -73,12 +79,13 @@ void cluster(string config_file, string data_file, string output_path) {
try {
gconf = GeometryConfiguration::from_legacy(data_file);
} catch (const OpenConfigurationFileException &ex) {
- printf("\nERROR: failed to open geometry configuration file.\n");
- printf("FILE: %s\n", ex.what());
+ logger->err("\nERROR: failed to open geometry configuration file.\n");
+ string message = "FILE: " + string(ex.what()) + "\n";
+ logger->err(message);
if (sconf) delete sconf;
exit(1);
}
- printf(" done.\n");
+ logger->log(" done.\n", LOG_INFO);
if (sconf->number_of_spheres == gconf->number_of_spheres) {
// Shortcuts to variables stored in configuration objects
int nsph = gconf->number_of_spheres;
@@ -278,9 +285,9 @@ void cluster(string config_file, string data_file, string output_path) {
tppoan.open(tppoan_name.c_str(), ios::out | ios::binary);
if (tppoan.is_open()) {
#ifdef USE_LAPACK
- printf("INFO: using LAPACK calls.\n");
+ logger->log("INFO: using LAPACK calls.\n", LOG_INFO);
#else
- printf("INFO: using fall-back lucin() calls.\n");
+ logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);
#endif
tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
@@ -299,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;
@@ -594,7 +601,7 @@ void cluster(string config_file, string data_file, string output_path) {
// 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, 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
@@ -613,7 +620,7 @@ void cluster(string config_file, string data_file, string output_path) {
tppoanp_2->close();
delete tppoanp_2;
delete[] gaps_2;
- for (int ti = 0; ti <2 -1; ti++) {
+ for (int ti = 0; ti < 2; ti++) {
delete[] tqse_2[ti];
delete[] tqce_2[ti];
delete[] tqcs_2[ti];
@@ -682,7 +689,8 @@ void cluster(string config_file, string data_file, string output_path) {
}
#pragma omp barrier
{
- printf("Closed thread-local output files of thread %d, syncing threads\n", myompthread);
+ string message = "Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
+ logger->log(message);
}
} // closes pragma omp parallel
#ifdef _OPENMP
@@ -692,7 +700,8 @@ void cluster(string config_file, string data_file, string output_path) {
for (int ri = 1; 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
string partial_file_name = output_path + "/c_OCLU_" + to_string(ri);
- printf("Copying ASCII output of thread %d of %d... ", ri, ompnumthreads - 1);
+ string message = "Copying ASCII output of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
+ logger->log(message);
FILE *partial_output = fopen(partial_file_name.c_str(), "r");
char c = fgetc(partial_output);
while (c != EOF) {
@@ -701,12 +710,10 @@ void cluster(string config_file, string data_file, string output_path) {
}
fclose(partial_output);
remove(partial_file_name.c_str());
- printf("done.\n");
- // }
- //for (int ri = 1; ri < ompnumthreads; ri++) {
- //string partial_file_name = output_path + "/c_TPPOAN_" + to_string(ri);
+ logger->log("done.\n");
partial_file_name = output_path + "/c_TPPOAN_" + to_string(ri);
- printf("Copying binary output of thread %d of %d... ", ri, ompnumthreads - 1);
+ message = "Copying binary output of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
+ logger->log(message);
fstream partial_tppoan;
partial_tppoan.open(partial_file_name.c_str(), ios::in | ios::binary);
partial_tppoan.seekg(0, ios::end);
@@ -718,14 +725,14 @@ void cluster(string config_file, string data_file, string output_path) {
partial_tppoan.close();
delete[] binary_buffer;
remove(partial_file_name.c_str());
- printf("done.\n");
+ logger->log("done.\n");
}
}
#endif
tppoanp->close();
delete tppoanp;
} else { // In case TPPOAN could not be opened. Should never happen.
- printf("\nERROR: failed to open TPPOAN file.\n");
+ logger->err("\nERROR: failed to open TPPOAN file.\n");
}
fclose(output);
// Clean memory
@@ -740,13 +747,13 @@ void cluster(string config_file, string data_file, string output_path) {
delete[] zpv;
delete c1;
delete c1ao;
+ delete c2;
delete c3;
delete c4;
delete c6;
delete c9;
delete[] am_vector;
delete[] am;
- //delete[] tam;
delete[] gaps;
for (int ti = 1; ti > -1; ti--) {
delete[] tqse[ti];
@@ -809,11 +816,12 @@ void cluster(string config_file, string data_file, string output_path) {
}
delete sconf;
delete gconf;
- printf("Finished: output written to %s.\n", (output_path + "/c_OCLU").c_str());
+ logger->log("Finished: output written to " + output_path + "/c_OCLU\n");
+ delete logger;
}
-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;
@@ -877,8 +885,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
// int ibf;
-
- printf("INFO: running scale iteration %d of %d...", jxi488, nxi);
+ logger->log("INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
int jaw = 1;
fprintf(output, "========== JXI =%3d ====================\n", jxi488);
double xi = sconf->scale_vec[jxi488 - 1];
@@ -1272,10 +1279,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(
@@ -1305,7 +1312,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.
@@ -1370,9 +1377,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++) {
@@ -1553,7 +1560,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
// delete[] am_vector;
// delete[] am;
- printf(" done.\n");
+ logger->log("INFO: finished scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
return jer;
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 4c0da316e7b0b43a583f054f62b1cc5d825221b9..90a911ebaf3c08e49e8d30b79070e46a27c34714 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -40,9 +40,18 @@ protected:
int nsph;
//! \brief Maximum order of field expansion.
int lm;
- //! \brief NLMMT. QUESTION: definition?
- int nlmmt;
+ //! \brief Contiguous RMI vector.
+ dcomplex *vec_rmi;
+ //! \brief Contiguous REI vector.
+ dcomplex *vec_rei;
+ //! \brief Contiguous W vector.
+ dcomplex *vec_w;
+ //! \brief Contiguous VINTS vector
+ dcomplex *vec_vints;
+
public:
+ //! \brief NLMMT = 2 * LM * (LM + 2).
+ int nlmmt;
//! \brief Number of configurations
int configurations;
//! \brief QUESTION: definition?
@@ -54,6 +63,8 @@ public:
//! \brief QUESTION: definition?
dcomplex *fsas;
//! \brief QUESTION: definition?
+ dcomplex *vint;
+ //! \brief QUESTION: definition?
dcomplex **vints;
//! \brief QUESTION: definition?
double *sscs;
@@ -94,11 +105,13 @@ public:
* \param iog: `int *` Vector of spherical units ID numbers.
*/
C1(int ns, int l_max, int *nshl, int *iog);
+
/*! \brief C1 instance constructor copying all contents from a preexisting template
*
* \param rhs: `C1` preexisting template.
*/
C1(const C1& rhs);
+
//! \brief C1 instance destroyer.
~C1();
};
@@ -107,12 +120,14 @@ public:
*
*/
class C2 {
+protected:
//! \brief Number of spheres.
int nsph;
//! \brief Number of required orders.
int nhspo;
//! \brief QUESTION: what is nl?
int nl;
+
public:
//! \brief QUESTION: definition?
dcomplex *ris;
@@ -133,6 +148,7 @@ public:
* \param npntts: `int`
*/
C2(int ns, int nl, int npnt, int npntts);
+
/*! \brief C2 instance constructor copying its contents from preexisting instance.
*
* \param rhs: `C2` object to copy contents from
@@ -163,6 +179,7 @@ public:
/*! \brief C3 instance constructor.
*/
C3();
+
/*! \brief C3 instance constructor copying its contents from a preexisting object.
*/
C3(const C3& rhs);
@@ -237,19 +254,6 @@ protected:
//! \brief QUESTION: definition?
int lm;
- /*! \brief Allocate the necessary common vectors depending on configuration.
- *
- * The size of the vectors and matrices defined in various common
- * blocks, and particularly in C1, depends on many settings of the
- * problem configuration, such as the number of spheres, the number
- * of layers the spheres are made of, the field expansion order and
- * others. This function collects the calculations needed to infer
- * the necessary amount of memory for these configurable elements,
- * thus making the class constructor more compact and easier to handle.
- *
- * \param c4: `C4 *` Pointer to a C4 structure.
- */
- //void allocate_vectors(C4 *c4);
public:
//! \brief QUESTION: definition?
dcomplex *vh;
@@ -264,12 +268,8 @@ public:
//! \brief QUESTION: definition?
dcomplex **am0m;
//! \brief QUESTION: definition?
- dcomplex *vint;
- //! \brief QUESTION: definition?
dcomplex *vintm;
//! \brief QUESTION: definition?
- dcomplex **vints;
- //! \brief QUESTION: definition?
dcomplex *vintt;
//! \brief QUESTION: definition?
dcomplex **fsac;
@@ -305,6 +305,7 @@ public:
* \param c4: `C4 *` Pointer to a C4 structure.
*/
C1_AddOns(C4 *c4);
+
/*! \brief C1_AddOns instance constructor copying contents from a preexisting object.
*
* \param rhs: `C1_AddOns` preexisting object to copy from.
@@ -319,7 +320,7 @@ public:
*/
class C6 {
public:
- //! \brief QUESTION: definition?
+ //! \brief LMTPO = 2 * LM + 1.
int lmtpo;
//! \brief QUESTION: definition?
double *rac3j;
@@ -329,6 +330,7 @@ public:
* \param lmtpo: `int` QUESTION: definition?
*/
C6(int lmtpo);
+
/*! \brief C6 instance constructor copying contents from preexisting object.
*
* \param lmtpo: `int` QUESTION: definition?
@@ -352,6 +354,7 @@ protected:
int nlem;
//! \brief QUESTION: definition?
int nlemt;
+
public:
//! \brief QUESTION: definition?
dcomplex **gis;
@@ -368,6 +371,7 @@ public:
* \param nlemt: `int` QUESTION: definition?
*/
C9(int ndi, int nlem, int ndit, int nlemt);
+
/*! \brief C9 instance constructor copying contents from preexisting object.
*
* \param rhs: `C9` preexisting object to copy from
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index 591fc0bad9d6b428e579896b447764a556115624..2efd8b2df9e8b1fd1c14777433e5adf6e7627c1e 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/include/List.h b/src/include/List.h
index f97ebc6004c0867ce552480c69e16a11d671bc07..6d2d74279e2858fe865d9e328fa4c41bb35f7134 100644
--- a/src/include/List.h
+++ b/src/include/List.h
@@ -73,6 +73,7 @@ public:
current = old->p_prev;
delete old;
}
+ delete current;
}
/*! \brief Append an element at the end of the List.
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 2e57fa71c3a55c7cffbd68a68b7bec5ccaa15b8e..4621148a45de3fcd4b0107d449dbff454fd66ed3 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -379,7 +379,7 @@ void tqr(
/*! \brief Calculate the single-centered inversion of the M-matrix.
*
- * This function computes the single-centered inverrted M-matrix appearing in Eq. (5.28)
+ * This function computes the single-centered inverted M-matrix appearing in Eq. (5.28)
* of Borghese, Denti & Saija (2007).
*
* \param am: `complex double **`
diff --git a/src/include/logging.h b/src/include/logging.h
new file mode 100644
index 0000000000000000000000000000000000000000..fac4d71d780848d35dd0050cda1f48d368b679c2
--- /dev/null
+++ b/src/include/logging.h
@@ -0,0 +1,98 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
+/*! \file logging.h
+ *
+ * \brief Definition of the logging system.
+ */
+#ifndef INCLUDE_LOGGING_H_
+#define INCLUDE_LOGGING_H_
+
+//! \brief Debug level logging (maximum verbosity).
+#define LOG_DEBG 0
+//! \brief Standard information level logging (default).
+#define LOG_INFO 1
+//! \brief Warning level logging (almost quiet).
+#define LOG_WARN 2
+//! \brief Error level logging (silent, unless in pain).
+#define LOG_ERRO 3
+
+//! \brief Macro to stringize code lines
+#define TOSTRING(ARG) string(#ARG)
+
+/*! \brief Logger class.
+ *
+ * Loggers are objects used to track the execution of a code, reporting activities
+ * such as function calls and parameter settings. They can be used to inform the
+ * user about the execution of operations at runtime (e.g. by printing messages to
+ * the terminal), as well as to record the execution history in appropriate log
+ * files. The `Logger` class offers an implementation of logging system complying
+ * with the requirements of the NP_TMcode project.
+ *
+ * The Logger class is designed to work with open files. It is a user responsibility
+ * to check that the required log files are properly opened before use, and closed
+ * thereafter, if they are not the standard `stdout` and `stderr` streams.
+ */
+class Logger {
+ protected:
+ //! \brief Pointer to error stream.
+ FILE *err_output;
+ //! \brief Pointer to logging stream.
+ FILE *log_output;
+ //! \brief Last logged message.
+ std::string *last_message;
+ //! \brief Threshold of logging level.
+ int log_threshold;
+ //! \brief Number of identical message repetitions.
+ long repetitions;
+
+ public:
+ /*! \brief Logger instance constructor.
+ *
+ * \param threshold: `int` Threshold of the messages to be included in log. Can
+ * be `LOG_DEBG` (log everything), `LOG_INFO` (give detailed information),
+ * `LOG_WARN` (log odd looking effects), or `LOG_ERRO` (print error messages,
+ * `always active). The default behaviour is `LOG_WARN`.
+ * \param logging_output: `FILE *` Pointer to an open output file for common messages
+ * (optional, default is `stdout`).
+ * \param error_output: `FILE *` Pointer to an open output file for error messages
+ * (optional, default is `stderr`).
+ */
+ Logger(int threshold, FILE *logging_output=stdout, FILE *error_output=stderr);
+
+ /*! \brief Logger instance destroyer.
+ */
+ ~Logger();
+
+ /*! \brief Print a message to the error output.
+ *
+ * \param message: `string` The message to be printed.
+ */
+ void err(std::string message);
+
+ /*! \brief Print a summary of recurrent messages and clear the stack.
+ *
+ * \param level: `int` The priority level (default is `LOG_DEBG = 0`).
+ */
+ void flush(int level=LOG_DEBG);
+
+ /*! \brief Print a message, depending on its logging level.
+ *
+ * \param message: `string` The message to be printed.
+ * \param level: `int` The priority level (default is `LOG_INFO = 1`).
+ */
+ void log(std::string message, int level=LOG_INFO);
+
+ /*! \brief Push a recurrent message to the message stack.
+ *
+ * When a long stream of identical messages is expected, it may be more
+ * convenient to put them in a stack. The stack is the error output stream,
+ * so that no memory is consumed. After the call stack is over, or is the
+ * message changes, the stack is flushed, meaning that a summary message
+ * is written to the logging output and the stack counter is reset to 0.
+ *
+ * \param message: `string` The message to be stacked.
+ */
+ void push(std::string message);
+};
+
+#endif
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index da4b62360dd72b2b2cac1355b762f4ab70d20d57..b10496381e7c340b8ae818d664ba8bd7673b4914 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -25,25 +25,30 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
nsph = ns;
lm = l_max;
+ vec_rmi = new dcomplex[nsph * lm]();
+ vec_rei = new dcomplex[nsph * lm]();
rmi = new dcomplex*[lm];
rei = new dcomplex*[lm];
for (int ri = 0; ri < lm; ri++) {
- rmi[ri] = new dcomplex[nsph]();
- rei[ri] = new dcomplex[nsph]();
+ rmi[ri] = &(vec_rmi[nsph * ri]);
+ rei[ri] = &(vec_rei[nsph * ri]);
}
+ vec_w = new dcomplex[4 * nlmmt]();
w = new dcomplex*[nlmmt];
- for (int wi = 0; wi < nlmmt; wi++) w[wi] = new dcomplex[4]();
+ for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);
configurations = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (_iog[ci - 1] >= ci) configurations++;
}
+ vint = new dcomplex[16]();
+ vec_vints = new dcomplex[nsph * 16]();
vints = new dcomplex*[nsph];
rc = new double*[nsph];
nshl = new int[nsph]();
iog = new int[nsph]();
int conf_index = 0;
for (int vi = 0; vi < nsph; vi++) {
- vints[vi] = new dcomplex[16]();
+ vints[vi] = &(vec_vints[vi * 16]);
iog[vi] = _iog[vi];
if (iog[vi] >= vi + 1) {
nshl[conf_index] = _nshl[conf_index];
@@ -77,33 +82,39 @@ C1::C1(const C1& rhs) {
nsph = rhs.nsph;
lm = rhs.lm;
+ vec_rmi = new dcomplex[nsph * lm]();
+ vec_rei = new dcomplex[nsph * lm]();
rmi = new dcomplex*[lm];
rei = new dcomplex*[lm];
for (int ri = 0; ri < lm; ri++) {
- rmi[ri] = new dcomplex[nsph]();
- rei[ri] = new dcomplex[nsph]();
+ rmi[ri] = &(vec_rmi[nsph * ri]);
+ rei[ri] = &(vec_rei[nsph * ri]);
/*! Copy the contents from the template */
for (int rj=0; rj<nsph; rj++) {
rmi[ri][rj] = rhs.rmi[ri][rj];
rei[ri][rj] = rhs.rei[ri][rj];
}
}
+ vec_w = new dcomplex[4 * nlmmt]();
w = new dcomplex*[nlmmt];
for (int wi = 0; wi < nlmmt; wi++) {
- w[wi] = new dcomplex[4]();
+ w[wi] = &(vec_w[4 * wi]);
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];
nshl = new int[nsph]();
- for (int ni=0; ni<nsph; ni++) nshl[ni] = rhs.nshl[ni];
+ for (int ni = 0; ni < nsph; ni++) nshl[ni] = rhs.nshl[ni];
iog = new int[nsph]();
- for (int ni=0; ni<nsph; ni++) iog[ni] = rhs.iog[ni];
+ for (int ni = 0; ni < nsph; ni++) iog[ni] = rhs.iog[ni];
int conf_index = 0;
for (int vi = 0; vi < nsph; vi++) {
- vints[vi] = new dcomplex[16]();
- for (int vj=0; vj<16; vj++) vints[vi][vj] = rhs.vints[vi][vj];
+ vints[vi] = &(vec_vints[16 * vi]);
+ for (int vj = 0; vj < 16; vj++) vints[vi][vj] = rhs.vints[vi][vj];
}
for (int ri=0; ri<configurations; ri++) {
rc[ri] = new double[nshl[ri]]();
@@ -147,17 +158,19 @@ C1::C1(const C1& rhs) {
}
C1::~C1() {
+ delete[] vec_rmi;
+ delete[] vec_rei;
delete[] rmi;
delete[] rei;
- for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
+ delete[] vec_w;
+ delete[] w;
int conf_index = 0;
for (int ci = 0; ci < configurations; ci++) {
delete[] rc[ci];
}
delete[] rc;
- for (int vi = nsph - 1; vi > - 1; vi--) {
- delete[] vints[vi];
- }
+ delete[] vint;
+ delete[] vec_vints;
delete[] vints;
for (int si = nsph - 1; si > -1; si--) {
delete[] sas[si][1];
@@ -193,18 +206,20 @@ C1_AddOns::C1_AddOns(C4 *c4) {
lm = c4->lm;
vh = new dcomplex[(nsph * nsph - 1) * litpo]();
vj0 = new dcomplex[nsph * lmtpo]();
- vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
+ vj = new dcomplex[1]();
vyhj = new dcomplex[(nsph * nsph - 1) * litpos]();
vyj0 = new dcomplex[nsph * lmtpos]();
am0m = new dcomplex*[nlemt];
for (int ai = 0; ai < nlemt; ai++) {
am0m[ai] = new dcomplex[nlemt]();
}
- vint = new dcomplex[16](); // QUESTION: is dimension 16 generally fixed?
vintm = new dcomplex[16]();
vintt = new dcomplex[16]();
- vints = new dcomplex*[nsph];
- for (int vi = 0; vi < nsph; vi++) vints[vi] = new dcomplex[16]();
+ // This calculates the size of v3j0
+ const int nv3j = c4->nv3j;
+ v3j0 = new double[nv3j]();
+ ind3j = new int*[lmpo];
+ for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
fsac = new dcomplex*[2];
sac = new dcomplex*[2];
fsacm = new dcomplex*[2];
@@ -217,10 +232,6 @@ C1_AddOns::C1_AddOns(C4 *c4) {
ecscp = new dcomplex[2]();
scscpm = new dcomplex[2]();
ecscpm = new dcomplex[2]();
- v3j0 = new double[nv3j]();
- ind3j = new int*[lmpo];
- for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[lm]();
- //allocate_vectors(c4);
sscs = new double[nsph]();
ecscm = new double[2]();
scscm = new double[2]();
@@ -257,19 +268,12 @@ 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];
+ for (int hi = 0; hi < 16; hi++) {
vintm[hi] = rhs.vintm[hi];
vintt[hi] = rhs.vintt[hi];
}
- vints = new dcomplex*[nsph];
- for (int vi = 0; vi < nsph; vi++) {
- vints[vi] = new dcomplex[16]();
- for (int hj=0; hj<16; hj++) vints[vi][hj] = rhs.vints[vi][hj];
- }
fsac = new dcomplex*[2];
sac = new dcomplex*[2];
fsacm = new dcomplex*[2];
@@ -288,16 +292,12 @@ C1_AddOns::C1_AddOns(const C1_AddOns& rhs) {
scscpm = new dcomplex[2]();
ecscpm = new dcomplex[2]();
v3j0 = new double[nv3j]();
- for (int hi=0; hi<nv3j; hi++) v3j0[hi] = rhs.v3j0[hi];
+ for (int hi = 0; hi < nv3j; hi++) v3j0[hi] = rhs.v3j0[hi];
ind3j = new int*[lmpo];
for (int ii = 0; ii < lmpo; ii++) {
ind3j[ii] = new int[lm]();
for (int ij=0; ij<lm; ij++) ind3j[ii][ij] = rhs.ind3j[ii][ij];
}
- // vyhj = new dcomplex[int vyhjsize = (nsph * nsph - 1) * litpos]();
- // for (int hi=0; hi<vyhjsize; hi++) vyhj[hi] = rhs.vyhj[hi]
- // vyj0 = new dcomplex[int vyj0size=lmtpos * nsph]();
- // for (int hi=0; hi<vyj0size; hi++) vyj0[hi] = rhs.vyj0[hi];
sscs = new double[nsph]();
for (int hi=0; hi<nsph; hi++) sscs[hi] = rhs.sscs[hi];
ecscm = new double[2]();
@@ -330,11 +330,8 @@ C1_AddOns::~C1_AddOns() {
delete[] am0m[ai];
}
delete[] am0m;
- delete[] vint;
delete[] vintm;
delete[] vintt;
- for (int vi = nsph - 1; vi > -1; vi--) delete[] vints[vi];
- delete[] vints;
for (int fi = 1; fi > -1; fi--) {
delete[] fsac[fi];
delete[] sac[fi];
@@ -353,21 +350,6 @@ C1_AddOns::~C1_AddOns() {
delete[] ecsc;
}
-// void C1_AddOns::allocate_vectors(C4 *c4) {
-// // This calculates the size of v3j0
-// //int lm = (c4->li > c4->le) ? c4->li : c4->le;
-// const int nv3j = c4->nv3j;
-// v3j0 = new double[nv3j]();
-// ind3j = new int*[lmpo];
-// for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
-// // This calculates the size of vyhj
-// // int ivy = (nsph * nsph - 1) * c4->litpos;
-// // vyhj = new dcomplex[ivy]();
-// // This calculates the size of vyj0
-// // ivy = c4->lmtpos * c4->nsph;
-// // vyj0 = new dcomplex[ivy]();
-// }
-
C2::C2(int ns, int _nl, int npnt, int npntts) {
nsph = ns;
int max_n = (npnt > npntts) ? npnt : npntts;
@@ -481,7 +463,7 @@ C6::C6(int _lmtpo) {
C6::C6(const C6& rhs) {
lmtpo = rhs.lmtpo;
rac3j = new double[lmtpo]();
- for (int ri=0; ri<lmtpo; ri++) rac3j[ri] = rhs.rac3j[ri];
+ for (int ri = 0; ri < lmtpo; ri++) rac3j[ri] = rhs.rac3j[ri];
}
C6::~C6() {
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index fa91d648c506b59a474105e039ec95e7114dfbbd..80fb5ecd92f49faa4b17d6febd2d155d65db5ee9 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -228,6 +228,7 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {
sc_ph_start, sc_ph_step, sc_ph_end,
_jwtm
);
+ delete[] file_lines;
return conf;
}
@@ -322,6 +323,7 @@ ScattererConfiguration::~ScattererConfiguration() {
for (int j = 0; j < number_of_spheres; j++) {
delete[] dc0_matrix[i][j];
}
+ delete[] dc0_matrix[i];
}
delete[] dc0_matrix;
for (int i = 0; i < configurations; i++) {
@@ -400,7 +402,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
variable_name = "XIV";
if (instpc < 1) { // The variable vector is explicitly defined.
double xi;
- List<double> xi_vector;
+ List<double> *xi_vector = new List<double>(1);
str_target = file_lines[++last_read_line];
re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
regex_search(str_target, m, re);
@@ -408,7 +410,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
str_number = regex_replace(str_number, regex("D"), "e");
str_number = regex_replace(str_number, regex("d"), "e");
xi = stod(str_number);
- xi_vector.set(0, xi);
+ xi_vector->set(0, xi);
for (int jxi310 = 1; jxi310 < nxi; jxi310++) {
str_target = file_lines[++last_read_line];
regex_search(str_target, m, re);
@@ -416,9 +418,10 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
str_number = regex_replace(str_number, regex("D"), "e");
str_number = regex_replace(str_number, regex("d"), "e");
xi = stod(str_number);
- xi_vector.append(xi);
+ xi_vector->append(xi);
}
- variable_vector = xi_vector.to_array();
+ variable_vector = xi_vector->to_array();
+ delete xi_vector;
} else { // instpc >= 1: the variable vector is defined in steps
double xi, xi_step;
str_target = file_lines[++last_read_line];
@@ -871,60 +874,60 @@ void ScattererConfiguration::write_binary(string file_name, string mode) {
void ScattererConfiguration::write_hdf5(string file_name) {
int ies = (use_external_sphere)? 1 : 0;
- List<string> rec_name_list(1);
- List<string> rec_type_list(1);
- List<void *> rec_ptr_list(1);
+ List<string> *rec_name_list = new List<string>(1);
+ List<string> *rec_type_list = new List<string>(1);
+ List<void *> *rec_ptr_list = new List<void *>(1);
string str_type, str_name;
int configurations = 0;
for (int ci = 1; ci <= number_of_spheres; ci++) {
if(iog_vec[ci - 1] >= ci) configurations++;
}
- rec_name_list.set(0, "NSPH");
- rec_type_list.set(0, "INT32_(1)");
- rec_ptr_list.set(0, &number_of_spheres);
- rec_name_list.append("IES");
- rec_type_list.append("INT32_(1)");
- rec_ptr_list.append(&ies);
- rec_name_list.append("IOGVEC");
+ rec_name_list->set(0, "NSPH");
+ rec_type_list->set(0, "INT32_(1)");
+ rec_ptr_list->set(0, &number_of_spheres);
+ rec_name_list->append("IES");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&ies);
+ rec_name_list->append("IOGVEC");
str_type = "INT32_(" + to_string(number_of_spheres) + ")";
- rec_type_list.append(str_type);
- rec_ptr_list.append(iog_vec);
- rec_name_list.append("EXDC");
- rec_type_list.append("FLOAT64_(1)");
- rec_ptr_list.append(&exdc);
- rec_name_list.append("WP");
- rec_type_list.append("FLOAT64_(1)");
- rec_ptr_list.append(&wp);
- rec_name_list.append("XIP");
- rec_type_list.append("FLOAT64_(1)");
- rec_ptr_list.append(&xip);
- rec_name_list.append("IDFC");
- rec_type_list.append("INT32_(1)");
- rec_ptr_list.append(&idfc);
- rec_name_list.append("NXI");
- rec_type_list.append("INT32_(1)");
- rec_ptr_list.append(&number_of_scales);
- rec_name_list.append("XIVEC");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(iog_vec);
+ rec_name_list->append("EXDC");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&exdc);
+ rec_name_list->append("WP");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&wp);
+ rec_name_list->append("XIP");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&xip);
+ rec_name_list->append("IDFC");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&idfc);
+ rec_name_list->append("NXI");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&number_of_scales);
+ rec_name_list->append("XIVEC");
str_type = "FLOAT64_(" + to_string(number_of_scales) + ")";
- rec_type_list.append(str_type);
- rec_ptr_list.append(scale_vec);
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(scale_vec);
for (int i115 = 1; i115 <= number_of_spheres; i115++) {
if (iog_vec[i115 - 1] < i115) continue;
str_name = "NSHL_" + to_string(i115);
- rec_name_list.append(str_name);
- rec_type_list.append("INT32_(1)");
- rec_ptr_list.append(&(nshl_vec[i115 - 1])); // was not from IOG
+ rec_name_list->append(str_name);
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&(nshl_vec[i115 - 1])); // was not from IOG
str_name = "ROS_" + to_string(i115);
- rec_name_list.append(str_name);
- rec_type_list.append("FLOAT64_(1)");
- rec_ptr_list.append(&(radii_of_spheres[i115 - 1])); // was not from IOG
+ rec_name_list->append(str_name);
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&(radii_of_spheres[i115 - 1])); // was not from IOG
int nsh = nshl_vec[i115 - 1]; // was not from IOG
if (i115 == 1) nsh += ies;
str_name = "RCF_" + to_string(i115); // was not from IOG
str_type = "FLOAT64_(" + to_string(nsh) + ")";
- rec_name_list.append(str_name);
- rec_type_list.append(str_type);
- rec_ptr_list.append(&(rcf[i115 - 1][0])); // was not from IOG
+ rec_name_list->append(str_name);
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(&(rcf[i115 - 1][0])); // was not from IOG
}
int dim3 = (idfc == 0) ? number_of_scales : 1;
@@ -949,25 +952,29 @@ void ScattererConfiguration::write_hdf5(string file_name) {
}
}
str_type = "FLOAT64_(" + to_string(dc0m_size) + ")";
- rec_name_list.append("DC0M");
- rec_type_list.append(str_type);
- rec_ptr_list.append(dc0m);
+ rec_name_list->append("DC0M");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(dc0m);
- string *rec_names = rec_name_list.to_array();
- string *rec_types = rec_type_list.to_array();
- void **rec_pointers = rec_ptr_list.to_array();
- const int rec_num = rec_name_list.length();
- FileSchema schema(rec_num, rec_types, rec_names);
- HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC);
+ string *rec_names = rec_name_list->to_array();
+ string *rec_types = rec_type_list->to_array();
+ void **rec_pointers = rec_ptr_list->to_array();
+ const int rec_num = rec_name_list->length();
+ FileSchema *schema = new FileSchema(rec_num, rec_types, rec_names);
+ HDFFile *hdf_file = HDFFile::from_schema(*schema, file_name, H5F_ACC_TRUNC);
for (int ri = 0; ri < rec_num; ri++)
hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
hdf_file->close();
// Clean memory
+ delete rec_name_list;
+ delete rec_type_list;
+ delete rec_ptr_list;
delete[] dc0m;
delete[] rec_names;
delete[] rec_types;
delete[] rec_pointers;
+ delete schema;
delete hdf_file;
}
diff --git a/src/libnptm/Makefile b/src/libnptm/Makefile
index 43b886ffe2745399f4e1389471c4050f1d0c4c48..de9085d0539347faef38aeb9f14489e58d02583e 100644
--- a/src/libnptm/Makefile
+++ b/src/libnptm/Makefile
@@ -19,11 +19,11 @@ endif
include ../make.inc
-CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tfrfme.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o $(OBJDIR)/lapack_calls.o $(OBJDIR)/algebraic.o $(OBJDIR)/types.o
+CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tfrfme.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o $(OBJDIR)/lapack_calls.o $(OBJDIR)/algebraic.o $(OBJDIR)/types.o $(OBJDIR)/logging.o
-CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tfrfme.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o $(DYNOBJDIR)/lapack_calls.o $(DYNOBJDIR)/algebraic.o $(DYNOBJDIR)/types.o
+CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tfrfme.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o $(DYNOBJDIR)/lapack_calls.o $(DYNOBJDIR)/algebraic.o $(DYNOBJDIR)/types.o $(DYNOBJDIR)/logging.o
-CXX_NPTM_DEBUG=$(OBJDIR)/Commons.g* $(OBJDIR)/Configuration.g* $(OBJDIR)/file_io.g* $(OBJDIR)/Parsers.g* $(OBJDIR)/sph_subs.g* $(OBJDIR)/clu_subs.g* $(OBJDIR)/tfrfme.g* $(OBJDIR)/tra_subs.g* $(OBJDIR)/TransitionMatrix.g* $(OBJDIR)/lapack_calls.g* $(OBJDIR)/algebraic.g* $(OBJDIR)/types.g* $(DYNOBJDIR)/Commons.g* $(DYNOBJDIR)/Configuration.g* $(DYNOBJDIR)/file_io.g* $(DYNOBJDIR)/Parsers.g* $(DYNOBJDIR)/sph_subs.g* $(DYNOBJDIR)/clu_subs.g* $(DYNOBJDIR)/tfrfme.g* $(DYNOBJDIR)/tra_subs.g* $(DYNOBJDIR)/TransitionMatrix.g* $(DYNOBJDIR)/lapack_calls.g* $(DYNOBJDIR)/algebraic.g* $(DYNOBJDIR)/types.g*
+CXX_NPTM_DEBUG=$(OBJDIR)/Commons.g* $(OBJDIR)/Configuration.g* $(OBJDIR)/file_io.g* $(OBJDIR)/Parsers.g* $(OBJDIR)/sph_subs.g* $(OBJDIR)/clu_subs.g* $(OBJDIR)/tfrfme.g* $(OBJDIR)/tra_subs.g* $(OBJDIR)/TransitionMatrix.g* $(OBJDIR)/lapack_calls.g* $(OBJDIR)/algebraic.g* $(OBJDIR)/types.g* $(OBJDIR)/logging.g* $(DYNOBJDIR)/Commons.g* $(DYNOBJDIR)/Configuration.g* $(DYNOBJDIR)/file_io.g* $(DYNOBJDIR)/Parsers.g* $(DYNOBJDIR)/sph_subs.g* $(DYNOBJDIR)/clu_subs.g* $(DYNOBJDIR)/tfrfme.g* $(DYNOBJDIR)/tra_subs.g* $(DYNOBJDIR)/TransitionMatrix.g* $(DYNOBJDIR)/lapack_calls.g* $(DYNOBJDIR)/algebraic.g* $(DYNOBJDIR)/types.g* $(DYNOBJDIR)/logging.g*
all: $(BUILDDIR_NPTM)/libnptm.a $(BUILDDIR_NPTM)/libnptm.so
diff --git a/src/libnptm/Parsers.cpp b/src/libnptm/Parsers.cpp
index 22ff56db05a692980df8f4140f97a1ba6c227c12..6c8326f7c4274258deb7f42bb16d47b050239bc9 100644
--- a/src/libnptm/Parsers.cpp
+++ b/src/libnptm/Parsers.cpp
@@ -21,19 +21,20 @@
std::string *load_file(std::string file_name, int *count = 0) {
std::fstream input_file(file_name.c_str(), std::ios::in);
- List<std::string> file_lines = List<std::string>();
+ List<std::string> *file_lines = new List<std::string>();
std::string line;
if (input_file.is_open()) {
getline(input_file, line);
- file_lines.set(0, line);
+ file_lines->set(0, line);
while (getline(input_file, line)) {
- file_lines.append(line);
+ file_lines->append(line);
}
input_file.close();
} else {
throw FILE_NOT_FOUND_ERROR;
}
- std::string *array_lines = file_lines.to_array();
- if (count != 0) *count = file_lines.length();
+ std::string *array_lines = file_lines->to_array();
+ if (count != 0) *count = file_lines->length();
+ delete file_lines;
return array_lines;
}
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index 744eae869b77731a46f8635b6deeddf3cbc318c9..2473d63602e9d589c412f8ffba01425ff0e6c6cb 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -1100,7 +1100,7 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
cc = dconjg(c1ao->sac[jpo1 - 1][ipo1 - 1]);
for (int ipo2 = 1; ipo2 <= 2; ipo2 ++) {
for (int jpo2 = 1; jpo2 <= 2; jpo2++) {
- c1ao->vint[i++] = c1ao->sac[jpo2 - 1][ipo2 - 1] * cc * cfsq;
+ c1->vint[i++] = c1ao->sac[jpo2 - 1][ipo2 - 1] * cc * cfsq;
} // jpo2 loop
} // ipo2 loop
} // jpo1 loop
@@ -1793,7 +1793,7 @@ void scr2(
for (int ipo2 = 1; ipo2 <= 2; ipo2++) {
for (int jpo2 = 1; jpo2 <= 2; jpo2++) {
j++;
- c1ao->vints[i24 - 1][j - 1] = c1->sas[i24 - 1][jpo2 - 1][ipo2 - 1] * cc * cfsq;
+ c1->vints[i24 - 1][j - 1] = c1->sas[i24 - 1][jpo2 - 1][ipo2 - 1] * cc * cfsq;
} // jpo2 loop
} // ipo2 loop
} // jpo1 loop
diff --git a/src/libnptm/logging.cpp b/src/libnptm/logging.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..67a586e467d4c0a5aa8f19a0197e4209dbe1355d
--- /dev/null
+++ b/src/libnptm/logging.cpp
@@ -0,0 +1,69 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
+/*! \file logging.cpp
+ *
+ * \brief Implementation of the logging system.
+ */
+#include <cstdio>
+#include <string>
+
+#ifndef INCLUDE_LOGGING_H_
+#include "../include/logging.h"
+#endif
+
+using namespace std;
+
+Logger::Logger(int threshold, FILE *logging_output, FILE *error_output) {
+ last_message = new string("");
+ log_threshold = threshold;
+ log_output = logging_output;
+ err_output = error_output;
+ repetitions = 0;
+}
+
+Logger::~Logger() {
+ delete last_message;
+}
+
+void Logger::err(std::string message) {
+ fprintf(err_output, "%s", message.c_str());
+ fflush(err_output);
+}
+
+void Logger::flush(int level) {
+ string summary = "\"" + *last_message + "\" issued " + to_string(repetitions);
+ if (repetitions == 1) summary += " time.\n";
+ else summary += " times.\n";
+ if (level == LOG_ERRO) err(summary);
+ else {
+ if (level >= log_threshold) {
+ fprintf(log_output, "%s", summary.c_str());
+ fflush(log_output);
+ }
+ }
+ delete last_message;
+ last_message = new string("");
+ repetitions = 0;
+}
+
+void Logger::log(std::string message, int level) {
+ if (level == LOG_ERRO) err(message);
+ else {
+ if (level >= log_threshold) {
+ fprintf(log_output, "%s", message.c_str());
+ fflush(log_output);
+ }
+ }
+}
+
+void Logger::push(std::string message) {
+ if (repetitions > 0) {
+ if (message.compare(*last_message) != 0) {
+ flush(LOG_DEBG);
+ }
+ }
+ log(message, LOG_DEBG);
+ delete last_message;
+ last_message = new string(message);
+ repetitions++;
+}
diff --git a/src/make.inc b/src/make.inc
index 962471f05863f79b4e84ace15e994e6e29fbc8a4..e854d59e55315a2a57c46c8e3c4e58d5d8678e54 100644
--- a/src/make.inc
+++ b/src/make.inc
@@ -82,6 +82,9 @@ override CXXFLAGS+= -DUSE_LAPACK -DLAPACK_ILP64
ifdef USE_MKL
override CXXFLAGS+= -DMKL_ILP64 -DUSE_MKL -I$(MKLROOT)/include
endif
+ifdef USE_OPENMP
+override CXXFLAGS+= -fopenmp
+endif
endif
endif
@@ -91,9 +94,12 @@ ifndef CXXLDFLAGS
ifndef HDF5_LIB
override HDF5_LIB=/usr/lib/x86_64-linux-gnu/hdf5/serial
endif
-override CXXLDFLAGS=-L/usr/lib64 -L$(HDF5_LIB) -lhdf5
+override CXXLDFLAGS=-L/usr/lib64 -L$(HDF5_LIB) -lhdf5 $(STATICFLAG)
ifdef USE_LAPACK
override CXXLDFLAGS+= $(LAPACK_LDFLAGS)
+ifdef USE_OPENMP
+override CXXLDFLAGS+= -lopenblas64
+endif
endif
override CXXLDFLAGS+= $(LDFLAGS)
endif
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index b4bf9652d3e0a9401e076f750b92615073750b97..94860567344395db2182f27b1019874f1460dc23 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -102,7 +102,6 @@ void sphere(string config_file, string data_file, string output_path) {
}
double frx = 0.0, fry = 0.0, frz = 0.0;
double cfmp, cfsp, sfmp, sfsp;
- dcomplex *vint = new dcomplex[16];
int jw;
int nsph = gconf->number_of_spheres;
C1 *c1 = new C1(nsph, gconf->l_max, sconf->nshl_vec, sconf->iog_vec);
@@ -509,8 +508,8 @@ void sphere(string config_file, string data_file, string output_path) {
output, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
frx, fry, frz
);
- for (int i225 = 0; i225 < 16; i225++) vint[i225] = c1->vints[ns226][i225];
- mmulc(vint, cmullr, cmul);
+ for (int i225 = 0; i225 < 16; i225++) c1->vint[i225] = c1->vints[ns226][i225];
+ mmulc(c1->vint, cmullr, cmul);
fprintf(output, " MULS\n ");
for (int imul = 0; imul < 4; imul++) {
for (int jmul = 0; jmul < 4; jmul++) {
@@ -528,9 +527,9 @@ void sphere(string config_file, string data_file, string output_path) {
else fprintf(output, "\n");
}
for (int vi = 0; vi < 16; vi++) {
- double value = real(vint[vi]);
+ double value = real(c1->vint[vi]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(vint[vi]);
+ value = imag(c1->vint[vi]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int imul = 0; imul < 4; imul++) {
@@ -577,7 +576,6 @@ void sphere(string config_file, string data_file, string output_path) {
delete[] upsmp;
delete[] unmp;
delete[] unsmp;
- delete[] vint;
delete[] argi;
delete[] args;
delete[] gaps;