diff --git a/src/include/TransitionMatrix.h b/src/include/TransitionMatrix.h
index 3966abdab0b48463eaefb256dd6d128fe84e4e75..54903a95d50bf3eb1882ed8fd6ae9ec471d0b734 100644
--- a/src/include/TransitionMatrix.h
+++ b/src/include/TransitionMatrix.h
@@ -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,
diff --git a/src/libnptm/TransitionMatrix.cpp b/src/libnptm/TransitionMatrix.cpp
index fc6c92af6c77dca0aabce667de72f1e51ef15fec..7468ec8dad6e693b3df88561b19e7d147d05c753 100644
--- a/src/libnptm/TransitionMatrix.cpp
+++ b/src/libnptm/TransitionMatrix.cpp
@@ -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;