From cfd0fa52fcd73f4de935a995084fa7eb636d680a Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Tue, 27 Feb 2024 20:22:12 +0100
Subject: [PATCH] Fix swapping progression on Swap1
---
src/include/tfrfme.h | 120 ++++++++++++---------------------------
src/include/tra_subs.h | 32 +++++------
src/libnptm/tfrfme.cpp | 111 ++++++++++++++++++++----------------
src/libnptm/tra_subs.cpp | 85 +++++++++++++--------------
src/trapping/cfrfme.cpp | 61 ++++++++++++--------
src/trapping/clffft.cpp | 19 ++++---
6 files changed, 201 insertions(+), 227 deletions(-)
diff --git a/src/include/tfrfme.h b/src/include/tfrfme.h
index 8e92f7af..6a6585c3 100644
--- a/src/include/tfrfme.h
+++ b/src/include/tfrfme.h
@@ -11,6 +11,8 @@
class Swap1 {
protected:
//! NLMMT = 2 * LM * (LM + 2)
+ int last_index;
+ int nkv;
int nlmmt;
//! QUESTION: definition?
@@ -46,13 +48,18 @@ public:
/*! \brief Swap1 instance constructor.
*
* \param lm: `int` Maximum field expansion order.
+ * \param _nkv: `int` Number of vector coordinates. QUESTION: correct?
*/
- Swap1(int lm);
+ Swap1(int lm, int _nkv);
/*! \brief Swap1 instance destroyer.
*/
~Swap1() { delete[] wk; }
+ /*! \brief Append an element at the end of the vector.
+ */
+ void append(std::complex<double> value) { wk[last_index++] = value; }
+
/*! \brief Load a Swap1 instance from binary file.
*
* \param file_name: `string` Name of the file.
@@ -62,27 +69,24 @@ public:
*/
static Swap1* from_binary(std::string file_name, std::string mode="LEGACY");
- /*! \brief Get an element from the WK vector.
- *
- * \param index: `int` Index of the desired element.
- * \return value: `complex<double>` The value of the requested element.
- */
- std::complex<double> get_element(int index) { return wk[index]; }
-
/*! \brief Calculate the necessary amount of memory to create a new instance.
*
* \param lm: `int` Maximum field expansion order.
+ * \param _nkv: `int` Number of vector coordinates. QUESTION: correct?
* \return size: `long` The necessary memory size in bytes.
*/
- static long get_memory_requirement(int lm);
+ static long get_memory_requirement(int lm, int _nkv);
- /*! \brief Set an element in the WK vector.
+ /*! \brief Get the pointer to the WK vector.
*
- * \param index: `int` Index of the desired element.
- * \param value: `complex<double>` The value to be stored in the vector.
+ * \return value: `complex<double> *` Memory address of the WK vector.
*/
- void set_element(int index, std::complex<double> value) { wk[index] = value; }
-
+ std::complex<double> *get_vector() { return wk; }
+
+ /*! \brief Bring the pointer to the next element at the start of vector.
+ */
+ void reset() { last_index = 0; }
+
/*! \brief Write a Swap1 instance to binary file.
*
* \param file_name: `string` Name of the file.
@@ -184,13 +188,11 @@ public:
*/
static Swap2* from_binary(std::string file_name, std::string mode="LEGACY");
- /*! \brief Get an element from the VKZM matrix.
+ /*! \brief Get the pointer to the VKZM matrix.
*
- * \param row: `int` Index of the desired element row.
- * \param col: `int` Index of the desired element column.
- * \return value: `double` The value of the requested element.
+ * \return value: `double **` Pointer to the VKZM matrix.
*/
- double get_matrix_element(int row, int col) { return vkzm[row][col]; }
+ double **get_matrix() { return vkzm; }
/*! \brief Calculate the necessary amount of memory to create a new instance.
*
@@ -206,21 +208,12 @@ public:
*/
double get_param(std::string param_name);
- /*! \brief Get an element from the VKV vector.
+ /*! \brief Get the pointer to the VKV vector.
*
- * \param index: `int` Index of the desired element.
- * \return value: `double` The value of the requested element.
+ * \return value: `double *` Pointer to the VKV vector.
*/
- double get_vector_element(int index) { return vkv[index]; }
+ double *get_vector() { return vkv; }
- /*! \brief Set an element in the VKZM matrix.
- *
- * \param row: `int` Index of the desired element row.
- * \param col: `int` Index of the desired element column.
- * \param value: `double` The value of the element to be stored.
- */
- void set_matrix_element(int row, int col, double value) { vkzm[row][col] = value; }
-
/*! \brief Set a parameter by its name and value.
*
* \param param_name: `string` Name of the parameter.
@@ -228,13 +221,6 @@ public:
*/
void set_param(std::string param_name, double value);
- /*! \brief Get an element from the VKV vector.
- *
- * \param index: `int` Index of the desired element.
- * \param value: `double` The value of element to be stored.
- */
- void set_vector_element(int index, double value) { vkv[index] = value; }
-
/*! \brief Write a Swap2 instance to binary file.
*
* \param file_name: `string` Name of the file.
@@ -355,13 +341,11 @@ public:
*/
static TFRFME* from_binary(std::string file_name, std::string mode="LEGACY");
- /*! \brief Get an element from the WSUM matrix.
+ /*! \brief Get the pointer to the WSUM matrix.
*
- * \param row: `int` Row index of the element to be read.
- * \param col: `int` Column index of the element to be read.
- * \return value: `complex<double>` The value of the requested element.
+ * \return value: `complex<double> **` Pointer to the WSUM matrix.
*/
- std::complex<double> get_matrix_element(int row, int col);
+ std::complex<double> **get_matrix() { return wsum; }
/*! \brief Calculate the necessary amount of memory to create a new instance.
*
@@ -385,35 +369,24 @@ public:
*/
double get_param(std::string param_name);
- /*! \brief Get the X coordinate of the computed point at given index.
+ /*! \brief Get the pointer to the X coordinates vector.
*
- * \param index: `int` Index of the requested point.
- * \return x: `double` X coordinate of the requested point.
+ * \return x: `double *` Pointer to X coordinates vector.
*/
- double get_x(int index) { return xv[index]; }
+ double *get_x() { return xv; }
- /*! \brief Get the Y coordinate of the computed point at given index.
+ /*! \brief Get the pointer to the Y coordinates vector.
*
- * \param index: `int` Index of the requested point.
- * \return y: `double` Y coordinate of the requested point.
+ * \return y: `double *` Pointer to Y coordinates vector.
*/
- double get_y(int index) { return yv[index]; }
+ double *get_y() { return yv; }
- /*! \brief Get the Z coordinate of the computed point at given index.
+ /*! \brief Get the pointer to the Z coordinates vector.
*
- * \param index: `int` Index of the requested point.
- * \return z: `double` Z coordinate of the requested point.
+ * \return z: `double *` Pointer to Z coordinates vector.
*/
- double get_z(int index) { return zv[index]; }
+ double *get_z() { return zv; }
- /*! \brief Set an element in the WSUM matrix.
- *
- * \param row: `int` Row index of the element to be read.
- * \param col: `int` Column index of the element to be read.
- * \param value: `complex<double>` The value to be placed in the matrix.
- */
- void set_matrix_element(int row, int col, std::complex<double> value);
-
/*! \brief Set a configuration parameter.
*
* \param param_name: `string` Name of the parameter.
@@ -421,27 +394,6 @@ public:
*/
void set_param(std::string param_name, double value);
- /*! \brief Set the X coordinate of the computed point at given index.
- *
- * \param index: `int` Index of the requested point.
- * \param value: `double` X coordinate of the requested point.
- */
- void set_x(int index, double value) { xv[index] = value; }
-
- /*! \brief Set the Y coordinate of the computed point at given index.
- *
- * \param index: `int` Index of the requested point.
- * \param value: `double` Y coordinate of the requested point.
- */
- void set_y(int index, double value) { yv[index] = value; }
-
- /*! \brief Set the Z coordinate of the computed point at given index.
- *
- * \param index: `int` Index of the requested point.
- * \param value: `double` Z coordinate of the requested point.
- */
- void set_z(int index, double value) { zv[index] = value; }
-
/*! \brief Write a trapping configuration instance to binary file.
*
* \param file_name: `string` Name of the file.
diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index 7038228c..c32aa9ae 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -90,19 +90,6 @@ void ffrf(
double *fffs, CIL *cil, CCR *ccr
);
-/*! C++ porting of FFRT
- *
- * \param ac: Vector of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
- * \param ffte: `double *`. QUESTION: definition?
- * \param ffts: `double *`. QUESTION: definition?
- * \param cil: `CIL *` Pointer to a CIL structure.
- */
-void ffrt(
- std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
- CIL *cil
-);
-
/*! C++ porting of FRFMER
*
* \param nkv: `int` QUESTION: definition?
@@ -119,12 +106,25 @@ void ffrt(
* \param tt1: `Swap1 *` Pointer to first swap object.
* \param tt2: `Swap2 *` Pointer to second swap object.
*/
-void frfmer(
+std::complex<double> *frfmer(
int nkv, double vkm, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
Swap1 *tt1, Swap2 *tt2
);
+/*! C++ porting of FFRT
+ *
+ * \param ac: Vector of complex. QUESTION: definition?
+ * \param ws: Vector of complex. QUESTION: definition?
+ * \param ffte: `double *`. QUESTION: definition?
+ * \param ffts: `double *`. QUESTION: definition?
+ * \param cil: `CIL *` Pointer to a CIL structure.
+ */
+void ffrt(
+ std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
+ CIL *cil
+);
+
/*! C++ porting of PWMALP
*
* \param w: Matrix of complex. QUESTION: definition?
@@ -166,7 +166,7 @@ void sampoa(
/*! C++ porting of WAMFF
*
- * \param tt1: `Swap1 *` Pointer to the swap object containing the WK vector.
+ * \param wk: `complex<double> *` QUESTION: definition?
* \param x: `double`
* \param y: `double`
* \param z: `double`
@@ -180,6 +180,6 @@ void sampoa(
* \param pmf: `double` QUESTION: definition?
*/
void wamff(
- Swap1 *tt1, double x, double y, double z, int lm, double apfafa,
+ std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
double tra, double spd, double rir, double ftcn, int lmode, double pmf
);
diff --git a/src/libnptm/tfrfme.cpp b/src/libnptm/tfrfme.cpp
index df4f435e..8933c657 100644
--- a/src/libnptm/tfrfme.cpp
+++ b/src/libnptm/tfrfme.cpp
@@ -28,9 +28,12 @@
using namespace std;
// >>> START OF Swap1 CLASS IMPLEMENTATION <<<
-Swap1::Swap1(int lm) {
+Swap1::Swap1(int lm, int _nkv) {
+ nkv = _nkv;
nlmmt = 2 * lm * (lm + 2);
- wk = new complex<double>[nlmmt]();
+ const int size = nkv * nkv * nlmmt;
+ wk = new complex<double>[size]();
+ last_index = 0;
}
Swap1* Swap1::from_binary(string file_name, string mode) {
@@ -53,20 +56,23 @@ Swap1* Swap1::from_hdf5(string file_name) {
herr_t status = hdf_file->get_status();
double *elements;
string str_type;
- int _nlmmt, lm, num_elements, index;
+ int _nlmmt, _nkv, lm, num_elements, index;
complex<double> value;
+ complex<double> *_wk = NULL;
if (status == 0) {
status = hdf_file->read("NLMMT", "INT32", &_nlmmt);
+ status = hdf_file->read("NKV", "INT32", &_nkv);
lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0);
- num_elements = 2 * _nlmmt;
- instance = new Swap1(lm);
+ num_elements = 2 * _nlmmt * _nkv * _nkv;
+ instance = new Swap1(lm, _nkv);
+ _wk = instance->get_vector();
elements = new double[num_elements]();
str_type = "FLOAT64_(" + to_string(num_elements) + ")";
status = hdf_file->read("WK", str_type, elements);
- for (int wi = 0; wi < _nlmmt; wi++) {
+ for (int wi = 0; wi < num_elements / 2; wi++) {
index = 2 * wi;
value = complex<double>(elements[index], elements[index + 1]);
- instance->set_element(wi, value);
+ _wk[wi] = value;
} // wi loop
delete[] elements;
status = hdf_file->close();
@@ -78,17 +84,21 @@ Swap1* Swap1::from_hdf5(string file_name) {
Swap1* Swap1::from_legacy(string file_name) {
fstream input;
Swap1 *instance = NULL;
- int _nlmmt, lm;
+ complex<double> *_wk = NULL;
+ int _nlmmt, _nkv, lm;
double rval, ival;
input.open(file_name.c_str(), ios::in | ios::binary);
if (input.is_open()) {
input.read(reinterpret_cast<char *>(&_nlmmt), sizeof(int));
lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0);
- instance = new Swap1(lm);
- for (int j = 0; j < _nlmmt; j++) {
+ input.read(reinterpret_cast<char *>(&_nkv), sizeof(int));
+ instance = new Swap1(lm, _nkv);
+ _wk = instance->get_vector();
+ int num_elements = _nlmmt * _nkv * _nkv;
+ for (int j = 0; j < num_elements; j++) {
input.read(reinterpret_cast<char *>(&rval), sizeof(double));
input.read(reinterpret_cast<char *>(&ival), sizeof(double));
- instance->set_element(j, complex<double>(rval, ival));
+ _wk[j] = complex<double>(rval, ival);
}
input.close();
} else {
@@ -97,9 +107,9 @@ Swap1* Swap1::from_legacy(string file_name) {
return instance;
}
-long Swap1::get_memory_requirement(int lm) {
- long size = (long)sizeof(int);
- size += (long)(sizeof(complex<double>) * 2 * lm * (lm + 2));
+long Swap1::get_memory_requirement(int lm, int _nkv) {
+ long size = (long)(3 * sizeof(int));
+ size += (long)(sizeof(complex<double>) * 2 * lm * (lm + 2) * _nkv * _nkv);
return size;
}
@@ -120,15 +130,18 @@ void Swap1::write_hdf5(string file_name) {
List<void *> rec_ptr_list(1);
herr_t status;
string str_type;
- int num_elements = 2 * nlmmt;
+ int num_elements = 2 * nlmmt * nkv * nkv;
rec_name_list.set(0, "NLMMT");
rec_type_list.set(0, "INT32_(1)");
rec_ptr_list.set(0, &nlmmt);
+ rec_name_list.append("NKV");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&nkv);
rec_name_list.append("WK");
str_type = "FLOAT64_(" + to_string(num_elements) + ")";
rec_type_list.append(str_type);
double *ptr_elements = new double[num_elements]();
- for (int wi = 0; wi < nlmmt; wi++) {
+ for (int wi = 0; wi < num_elements / 2; wi++) {
ptr_elements[2 * wi] = wk[wi].real();
ptr_elements[2 * wi + 1] = wk[wi].imag();
}
@@ -156,8 +169,10 @@ void Swap1::write_legacy(string file_name) {
double rval, ival;
output.open(file_name.c_str(), ios::out | ios::binary);
if (output.is_open()) {
+ int num_elements = nlmmt * nkv * nkv;
output.write(reinterpret_cast<char *>(&nlmmt), sizeof(int));
- for (int j = 0; j < nlmmt; j++) {
+ output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
+ for (int j = 0; j < num_elements; j++) {
rval = wk[j].real();
ival = wk[j].imag();
output.write(reinterpret_cast<char *>(&rval), sizeof(double));
@@ -173,7 +188,11 @@ bool Swap1::operator ==(Swap1 &other) {
if (nlmmt != other.nlmmt) {
return false;
}
- for (int i = 0; i < nlmmt; i++) {
+ if (nkv != other.nkv) {
+ return false;
+ }
+ const int num_elements = nlmmt * nkv * nkv;
+ for (int i = 0; i < num_elements; i++) {
if (wk[i] != other.wk[i]) {
return false;
}
@@ -260,19 +279,23 @@ Swap2* Swap2::from_legacy(string file_name) {
fstream input;
Swap2 *instance = NULL;
int _nkv, _nlmmt, _nrvc;
+ double **_vkzm = NULL;
+ double *_vkv = NULL;
double value;
input.open(file_name.c_str(), ios::in | ios::binary);
if (input.is_open()) {
input.read(reinterpret_cast<char *>(&_nkv), sizeof(int));
instance = new Swap2(_nkv);
+ _vkzm = instance->get_matrix();
+ _vkv = instance->get_vector();
for (int vj = 0; vj < _nkv; vj++) {
input.read(reinterpret_cast<char *>(&value), sizeof(double));
- instance->set_vector_element(vj, value);
+ _vkv[vj] = value;
}
for (int mi = 0; mi < _nkv; mi++) {
for (int mj = 0; mj < _nkv; mj++) {
input.read(reinterpret_cast<char *>(&value), sizeof(double));
- instance->set_matrix_element(mi, mj, value);
+ _vkzm[mi][mj] = value;
}
}
input.read(reinterpret_cast<char *>(&value), sizeof(double));
@@ -589,10 +612,11 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
unsigned int flags = H5F_ACC_RDONLY;
HDFFile *hdf_file = new HDFFile(file_name, flags);
herr_t status = hdf_file->get_status();
- double *coord_vec, *elements;
+ double *elements;
string str_type;
int _nlmmt, _nrvc, num_elements;
complex<double> value;
+ complex<double> **_wsum = NULL;
if (status == 0) {
int lmode, lm, nkv, nxv, nyv, nzv;
double vk, exri, an, ff, tra, spd, frsh, exril;
@@ -611,6 +635,7 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
status = hdf_file->read("FRSH", "FLOAT64", &frsh);
status = hdf_file->read("EXRIL", "FLOAT64", &exril);
instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
+ _wsum = instance->get_matrix();
instance->set_param("vk", vk);
instance->set_param("exri", exri);
instance->set_param("an", an);
@@ -621,24 +646,12 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
instance->set_param("exril", exril);
// Vectors and matrix need to be copied element-wise and
// subsequently deleted.
- coord_vec = new double[nxv]();
str_type = "FLOAT64_(" + to_string(nxv) + ")";
- status = hdf_file->read("XVEC", str_type, coord_vec);
- for (int xi = 0; xi < nxv; xi++)
- instance->set_x(xi, coord_vec[xi]);
- delete[] coord_vec;
- coord_vec = new double[nyv]();
+ status = hdf_file->read("XVEC", str_type, instance->xv);
str_type = "FLOAT64_(" + to_string(nyv) + ")";
- status = hdf_file->read("YVEC", str_type, coord_vec);
- for (int yi = 0; yi < nyv; yi++)
- instance->set_y(yi, coord_vec[yi]);
- delete[] coord_vec;
- coord_vec = new double[nzv]();
+ status = hdf_file->read("YVEC", str_type, instance->yv);
str_type = "FLOAT64_(" + to_string(nzv) + ")";
- status = hdf_file->read("ZVEC", str_type, coord_vec);
- for (int zi = 0; zi < nzv; zi++)
- instance->set_z(zi, coord_vec[zi]);
- delete[] coord_vec;
+ status = hdf_file->read("ZVEC", str_type, instance->zv);
_nlmmt = 2 * lm * (lm + 2);
_nrvc = nxv * nyv * nzv;
num_elements = 2 * _nlmmt * _nrvc;
@@ -649,7 +662,7 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
for (int wj = 0; wj < _nrvc; wj++) {
int index = (2 * _nrvc) * wi + 2 * wj;
value = complex<double>(elements[index], elements[index + 1]);
- instance->set_matrix_element(wi, wj, value);
+ _wsum[wi][wj] = value;
} // wj loop
} // wi loop
delete[] elements;
@@ -662,6 +675,8 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
TFRFME* TFRFME::from_legacy(string file_name) {
fstream input;
TFRFME *instance = NULL;
+ complex<double> **_wsum = NULL;
+ double *coord_vec = NULL;
input.open(file_name.c_str(), ios::in | ios::binary);
if (input.is_open()) {
int lmode, lm, nkv, nxv, nyv, nzv;
@@ -682,6 +697,7 @@ TFRFME* TFRFME::from_legacy(string file_name) {
input.read(reinterpret_cast<char *>(&frsh), sizeof(double));
input.read(reinterpret_cast<char *>(&exril), sizeof(double));
instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
+ _wsum = instance->get_matrix();
instance->set_param("vk", vk);
instance->set_param("exri", exri);
instance->set_param("an", an);
@@ -690,17 +706,20 @@ TFRFME* TFRFME::from_legacy(string file_name) {
instance->set_param("spd", spd);
instance->set_param("frsh", frsh);
instance->set_param("exril", exril);
+ coord_vec = instance->get_x();
for (int xi = 0; xi < nxv; xi++) {
input.read(reinterpret_cast<char *>(&dval), sizeof(double));
- instance->set_x(xi, dval);
+ coord_vec[xi] = dval;
}
+ coord_vec = instance->get_y();
for (int yi = 0; yi < nyv; yi++) {
input.read(reinterpret_cast<char *>(&dval), sizeof(double));
- instance->set_y(yi, dval);
+ coord_vec[yi] = dval;
}
+ coord_vec = instance->get_z();
for (int zi = 0; zi < nzv; zi++) {
input.read(reinterpret_cast<char *>(&dval), sizeof(double));
- instance->set_z(zi, dval);
+ coord_vec[zi] = dval;
}
int _nlmmt = 2 * lm * (lm + 2);
int _nrvc = nxv * nyv * nzv;
@@ -709,7 +728,7 @@ TFRFME* TFRFME::from_legacy(string file_name) {
input.read(reinterpret_cast<char *>(&rval), sizeof(double));
input.read(reinterpret_cast<char *>(&ival), sizeof(double));
complex<double> value(rval, ival);
- instance->set_matrix_element(wi, wj, value);
+ _wsum[wi][wj] = value;
} // wj loop
} // wi loop
input.close();
@@ -719,10 +738,6 @@ TFRFME* TFRFME::from_legacy(string file_name) {
return instance;
}
-std::complex<double> TFRFME::get_matrix_element(int row, int col) {
- return wsum[row][col];
-}
-
long TFRFME::get_memory_requirement(
int _lmode, int _lm, int _nkv, int _nxv,
int _nyv, int _nzv
@@ -758,10 +773,6 @@ double TFRFME::get_param(string param_name) {
return value;
}
-void TFRFME::set_matrix_element(int row, int col, complex<double> value) {
- wsum[row][col] = value;
-}
-
void TFRFME::set_param(string param_name, double value) {
if (param_name.compare("vk") == 0) vk = value;
else if (param_name.compare("exri") == 0) exri = value;
@@ -814,7 +825,7 @@ void TFRFME::write_hdf5(string file_name) {
rec_ptr_list.append(&nzv);
rec_name_list.append("VK");
rec_type_list.append("FLOAT64_(1)");
- rec_ptr_list.append(&lm);
+ rec_ptr_list.append(&vk);
rec_name_list.append("EXRI");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&exri);
diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp
index b87ea03d..c9a490c4 100644
--- a/src/libnptm/tra_subs.cpp
+++ b/src/libnptm/tra_subs.cpp
@@ -244,7 +244,7 @@ void ffrt(
delete[] ctqcs;
}
-void frfmer(
+complex<double> *frfmer(
int nkv, double vkm, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
Swap1 *tt1, Swap2 *tt2
@@ -252,39 +252,32 @@ void frfmer(
const int nlemt = le * (le + 2) * 2;
const complex<double> cc0(0.0, 0.0);
double sq = vkm * vkm;
+ double *_vkv = tt2->get_vector();
+ double **_vkzm = tt2->get_matrix();
+ complex<double> *wk = new complex<double>[nlemt]();
for (int jy90 = 0; jy90 < nkv; jy90++) {
- double vky = tt2->get_vector_element(jy90);
+ double vky = _vkv[jy90];
double sqy = vky * vky;
for (int jx80 = 0; jx80 < nkv; jx80++) {
- double vkx = tt2->get_vector_element(jx80);
+ double vkx = _vkv[jx80];
double sqx = vkx * vkx;
double sqn = sqx + sqy;
double vkn = sqrt(sqn);
if (vkn <= vknmx) {
double vkz = sqrt(sq - sqn);
- wamff(tt1, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf);
- /*for (int j = 0; j < nlemt; j++) {
- double vreal = wk[j].real();
- double vimag = wk[j].imag();
- tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- }*/
+ wamff(wk, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf);
+ for (int j = 0; j < nlemt; j++) tt1->append(wk[j]);
//tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
- tt2->set_matrix_element(jx80, jy90, vkz);
+ _vkzm[jx80][jy90] = vkz;
} else { // label 50
- for (int j = 0; j < nlemt; j++) {
- //double vreal = 0.0;
- //double vimag = 0.0;
- //tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- //tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- tt1->set_element(j, cc0);
- }
+ for (int j = 0; j < nlemt; j++) tt1->append(cc0);
//double vkz = 0.0;
//tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
- tt2->set_matrix_element(jx80, jy90, 0.0);
+ _vkzm[jx80][jy90] = 0.0;
}
} // jx80 loop
} // jy90 loop
+ return wk;
}
void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw) {
@@ -369,7 +362,7 @@ void sampoa(
}
void wamff(
- Swap1 *tt1, double x, double y, double z, int lm, double apfafa,
+ complex<double> *wk, double x, double y, double z, int lm, double apfafa,
double tra, double spd, double rir, double ftcn, int lmode, double pmf
) {
const int nlmm = lm * (lm + 2);
@@ -389,27 +382,29 @@ void wamff(
bool onx = (y == 0.0);
bool ony = (x == 0.0);
bool onz = (onx && ony);
- if (!(onz && onx && ony)) {
- rho = sqrt(x * x + y * y);
- cph = x / rho;
- sph = y / rho;
- } else {
- if (onz) { // label 10
- cph = 1.0;
- sph = 0.0;
- } else {
- if (onx) { // label 12
- rho = sqrt(x * x);
- cph = (x < 0.0)? -1.0 : 1.0;
- sph = 0.0;
- } else {
- if (ony) { // label 13
- rho = sqrt(y * y);
- cph = 0.0;
- sph = (y < 0.0)? -1.0: 1.0;
- }
+ if (!onz) {
+ if (!onx) {
+ if (!ony) {
+ rho = sqrt(x * x + y * y);
+ cph = x / rho;
+ sph = y / rho;
+ // goes to 15
+ } else { // label 13
+ rho = (y < 0.0) ? -y : y;
+ cph = 0.0;
+ sph = (y < 0.0) ? -1.0 : 1.0;
+ // goes to 15
}
+ } else { // label 12
+ rho = (x < 0.0) ? -x : x;
+ cph = (x < 0.0) ? -1.0 : 1.0;
+ sph = 0.0;
+ // goes to 15
}
+ } else { // label 10
+ cph = 1.0;
+ sph = 0.0;
+ // goes to 15
}
// label 15
if (z == 0.0) {
@@ -464,13 +459,13 @@ void wamff(
} else if (lmode == 2) { // label 50
ftc1 = 2.0 * cthl / (cthl * rir + cth);
cfam *= (sth * pmf * ftc1);
- for (int i52 = 0; i52 < nlmmt; i52++) tt1->set_element(i52, w[i52][0] * cfam);
+ for (int i52 = 0; i52 < nlmmt; i52++) wk[i52] = w[i52][0] * cfam;
// returns
skip62 = true;
} else if (lmode == 3) { // label 53
ftc2 = 2.0 * cthl / (cthl + cth * rir);
cfam *= (sth * pmf * ftc2);
- for (int i55 = 0; i55 < nlmmt; i55++) tt1->set_element(i55, w[i55][1] * cfam);
+ for (int i55 = 0; i55 < nlmmt; i55++) wk[i55] = w[i55][1] * cfam;
// returns
skip62 = true;
}
@@ -485,7 +480,7 @@ void wamff(
if (lmode == 0) {
double f1 = cph * fam;
double f2 = -sph * fam;
- for (int i58 = 0; i58 < nlmmt; i58++) tt1->set_element(i58, w[i58][0] * f1 + w[i58][1] * f2);
+ for (int i58 = 0; i58 < nlmmt; i58++) wk[i58] = w[i58][0] * f1 + w[i58][1] * f2;
// returns
skip62 = true;
} else if (lmode == 1) { // label 60
@@ -496,19 +491,19 @@ void wamff(
skip62 = false;
} else if (lmode == 2) { // label 65
fam *= (pmf * sth);
- for (int i67 = 0; i67 < nlmmt; i67++) tt1->set_element(i67, w[i67][0] * fam);
+ for (int i67 = 0; i67 < nlmmt; i67++) wk[i67] = w[i67][0] * fam;
// returns
skip62 = true;
} else if (lmode == 3) { // label 68
fam *= (pmf * sth);
- for (int i70 = 0; i70 < nlmmt; i70++) tt1->set_element(i70, w[i70][1] * fam);
+ for (int i70 = 0; i70 < nlmmt; i70++) wk[i70] = w[i70][1] * fam;
// returns
skip62 = true;
}
}
if (!skip62) {
if (lmode == 0 || lmode == 1) { // label 62
- for (int i63 = 0; i63 < nlmmt; i63++) tt1->set_element(i63, w[i63][0] * cf1 + w[i63][1] * cf2);
+ for (int i63 = 0; i63 < nlmmt; i63++) wk[i63] = w[i63][0] * cf1 + w[i63][1] * cf2;
}
}
}
diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
index 79c68e35..d96a4fad 100644
--- a/src/trapping/cfrfme.cpp
+++ b/src/trapping/cfrfme.cpp
@@ -48,6 +48,7 @@ void frfme(string data_file, string output_path) {
char namef[7];
char more;
complex<double> **w = NULL;
+ complex<double> *wk = NULL;
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
int line_count = 0, last_read_line = 0;
@@ -228,13 +229,15 @@ void frfme(string data_file, string output_path) {
long swap1_size, swap2_size, tfrfme_size;
double size_mb;
printf("INFO: calculation memory requirements\n");
- swap1_size = Swap1::get_memory_requirement(lm);
+ swap1_size = Swap1::get_memory_requirement(lm, nkv);
size_mb = 1.0 * swap1_size / 1024.0 / 1024.0;
printf("Swap 1: %.2lg MB\n", size_mb);
swap2_size = Swap2::get_memory_requirement(nkv);
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();
// End of array initialization
double vkm = vk * exri;
vknmx = vk * an;
@@ -258,23 +261,27 @@ void frfme(string data_file, string output_path) {
size_mb = 1.0 * (swap1_size + swap2_size + tfrfme_size) / 1024.0 / 1024.0;
printf("TOTAL: %.2lg MB\n", size_mb);
tfrfme = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
+ double *_xv = tfrfme->get_x();
+ double *_yv = tfrfme->get_y();
+ double *_zv = tfrfme->get_z();
+ complex<double> **_wsum = tfrfme->get_matrix();
for (int i24 = nxshpo; i24 <= nxs; i24++) {
- tfrfme->set_x(i24, tfrfme->get_x(i24 - 1) + delxyz);
- tfrfme->set_x(nxv - i24 - 1, -tfrfme->get_x(i24));
+ _xv[i24] = _xv[i24 - 1] + delxyz;
+ _xv[nxv - i24 - 1] = -_xv[i24];
} // i24 loop
for (int i25 = nyshpo; i25 <= nys; i25++) {
- tfrfme->set_y(i25, tfrfme->get_y(i25 - 1) + delxyz);
- tfrfme->set_y(nyv - i25 - 1, -tfrfme->get_y(i25));
+ _yv[i25] = _yv[i25 - 1] + delxyz;
+ _yv[nyv - i25 - 1] = -_yv[i25];
} // i25 loop
for (int i27 = nzshpo; i27 <= nzs; i27++) {
- tfrfme->set_z(i27, tfrfme->get_z(i27 - 1) + delxyz);
- tfrfme->set_z(nzv - i27 - 1, -tfrfme->get_z(i27));
+ _zv[i27] = _zv[i27 - 1] + delxyz;
+ _zv[nzv - i27 - 1] = -_zv[i27];
} // i27 loop
int nrvc = nxv * nyv * nzv;
int nkshpo = nksh + 1;
for (int i28 = nkshpo; i28 <= nks; i28++) {
- tt2->set_vector_element(i28, tt2->get_vector_element(i28 - 1) + delk);
- tt2->set_vector_element(nkv - i28 - 1, -tt2->get_vector_element(i28));
+ _vkv[i28] = _vkv[i28 - 1] + delk;
+ _vkv[nkv - i28 - 1] = -_vkv[i28];
} // i28 loop
if (tfrfme != NULL) {
tfrfme->set_param("vk", vk);
@@ -289,13 +296,14 @@ void frfme(string data_file, string output_path) {
string temp_name1 = output_path + "/c_TEMPTAPE1.hd5";
string temp_name2 = output_path + "/c_TEMPTAPE2.hd5";
//temptape1.open(temp_name1.c_str(), ios::out | ios::binary);
- tt1 = new Swap1(lm);
+ 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));*/
- frfmer(nkv, vkm, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, tt1, tt2);
- //temptape1.close();
+ 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");
+ //delete tt1;
tt2->set_param("apfafa", apfafa);
tt2->set_param("pmf", pmf);
tt2->set_param("spd", spd);
@@ -325,6 +333,8 @@ void frfme(string data_file, string output_path) {
} // jy40 loop
temptape2.close();*/
for (int j80 = jlmf - 1; j80 < jlml; j80++) {
+ complex<double> *tt1_wk = tt1->get_vector();
+ int wk_index = 0;
// w matrix
if (w != NULL) {
for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
@@ -343,38 +353,39 @@ void frfme(string data_file, string output_path) {
temptape1.read(reinterpret_cast<char *>(&vimag), sizeof(double));
wk[i] = complex<double>(vreal, vimag);
}*/
- w[jx50][jy50] = tt1->get_element(j80);
+ for (int wi = 0; wi < nlmmt; wi++) wk[wi] = tt1_wk[wk_index++];
+ w[jx50][jy50] = wk[j80];
} // jx50
} // jy50 loop
//temptape1.close();
int ixyz = 0;
- for (int wj = 0; wj < nrvc; wj++) tfrfme->set_matrix_element(j80, wj, cc0);
+ for (int wj = 0; wj < nrvc; wj++) _wsum[j80][wj] = cc0;
for (int iz75 = 0; iz75 < nzv; iz75++) {
- double z = tfrfme->get_z(iz75) + frsh;
+ double z = _zv[iz75] + frsh;
for (int iy70 = 0; iy70 < nyv; iy70++) {
- double y = tfrfme->get_y(iy70);
+ double y = _yv[iy70];
for (int ix65 = 0; ix65 < nxv; ix65++) {
- double x = tfrfme->get_x(ix65);
+ double x = _xv[ix65];
ixyz++;
complex<double> sumy = cc0;
for (int jy60 = 0; jy60 < nkv; jy60++) {
- double vky = tt2->get_vector_element(jy60);
- double vkx = tt2->get_vector_element(nkv - 1);
- double vkzf = tt2->get_matrix_element(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 = tt2->get_matrix_element(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 = tt2->get_vector_element(jx55 - 1);
- double vkz = tt2->get_matrix_element(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
- tfrfme->set_matrix_element(j80, ixyz - 1, sumy * delks);
+ _wsum[j80][ixyz - 1] = sumy * delks;
} // ix65 loop
} // iy70 loop
} // iz75 loop
@@ -425,7 +436,7 @@ void frfme(string data_file, string output_path) {
for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
delete[] w;
}
- //if (wk != NULL) delete[] wk;
+ if (wk != NULL) delete[] wk;
if (tt1 != NULL) delete tt1;
printf("FRFME: Done.\n");
}
diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp
index d44960b0..e2c0793e 100644
--- a/src/trapping/clffft.cpp
+++ b/src/trapping/clffft.cpp
@@ -48,6 +48,7 @@ void lffft(string data_file, string output_path) {
complex<double> **amd = NULL;
int **indam = NULL;
complex<double> *tmsm = NULL, *tmse = NULL, **tms = NULL;
+ complex<double> **_wsum = NULL;
int jft, jss, jtw;
int is, le, nvam = 0;
double vks, exris;
@@ -159,6 +160,10 @@ void lffft(string data_file, string output_path) {
nxv = (int)tfrfme->get_param("nxv");
nyv = (int)tfrfme->get_param("nyv");
nzv = (int)tfrfme->get_param("nzv");
+ _wsum = tfrfme->get_matrix();
+ double *_xv = tfrfme->get_x();
+ double *_yv = tfrfme->get_y();
+ double *_zv = tfrfme->get_z();
if (lm >= le) {
vk = tfrfme->get_param("vk");
exri = tfrfme->get_param("exri");
@@ -203,15 +208,15 @@ void lffft(string data_file, string output_path) {
tlfff.write(reinterpret_cast<char *>(&frsh), sizeof(double));
tlfff.write(reinterpret_cast<char *>(&exril), sizeof(double));
for (int i = 0; i < nxv; i++) {
- double x = tfrfme->get_x(i);
+ double x = _xv[i];
tlfff.write(reinterpret_cast<char *>(&x), sizeof(double));
}
for (int i = 0; i < nyv; i++) {
- double y = tfrfme->get_y(i);
+ double y = _yv[i];
tlfff.write(reinterpret_cast<char *>(&y), sizeof(double));
}
for (int i = 0; i < nzv; i++) {
- double z = tfrfme->get_z(i);
+ double z = _zv[i];
tlfff.write(reinterpret_cast<char *>(&z), sizeof(double));
}
if (jft < 0) goto160 = true;
@@ -241,15 +246,15 @@ void lffft(string data_file, string output_path) {
tlfft.write(reinterpret_cast<char *>(&frsh), sizeof(double));
tlfft.write(reinterpret_cast<char *>(&exril), sizeof(double));
for (int i = 0; i < nxv; i++) {
- double x = tfrfme->get_x(i);
+ double x = _xv[i];
tlfft.write(reinterpret_cast<char *>(&x), sizeof(double));
}
for (int i = 0; i < nyv; i++) {
- double y = tfrfme->get_y(i);
+ double y = _yv[i];
tlfft.write(reinterpret_cast<char *>(&y), sizeof(double));
}
for (int i = 0; i < nzv; i++) {
- double z = tfrfme->get_z(i);
+ double z = _zv[i];
tlfft.write(reinterpret_cast<char *>(&z), sizeof(double));
}
} else { // Should never happen.
@@ -277,7 +282,7 @@ void lffft(string data_file, string output_path) {
//binary_input.read(reinterpret_cast<char *>(&vimag), sizeof(double));
int row = i;
int col = (nzv * iz475) + (nyv * iy475) + ix475;
- complex<double> value = tfrfme->get_matrix_element(row, col);
+ complex<double> value = _wsum[row][col];
if (lm <= le) {
//ws[i] = complex<double>(vreal, vimag);
ws[i] = value;
--
GitLab