diff --git a/src/include/tfrfme.h b/src/include/tfrfme.h
index 6a6585c3431a0760ff1e01a9ca597dab9619703b..ba2954a742b5458f8cabb1ca64729befef232e20 100644
--- a/src/include/tfrfme.h
+++ b/src/include/tfrfme.h
@@ -10,9 +10,11 @@
*/
class Swap1 {
protected:
- //! NLMMT = 2 * LM * (LM + 2)
+ //! Index of the last element to be filled.
int last_index;
+ //! Number of vector coordinates. QUESTION: correct?
int nkv;
+ //! NLMMT = 2 * LM * (LM + 2)
int nlmmt;
//! QUESTION: definition?
@@ -57,6 +59,8 @@ public:
~Swap1() { delete[] wk; }
/*! \brief Append an element at the end of the vector.
+ *
+ * \param value: `complex<double>` The value to be added to the vector.
*/
void append(std::complex<double> value) { wk[last_index++] = value; }
@@ -109,7 +113,11 @@ public:
*/
class Swap2 {
protected:
- //! Number of vector coordinates
+ //! Index of the last vector element to be filled.
+ int last_vector;
+ //! Index of the last matrix element to be filled.
+ int last_matrix;
+ //! Number of vector coordinates. QUESTION: correct?
int nkv;
//! QUESTION: definition?
double *vkv;
@@ -178,7 +186,7 @@ public:
/*! \brief Swap2 instance destroyer.
*/
~Swap2();
-
+
/*! \brief Load a Swap2 instance from binary file.
*
* \param file_name: `string` Name of the file.
@@ -214,6 +222,26 @@ public:
*/
double *get_vector() { return vkv; }
+ /*! \brief Append an element at the end of the matrix.
+ *
+ * \param value: `double` The value to be pushed in the matrix.
+ */
+ void push_matrix(double value);
+
+ /*! \brief Append an element at the end of the vector.
+ *
+ * \param value: `double` The value to be pushed in the vector.
+ */
+ void push_vector(double value) { vkv[last_vector++] = value; }
+
+ /*! \brief Bring the matrix pointer to the start of the array.
+ */
+ void reset_matrix() { last_matrix = 0; }
+
+ /*! \brief Bring the vector pointer to the start of the array.
+ */
+ void reset_vector() { last_vector = 0; }
+
/*! \brief Set a parameter by its name and value.
*
* \param param_name: `string` Name of the parameter.
diff --git a/src/libnptm/tfrfme.cpp b/src/libnptm/tfrfme.cpp
index 8933c6577b5857b775bcc28a28ba95596cc1897b..b3478d3cc34850289b6a583a6db78fbd0ae0f663 100644
--- a/src/libnptm/tfrfme.cpp
+++ b/src/libnptm/tfrfme.cpp
@@ -207,6 +207,8 @@ Swap2::Swap2(int _nkv) {
vkv = new double[nkv]();
vkzm = new double*[nkv];
for (int vi = 0; vi < nkv; vi++) vkzm[vi] = new double[nkv]();
+ last_vector = 0;
+ last_matrix = 0;
}
Swap2::~Swap2() {
@@ -332,7 +334,7 @@ Swap2* Swap2::from_legacy(string file_name) {
}
long Swap2::get_memory_requirement(int _nkv) {
- long size = (long)(3 * sizeof(int) + 11 * sizeof(double));
+ long size = (long)(5 * sizeof(int) + 11 * sizeof(double));
size += (long)(sizeof(double) * _nkv * (_nkv + 1));
return size;
}
@@ -360,6 +362,13 @@ double Swap2::get_param(string param_name) {
return value;
}
+void Swap2::push_matrix(double value) {
+ int col = last_matrix % (nkv - 1);
+ int row = last_matrix - (nkv * row);
+ vkzm[row][col] = value;
+ last_matrix++;
+}
+
void Swap2::set_param(string param_name, double value) {
if (param_name.compare("nkv") == 0) nkv = (int)value;
else if (param_name.compare("apfafa") == 0) apfafa = value;
@@ -644,8 +653,6 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
instance->set_param("spd", spd);
instance->set_param("frsh", frsh);
instance->set_param("exril", exril);
- // Vectors and matrix need to be copied element-wise and
- // subsequently deleted.
str_type = "FLOAT64_(" + to_string(nxv) + ")";
status = hdf_file->read("XVEC", str_type, instance->xv);
str_type = "FLOAT64_(" + to_string(nyv) + ")";
@@ -658,13 +665,14 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
elements = new double[num_elements]();
str_type = "FLOAT64_(" + to_string(num_elements) + ")";
status = hdf_file->read("WSUM", str_type, elements);
- for (int wi = 0; wi < _nlmmt; wi++) {
- for (int wj = 0; wj < _nrvc; wj++) {
- int index = (2 * _nrvc) * wi + 2 * wj;
+ int index = 0;
+ for (int wj = 0; wj < _nrvc; wj++) {
+ for (int wi = 0; wi < _nlmmt; wi++) {
value = complex<double>(elements[index], elements[index + 1]);
_wsum[wi][wj] = value;
- } // wj loop
- } // wi loop
+ index += 2;
+ } // wi loop
+ } // wj loop
delete[] elements;
status = hdf_file->close();
delete hdf_file;
@@ -723,14 +731,14 @@ TFRFME* TFRFME::from_legacy(string file_name) {
}
int _nlmmt = 2 * lm * (lm + 2);
int _nrvc = nxv * nyv * nzv;
- for (int wi = 0; wi < _nlmmt; wi++) {
- for (int wj = 0; wj < _nrvc; wj++) {
+ for (int wj = 0; wj < _nrvc; wj++) {
+ for (int wi = 0; wi < _nlmmt; wi++) {
input.read(reinterpret_cast<char *>(&rval), sizeof(double));
input.read(reinterpret_cast<char *>(&ival), sizeof(double));
complex<double> value(rval, ival);
_wsum[wi][wj] = value;
- } // wj loop
- } // wi loop
+ } // wi loop
+ } // wj loop
input.close();
} else {
printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
@@ -860,19 +868,19 @@ void TFRFME::write_hdf5(string file_name) {
rec_type_list.append(str_type);
rec_ptr_list.append(zv);
rec_name_list.append("WSUM");
- str_type = "FLOAT64_(" + to_string(nlmmt) + "," + to_string(2 * nrvc) + ")";
- rec_type_list.append(str_type);
- // The (N x M) matrix of complex is converted to a (N x 2M) matrix of double,
- // where REAL(E_i,j) -> E_i,(2 j) and IMAG(E_i,j) -> E_i,(2 j + 1)
int num_elements = 2 * nlmmt * nrvc;
+ str_type = "FLOAT64_(" + to_string(num_elements) + ")";
+ rec_type_list.append(str_type);
+ // The (N x M) matrix of complex is converted to a vector of double
+ // with length (2 * N * M)
double *ptr_elements = new double[num_elements]();
- for (int wi = 0; wi < nlmmt; wi++) {
- for (int wj = 0; wj < nrvc; wj++) {
- int index = (2 * nrvc) * wi + 2 * wj;
- ptr_elements[index] = wsum[wi][wj].real();
- ptr_elements[index + 1] = wsum[wi][wj].imag();
- } // wj loop
- } // wi loop
+ int index = 0;
+ for (int wj = 0; wj < nrvc; wj++) {
+ for (int wi = 0; wi < nlmmt; wi++) {
+ ptr_elements[index++] = wsum[wi][wj].real();
+ ptr_elements[index++] = wsum[wi][wj].imag();
+ } // wi loop
+ } // wj loop
rec_ptr_list.append(ptr_elements);
string *rec_names = rec_name_list.to_array();
@@ -916,14 +924,14 @@ void TFRFME::write_legacy(string file_name) {
output.write(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));
for (int zi = 0; zi < nzv; zi++)
output.write(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
- for (int wi = 0; wi < nlmmt; wi++) {
- for (int wj = 0; wj < nrvc; wj++) {
+ for (int wj = 0; wj < nrvc; wj++) {
+ for (int wi = 0; wi < nlmmt; wi++) {
double rval = wsum[wi][wj].real();
double ival = wsum[wi][wj].imag();
output.write(reinterpret_cast<char *>(&rval), sizeof(double));
output.write(reinterpret_cast<char *>(&ival), sizeof(double));
- } // wj loop
- } // wi loop
+ } // wi loop
+ } // wj loop
output.close();
} else { // Should never occur
printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
index d96a4fad6c2733531fc046c01c897759e48431dc..c4866e1f65073bcc2f61ec54d7a23a125ac315ca 100644
--- a/src/trapping/cfrfme.cpp
+++ b/src/trapping/cfrfme.cpp
@@ -68,6 +68,8 @@ void frfme(string data_file, string output_path) {
double apfafa = 0.0, pmf = 0.0, spd = 0.0, rir = 0.0, ftcn = 0.0, fshmx = 0.0;
double vxyzmx = 0.0, delxyz = 0.0, vknmx = 0.0, delk = 0.0, delks = 0.0;
double frsh = 0.0, exril = 0.0;
+ double **vkzm = NULL;
+ double *vkv = NULL;
int nlmmt = 0, nrvc = 0;
// Vector size variables
int wsum_size;
@@ -90,19 +92,11 @@ void frfme(string data_file, string output_path) {
spd = tfrfme->get_param("spd");
frsh = tfrfme->get_param("frsh");
exril = tfrfme->get_param("exril");
- //fstream temptape2;
string tempname2 = output_path + "/c_TEMPTAPE2.hd5";
- //temptape2.open(tempname2.c_str(), ios::in | ios::binary);
if (tt2 == NULL) tt2 = Swap2::from_binary(tempname2, "HDF5");
if (tt2 != NULL) {
- /*for (int jx = 0; jx < nkv; jx++) temptape2.read(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
- vkzm = new double*[nkv];
- for (int vki = 0; vki < nkv; vki++) vkzm[vki] = new double[nkv]();
- for (int jy10 = 0; jy10 < nkv; jy10++) {
- for (int jx10 = 0; jx10 < nkv; jx10++) {
- temptape2.read(reinterpret_cast<char *>(&(vkzm[jx10][jy10])), sizeof(double));
- } //jx10 loop
- } // jy10 loop*/
+ vkv = tt2->get_vector();
+ vkzm = tt2->get_matrix();
apfafa = tt2->get_param("apfafa");
pmf = tt2->get_param("pmf");
spd = tt2->get_param("spd");
@@ -223,9 +217,6 @@ void frfme(string data_file, string output_path) {
nks = nksh * 2;
nkv = nks + 1;
// Array initialization
- /*vkv = new double[nkv]();
- vkzm = new double*[nkv];
- for (int vi = 0; vi < nkv; vi++) vkzm[vi] = new double[nkv];*/
long swap1_size, swap2_size, tfrfme_size;
double size_mb;
printf("INFO: calculation memory requirements\n");
@@ -236,8 +227,8 @@ void frfme(string data_file, string output_path) {
size_mb = 1.0 * swap2_size / 1024.0 / 1024.0;
printf("Swap 2: %.2lg MB\n", size_mb);
tt2 = new Swap2(nkv);
- double *_vkv = tt2->get_vector();
- double **_vkzm = tt2->get_matrix();
+ vkv = tt2->get_vector();
+ vkzm = tt2->get_matrix();
// End of array initialization
double vkm = vk * exri;
vknmx = vk * an;
@@ -280,8 +271,8 @@ void frfme(string data_file, string output_path) {
int nrvc = nxv * nyv * nzv;
int nkshpo = nksh + 1;
for (int i28 = nkshpo; i28 <= nks; i28++) {
- _vkv[i28] = _vkv[i28 - 1] + delk;
- _vkv[nkv - i28 - 1] = -_vkv[i28];
+ vkv[i28] = vkv[i28 - 1] + delk;
+ vkv[nkv - i28 - 1] = -vkv[i28];
} // i28 loop
if (tfrfme != NULL) {
tfrfme->set_param("vk", vk);
@@ -297,9 +288,6 @@ void frfme(string data_file, string output_path) {
string temp_name2 = output_path + "/c_TEMPTAPE2.hd5";
//temptape1.open(temp_name1.c_str(), ios::out | ios::binary);
tt1 = new Swap1(lm, nkv);
- /*temptape2.open(temp_name2.c_str(), ios::out | ios::binary);
- for (int jx = 0; jx < nkv; jx++)
- temptape2.write(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));*/
if (wk != NULL) delete[] wk;
wk = frfmer(nkv, vkm, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, tt1, tt2);
tt1->write_binary(temp_name1, "HDF5");
@@ -318,21 +306,7 @@ void frfme(string data_file, string output_path) {
tt2->set_param("nlmmt", 1.0 * nlmmt);
tt2->set_param("nrvc", 1.0 * nrvc);
tt2->write_binary(temp_name2, "HDF5");
- //temptape2.open("c_TEMPTAPE2", ios::in | ios::binary);
- /*for (int jx = 0; jx < nkv; jx++) {
- double value = 0.0;
- temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));
- vkv[jx] = value;
- }
- for (int jy40 = 0; jy40 < nkv; jy40++) {
- for (int jx40 = 0; jx40 < nkv; jx40++) {
- double value = 0.0;
- temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));
- vkzm[jx40][jy40] = value;
- }
- } // jy40 loop
- temptape2.close();*/
- for (int j80 = jlmf - 1; j80 < jlml; j80++) {
+ for (int j80 = jlmf; j80 <= jlml; j80++) {
complex<double> *tt1_wk = tt1->get_vector();
int wk_index = 0;
// w matrix
@@ -342,24 +316,14 @@ void frfme(string data_file, string output_path) {
}
w = new complex<double>*[nkv];
for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv]();
- //if (wk != NULL) delete[] wk;
- //wk = new complex<double>[nlmmt]();
- //temptape1.open(temp_name1.c_str(), ios::in | ios::binary);
for (int jy50 = 0; jy50 < nkv; jy50++) {
for (int jx50 = 0; jx50 < nkv; jx50++) {
- /*for (int i = 0; i < nlmmt; i++) {
- double vreal, vimag;
- temptape1.read(reinterpret_cast<char *>(&vreal), sizeof(double));
- temptape1.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- wk[i] = complex<double>(vreal, vimag);
- }*/
for (int wi = 0; wi < nlmmt; wi++) wk[wi] = tt1_wk[wk_index++];
- w[jx50][jy50] = wk[j80];
+ w[jx50][jy50] = wk[j80 - 1];
} // jx50
} // jy50 loop
- //temptape1.close();
int ixyz = 0;
- for (int wj = 0; wj < nrvc; wj++) _wsum[j80][wj] = cc0;
+ for (int wj = 0; wj < nrvc; wj++) _wsum[j80 - 1][wj] = cc0;
for (int iz75 = 0; iz75 < nzv; iz75++) {
double z = _zv[iz75] + frsh;
for (int iy70 = 0; iy70 < nyv; iy70++) {
@@ -369,43 +333,27 @@ void frfme(string data_file, string output_path) {
ixyz++;
complex<double> sumy = cc0;
for (int jy60 = 0; jy60 < nkv; jy60++) {
- double vky = _vkv[jy60];
- double vkx = _vkv[nkv - 1];
- double vkzf = _vkzm[0][jy60];
+ double vky = vkv[jy60];
+ double vkx = vkv[nkv - 1];
+ double vkzf = vkzm[0][jy60];
complex<double> phasf = exp(uim * (-vkx * x + vky * y + vkzf * z));
- double vkzl = _vkzm[nkv - 1][jy60];
+ double vkzl = vkzm[nkv - 1][jy60];
complex<double> phasl = exp(uim * (vkx * x + vky * y + vkzl * z));
complex<double> sumx = 0.5 * (w[0][jy60] * phasf + w[nkv - 1][jy60] * phasl);
for (int jx55 = 2; jx55 <= nks; jx55++) {
- vkx = _vkv[jx55 - 1];
- double vkz = _vkzm[jx55 - 1][jy60];
+ vkx = vkv[jx55 - 1];
+ double vkz = vkzm[jx55 - 1][jy60];
complex<double> phas = exp(uim * (vkx * x + vky * y + vkz * z));
sumx += (w[jx55 - 1][jy60] * phas);
} // jx55 loop
if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;
sumy += sumx;
} // jy60 loop
- _wsum[j80][ixyz - 1] = sumy * delks;
+ _wsum[j80 - 1][ixyz - 1] = sumy * delks;
} // ix65 loop
} // iy70 loop
} // iz75 loop
} // j80 loop
- if (jlmf != 1) {
- lmode = (int)tfrfme->get_param("lmode");
- lm = (int)tfrfme->get_param("lm");
- nkv = (int)tfrfme->get_param("nkv");
- nxv = (int)tfrfme->get_param("nxv");
- nyv = (int)tfrfme->get_param("nyv");
- nzv = (int)tfrfme->get_param("nzv");
- vk = tfrfme->get_param("vk");
- exri = tfrfme->get_param("exri");
- an = tfrfme->get_param("an");
- ff = tfrfme->get_param("ff");
- tra = tfrfme->get_param("tra");
- spd = tfrfme->get_param("spd");
- frsh = tfrfme->get_param("frsh");
- exril = tfrfme->get_param("exril");
- }
// label 88
tfrfme->write_binary(tfrfme_name, "HDF5");
string output_name = output_path + "/c_OFRFME";
diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp
index e2c0793e6b52a2c83b2067cb6a97dc9d572dbf65..7e3619a252ada2d7102a4d88b3b2a2ba97e60f19 100644
--- a/src/trapping/clffft.cpp
+++ b/src/trapping/clffft.cpp
@@ -281,7 +281,7 @@ void lffft(string data_file, string output_path) {
//binary_input.read(reinterpret_cast<char *>(&vreal), sizeof(double));
//binary_input.read(reinterpret_cast<char *>(&vimag), sizeof(double));
int row = i;
- int col = (nzv * iz475) + (nyv * iy475) + ix475;
+ int col = (nyv * nxv * iz475) + (nxv * iy475) + ix475;
complex<double> value = _wsum[row][col];
if (lm <= le) {
//ws[i] = complex<double>(vreal, vimag);