Commit 67e9e04a authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Implement static functions to write TM data for cluster

parent d614f766
Loading
Loading
Loading
Loading
+69 −3
Original line number Diff line number Diff line
@@ -48,19 +48,58 @@ class TransitionMatrix {

  /*! \brief Write the Transition Matrix to HDF5 binary output.
   *
   * This function takes care of the specific task of building a transition
   * matrix memory data structure from a binary input file formatted according
   * to the structure used by the original FORTRAN code.
   * This function takes care of the specific task of writing the transition
   * matrix memory data structure to a binary output file formatted according
   * to the HDF5 standard.
   *
   * \param file_name: `string` Name of the binary configuration data file.
   */
  void write_hdf5(std::string file_name);
  
  /*! \brief Write transition matrix data to HDF5 binary output.
   *
   * This function takes care of the specific task of writing the transition
   * matrix memory data structure to a binary output file formatted according
   * to the HDF5 standard. It is designed to work for the case of clusters of
   * spheres.
   *
   * \param file_name: `string` Name of the binary configuration data file.
   * \param _nlemt: `np_int` Size of the matrix (2 * LE * (LE + 2)).
   * \param _lm: `int` Maximum field expansion order.
   * \param _vk: `double` Wave number in scale units.
   * \param _exri: `double` External medium refractive index.
   * \param _am0m: `complex double **`
   */
  static void write_hdf5(
			 std::string file_name, np_int _nlemt, int _lm, double _vk,
			 double _exri, dcomplex **_am0m
  );
  
  /*! \brief Write the Transition Matrix to legacy binary output.
   *
   * \param file_name: `string` Name of the binary configuration data file.
   */
  void write_legacy(std::string file_name);

  /*! \brief Write transition matrix data to HDF5 binary output.
   *
   * This function takes care of the specific task of writing the transition
   * matrix memory data structure to a binary output file formatted according
   * to the format used by the legacy FORTRAN code. It is designed to work for
   * the case of clusters of spheres.
   *
   * \param file_name: `string` Name of the binary configuration data file.
   * \param _nlemt: `np_int` Size of the matrix (2 * LE * (LE + 2)).
   * \param _lm: `int` Maximum field expansion order.
   * \param _vk: `double` Wave number in scale units.
   * \param _exri: `double` External medium refractive index.
   * \param _am0m: `complex double **`
   */
  static void write_legacy(
			   std::string file_name, np_int _nlemt, int _lm, double _vk,
			   double _exri, dcomplex **_am0m
  );
  
 public:
  /*! \brief Default Transition Matrix instance constructor.
   *
@@ -142,6 +181,33 @@ class TransitionMatrix {
   */
  void write_binary(std::string file_name, std::string mode="LEGACY");
  
  /*! \brief Write a Transition Matrix to a binary file without instanciating it.
   *
   * Transition Matrix data can take a large amount of memory. For such reason, attempts
   * to create TransitionMatrix instances only for writing purposes can create
   * unnecessary resource consumption and computing time to duplicate the data into
   * the output buffer. This function offers output to file as a static method. It
   * takes the arguments of a constructor together with the usual arguments to specify
   * the output file name and format, to write the required data directly to a file,
   * without creating a new TransitionMatrix instance. The implementation works for
   * TransitionMatrix objects built for the CLUSTER case. It belongs to the public class
   * interface and it calls the proper versions of `write_legacy()` and `write_hdf5()`,
   * depending on the requested output format.
   * 
   * \param file_name: `string` Name of the file to be written.
   * \param _nlemt: `np_int` Size of the matrix (2 * LE * (LE + 2)).
   * \param _lm: `int` Maximum field expansion order.
   * \param _vk: `double` Wave number in scale units.
   * \param _exri: `double` External medium refractive index.
   * \param _am0m: `complex double **`
   * \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"] . Optional
   * (default is "LEGACY").
   */
  static void write_binary(
			   std::string file_name, np_int _nlemt, int _lm, double _vk,
			   double _exri, dcomplex **_am0m, std::string mode="LEGACY"
  );
  
  /*! \brief Test whether two instances of TransitionMatrix are equal.
   *
   * Transition matrices can be the result of a calculation or of a file input operation,
+91 −11
Original line number Diff line number Diff line
@@ -218,6 +218,20 @@ void TransitionMatrix::write_binary(string file_name, string mode) {
  }
}

void TransitionMatrix::write_binary(
				   std::string file_name, np_int _nlemt, int _lm, double _vk,
				   double _exri, dcomplex **_am0m, std::string mode
) {
  if (mode.compare("LEGACY") == 0) {
    write_legacy(file_name, _nlemt, _lm, _vk, _exri, _am0m);
  } else if (mode.compare("HDF5") == 0) {
    write_hdf5(file_name, _nlemt, _lm, _vk, _exri, _am0m);
  } else {
    string message = "Unknown format mode: \"" + mode + "\"";
    throw UnrecognizedFormatException(message);
  }
}

void TransitionMatrix::write_hdf5(string file_name) {
  if (is == 1 || is == 1111) {
    List<string> rec_name_list(1);
@@ -237,17 +251,10 @@ void TransitionMatrix::write_hdf5(string file_name) {
    rec_type_list.append("FLOAT64_(1)");
    rec_ptr_list.append(&exri);
    rec_name_list.append("ELEMENTS");
    str_type = "FLOAT64_(" + to_string(shape[0]) + "," + to_string(2 * shape[1]) + ")";
    str_type = "COMPLEX128_(" + to_string(shape[0]) + "," + to_string(shape[1]) + ")";
    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 * shape[0] * shape[1];
    double *ptr_elements = new double[num_elements]();
    for (int ei = 0; ei < num_elements / 2; ei++) {
      ptr_elements[2 * ei] = real(elements[ei]);
      ptr_elements[2 * ei + 1] = imag(elements[ei]);
    }
    rec_ptr_list.append(ptr_elements);
    dcomplex *p_first = elements;
    rec_ptr_list.append(p_first);
    if (is == 1111) {
      rec_name_list.append("RADIUS");
      rec_type_list.append("FLOAT64_(1)");
@@ -264,7 +271,7 @@ void TransitionMatrix::write_hdf5(string file_name) {
      hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
    hdf_file->close();
    
    delete[] ptr_elements;
    p_first = NULL;
    delete[] rec_names;
    delete[] rec_types;
    delete[] rec_pointers;
@@ -275,6 +282,52 @@ void TransitionMatrix::write_hdf5(string file_name) {
  }
}

void TransitionMatrix::write_hdf5(
				  std::string file_name, np_int _nlemt, int _lm, double _vk,
				  double _exri, dcomplex **_am0m
) {
  int is = 1;
  List<string> rec_name_list(1);
  List<string> rec_type_list(1);
  List<void *> rec_ptr_list(1);
  string str_type, str_name;
  rec_name_list.set(0, "IS");
  rec_type_list.set(0, "INT32_(1)");
  rec_ptr_list.set(0, &is);
  rec_name_list.append("L_MAX");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&_lm);
  rec_name_list.append("VK");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&_vk);
  rec_name_list.append("EXRI");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&_exri);
  rec_name_list.append("ELEMENTS");
  str_type = "COMPLEX128_(" + to_string(_nlemt) + "," + to_string(_nlemt) + ")";
  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)
  dcomplex *p_first = _am0m[0];
  rec_ptr_list.append(p_first);
  
  string *rec_names = rec_name_list.to_array();
  string *rec_types = rec_type_list.to_array();
  void **rec_pointers = rec_ptr_list.to_array();
  const int rec_num = rec_name_list.length();
  FileSchema schema(rec_num, rec_types, rec_names);
  HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC);
  for (int ri = 0; ri < rec_num; ri++)
    hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
  hdf_file->close();

  p_first = NULL;
  delete[] rec_names;
  delete[] rec_types;
  delete[] rec_pointers;
  delete hdf_file;
}

void TransitionMatrix::write_legacy(string file_name) {
  fstream ttms;
  if (is == 1111 || is == 1) {
@@ -308,6 +361,33 @@ void TransitionMatrix::write_legacy(string file_name) {
  }
}

void TransitionMatrix::write_legacy(
				    std::string file_name, np_int _nlemt, int _lm, double _vk,
				    double _exri, dcomplex **_am0m
) {
  fstream ttms;
  int is = 1;
  ttms.open(file_name, ios::binary | ios::out);
  if (ttms.is_open()) {
    ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
    ttms.write(reinterpret_cast<char *>(&_lm), sizeof(int));
    ttms.write(reinterpret_cast<char *>(&_vk), sizeof(double));
    ttms.write(reinterpret_cast<char *>(&_exri), sizeof(double));
    double rval, ival;
    for (np_int ei = 0; ei < _nlemt; ei++) {
      for (np_int ej = 0; ej < _nlemt; ej++) {
	rval = real(_am0m[ei][ej]);
	ival = imag(_am0m[ei][ej]);
	ttms.write(reinterpret_cast<char *>(&rval), sizeof(double));
	ttms.write(reinterpret_cast<char *>(&ival), sizeof(double));
      }
    }
    ttms.close();
  } else { // Failed to open output file. Should never happen.
    printf("ERROR: could not open Transition Matrix file for writing.\n");
  }
}

bool TransitionMatrix::operator ==(TransitionMatrix &other) {
  if (is != other.is) {
    return false;