Commit 76fbd1bb authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Fix binary I/O to TFRFME file

parent cfd0fa52

src/include/tfrfme.h

+31 −3
Original line number Diff line number Diff line
@@ -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;
@@ -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.

src/libnptm/tfrfme.cpp

+35 −27
Original line number Diff line number Diff line
@@ -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++) {

    int index = 0;

    for (int wj = 0; wj < _nrvc; wj++) {

	int index = (2 * _nrvc) * wi + 2 * wj;

      for (int wi = 0; wi < _nlmmt; wi++) {

	value = complex<double>(elements[index], elements[index + 1]);

	_wsum[wi][wj] = value;

      } // wj 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 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

    } // 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++) {

  int index = 0;

  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

    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 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

    } // wj loop

    output.close();

  } else { // Should never occur

    printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());

src/trapping/cfrfme.cpp

+18 −70
Original line number Diff line number Diff line
@@ -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";

src/trapping/clffft.cpp

+1 −1
Original line number Diff line number Diff line
@@ -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);