Skip to content
Commits on Source (5)
Hide whitespace changes
Inline Side-by-side
src/cluster/cluster.cpp
View file @ f88ad025
......@@ -917,10 +917,9 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
if (jxi488 == jwtm) {
int nlemt = 2 * c4->nlem;
string ttms_name = output_path + "/c_TTMS.hd5";
TransitionMatrix ttms(nlemt, lm, vk, exri, c1ao->am0m);
ttms.write_binary(ttms_name, "HDF5");
TransitionMatrix::write_binary(ttms_name, nlemt, lm, vk, exri, c1ao->am0m, "HDF5");
ttms_name = output_path + "/c_TTMS";
ttms.write_binary(ttms_name);
TransitionMatrix::write_binary(ttms_name, nlemt, lm, vk, exri, c1ao->am0m);
}
}
// label 156: continue from here
......
src/include/TransitionMatrix.h
View file @ f88ad025
......@@ -48,19 +48,98 @@ 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 without a pre-existing instance. It is designed to
* work for the case of a cluster 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 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 without a pre-existing instance. It is designed to
* work for the case of a single sphere.
*
* \param file_name: `string` Name of the binary configuration data file.
* \param _lm: `int` Maximum field expansion order.
* \param _vk: `double` Wave number in scale units.
* \param _exri: `double` External medium refractive index.
* \param _rmi: `complex double **`
* \param _rei: `complex double **`
* \param _sphere_radius: `double` Radius of the sphere.
*/
static void write_hdf5(
std::string file_name, int _lm, double _vk, double _exri,
dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
);
/*! \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 binary output using legacy format.
*
* 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
);
/*! \brief Write transition matrix data to binary output using legacy format.
*
* 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 original code structure without a pre-existing instance. It is designed
* to work for the case of a single sphere.
*
* \param file_name: `string` Name of the binary configuration data file.
* \param _lm: `int` Maximum field expansion order.
* \param _vk: `double` Wave number in scale units.
* \param _exri: `double` External medium refractive index.
* \param _rmi: `complex double **`
* \param _rei: `complex double **`
* \param _sphere_radius: `double` Radius of the sphere.
*/
static void write_legacy(
std::string file_name, int _lm, double _vk, double _exri,
dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
);
public:
/*! \brief Default Transition Matrix instance constructor.
*
......@@ -142,6 +221,62 @@ class TransitionMatrix {
*/
void write_binary(std::string file_name, std::string mode="LEGACY");
/*! \brief Write a cluster 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 Write a single sphere 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 single sphere 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 _lm: `int` Maximum field expansion order.
* \param _vk: `double` Wave number in scale units.
* \param _exri: `double` External medium refractive index.
* \param _rmi: `complex double **`
* \param _rei: `complex double **`
* \param _sphere_radius: `double` Radius of the sphere.
* \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"] . Optional
* (default is "LEGACY").
*/
static void write_binary(
std::string file_name, int _lm, double _vk, double _exri,
dcomplex **_rmi, dcomplex **_rei, double _sphere_radius,
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,
......
src/libnptm/TransitionMatrix.cpp
View file @ f88ad025
......@@ -94,7 +94,7 @@ TransitionMatrix::TransitionMatrix(
}
}
TransitionMatrix* TransitionMatrix::from_binary(string file_name, string mode) {
TransitionMatrix* TransitionMatrix::from_binary(std::string file_name, std::string mode) {
TransitionMatrix *tm = NULL;
if (mode.compare("LEGACY") == 0) {
tm = TransitionMatrix::from_legacy(file_name);
......@@ -107,7 +107,7 @@ TransitionMatrix* TransitionMatrix::from_binary(string file_name, string mode) {
return tm;
}
TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
TransitionMatrix* TransitionMatrix::from_hdf5(std::string file_name) {
TransitionMatrix *tm = NULL;
unsigned int flags = H5F_ACC_RDONLY;
HDFFile *hdf_file = new HDFFile(file_name, flags);
......@@ -160,7 +160,7 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
return tm;
}
TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
TransitionMatrix* TransitionMatrix::from_legacy(std::string file_name) {
fstream ttms;
TransitionMatrix *tm = NULL;
ttms.open(file_name, ios::binary | ios::in);
......@@ -207,7 +207,7 @@ TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
return tm;
}
void TransitionMatrix::write_binary(string file_name, string mode) {
void TransitionMatrix::write_binary(std::string file_name, std::string mode) {
if (mode.compare("LEGACY") == 0) {
write_legacy(file_name);
} else if (mode.compare("HDF5") == 0) {
......@@ -218,7 +218,36 @@ void TransitionMatrix::write_binary(string file_name, string mode) {
}
}
void TransitionMatrix::write_hdf5(string file_name) {
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_binary(
std::string file_name, int _lm, double _vk, double _exri,
dcomplex **_rmi, dcomplex **_rei, double _sphere_radius,
std::string mode
) {
if (mode.compare("LEGACY") == 0) {
write_legacy(file_name, _lm, _vk, _exri, _rmi, _rei, _sphere_radius);
} else if (mode.compare("HDF5") == 0) {
write_hdf5(file_name, _lm, _vk, _exri, _rmi, _rei, _sphere_radius);
} else {
string message = "Unknown format mode: \"" + mode + "\"";
throw UnrecognizedFormatException(message);
}
}
void TransitionMatrix::write_hdf5(std::string file_name) {
if (is == 1 || is == 1111) {
List<string> rec_name_list(1);
List<string> rec_type_list(1);
......@@ -237,17 +266,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 +286,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,7 +297,104 @@ void TransitionMatrix::write_hdf5(string file_name) {
}
}
void TransitionMatrix::write_legacy(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_hdf5(
std::string file_name, int _lm, double _vk, double _exri,
dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
) {
int is = 1111;
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);
dcomplex *_elements = new dcomplex[2 * _lm]();
for (int ei = 0; ei < _lm; ei++) {
_elements[2 * ei] = -1.0 / _rmi[ei][0];
_elements[2 * ei + 1] = -1.0 / _rei[ei][0];
}
rec_name_list.append("ELEMENTS");
str_type = "COMPLEX128_(" + to_string(_lm) + "," + to_string(2) + ")";
rec_type_list.append(str_type);
rec_ptr_list.append(_elements);
rec_name_list.append("RADIUS");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&_sphere_radius);
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();
delete[] _elements;
delete[] rec_names;
delete[] rec_types;
delete[] rec_pointers;
delete hdf_file;
}
void TransitionMatrix::write_legacy(std::string file_name) {
fstream ttms;
if (is == 1111 || is == 1) {
ttms.open(file_name, ios::binary | ios::out);
......@@ -308,6 +427,66 @@ 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");
}
}
void TransitionMatrix::write_legacy(
std::string file_name, int _lm, double _vk, double _exri,
dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
) {
fstream ttms;
int is = 1111;
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;
dcomplex element;
for (int ei = 0; ei < _lm; ei++) {
element = -1.0 / _rmi[ei][0];
rval = real(element);
ival = imag(element);
ttms.write(reinterpret_cast<char *>(&rval), sizeof(double));
ttms.write(reinterpret_cast<char *>(&ival), sizeof(double));
element = -1.0 / _rei[ei][0];
rval = real(element);
ival = imag(element);
ttms.write(reinterpret_cast<char *>(&rval), sizeof(double));
ttms.write(reinterpret_cast<char *>(&ival), sizeof(double));
}
ttms.write(reinterpret_cast<char *>(&_sphere_radius), 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;
......
src/libnptm/file_io.cpp
View file @ f88ad025
......@@ -81,7 +81,7 @@ HDFFile* HDFFile::from_schema(
herr_t status;
string *rec_types = schema.get_record_types();
string *rec_names = schema.get_record_names();
string known_types[] = {"INT32", "FLOAT64"};
string known_types[] = {"INT32", "FLOAT64", "COMPLEX128"};
int rec_num = schema.get_record_number();
regex re;
smatch m;
......@@ -99,8 +99,9 @@ HDFFile* HDFFile::from_schema(
str_target = m.suffix().str();
if (type_index == 1) data_type = H5Tcopy(H5T_NATIVE_INT);
else if (type_index == 2) data_type = H5Tcopy(H5T_NATIVE_DOUBLE);
else if (type_index == 3) data_type = H5Tcopy(H5T_NATIVE_DOUBLE);
}
if (type_index == 2) break;
if (type_index == 3) break;
}
if (found_type) {
re = regex("[0-9]+");
......@@ -119,6 +120,9 @@ HDFFile* HDFFile::from_schema(
max_dims[ti] = dim;
str_target = m.suffix().str();
}
int multiplier = (type_index == 3) ? 2 : 1;
dims[rank - 1] *= multiplier;
max_dims[rank - 1] *= multiplier;
hid_t dataspace_id = H5Screate_simple(rank, dims, max_dims);
hid_t dataset_id = H5Dcreate(
file_id, rec_names[ri].c_str(), data_type, dataspace_id, H5P_DEFAULT,
......@@ -183,7 +187,7 @@ herr_t HDFFile::write(
hid_t mem_space_id, hid_t file_space_id, hid_t dapl_id,
hid_t dxpl_id
) {
string known_types[] = {"INT32", "FLOAT64"};
string known_types[] = {"INT32", "FLOAT64", "COMPLEX128"};
regex re;
smatch m;
bool found_type = false;
......@@ -191,7 +195,7 @@ herr_t HDFFile::write(
while (!found_type) {
re = regex(known_types[type_index++]);
found_type = regex_search(data_type, m, re);
if (type_index == 2) break;
if (type_index == 3) break;
}
if (found_type) {
hid_t dataset_id = H5Dopen2(file_id, dataset_name.c_str(), dapl_id);
......@@ -201,6 +205,8 @@ herr_t HDFFile::write(
mem_type_id = H5T_NATIVE_INT; break;
case 2:
mem_type_id = H5T_NATIVE_DOUBLE; break;
case 3:
mem_type_id = H5T_NATIVE_DOUBLE; break;
default:
throw UnrecognizedParameterException("unrecognized data type \"" + data_type + "\"");
}
......
src/sphere/sphere.cpp
View file @ f88ad025
......@@ -369,11 +369,16 @@ void sphere(string config_file, string data_file, string output_path) {
} // i132
if (idfc >= 0 and nsph == 1 and jxi == jwtm) {
// This is the condition that writes the transition matrix to output.
TransitionMatrix ttms(l_max, vk, exri, c1->rmi, c1->rei, sconf->get_radius(0));
string ttms_name = output_path + "/c_TTMS.hd5";
ttms.write_binary(ttms_name, "HDF5");
TransitionMatrix::write_binary(
ttms_name, l_max, vk, exri, c1->rmi, c1->rei,
sconf->get_radius(0), "HDF5"
);
ttms_name = output_path + "/c_TTMS";
ttms.write_binary(ttms_name);
TransitionMatrix::write_binary(
ttms_name, l_max, vk, exri, c1->rmi, c1->rei,
sconf->get_radius(0)
);
}
double cs0 = 0.25 * vk * vk * vk / half_pi;
sscr0(tfsas, nsph, l_max, vk, exri, c1);
......