diff --git a/build/README.md b/build/README.md
index 26ca5880008cb1d4adcfddfb3a5e6c144f5c4a55..554c0d68035e9e089330df47a8493d58edb53a2f 100644
--- a/build/README.md
+++ b/build/README.md
@@ -111,4 +111,3 @@ where the arguments must be valid paths to binary transition matrix files saved
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-
\ No newline at end of file
diff --git a/containers/docker/README.md b/containers/docker/README.md
index 46e5d752af612bf4850f0e5d65cadafb01c27140..d7f8e8f9b6f0681ab3513db167b6b26f136333c1 100644
--- a/containers/docker/README.md
+++ b/containers/docker/README.md
@@ -26,4 +26,3 @@ where `<image name>` is either `np-tmcode` or `np-tmcode-run`. One may also add
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-
\ No newline at end of file
diff --git a/containers/singularity/README.md b/containers/singularity/README.md
index 0f400330c55328f9621059db4cda984e1ba59cae..fad127007ba4d2b5d6dbfcb2f093ecab6f8b176c 100644
--- a/containers/singularity/README.md
+++ b/containers/singularity/README.md
@@ -33,4 +33,3 @@ where `<full path to image name>` is the name of the sif image, including full o
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-
\ No newline at end of file
diff --git a/doc/src/README.md b/doc/src/README.md
index 55f445eb0cbcd88e5d5ace491628e0394c034be5..28f258ffdaae2cafb173ebc798fe974d1618c1bc 100644
--- a/doc/src/README.md
+++ b/doc/src/README.md
@@ -29,4 +29,3 @@ Alternatively, you can use `make` from the project `src` folder (named `np_tmcod
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-
\ No newline at end of file
diff --git a/src/README.md b/src/README.md
index a665a1d843acfe397daa639d56892aa34f173646..fe3cbeef3f96417a1a708c6006a7923c7b2eb58e 100644
--- a/src/README.md
+++ b/src/README.md
@@ -51,4 +51,3 @@ In all cases, build commands executed through `make` will output the object file
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-
\ No newline at end of file
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index faabf2163db9b6902b868877742d5115144ec6d3..b35fe1d9b996816382e413089321f30548e6880c 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -30,6 +30,12 @@
#include <mpi.h>
#endif
#endif
+#ifdef USE_NVTX
+#include <nvtx3/nvToolsExt.h>
+#endif
+#ifdef USE_MAGMA
+#include <cuda_runtime.h>
+#endif
#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
@@ -86,9 +92,27 @@ void cluster(const string& config_file, const string& data_file, const string& o
string timing_name = output_path + "/c_timing_mpi"+ to_string(mpidata->rank) +".log";
FILE *timing_file = fopen(timing_name.c_str(), "w");
Logger *time_logger = new Logger(LOG_DEBG, timing_file);
- Logger *logger = new Logger(LOG_INFO);
+ Logger *logger = new Logger(LOG_DEBG);
+#ifdef USE_MAGMA
+ int device_count;
+ cudaGetDeviceCount(&device_count);
+ logger->log("DEBUG: Proc-" + to_string(mpidata->rank) + " found " + to_string(device_count) + " CUDA devices.\n", LOG_DEBG);
+ logger->log("INFO: Process " + to_string(mpidata->rank) + " initializes MAGMA.\n");
+ magma_int_t magma_result = magma_init();
+ if (magma_result != MAGMA_SUCCESS) {
+ logger->err("ERROR: Process " + to_string(mpidata->rank) + " failed to initilize MAGMA.\n");
+ logger->err("PROC-" + to_string(mpidata->rank) + ": MAGMA error code " + to_string(magma_result) + "\n");
+ fclose(timing_file);
+ delete time_logger;
+ delete logger;
+ return;
+ }
+#endif
// the following only happens on MPI process 0
if (mpidata->rank == 0) {
+#ifdef USE_NVTX
+ nvtxRangePush("Set up");
+#endif
logger->log("INFO: making legacy configuration...", LOG_INFO);
ScattererConfiguration *sconf = NULL;
try {
@@ -97,7 +121,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
logger->err("\nERROR: failed to open scatterer configuration file.\n");
string message = "FILE: " + string(ex.what()) + "\n";
logger->err(message);
- exit(1);
+ fclose(timing_file);
+ delete time_logger;
+ delete logger;
+ return;
}
sconf->write_formatted(output_path + "/c_OEDFB");
sconf->write_binary(output_path + "/c_TEDF");
@@ -110,9 +137,15 @@ void cluster(const string& config_file, const string& data_file, const string& o
string message = "FILE: " + string(ex.what()) + "\n";
logger->err(message);
if (sconf) delete sconf;
- exit(1);
+ fclose(timing_file);
+ delete time_logger;
+ delete logger;
+ return;
}
logger->log(" done.\n", LOG_INFO);
+#ifdef USE_NVTX
+ nvtxRangePop();
+#endif
int s_nsph = sconf->number_of_spheres;
int nsph = gconf->number_of_spheres;
if (s_nsph == nsph) {
@@ -167,7 +200,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
string tppoan_name = output_path + "/c_TPPOAN";
tppoan.open(tppoan_name.c_str(), ios::out | ios::binary);
if (tppoan.is_open()) {
-#ifdef USE_LAPACK
+#ifdef USE_MAGMA
+ logger->log("INFO: using MAGMA calls.\n", LOG_INFO);
+#elif defined USE_LAPACK
logger->log("INFO: using LAPACK calls.\n", LOG_INFO);
#else
logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);
@@ -196,7 +231,13 @@ void cluster(const string& config_file, const string& data_file, const string& o
// do the first iteration on jxi488 separately, since it seems to be different from the others
int jxi488 = 1;
chrono::time_point<chrono::high_resolution_clock> start_iter_1 = chrono::high_resolution_clock::now();
+#ifdef USE_NVTX
+ nvtxRangePush("First iteration");
+#endif
int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, output, output_path, tppoan);
+#ifdef USE_NVTX
+ nvtxRangePop();
+#endif
chrono::time_point<chrono::high_resolution_clock> end_iter_1 = chrono::high_resolution_clock::now();
elapsed = start_iter_1 - t_start;
string message = "INFO: Calculation setup took " + to_string(elapsed.count()) + "s.\n";
@@ -206,6 +247,19 @@ void cluster(const string& config_file, const string& data_file, const string& o
message = "INFO: First iteration took " + to_string(elapsed.count()) + "s.\n";
logger->log(message);
time_logger->log(message);
+ if (jer != 0) {
+ // First loop failed. Halt the calculation.
+ tppoan.close();
+ fclose(timing_file);
+ fclose(output);
+ delete p_scattering_angles;
+ delete cid;
+ delete logger;
+ delete time_logger;
+ delete sconf;
+ delete gconf;
+ return;
+ }
// here go the calls that send data to be duplicated on other MPI processes from process 0 to others, using MPI broadcasts, but only if MPI is actually used
#ifdef MPI_VERSION
@@ -219,6 +273,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
// 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;
+#ifdef USE_NVTX
+ nvtxRangePush("Parallel loop");
+#endif
#pragma omp parallel
{
// Create and initialise this variable here, so that if OpenMP is enabled it is local to the thread, and if OpenMP is not enabled it has a well-defiled value anyway
@@ -247,7 +304,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
fstream &tppoan_2 = *tppoanp_2;
// make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads
#pragma omp barrier
- if (myompthread==0) logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
+ if (myompthread==0) {
+ logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
+ }
// ok, now I can actually start the parallel calculations
#pragma omp for
for (jxi488 = cid_2->firstxi; jxi488 <= cid_2->lastxi; jxi488++) {
@@ -268,6 +327,11 @@ void cluster(const string& config_file, const string& data_file, const string& o
logger->log(message);
}
} // closes pragma omp parallel
+#ifdef USE_NVTX
+ nvtxRangePop();
+
+ nvtxRangePush("Output concatenation");
+#endif
#ifdef _OPENMP
#pragma omp barrier
{
@@ -340,6 +404,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
}
}
}
+#endif
+#ifdef USE_NVTX
+ nvtxRangePop();
#endif
tppoanp->close();
delete tppoanp;
@@ -510,6 +577,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
delete gconf;
}
+#endif
+#ifdef USE_MAGMA
+ logger->log("INFO: Process " + to_string(mpidata->rank) + " finalizes MAGMA.\n");
+ magma_finalize();
#endif
fclose(timing_file);
delete time_logger;
@@ -520,7 +591,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
{
int nxi = sconf->number_of_scales;
string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";
- Logger *logger = new Logger(LOG_INFO);
+ Logger *logger = new Logger(LOG_DEBG);
logger->log(message);
chrono::duration<double> elapsed;
chrono::time_point<chrono::high_resolution_clock> interval_start, interval_end;
@@ -542,7 +613,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
int jwtm = gconf->jwtm;
np_int ndit = 2 * nsph * cid->c4->nlim;
int isq, ibf;
-
+
+#ifdef USE_NVTX
+ nvtxRangePush("Prepare matrix calculation");
+#endif
fprintf(output, "========== JXI =%3d ====================\n", jxi488);
double xi = sconf->get_scale(jxi488 - 1);
double exdc = sconf->exdc;
@@ -599,23 +673,43 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
//break;
}
} // i132 loop
+#ifdef USE_NVTX
+ nvtxRangePop();
+#endif
interval_start = chrono::high_resolution_clock::now();
+#ifdef USE_NVTX
+ nvtxRangePush("Calculate inverted matrix");
+#endif
cms(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6);
+#ifdef USE_NVTX
+ nvtxRangePop();
+#endif
interval_end = chrono::high_resolution_clock::now();
elapsed = interval_end - interval_start;
message = "INFO: matrix calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
logger->log(message);
interval_start = chrono::high_resolution_clock::now();
+#ifdef USE_NVTX
+ nvtxRangePush("Invert the matrix");
+#endif
invert_matrix(cid->am, ndit, jer, mxndm);
+#ifdef USE_NVTX
+ nvtxRangePop();
+#endif
interval_end = chrono::high_resolution_clock::now();
elapsed = interval_end - interval_start;
message = "INFO: matrix inversion for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
logger->log(message);
if (jer != 0) {
+ message = "ERROR: matrix inversion ended with error code " + to_string(jer) + ".\n";
+ logger->err(message);
return jer;
// break; // jxi488 loop: goes to memory clean
}
interval_start = chrono::high_resolution_clock::now();
+#ifdef USE_NVTX
+ nvtxRangePush("Average calculation");
+#endif
ztm(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6, cid->c9);
if (idfc >= 0) {
if (jxi488 == jwtm) {
@@ -683,11 +777,17 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
tppoan.write(reinterpret_cast<char *>(&(cid->vk)), sizeof(double));
pcrsm0(cid->vk, exri, inpol, cid->c1, cid->c1ao, cid->c4);
apcra(cid->zpv, cid->c4->le, cid->c1ao->am0m, inpol, sqk, cid->gapm, cid->gappm);
+#ifdef USE_NVTX
+ nvtxRangePop();
+#endif
interval_end = chrono::high_resolution_clock::now();
elapsed = interval_end - interval_start;
message = "INFO: average calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
logger->log(message);
interval_start = chrono::high_resolution_clock::now();
+#ifdef USE_NVTX
+ nvtxRangePush("Angle loop");
+#endif
double th = sa->th;
for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable?
double ph = sa->ph;
@@ -1172,6 +1272,9 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
} // jph484 loop
th += sa->thstp;
} // jth486 loop
+#ifdef USE_NVTX
+ nvtxRangePop();
+#endif
interval_end = chrono::high_resolution_clock::now();
elapsed = interval_end - interval_start;
message = "INFO: angle loop for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
diff --git a/src/include/magma_calls.h b/src/include/magma_calls.h
new file mode 100644
index 0000000000000000000000000000000000000000..4ee244c52cb5e1da253e83e1f684d59dd10db2cc
--- /dev/null
+++ b/src/include/magma_calls.h
@@ -0,0 +1,36 @@
+/* Copyright 2004 INAF - Osservatorio Astronomico di Cagliari
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+ */
+
+/*! \file magma_calls.h
+ *
+ * \brief C++ interface to MAGMA calls.
+ *
+ */
+
+#ifndef INCLUDE_MAGMA_CALLS_H_
+#define INCLUDE_MAGMA_CALLS_H_
+
+/*! \brief Invert a complex matrix with double precision elements.
+ *
+ * Use LAPACKE64 to perform an in-place matrix inversion for a complex
+ * matrix with double precision elements.
+ *
+ * \param mat: Matrix of complex. The matrix to be inverted.
+ * \param n: `np_int` The number of rows and columns of the [n x n] matrix.
+ * \param jer: `int &` Reference to an integer return flag.
+ */
+void magma_zinvert(dcomplex **mat, np_int n, int &jer);
+
+#endif
diff --git a/src/include/types.h b/src/include/types.h
index e258c75e8bfa6a448de8f5a35dddfbbef5f2e170..0ac02c192f1787e065af58b4553150c5583e3ecb 100644
--- a/src/include/types.h
+++ b/src/include/types.h
@@ -25,24 +25,41 @@
typedef __complex__ double dcomplex;
-#ifdef USE_LAPACK
#ifdef USE_MKL
+#ifdef USE_ILP64
#ifndef MKL_INT
#define MKL_INT int64_t
-#endif
+#endif // MKL_INT
+#else
+#ifndef MKL_INT
+#define MKL_INT int32_t
+#endif // MKL_INT
+#endif // USE_ILP64
+#endif // USE_MKL
+
+#ifdef USE_LAPACK
+#ifdef USE_MKL
#include <mkl_lapacke.h>
#else
#include <lapacke.h>
-#endif
+#endif // USE_MKL
+#endif // USE_LAPACK
+
+#ifdef USE_MAGMA
+#include <magma_v2.h>
#endif
#ifndef np_int
#ifdef lapack_int
#define np_int lapack_int
#else
+#ifdef USE_ILP64
#define np_int int64_t
-#endif
-#endif
+#else
+#define np_int int32_t
+#endif // USE_ILP64
+#endif // lapack_int
+#endif // np_int
/*! \brief Get the real part of a complex number.
*
diff --git a/src/libnptm/Makefile b/src/libnptm/Makefile
index de9085d0539347faef38aeb9f14489e58d02583e..b02e5f2e9f27e8a77df5db27edd5c55713984c30 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 $(OBJDIR)/logging.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)/magma_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 $(DYNOBJDIR)/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)/magma_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* $(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*
+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)/magma_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)/magma_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/algebraic.cpp b/src/libnptm/algebraic.cpp
index 637091765d2fe7d5c1de7d6d332de69ce208a1db..53b7dec9d6904e01de8c600c725e0ea3c90475ba 100644
--- a/src/libnptm/algebraic.cpp
+++ b/src/libnptm/algebraic.cpp
@@ -27,6 +27,12 @@
#endif
#endif
+#ifdef USE_MAGMA
+#ifndef INCLUDE_MAGMA_CALLS_H_
+#include "../include/magma_calls.h"
+#endif
+#endif
+
#ifndef INCLUDE_ALGEBRAIC_H_
#include "../include/algebraic.h"
#endif
@@ -39,7 +45,9 @@ using namespace std;
void invert_matrix(dcomplex **mat, np_int size, int &ier, np_int max_size) {
ier = 0;
-#ifdef USE_LAPACK
+#ifdef USE_MAGMA
+ magma_zinvert(mat, size, ier);
+#elif defined USE_LAPACK
zinvert(mat, size, ier);
#else
lucin(mat, max_size, size, ier);
diff --git a/src/libnptm/lapack_calls.cpp b/src/libnptm/lapack_calls.cpp
index 718b1ecbdeb14b79849ba9d6cb20d201e0be8178..a9e39cd549eabe9cb8c18ab85a89743b758bc867 100644
--- a/src/libnptm/lapack_calls.cpp
+++ b/src/libnptm/lapack_calls.cpp
@@ -21,19 +21,20 @@
#include "../include/types.h"
#endif
+/*
#ifdef USE_LAPACK
#ifdef USE_MKL
#include <mkl_lapacke.h>
#else
#include <lapacke.h>
#endif
+*/
+#ifdef USE_LAPACK
#ifndef INCLUDE_LAPACK_CALLS_H_
#include "../include/lapack_calls.h"
#endif
-#endif
-#ifdef USE_LAPACK
void zinvert(dcomplex **mat, np_int n, int &jer) {
jer = 0;
dcomplex *arr = &(mat[0][0]);
diff --git a/src/libnptm/logging.cpp b/src/libnptm/logging.cpp
index 7565e91305f9f8e00d75714a5b8fc6279a9b50e1..2f7526e58f6da76787b3bf4a0d4789597f428d21 100644
--- a/src/libnptm/logging.cpp
+++ b/src/libnptm/logging.cpp
@@ -39,56 +39,44 @@ Logger::~Logger() {
}
void Logger::err(const std::string& message) {
-//#pragma omp critical
- {
- fprintf(err_output, "%s", message.c_str());
- fflush(err_output);
- }
+ fprintf(err_output, "%s", message.c_str());
+ fflush(err_output);
}
void Logger::flush(int level) {
-//#pragma omp critical
- {
- 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);
- }
+ 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;
}
+ delete last_message;
+ last_message = new string("");
+ repetitions = 0;
}
void Logger::log(const std::string& message, int level) {
-//#pragma omp critical
- {
- if (level == LOG_ERRO) err(message);
- else {
- if (level >= log_threshold) {
- fprintf(log_output, "%s", message.c_str());
- fflush(log_output);
- }
+ if (level == LOG_ERRO) err(message);
+ else {
+ if (level >= log_threshold) {
+ fprintf(log_output, "%s", message.c_str());
+ fflush(log_output);
}
}
}
void Logger::push(const std::string& message) {
-//#pragma omp critical
- {
- if (repetitions > 0) {
- if (message.compare(*last_message) != 0) {
- flush(LOG_DEBG);
- }
+ 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++;
}
+ log(message, LOG_DEBG);
+ delete last_message;
+ last_message = new string(message);
+ repetitions++;
}
diff --git a/src/libnptm/magma_calls.cpp b/src/libnptm/magma_calls.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ab83103a4c0c5ddc469245291f6e6e6974e27b54
--- /dev/null
+++ b/src/libnptm/magma_calls.cpp
@@ -0,0 +1,59 @@
+/* Copyright 2004 INAF - Osservatorio Astronomico di Cagliari
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+ */
+
+/*! \file magma_calls.cpp
+ *
+ * \brief Implementation of the interface with MAGMA libraries.
+ */
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
+#ifdef USE_MAGMA
+#ifndef INCLUDE_MAGMA_CALLS_H_
+#include "../include/magma_calls.h"
+#endif
+
+void magma_zinvert(dcomplex **mat, np_int n, int &jer) {
+ // magma_int_t result = magma_init();
+ magma_int_t err = MAGMA_SUCCESS;
+ magma_queue_t queue = NULL;
+ magma_device_t dev = 0;
+ magma_queue_create(dev, &queue);
+ magmaDoubleComplex *dwork; // dwork - workspace
+ magma_int_t ldwork; // size of dwork
+ magma_int_t *piv , info; // piv - array of indices of inter -
+ magma_int_t m = (magma_int_t)n; // changed rows; a - mxm matrix
+ magma_int_t mm = m * m; // size of a, r, c
+ magmaDoubleComplex *a = (magmaDoubleComplex *)&(mat[0][0]); // a - mxm matrix on the host
+ magmaDoubleComplex *d_a; // d_a - mxm matrix a on the device
+ ldwork = m * magma_get_zgetri_nb(m); // optimal block size
+ // allocate matrices
+ err = magma_zmalloc(&d_a, mm); // device memory for a
+ err = magma_zmalloc(&dwork, ldwork); // dev. mem. for ldwork
+ piv = new magma_int_t[m]; // host mem.
+ magma_zsetmatrix(m, m, a, m, d_a , m, queue); // copy a -> d_a
+
+ magma_zgetrf_gpu(m, m, d_a, m, piv, &info);
+ magma_zgetri_gpu(m, d_a, m, piv, dwork, ldwork, &info);
+
+ magma_zgetmatrix( m, m, d_a , m, a, m, queue); // copy d_a -> a
+ delete[] piv; // free host memory
+ magma_free(d_a); // free device memory
+ magma_queue_destroy(queue); // destroy queue
+ // result = magma_finalize();
+ jer = (int)err;
+}
+#endif
diff --git a/src/make.inc b/src/make.inc
index 28cee0656e6208962ef313de67cad48cb6ff9580..6f65edda86bf6399edff50f2172a70007a2d7e5d 100644
--- a/src/make.inc
+++ b/src/make.inc
@@ -1,117 +1,190 @@
# FC defines the fortran compiler to use. If undefined, GNU Make tries to use `f77`
ifndef FC
-override FC=gfortran
+ override FC=gfortran
endif
# FCFLAGS defines the compilation options for the fortran compiler
ifndef FCFLAGS
-override FCFLAGS=-std=legacy -O3
+ override FCFLAGS=-std=legacy -O3
endif
# AR defines the command to create static library files
ifndef AR
-override AR=ar
+ override AR=ar
endif
# ARFLAGS defines the flags for AR to create static library files
ifndef ARFLAGS
-override ARFLAGS=-rs
+ override ARFLAGS=-rs
endif
# PICFLAGS defines the additional flags for the c++ compiler to create objects suitable for shared library creation
ifndef PICFLAGS
-override PICFLAGS=-fPIC
+ override PICFLAGS=-fPIC
endif
# LDFLAGS defines the default linker flags
ifndef LDFLAGS
-override LDFLAGS=
+ override LDFLAGS=
endif
# CXX defines the default C++ compiler to use. If undefined, GNU Make tries to use g++
ifndef CXX
-ifdef USE_MPI
-override CXX=mpicxx
-else
-override CXX=g++
-endif
+ ifdef USE_MPI
+ override CXX=mpicxx
+ else
+ override CXX=g++
+ endif
endif
ifdef USE_MPI
-override MPI_CXXFLAGS+=-DUSE_MPI
+ override MPI_CXXFLAGS+=-DUSE_MPI
endif
# HDF5_INCLUDE defines the default path to the HDF5 include files to use
ifndef HDF5_INCLUDE
-override HDF5_INCLUDE=/usr/include/hdf5/serial
+ override HDF5_INCLUDE=/usr/include/hdf5/serial
endif
# define (outside) USE_LAPACK for lapacke support, LAPACK_ILP64 for ilp64 interface, MKL_ILP64 the same if using MKL implementation
ifdef USE_LAPACK
-ifndef LAPACK_ILP64
-override LAPACK_ILP64=1
-endif
+# define (outside) USE_ILP64 for long long int support in lapack/mkl/magma interfaces
+ ifdef USE_ILP64
+ ifndef LAPACK_ILP64
+ override LAPACK_ILP64=1
+ endif #LAPACK_ILP64
+ endif
# define (outside) USE_MKL to use the MKL implementation of lapacke
-ifdef USE_MKL
-ifndef MKL_ILP64
-override MKL_ILP64=1
-endif
-ifndef LAPACK_INCLUDE
+ ifdef USE_MKL
+# define (outside) USE_ILP64 for long long int support in lapack/mkl/magma interfaces
+ ifdef USE_ILP64
+ ifndef MKL_ILP64
+ override MKL_ILP64=1
+ endif #MKL_ILP64
+ endif #USE_ILP64
+ ifndef LAPACK_INCLUDE
# this is for the MKL implementation
-override LAPACK_INCLUDE=$(MKLROOT)/include
-endif
-ifndef LAPACK_LDFLAGS
+ override LAPACK_INCLUDE=$(MKLROOT)/include
+ endif #LAPACK_INCLUDE
+ ifndef LAPACK_LDFLAGS
# this is for the MKL implementation
-override LAPACK_LDFLAGS=-L$(MKLROOT)/lib -Wl,--no-as-needed -lmkl_intel_ilp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl
-endif
-# the next else refers to USE_MKL
+# define (outside) USE_ILP64 for long long int support in lapack/mkl/magma interfaces
+ ifdef USE_ILP64
+ override LAPACK_LDFLAGS=-L$(MKLROOT)/lib -Wl,--no-as-needed -lmkl_intel_ilp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl
else
-ifndef LAPACK_INCLUDE
+ override LAPACK_LDFLAGS=-L$(MKLROOT)/lib -Wl,--no-as-needed -lmkl_intel_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl
+ endif #USE_ILP64
+ endif #LAPACK_LDFLAGS
+# the next else refers to USE_MKL
+ else #this is for when USE_MKL is _not_ defined
+ ifndef LAPACK_INCLUDE
# this is for standard "vanilla" lapacke64
-override LAPACK_INCLUDE=/usr/include
-endif
-ifndef LAPACK_LDFLAGS
+ override LAPACK_INCLUDE=/usr/include
+ endif # LAPACK_INCLUDE
+ ifndef LAPACK_LDFLAGS
+ ifdef USE_ILP64
# this is for standard "vanilla" lapacke64
-override LAPACK_LDFLAGS=-llapacke64
-endif
+ override LAPACK_LDFLAGS=-llapacke64
+ else
+ override LAPACK_LDFLAGS=-llapacke
+ endif #USE_ILP64
+ endif #LAPACK_LDFLAGS
# the next endif is for USE_MKL
-endif
+ endif
#the next endif is for USE_LAPACK
endif
+# define (outside) USE_MAGMA for magma support
+ifdef USE_MAGMA
+ ifndef MAGMA_LDFLAGS
+ ifdef MAGMA_LIB
+ override MAGMA_LDFLAGS= -L$(MAGMA_LIB)
+ endif
+ ifdef CUDA_HOME
+ override MAGMA_LDFLAGS+= -L$(CUDA_HOME)/lib64
+ endif
+ override MAGMA_LDFLAGS+= -lmagma -lcudart
+#the next endif is for MAGMA_LDFLAGS
+ endif
+#the next endif is for USE_MAGMA
+endif
+
+# define (outside) USE_NVTX for NVIDIA profiling
+ifdef USE_NVTX
+ ifndef NVTX_CXXFLAGS
+ override NVTX_CXXFLAGS= -DUSE_NVTX
+ ifdef CUDA_HOME
+ override NVTX_CXXFLAGS+= -I$(CUDA_HOME)/include
+# closes CUDA_HOME
+ endif
+# closes NVTX_CXXFLAGS
+ endif
+# closes USE_NVTX
+endif
+
# CXXFLAGS defines the default compilation options for the C++ compiler
ifndef CXXFLAGS
-override CXXFLAGS=-O3 -ggdb -pg -coverage -I$(HDF5_INCLUDE) $(MPI_CXXFLAGS)
-ifdef USE_OPENMP
-override CXXFLAGS+= -fopenmp
-endif
-ifdef USE_LAPACK
-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
+ override CXXFLAGS=-O3 -ggdb -pg -coverage -I$(HDF5_INCLUDE) $(MPI_CXXFLAGS) $(NVTX_CXXFLAGS)
+ ifdef USE_OPENMP
+ override CXXFLAGS+= -fopenmp
+# closes USE_OPENMP
+ endif
+ ifdef USE_ILP64
+ override CXXFLAGS+= -DUSE_ILP64
+ endif
+ ifdef USE_LAPACK
+ override CXXFLAGS+= -DUSE_LAPACK
+ ifdef USE_ILP64
+ override CXXFLAGS+= -DLAPACK_ILP64
+ endif
+# closes USE_LAPACK
+ endif
+ ifdef USE_MKL
+ override CXXFLAGS+= -DUSE_MKL -I$(MKLROOT)/include
+ ifdef USE_ILP64
+ override CXXFLAGS+= -DMKL_ILP64
+ endif
+# closes USE_MKL
+ endif
+ ifdef USE_OPENMP
+ override CXXFLAGS+= -fopenmp
+# closes USE_OPENMP
+ endif
+ ifdef USE_MAGMA
+ override CXXFLAGS+= -DUSE_MAGMA
+ ifdef CUDA_HOME
+ override CXXFLAGS+= -I$(CUDA_HOME)/include
+ endif
+ ifdef MAGMA_INCLUDE
+ override CXXFLAGS+= -I$(MAGMA_INCLUDE)
+ endif
+ ifdef USE_ILP64
+ override CXXFLAGS+= -DMAGMA_ILP64
+ endif
+# closes USE_MAGMA
+ endif
+#closes CXXFLAGS
endif
# HDF5_LIB defines the default path to the HDF5 libraries to use
# CXXLDFLAGS defines the default linker flags to use for C++ codes
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 $(STATICFLAG)
-ifdef USE_LAPACK
-override CXXLDFLAGS+= $(LAPACK_LDFLAGS)
-endif
-override CXXLDFLAGS+= $(LDFLAGS)
+ ifndef HDF5_LIB
+ override HDF5_LIB=/usr/lib/x86_64-linux-gnu/hdf5/serial
+ endif
+ override CXXLDFLAGS=-L/usr/lib64 -L$(HDF5_LIB) -lhdf5 $(STATICFLAG)
+ ifdef USE_LAPACK
+ override CXXLDFLAGS+= $(LAPACK_LDFLAGS)
+ endif
+ ifdef USE_MAGMA
+ override CXXLDFLAGS+= $(MAGMA_LDFLAGS)
+ endif
+ override CXXLDFLAGS+= $(LDFLAGS)
endif
#SOFLAGS defines the additional flags for the c++ compiler to create a shared object file
ifndef SOFLAGS
-override SOFLAGS=-shared
+ override SOFLAGS=-shared
endif
%.o : %.f
diff --git a/src/scripts/README.md b/src/scripts/README.md
index 2bda93b8fd8a238cfc1760f977b099db5207c6f4..de3fead8ccbe0ec983758d66a2742f14cc3625a3 100644
--- a/src/scripts/README.md
+++ b/src/scripts/README.md
@@ -27,6 +27,14 @@ or just:
will print a help screen, giving a brief explanation of all the possible options.
+## Estimating the execution time from a terminal log
+
+Performance estimates can be obtained from the code logging system, assuming the uses chose to save the terminal output to some log file. To obtain the time spent by the code in performing a specific operation, the syntax is:
+
+ > $PATH_TO_SCRIPT/pytiming.py --logname=LOG_FILE [--filter=FILTER --threads=NUM_THREADS]
+
+where `LOG_FILE` must be the name of a file containing the output that would normally go to terminal, `FILTER` must be the starting character sequence of the log line identifying the operation that should be taken into account, and `NUM_THREADS` is the number of processes that were used to perform the whole calculation loop. In case no filter is given, the script provides an estimate of the total amount of time spent in doing the calculation. This estimate, however, is known not to be reliable, because it ignores thread concurrency effects. A more accurate estimate of the total time spent in executing the code is always saved in a file named `c_timing.log`.
+
# License
Copyright 2004 INAF - Osservatorio Astronomico di Cagliari
@@ -42,4 +50,3 @@ will print a help screen, giving a brief explanation of all the possible options
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-
\ No newline at end of file
diff --git a/src/scripts/pycompare.py b/src/scripts/pycompare.py
index 8e2fd8a93ac904a95dace2df1a3d4800db8892fc..c1e6f639e881c5816ce2a7192b33e75a4fcf025e 100755
--- a/src/scripts/pycompare.py
+++ b/src/scripts/pycompare.py
@@ -19,6 +19,7 @@
# The script execution requires python3.
import re
+import os
from math import log10
from sys import argv
@@ -99,51 +100,98 @@ def compare_files(config):
fortran_file = open(config['fortran_file_name'], 'r')
c_file = open(config['c_file_name'], 'r')
l_file = None
- f_lines = fortran_file.readlines()
- c_lines = c_file.readlines()
- fortran_file.close()
- c_file.close()
- if (len(f_lines) == len(c_lines)):
+ f_lines = []
+ c_lines = []
+ line_count = 0
+ if (not config['linewise']):
+ f_lines = fortran_file.readlines()
+ c_lines = c_file.readlines()
line_count = len(f_lines)
- num_len = 1
- if (line_count > 0):
- num_len = max(4, int(log10(line_count)) + 1)
- if (config['log_html']):
- l_file = open(config['html_output'], 'w')
- l_file.write("<!DOCTYPE html>\n")
- l_file.write("<html xmnls=\"http://www.w3.org/1999/xhtml\">\n")
- l_file.write(" <header>\n")
- l_file.write(
- " <h1>Comparison between {0:s} and {1:s}</h1>\n".format(
- config['fortran_file_name'], config['c_file_name']
- )
+ fortran_file.close()
+ c_file.close()
+ else: # line-wise comparison mode
+ f_lines = [fortran_file.readline()]
+ c_lines = [c_file.readline()]
+ print("INFO: using line-wise mode")
+ print("INFO: counting result lines...")
+ while (f_lines[0] != ''):
+ if (c_lines[0] != ''):
+ line_count += 1
+ else:
+ print("ERROR: C++ file is shorter than FORTRAN file.")
+ fortran_file.close()
+ c_file.close()
+ mismatch_count['errors'] = line_count
+ return mismatch_count
+ f_lines[0] = fortran_file.readline()
+ c_lines[0] = c_file.readline()
+ if (c_lines[0] != ''):
+ print("ERROR: C++ file is longer than FORTRAN file.")
+ fortran_file.close()
+ c_file.close()
+ mismatch_count['errors'] = line_count
+ return mismatch_count
+ fortran_file.close()
+ c_file.close()
+ print("INFO: the output files have %d lines"%line_count)
+ fortran_file = open(config['fortran_file_name'], 'r')
+ c_file = open(config['c_file_name'], 'r')
+ num_read_lines = 0
+ # LOG FILE INITIALIZATION #
+ if (config['log_html']):
+ l_file = open(config['html_output'], 'w')
+ l_file.write("<!DOCTYPE html>\n")
+ l_file.write("<html xmnls=\"http://www.w3.org/1999/xhtml\">\n")
+ l_file.write(" <header>\n")
+ l_file.write(
+ " <h1>Comparison between {0:s} and {1:s}</h1>\n".format(
+ config['fortran_file_name'], config['c_file_name']
)
- l_file.write(" </header>\n")
- l_file.write(" <body>\n")
- l_file.write(" <div>Numeric noise is marked <span style=\"font-weight: bold; color: rgb(0,185,0)\">"
- + "GREEN</span>, warnings are marked <span style=\"font-weight: bold; color: rgb(0,0,255)\">"
- + "BLUE</span> and errors are marked <span style=\"font-weight: bold; color: rgb(255,0,0)\">"
- + "RED</span>.</div>\n")
- for li in range(line_count):
- line_result = compare_lines(f_lines[li], c_lines[li], config, li + 1, num_len, l_file)
- mismatch_count['errors'] += line_result[0]
- mismatch_count['warnings'] += line_result[1]
- mismatch_count['noisy'] += line_result[2]
- if (mismatch_count['errors'] > 0 and not config['check_all']):
- print("INFO: mismatch found at line %d"%(li + 1))
- break
- if l_file is not None:
- l_file.write(" </body>\n")
- l_file.write("</html>\n")
- l_file.close()
- else:
- mismatch_count['errors'] = len(c_lines)
- print("ERROR: {0:s} and {1:s} have different numbers of lines!".format(
- config['fortran_file_name'], config['c_file_name']
- ))
- if (config['log_html']):
- print("Different file sizes. No log produced.")
- config['log_html'] = False
+ )
+ l_file.write(" </header>\n")
+ l_file.write(" <body>\n")
+ l_file.write(" <div>Numeric noise is marked <span style=\"font-weight: bold; color: rgb(0,185,0)\">"
+ + "GREEN</span>, warnings are marked <span style=\"font-weight: bold; color: rgb(0,0,255)\">"
+ + "BLUE</span> and errors are marked <span style=\"font-weight: bold; color: rgb(255,0,0)\">"
+ + "RED</span>.</div>\n")
+ # END LOG FILE INITIALIZATION #
+ line_loop = True
+ num_len = 1
+ if (line_count > 0):
+ num_len = max(4, int(log10(line_count)) + 1)
+ print("INFO: checking file contents...")
+ while (line_loop):
+ if (not config['linewise']):
+ line_loop = False
+ else:
+ f_lines = [fortran_file.readline()]
+ c_lines = [c_file.readline()]
+ num_read_lines += 1
+ # Start here the comparison loop
+ if (len(f_lines) == len(c_lines)):
+ for li in range(len(f_lines)):
+ line_result = compare_lines(f_lines[li], c_lines[li], config, li + 1, num_len, l_file)
+ mismatch_count['errors'] += line_result[0]
+ mismatch_count['warnings'] += line_result[1]
+ mismatch_count['noisy'] += line_result[2]
+ if (mismatch_count['errors'] > 0 and not config['check_all']):
+ print("INFO: mismatch found at line %d"%(li + 1))
+ break
+ else:
+ mismatch_count['errors'] = len(c_lines)
+ print("ERROR: {0:s} and {1:s} have different numbers of lines!".format(
+ config['fortran_file_name'], config['c_file_name']
+ ))
+ if (config['log_html']):
+ print("Different file sizes. No log produced.")
+ config['log_html'] = False
+ if (num_read_lines >= line_count):
+ line_loop = False
+ #End line loop
+ if l_file is not None:
+ l_file.write(" </body>\n")
+ l_file.write("</html>\n")
+ l_file.close()
return mismatch_count
## \brief Perform the comparison of two file lines.
@@ -381,6 +429,7 @@ def parse_arguments():
'fortran_file_name': '',
'c_file_name': '',
'full_log': False,
+ 'linewise': False,
'log_html': False,
'html_output': 'pycompare.html',
'warning_threshold': 0.005,
@@ -403,6 +452,8 @@ def parse_arguments():
config['warning_threshold'] = float(split_arg[1])
elif (arg.startswith("--help")):
config['help_mode'] = True
+ elif (arg.startswith("--linewise")):
+ config['linewise'] = True
elif (arg.startswith("--quick")):
config['check_all'] = False
else:
@@ -424,8 +475,9 @@ def print_help():
print("--full Print all lines to log file (default prints only mismatches).")
print("--help Print this help and exit.")
print("--html[=OPT_OUTPUT_NAME] Enable logging to HTML file (default logs to \"pycompare.html\").")
+ print("--linewise Load only one line at a time. Useful to compare big files (false by default).")
print("--quick Stop on first mismatch (default is to perform a full check).")
- print("--warn Set a fractional threshold for numeric warning (default=0.005).")
+ print("--warn Set a fractional threshold for numeric warning (default = 0.005).")
print(" ")
## \brief Add summary information to the HTML log file
@@ -440,20 +492,25 @@ def print_help():
# \param noisy: `int` The number of noisy values detected by the comparison.
def reformat_log(config, errors, warnings, noisy):
log_file = open(config['html_output'], 'r')
- log_lines = log_file.readlines()
- log_file.close()
- log_file = open(config['html_output'], 'w')
- for i in range(7): log_file.write(log_lines[i] + "\n")
+ new_file = open("PYCOMPARE_TEMPORARY_LOG.html", 'w')
+ for hi in range(7):
+ log_line = log_file.readline()
+ new_file.write(log_line)
str_errors = "error" if errors == 1 else "errors"
str_warnings = "warning" if warnings == 1 else "warnings"
str_noisy = "noisy value" if noisy == 1 else "noisy values"
summary = " <div>Comparison yielded %d %s"%(errors, str_errors)
summary = summary + ", %d %s"%(warnings, str_warnings)
summary = summary + " and %d %s.</div>\n"%(noisy, str_noisy)
- log_file.write(summary)
- for i in range(7, len(log_lines)): log_file.write(log_lines[i] + "\n")
+ new_file.write(summary)
+ log_line = log_file.readline()
+ while (log_line != ''):
+ new_file.write(log_line)
+ log_line = log_file.readline()
log_file.close()
-
+ new_file.close()
+ os.remove(config['html_output'])
+ os.rename("PYCOMPARE_TEMPORARY_LOG.html", config['html_output'])
# ### PROGRAM EXECUTION ###
## \cond
diff --git a/src/scripts/pytiming.py b/src/scripts/pytiming.py
new file mode 100755
index 0000000000000000000000000000000000000000..3887702e088249715c428cad49ec2cc24a97f1cb
--- /dev/null
+++ b/src/scripts/pytiming.py
@@ -0,0 +1,123 @@
+#!/bin/python3
+
+## @package pycompare
+# \brief Script to calculate the execution time of logged operations
+#
+# Comparing the numeric output can be rendered hard by the amount of information
+# contained in a typical output file and the necessity to determine whether a
+# difference is actually significant or just caused by numeric noise hitting
+# negligible values. The task of `pycompare.py` is to compare two output files, in
+# the assumption that they were written by the FORTRAN and the C++ versions of
+# the code and to flag all the possible inconsistencies according to various
+# severity levels (namely: NOISE, WARNING, and ERROR).
+#
+# After execution, the script returns an exit code, which is set to 0, if no
+# error-level inconsistencies were found, or 1 otherwise. This can be used by
+# subsequent system calls to set up a testing suite checking whether the code
+# is able to reproduce legacy results.
+#
+# The script execution requires python3.
+
+import re
+
+from sys import argv
+
+## \cond
+time_reg = re.compile(r'[0-9]+\.[0-9]+s')
+## \endcond
+
+## \brief Main execution code
+#
+# `main()` is the function that handles the creation of the script configuration
+# and the execution of the calculation. It returns 0 on successful completion.
+#
+# \returns exit_code: `int` 0 on successful completion.
+def main():
+ config = parse_arguments()
+ exit_code = 0
+ if config['help_mode'] or len(argv) == 1:
+ config['help_mode'] = True
+ print_help()
+ else:
+ if config['log_name'] is None:
+ exit_code = 1
+ else:
+ operation_time = get_time_from_log(config)
+ print("Calculation took %fs."%operation_time)
+ return exit_code
+
+## \brief Parse a log file and extract the time.
+#
+# \param config: `dict` A dictionary containing the script configuration.
+#
+# \returns operation_time: `float` The time of the requested operation in seconds.
+def get_time_from_log(config):
+ op_time = 0.0
+ log_file = open(config['log_name'], 'r')
+ file_lines = log_file.readlines()
+ log_file.close()
+ for li in range(len(file_lines)):
+ str_line = file_lines[li]
+ if (config['filter'] == "" or str_line.startswith(config['filter'])):
+ time_iters = time_reg.finditer(str_line)
+ time_groups = []
+ for ti in time_iters:
+ time_groups.append(ti.group())
+ if len(time_groups) == 1:
+ op_time += float(time_groups[0][:-1])
+ if config['threads'] > 1:
+ op_time /= config['threads']
+ return op_time
+
+## \brief Parse the command line arguments.
+#
+# The script behaviour can be modified through a set of mandatory and optional
+# arguments. The only mandatory argument is the name of the log file to be
+# parsed. Additional optional arguments are an operation filter, which should
+# be the starting sequence of the log strings to pe included in the timing
+# calculation and the number of threads used during code execution.
+#
+# \returns config: `dict` A dictionary containing the script configuration.
+def parse_arguments():
+ config = {
+ 'log_name': None,
+ 'help_mode': False,
+ 'filter': "",
+ 'threads': 1,
+ }
+ for arg in argv[1:]:
+ split_arg = arg.split("=")
+ if (arg.startswith("--logname")):
+ config['log_name'] = split_arg[1]
+ elif (arg.startswith("--filter")):
+ config['filter'] = split_arg[1]
+ elif (arg.startswith("--threads")):
+ config['threads'] = int(split_arg[1])
+ elif (arg.startswith("--help")):
+ config['help_mode'] = True
+ else:
+ raise Exception("Unrecognized argument \'{0:s}\'".format(arg))
+ return config
+
+## \brief Print a command-line help summary.
+def print_help():
+ print(" ")
+ print("*** PYTIMING ***")
+ print(" ")
+ print("Get the amount of time spent in calculation.")
+ print(" ")
+ print("Usage: \"./pytiming.py OPTIONS\" ")
+ print(" ")
+ print("Valid options are: ")
+ print("--logname=TIMING_LOG File containing log of timing (mandatory).")
+ print("--filter=FILTER Start of the log lines to be accounted for (optional).")
+ print("--help Print this help and exit.")
+ print("--threads=NUM_THREADS Number of threads or processes used in calculation (optional).")
+ print(" ")
+
+
+# ### PROGRAM EXECUTION ###
+## \cond
+res = main()
+## \endcond
+exit(res)
diff --git a/test_data/README.md b/test_data/README.md
index 7c60da1433a04d4f4174e1a7f64431f42a23acdd..5926500edd6b7edbbbc4c666a4ad6abb769cddb8 100644
--- a/test_data/README.md
+++ b/test_data/README.md
@@ -75,4 +75,3 @@ were the different lines have the following roles:
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-
\ No newline at end of file
diff --git a/test_data/cluster/case_48_long/DCLU b/test_data/cluster/case_48_long/DCLU
new file mode 100644
index 0000000000000000000000000000000000000000..7b63fd42a24fa5a468fdee8d44f7bd9914d2efbf
--- /dev/null
+++ b/test_data/cluster/case_48_long/DCLU
@@ -0,0 +1,53 @@
+ 48 10 10 5760 1 149 300 0 0
+ 0.0000000e+00 1.0000000e-08 0.0000000e+00
+ 0.0000000e+00 -1.0000000e-08 0.0000000e+00
+ 8.5355339e-09 5.0000000e-09 0.0000000e+00
+ 8.5355339e-09 -5.0000000e-09 0.0000000e+00
+ -8.5355339e-09 5.0000000e-09 0.0000000e+00
+ -8.5355339e-09 -5.0000000e-09 0.0000000e+00
+ 0.0000000e+00 1.0000000e-08 1.0000000e-08
+ 0.0000000e+00 -1.0000000e-08 1.0000000e-08
+ 8.5355339e-09 5.0000000e-09 1.0000000e-08
+ 8.5355339e-09 -5.0000000e-09 1.0000000e-08
+ -8.5355339e-09 5.0000000e-09 1.0000000e-08
+ -8.5355339e-09 -5.0000000e-09 1.0000000e-08
+ 0.0000000e+00 1.0000000e-08 2.0000000e-08
+ 0.0000000e+00 -1.0000000e-08 2.0000000e-08
+ 8.5355339e-09 5.0000000e-09 2.0000000e-08
+ 8.5355339e-09 -5.0000000e-09 2.0000000e-08
+ -8.5355339e-09 5.0000000e-09 2.0000000e-08
+ -8.5355339e-09 -5.0000000e-09 2.0000000e-08
+ 0.0000000e+00 1.0000000e-08 3.0000000e-08
+ 0.0000000e+00 -1.0000000e-08 3.0000000e-08
+ 8.5355339e-09 5.0000000e-09 3.0000000e-08
+ 8.5355339e-09 -5.0000000e-09 3.0000000e-08
+ -8.5355339e-09 5.0000000e-09 3.0000000e-08
+ -8.5355339e-09 -5.0000000e-09 3.0000000e-08
+ -6.4278761e-09 7.6604444e-09 4.0000000e-08
+ 6.4278761e-09 -7.6604444e-09 4.0000000e-08
+ 3.3246603e-09 9.3167577e-09 4.0000000e-08
+ 9.7525364e-09 1.6563132e-09 4.0000000e-08
+ -9.7525364e-09 -1.6563132e-09 4.0000000e-08
+ -3.3246603e-09 -9.3167577e-09 4.0000000e-08
+ -8.6602540e-09 5.0000000e-09 5.0000000e-08
+ 8.6602540e-09 -5.0000000e-09 5.0000000e-08
+ -6.2360066e-11 9.8919892e-09 5.0000000e-08
+ 8.5978940e-09 4.8919892e-09 5.0000000e-08
+ -8.5978940e-09 -4.8919892e-09 5.0000000e-08
+ 6.2360066e-11 -9.8919892e-09 5.0000000e-08
+ -5.0000000e-09 8.6602540e-09 6.0000000e-08
+ 5.0000000e-09 -8.6602540e-09 6.0000000e-08
+ 4.8919892e-09 8.5978940e-09 6.0000000e-08
+ 9.8919892e-09 -6.2360066e-11 6.0000000e-08
+ -9.8919892e-09 6.2360066e-11 6.0000000e-08
+ -4.8919892e-09 -8.5978940e-09 6.0000000e-08
+ -2.5881905e-09 9.6592583e-09 7.0000000e-08
+ 2.5881905e-09 -9.6592583e-09 7.0000000e-08
+ 6.9505974e-09 7.0387879e-09 7.0000000e-08
+ 9.5387879e-09 -2.6204704e-09 7.0000000e-08
+ -9.5387879e-09 2.6204704e-09 7.0000000e-08
+ -6.9505974e-09 -7.0387879e-09 7.0000000e-08
+ 0.00D01 0.00D00 0.00D01 180.00D00 0.00D01 0.00D01
+ 0.00D01 0.00D00 0.00D01 0.00D00 0.00D01 0.00D01
+ 1
+0
diff --git a/test_data/cluster/case_48_long/DEDFB b/test_data/cluster/case_48_long/DEDFB
new file mode 100644
index 0000000000000000000000000000000000000000..ac4071a84011eead571b865749f22b82d0393435
--- /dev/null
+++ b/test_data/cluster/case_48_long/DEDFB
@@ -0,0 +1,370 @@
+ 48 0
+ 1.000D00 3.0000000D08 1.000000D00 0 181 0 3
+4.2000000e-07
+4.2100000e-07
+4.2200000e-07
+4.2300000e-07
+4.2400000e-07
+4.2500000e-07
+4.2600000e-07
+4.2700000e-07
+4.2800000e-07
+4.2900000e-07
+4.3000000e-07
+4.3100000e-07
+4.3200000e-07
+4.3300000e-07
+4.3400000e-07
+4.3500000e-07
+4.3600000e-07
+4.3700000e-07
+4.3800000e-07
+4.3900000e-07
+4.4000000e-07
+4.4100000e-07
+4.4200000e-07
+4.4300000e-07
+4.4400000e-07
+4.4500000e-07
+4.4600000e-07
+4.4700000e-07
+4.4800000e-07
+4.4900000e-07
+4.5000000e-07
+4.5100000e-07
+4.5200000e-07
+4.5300000e-07
+4.5400000e-07
+4.5500000e-07
+4.5600000e-07
+4.5700000e-07
+4.5800000e-07
+4.5900000e-07
+4.6000000e-07
+4.6100000e-07
+4.6200000e-07
+4.6300000e-07
+4.6400000e-07
+4.6500000e-07
+4.6600000e-07
+4.6700000e-07
+4.6800000e-07
+4.6900000e-07
+4.7000000e-07
+4.7100000e-07
+4.7200000e-07
+4.7300000e-07
+4.7400000e-07
+4.7500000e-07
+4.7600000e-07
+4.7700000e-07
+4.7800000e-07
+4.7900000e-07
+4.8000000e-07
+4.8100000e-07
+4.8200000e-07
+4.8300000e-07
+4.8400000e-07
+4.8500000e-07
+4.8600000e-07
+4.8700000e-07
+4.8800000e-07
+4.8900000e-07
+4.9000000e-07
+4.9100000e-07
+4.9200000e-07
+4.9300000e-07
+4.9400000e-07
+4.9500000e-07
+4.9600000e-07
+4.9700000e-07
+4.9800000e-07
+4.9900000e-07
+5.0000000e-07
+5.0100000e-07
+5.0200000e-07
+5.0300000e-07
+5.0400000e-07
+5.0500000e-07
+5.0600000e-07
+5.0700000e-07
+5.0800000e-07
+5.0900000e-07
+5.1000000e-07
+5.1100000e-07
+5.1200000e-07
+5.1300000e-07
+5.1400000e-07
+5.1500000e-07
+5.1600000e-07
+5.1700000e-07
+5.1800000e-07
+5.1900000e-07
+5.2000000e-07
+5.2100000e-07
+5.2200000e-07
+5.2300000e-07
+5.2400000e-07
+5.2500000e-07
+5.2600000e-07
+5.2700000e-07
+5.2800000e-07
+5.2900000e-07
+5.3000000e-07
+5.3100000e-07
+5.3200000e-07
+5.3300000e-07
+5.3400000e-07
+5.3500000e-07
+5.3600000e-07
+5.3700000e-07
+5.3800000e-07
+5.3900000e-07
+5.4000000e-07
+5.4100000e-07
+5.4200000e-07
+5.4300000e-07
+5.4400000e-07
+5.4500000e-07
+5.4600000e-07
+5.4700000e-07
+5.4800000e-07
+5.4900000e-07
+5.5000000e-07
+5.5100000e-07
+5.5200000e-07
+5.5300000e-07
+5.5400000e-07
+5.5500000e-07
+5.5600000e-07
+5.5700000e-07
+5.5800000e-07
+5.5900000e-07
+5.6000000e-07
+5.6100000e-07
+5.6200000e-07
+5.6300000e-07
+5.6400000e-07
+5.6500000e-07
+5.6600000e-07
+5.6700000e-07
+5.6800000e-07
+5.6900000e-07
+5.7000000e-07
+5.7100000e-07
+5.7200000e-07
+5.7300000e-07
+5.7400000e-07
+5.7500000e-07
+5.7600000e-07
+5.7700000e-07
+5.7800000e-07
+5.7900000e-07
+5.8000000e-07
+5.8100000e-07
+5.8200000e-07
+5.8300000e-07
+5.8400000e-07
+5.8500000e-07
+5.8600000e-07
+5.8700000e-07
+5.8800000e-07
+5.8900000e-07
+5.9000000e-07
+5.9100000e-07
+5.9200000e-07
+5.9300000e-07
+5.9400000e-07
+5.9500000e-07
+5.9600000e-07
+5.9700000e-07
+5.9800000e-07
+5.9900000e-07
+6.0000000e-07
+ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 5.00d-9
+ 1.000000D00
+ ( 0.6965478043675424 , 0.012177869466280656 )
+ ( 0.6910987175052832 , 0.012554064506388863 )
+ ( 0.6855360485660978 , 0.012942754740189038 )
+ ( 0.6798066752951744 , 0.013350374401007377 )
+ ( 0.6738920152444744 , 0.013779653436464834 )
+ ( 0.6678508686371436 , 0.01422354911340056 )
+ ( 0.6616505736901351 , 0.014686154187061628 )
+ ( 0.6552091880683985 , 0.015178506717677436 )
+ ( 0.648626795673544 , 0.01568797851280193 )
+ ( 0.6418846479452085 , 0.0162169277351358 )
+ ( 0.6348445917724804 , 0.016784713760307202 )
+ ( 0.6276474807258959 , 0.0173725368046002 )
+ ( 0.6202769786636171 , 0.01798241937365526 )
+ ( 0.6125531613000286 , 0.01864103719769382 )
+ ( 0.6046548955664187 , 0.01932301612058496 )
+ ( 0.5965489649813588 , 0.02003289743754264 )
+ ( 0.5880407157515819 , 0.02080162282741932 )
+ ( 0.5793396913640244 , 0.02159734076893952 )
+ ( 0.5703638195559585 , 0.0224326248886312 )
+ ( 0.5609512658546029 , 0.0233357270125777 )
+ ( 0.5513282106958756 , 0.024269375252460982 )
+ ( 0.5413110764950205 , 0.025264669694510276 )
+ ( 0.5308516974311571 , 0.026332713845906165 )
+ ( 0.520167523655191 , 0.0274338031711828 )
+ ( 0.5088858212296016 , 0.0286376825595599 )
+ ( 0.4972110480736359 , 0.02990890485111186 )
+ ( 0.48526191306235184 , 0.03122138600075616 )
+ ( 0.47246314200315426 , 0.03269528657253276 )
+ ( 0.45937630121857237 , 0.0342162843063422 )
+ ( 0.4457057614275002 , 0.03584434633527552 )
+ ( 0.4312635907829137 , 0.03762906608903808 )
+ ( 0.41654518905862686 , 0.039451835640326824 )
+ ( 0.40068237774761795 , 0.0415236061440668 )
+ ( 0.38431490386773487 , 0.0436942853071663 )
+ ( 0.3672967795844664 , 0.045987918531406206 )
+ ( 0.3490036098881765 , 0.0485936998775316 )
+ ( 0.3304185923407915 , 0.05121856408811032 )
+ ( 0.3101531206636717 , 0.054283277567830475 )
+ ( 0.2891722918455745 , 0.0575074338345765 )
+ ( 0.2672177912261406 , 0.06092302679777248 )
+ ( 0.24336337771698055 , 0.06494756656796982 )
+ ( 0.2193964481088083 , 0.06876714714786752 )
+ ( 0.1922129574569899 , 0.07382533351429672 )
+ ( 0.1646196822471307 , 0.07883893467287696 )
+ ( 0.1347899789733728 , 0.0844868318367246 )
+ ( 0.102904358853672 , 0.0914416248339538 )
+ ( 0.06995886457571132 , 0.0972166291339016 )
+ ( 0.0336573561886173 , 0.1069703394877764 )
+ ( -0.004289272609839996 , 0.1129330599177 )
+ ( -0.04543704054779692 , 0.124436594734044 )
+ ( -0.09055229411010951 , 0.1334794772997432 )
+ ( -0.13906545545306717 , 0.14590433023743538 )
+ ( -0.1921634783990852 , 0.1592356281203136 )
+ ( -0.2499532167080883 , 0.1745908033846756 )
+ ( -0.31306886780850723 , 0.1930459138035904 )
+ ( -0.3821948764520447 , 0.2131302103410696 )
+ ( -0.45880969038541675 , 0.23864019562748257 )
+ ( -0.5426314607294217 , 0.26629355627414564 )
+ ( -0.6374452217842125 , 0.30196385436707 )
+ ( -0.7417742639164767 , 0.3421508090546344 )
+ ( -0.8608016627425203 , 0.39317190972679195 )
+ ( -0.9941122583780242 , 0.45411611308527194 )
+ ( -1.1465285813934336 , 0.530682040348918 )
+ ( -1.3215857509053397 , 0.628492182304824 )
+ ( -1.5210890331876734 , 0.7503518904666359 )
+ ( -1.7564491751987403 , 0.91926912101532 )
+ ( -2.022101522140445 , 1.128179080660086 )
+ ( -2.3379525273128863 , 1.4481519007414438 )
+ ( -2.6860300009391516 , 1.840934433099552 )
+ ( -3.0538312819720566 , 2.511362127727128 )
+ ( -3.4248475990483094 , 3.323750446248 )
+ ( -3.429351678875 , 4.74026629667472 )
+ ( -3.24873175766853 , 6.37239096969004 )
+ ( -1.38445729707063 , 7.90050270168416 )
+ ( 0.9202337980964693 , 9.20125829609504 )
+ ( 3.1923614071008903 , 8.0083294505548 )
+ ( 5.200868270044681 , 6.48976807543376 )
+ ( 5.389481015208544 , 4.910009396573518 )
+ ( 5.439614157410399 , 3.4247324884147194 )
+ ( 5.089834070019386 , 2.644146581021736 )
+ ( 4.727701207513338 , 1.922506424682774 )
+ ( 4.400433152414912 , 1.5484027457523057 )
+ ( 4.078921889689015 , 1.1989432696703282 )
+ ( 3.8294402542062276 , 0.9984297530452502 )
+ ( 3.5874365470640392 , 0.811931484854084 )
+ ( 3.399393394287642 , 0.69366257557086 )
+ ( 3.2188783067528077 , 0.584722040903896 )
+ ( 3.0742867218487673 , 0.50950169214669 )
+ ( 2.9369964846928385 , 0.441103992398784 )
+ ( 2.8230546982878777 , 0.39037046272590004 )
+ ( 2.716127917923021 , 0.34492524419459203 )
+ ( 2.62424398991363 , 0.30908449982984326 )
+ ( 2.5390775935854393 , 0.2775023586764044 )
+ ( 2.463461501702361 , 0.25122342232074646 )
+ ( 2.3942707614986856 , 0.22846730229542658 )
+ ( 2.3309393582835285 , 0.20860500019949857 )
+ ( 2.2737507921562057 , 0.19171562652773358 )
+ ( 2.2198969789968093 , 0.17632021438227 )
+ ( 2.171917054188519 , 0.16347309989074899 )
+ ( 2.1255165589415856 , 0.1512847447725438 )
+ ( 2.0847359367724327 , 0.141308495240978 )
+ ( 2.0443379806635686 , 0.131498891134892 )
+ ( 2.009236483354414 , 0.12359982880898822 )
+ ( 1.9744722111439765 , 0.11583354134820201 )
+ ( 1.9431887150450162 , 0.1092321656427364 )
+ ( 1.9129734054688714 , 0.102992148796463 )
+ ( 1.8848880638244556 , 0.0974189976962562 )
+ ( 1.8583927735848502 , 0.09234098650103119 )
+ ( 1.8330116877979346 , 0.0875928367580112 )
+ ( 1.8095927576758728 , 0.0834145464626946 )
+ ( 1.7865175100024524 , 0.07933599887988799 )
+ ( 1.765667550947678 , 0.07586478955595401 )
+ ( 1.7449383965552594 , 0.072429069032981 )
+ ( 1.7258883416491027 , 0.069427504243433 )
+ ( 1.707302540912686 , 0.06654678946807979 )
+ ( 1.6896620742411543 , 0.0638992661850522 )
+ ( 1.6728981160273444 , 0.061466594363694796 )
+ ( 1.6565019457552344 , 0.059121438346124995 )
+ ( 1.6412967952834274 , 0.05705443293865 )
+ ( 1.6261610412894691 , 0.055004043656375 )
+ ( 1.6121408913824589 , 0.053203353123279604 )
+ ( 1.598336906641767 , 0.051450999618952796 )
+ ( 1.5851292009039528 , 0.049827203550821195 )
+ ( 1.5724789224263418 , 0.0483230991350568 )
+ ( 1.5600064749339058 , 0.04685598599635799 )
+ ( 1.5483605431614902 , 0.045560510048340004 )
+ ( 1.5367577070391734 , 0.044273451758516394 )
+ ( 1.5257870965257398 , 0.0431137453907016 )
+ ( 1.515056657159423 , 0.0420018582725852 )
+ ( 1.5045907261931644 , 0.0409416254215892 )
+ ( 1.4946255261901655 , 0.039979484469390605 )
+ ( 1.4846932765288385 , 0.03902284814410801 )
+ ( 1.4753341919689693 , 0.038176572535969996 )
+ ( 1.4660717475358582 , 0.0373479309580838 )
+ ( 1.4570681572227873 , 0.0365678477170412 )
+ ( 1.4483969715085636 , 0.0358504110436288 )
+ ( 1.4397515780403363 , 0.0351366138118164 )
+ ( 1.4315780201280788 , 0.0345116881247668 )
+ ( 1.4234529070163335 , 0.03389458780158 )
+ ( 1.4155346474990393 , 0.033315704373481594 )
+ ( 1.4078702633012197 , 0.032783989324985004 )
+ ( 1.4002268775423472 , 0.0322546777260286 )
+ ( 1.3929392126217983 , 0.0317928404762148 )
+ ( 1.3857017802242346 , 0.03133907680392 )
+ ( 1.378601279830167 , 0.030910624492512598 )
+ ( 1.3717242329561528 , 0.0305247366302792 )
+ ( 1.364864362203135 , 0.0301403948868218 )
+ ( 1.3582461227329714 , 0.029803725562952 )
+ ( 1.3517035905684374 , 0.029480570775798003 )
+ ( 1.3452249869378365 , 0.029168989969302 )
+ ( 1.338962837935282 , 0.028901850818829997 )
+ ( 1.332715291923563 , 0.028635644511242998 )
+ ( 1.3266066854658303 , 0.028398878037701202 )
+ ( 1.3206049765829333 , 0.0281842302309696 )
+ ( 1.3146168718636133 , 0.0279702649333492 )
+ ( 1.3088167432733353 , 0.0278003620538208 )
+ ( 1.3030431010506038 , 0.027634395929356997 )
+ ( 1.2973235664587244 , 0.027480307625799998 )
+ ( 1.2917362976083313 , 0.027359763732548004 )
+ ( 1.2861613130947016 , 0.0272394871207368 )
+ ( 1.2806744162823347 , 0.0271435513663328 )
+ ( 1.2752591352477964 , 0.027066605718989202 )
+ ( 1.269855101683756 , 0.026989740393865598 )
+ ( 1.264557316416183 , 0.0269481553868072 )
+ ( 1.2592882655327051 , 0.0269131755949056 )
+ ( 1.2540389067778386 , 0.0268820241036404 )
+ ( 1.2488816988853564 , 0.026888250901969 )
+ ( 1.2437351118388418 , 0.026894223871653602 )
+ ( 1.238621040074352 , 0.026913335711270397 )
+ ( 1.2335638459900171 , 0.026960249086124002 )
+ ( 1.2285169871456798 , 0.0270067430837876 )
+ ( 1.223505221790627 , 0.027074932463126405 )
+ ( 1.2185259866666864 , 0.0271626772847928 )
+ ( 1.2135568879779164 , 0.02724983915602 )
+ ( 1.2086162750294718 , 0.0273667804999542 )
+ ( 1.203693469750129 , 0.027496092484252 )
+ ( 1.1987806863157806 , 0.0276246509562568 )
+ ( 1.1938825056702782 , 0.027790829566725603 )
+ ( 1.188995328487824 , 0.027963064238444003 )
+ ( 1.1841179226025922 , 0.028134356112973602 )
+ ( 1.1792351783859978 , 0.0283510332135992 )
+ 0
\ No newline at end of file
diff --git a/test_data/cluster/case_48_short/DCLU b/test_data/cluster/case_48_short/DCLU
new file mode 100644
index 0000000000000000000000000000000000000000..7b63fd42a24fa5a468fdee8d44f7bd9914d2efbf
--- /dev/null
+++ b/test_data/cluster/case_48_short/DCLU
@@ -0,0 +1,53 @@
+ 48 10 10 5760 1 149 300 0 0
+ 0.0000000e+00 1.0000000e-08 0.0000000e+00
+ 0.0000000e+00 -1.0000000e-08 0.0000000e+00
+ 8.5355339e-09 5.0000000e-09 0.0000000e+00
+ 8.5355339e-09 -5.0000000e-09 0.0000000e+00
+ -8.5355339e-09 5.0000000e-09 0.0000000e+00
+ -8.5355339e-09 -5.0000000e-09 0.0000000e+00
+ 0.0000000e+00 1.0000000e-08 1.0000000e-08
+ 0.0000000e+00 -1.0000000e-08 1.0000000e-08
+ 8.5355339e-09 5.0000000e-09 1.0000000e-08
+ 8.5355339e-09 -5.0000000e-09 1.0000000e-08
+ -8.5355339e-09 5.0000000e-09 1.0000000e-08
+ -8.5355339e-09 -5.0000000e-09 1.0000000e-08
+ 0.0000000e+00 1.0000000e-08 2.0000000e-08
+ 0.0000000e+00 -1.0000000e-08 2.0000000e-08
+ 8.5355339e-09 5.0000000e-09 2.0000000e-08
+ 8.5355339e-09 -5.0000000e-09 2.0000000e-08
+ -8.5355339e-09 5.0000000e-09 2.0000000e-08
+ -8.5355339e-09 -5.0000000e-09 2.0000000e-08
+ 0.0000000e+00 1.0000000e-08 3.0000000e-08
+ 0.0000000e+00 -1.0000000e-08 3.0000000e-08
+ 8.5355339e-09 5.0000000e-09 3.0000000e-08
+ 8.5355339e-09 -5.0000000e-09 3.0000000e-08
+ -8.5355339e-09 5.0000000e-09 3.0000000e-08
+ -8.5355339e-09 -5.0000000e-09 3.0000000e-08
+ -6.4278761e-09 7.6604444e-09 4.0000000e-08
+ 6.4278761e-09 -7.6604444e-09 4.0000000e-08
+ 3.3246603e-09 9.3167577e-09 4.0000000e-08
+ 9.7525364e-09 1.6563132e-09 4.0000000e-08
+ -9.7525364e-09 -1.6563132e-09 4.0000000e-08
+ -3.3246603e-09 -9.3167577e-09 4.0000000e-08
+ -8.6602540e-09 5.0000000e-09 5.0000000e-08
+ 8.6602540e-09 -5.0000000e-09 5.0000000e-08
+ -6.2360066e-11 9.8919892e-09 5.0000000e-08
+ 8.5978940e-09 4.8919892e-09 5.0000000e-08
+ -8.5978940e-09 -4.8919892e-09 5.0000000e-08
+ 6.2360066e-11 -9.8919892e-09 5.0000000e-08
+ -5.0000000e-09 8.6602540e-09 6.0000000e-08
+ 5.0000000e-09 -8.6602540e-09 6.0000000e-08
+ 4.8919892e-09 8.5978940e-09 6.0000000e-08
+ 9.8919892e-09 -6.2360066e-11 6.0000000e-08
+ -9.8919892e-09 6.2360066e-11 6.0000000e-08
+ -4.8919892e-09 -8.5978940e-09 6.0000000e-08
+ -2.5881905e-09 9.6592583e-09 7.0000000e-08
+ 2.5881905e-09 -9.6592583e-09 7.0000000e-08
+ 6.9505974e-09 7.0387879e-09 7.0000000e-08
+ 9.5387879e-09 -2.6204704e-09 7.0000000e-08
+ -9.5387879e-09 2.6204704e-09 7.0000000e-08
+ -6.9505974e-09 -7.0387879e-09 7.0000000e-08
+ 0.00D01 0.00D00 0.00D01 180.00D00 0.00D01 0.00D01
+ 0.00D01 0.00D00 0.00D01 0.00D00 0.00D01 0.00D01
+ 1
+0
diff --git a/test_data/cluster/case_48_short/DEDFB b/test_data/cluster/case_48_short/DEDFB
new file mode 100644
index 0000000000000000000000000000000000000000..1484cd7f32aedb538155051fa3be1e7eaff535e5
--- /dev/null
+++ b/test_data/cluster/case_48_short/DEDFB
@@ -0,0 +1,12 @@
+ 48 0
+ 1.000D00 3.0000000D08 1.000000D00 0 2 0 3
+4.2000000e-07
+4.2100000e-07
+ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 5.00d-9
+ 1.000000D00
+ ( 0.6965478043675424 , 0.012177869466280656 )
+ ( 0.6910987175052832 , 0.012554064506388863 )
+ 0
\ No newline at end of file