Something went wrong on our end
Select Git revision
algebraic.cpp
-
Giovanni La Mura authoredGiovanni La Mura authored
tfrfme.cpp 33.94 KiB
/* Copyright 2004 INAF - Osservatorio Astronomico di Cagliari
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*! \file tfrfme.cpp
*
* \brief Implementation of the trapping calculation objects.
*/
#include <exception>
#include <fstream>
#include <hdf5.h>
#include <string>
#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
#endif
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
#endif
#ifndef INCLUDE_TFRFME_H_
#include "../include/tfrfme.h"
#endif
#ifndef INCLUDE_FILE_IO_H_
#include "../include/file_io.h"
#endif
using namespace std;
// >>> START OF Swap1 CLASS IMPLEMENTATION <<<
Swap1::Swap1(int lm, int _nkv) {
nkv = _nkv;
nlmmt = 2 * lm * (lm + 2);
const int size = nkv * nkv * nlmmt;
wk = new dcomplex[size]();
last_index = 0;
}
Swap1* Swap1::from_binary(const std::string& file_name, const std::string& mode) {
Swap1 *instance = NULL;
if (mode.compare("LEGACY") == 0) {
instance = from_legacy(file_name);
} else if (mode.compare("HDF5") == 0) {
instance = from_hdf5(file_name);
} else {
string message = "Unknown format mode: \"" + mode + "\"";
throw UnrecognizedFormatException(message);
}
return instance;
}
Swap1* Swap1::from_hdf5(const std::string& file_name) { Swap1 *instance = NULL;
unsigned int flags = H5F_ACC_RDONLY;
HDFFile *hdf_file = new HDFFile(file_name, flags);
herr_t status = hdf_file->get_status();
double *elements;
string str_type;
int _nlmmt, _nkv, lm, num_elements, index;
dcomplex value;
dcomplex *_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 * _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 < num_elements / 2; wi++) {
index = 2 * wi;
value = elements[index] + elements[index + 1] * I;
_wk[wi] = value;
} // wi loop
delete[] elements;
status = hdf_file->close();
delete hdf_file;
}
return instance;
}
Swap1* Swap1::from_legacy(const std::string& file_name) {
fstream input;
Swap1 *instance = NULL;
dcomplex *_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);
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));
_wk[j] = rval + ival * I;
}
input.close();
} else {
printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
}
return instance;
}
long Swap1::get_memory_requirement(int lm, int _nkv) {
long size = (long)(3 * sizeof(int));
size += (long)(sizeof(dcomplex) * 2 * lm * (lm + 2) * _nkv * _nkv);
return size;
}
void Swap1::write_binary(const std::string& file_name, const std::string& mode) {
if (mode.compare("LEGACY") == 0) {
write_legacy(file_name);
} else if (mode.compare("HDF5") == 0) {
write_hdf5(file_name);
} else {
string message = "Unknown format mode: \"" + mode + "\""; throw UnrecognizedFormatException(message);
}
}
void Swap1::write_hdf5(const std::string& file_name) {
List<string> rec_name_list(1);
List<string> rec_type_list(1);
List<void *> rec_ptr_list(1);
herr_t status;
string str_type;
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 < num_elements / 2; wi++) {
ptr_elements[2 * wi] = real(wk[wi]);
ptr_elements[2 * wi + 1] = imag(wk[wi]);
}
rec_ptr_list.append(ptr_elements);
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++)
status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
status = hdf_file->close();
delete[] ptr_elements;
delete[] rec_names;
delete[] rec_types;
delete[] rec_pointers;
delete hdf_file;
}
void Swap1::write_legacy(const std::string& file_name) {
fstream output;
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));
output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
for (int j = 0; j < num_elements; j++) {
rval = real(wk[j]);
ival = imag(wk[j]);
output.write(reinterpret_cast<char *>(&rval), sizeof(double));
output.write(reinterpret_cast<char *>(&ival), sizeof(double));
}
output.close();
} else { // Should never occur
printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
}
}
bool Swap1::operator ==(Swap1 &other) {
if (nlmmt != other.nlmmt) {
return false;
}
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;
}
}
return true;
}
// >>> END OF Swap1 CLASS IMPLEMENTATION <<<
// >>> START OF Swap2 CLASS IMPLEMENTATION <<<
Swap2::Swap2(int _nkv) {
nkv = _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() {
delete[] vkv;
for (int vi = nkv - 1; vi > -1; vi--) delete[] vkzm[vi];
delete[] vkzm;
}
Swap2* Swap2::from_binary(const std::string& file_name, const std::string& mode) {
Swap2 *instance = NULL;
if (mode.compare("LEGACY") == 0) {
instance = from_legacy(file_name);
} else if (mode.compare("HDF5") == 0) {
instance = from_hdf5(file_name);
} else {
string message = "Unknown format mode: \"" + mode + "\"";
throw UnrecognizedFormatException(message);
}
return instance;
}
Swap2* Swap2::from_hdf5(const std::string& file_name) {
Swap2 *instance = NULL;
unsigned int flags = H5F_ACC_RDONLY;
HDFFile *hdf_file = new HDFFile(file_name, flags);
herr_t status = hdf_file->get_status();
string str_type;
int _nkv, _nlmmt, _nrvc;
double value;
if (status == 0) {
status = hdf_file->read("NKV", "INT32", &_nkv);
instance = new Swap2(_nkv);
str_type = "FLOAT64_(" + to_string(_nkv) + ")";
status = hdf_file->read("VKV", str_type, instance->vkv);
str_type = "FLOAT64_(" + to_string(_nkv) + "," + to_string(_nkv) + ")";
status = hdf_file->read("VKZM", str_type, instance->vkzm);
status = hdf_file->read("APFAFA", "FLOAT64", &value);
instance->set_param("apfafa", value);
status = hdf_file->read("PMF", "FLOAT64", &value);
instance->set_param("pmf", value);
status = hdf_file->read("SPD", "FLOAT64", &value);
instance->set_param("spd", value);
status = hdf_file->read("RIR", "FLOAT64", &value);
instance->set_param("rir", value);
status = hdf_file->read("FTCN", "FLOAT64", &value);
instance->set_param("ftcn", value);
status = hdf_file->read("FSHMX", "FLOAT64", &value);
instance->set_param("fshmx", value);
status = hdf_file->read("VXYZMX", "FLOAT64", &value);
instance->set_param("vxyzmx", value);
status = hdf_file->read("DELXYZ", "FLOAT64", &value); instance->set_param("delxyz", value);
status = hdf_file->read("VKNMX", "FLOAT64", &value);
instance->set_param("vknmx", value);
status = hdf_file->read("DELK", "FLOAT64", &value);
instance->set_param("delk", value);
status = hdf_file->read("DELKS", "FLOAT64", &value);
instance->set_param("delks", value);
status = hdf_file->read("NLMMT", "INT32", &_nlmmt);
instance->set_param("nlmmt", 1.0 * _nlmmt);
status = hdf_file->read("NRVC", "INT32", &_nrvc);
instance->set_param("nrvc", 1.0 * _nrvc);
status = hdf_file->close();
delete hdf_file;
}
return instance;
}
Swap2* Swap2::from_legacy(const std::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));
_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));
_vkzm[mi][mj] = value;
}
}
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("apfafa", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("pmf", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("spd", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("rir", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("ftcn", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("fshmx", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("vxyzmx", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("delxyz", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("vknmx", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("delk", value);
input.read(reinterpret_cast<char *>(&value), sizeof(double));
instance->set_param("delks", value);
input.read(reinterpret_cast<char *>(&_nlmmt), sizeof(int));
instance->set_param("nlmmt", 1.0 * _nlmmt);
input.read(reinterpret_cast<char *>(&_nrvc), sizeof(int));
instance->set_param("nrvc", 1.0 * _nrvc);
input.close();
} else {
printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
} return instance;
}
long Swap2::get_memory_requirement(int _nkv) {
long size = (long)(5 * sizeof(int) + 11 * sizeof(double));
size += (long)(sizeof(double) * _nkv * (_nkv + 1));
return size;
}
double Swap2::get_param(const std::string& param_name) {
double value;
if (param_name.compare("nkv") == 0) value = 1.0 * nkv;
else if (param_name.compare("apfafa") == 0) value = apfafa;
else if (param_name.compare("pmf") == 0) value = pmf;
else if (param_name.compare("spd") == 0) value = spd;
else if (param_name.compare("rir") == 0) value = rir;
else if (param_name.compare("ftcn") == 0) value = ftcn;
else if (param_name.compare("fshmx") == 0) value = fshmx;
else if (param_name.compare("vxyzmx") == 0) value = vxyzmx;
else if (param_name.compare("delxyz") == 0) value = delxyz;
else if (param_name.compare("vknmx") == 0) value = vknmx;
else if (param_name.compare("delk") == 0) value = delk;
else if (param_name.compare("delks") == 0) value = delks;
else if (param_name.compare("nlmmt") == 0) value = 1.0 * nlmmt;
else if (param_name.compare("nrvc") == 0) value = 1.0 * nrvc;
else {
string message = "Unrecognized parameter name \"" + param_name + "\"";
throw UnrecognizedParameterException(message);
}
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(const std::string& param_name, double value) {
if (param_name.compare("nkv") == 0) nkv = (int)value;
else if (param_name.compare("apfafa") == 0) apfafa = value;
else if (param_name.compare("pmf") == 0) pmf = value;
else if (param_name.compare("spd") == 0) spd = value;
else if (param_name.compare("rir") == 0) rir = value;
else if (param_name.compare("ftcn") == 0) ftcn = value;
else if (param_name.compare("fshmx") == 0) fshmx = value;
else if (param_name.compare("vxyzmx") == 0) vxyzmx = value;
else if (param_name.compare("delxyz") == 0) delxyz = value;
else if (param_name.compare("vknmx") == 0) vknmx = value;
else if (param_name.compare("delk") == 0) delk = value;
else if (param_name.compare("delks") == 0) delks = value;
else if (param_name.compare("nlmmt") == 0) nlmmt = (int)value;
else if (param_name.compare("nrvc") == 0) nrvc = (int)value;
else {
string message = "Unrecognized parameter name \"" + param_name + "\"";
throw UnrecognizedParameterException(message);
}
}
void Swap2::write_binary(const std::string& file_name, const std::string& mode) {
if (mode.compare("LEGACY") == 0) {
write_legacy(file_name);
} else if (mode.compare("HDF5") == 0) {
write_hdf5(file_name);
} else {
string message = "Unknown format mode: \"" + mode + "\"";
throw UnrecognizedFormatException(message);
}
}
void Swap2::write_hdf5(const std::string& file_name) {
List<string> rec_name_list(1);
List<string> rec_type_list(1);
List<void *> rec_ptr_list(1);
herr_t status;
string str_type;
rec_name_list.set(0, "NKV");
rec_type_list.set(0, "INT32_(1)");
rec_ptr_list.set(0, &nkv);
rec_name_list.append("VKV");
str_type = "FLOAT64_(" + to_string(nkv) + ")";
rec_type_list.append(str_type);
rec_ptr_list.append(vkv);
rec_name_list.append("VKZM");
str_type = "FLOAT64_(" + to_string(nkv) + "," + to_string(nkv) + ")";
rec_type_list.append(str_type);
rec_ptr_list.append(vkzm);
rec_name_list.append("APFAFA");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&apfafa);
rec_name_list.append("PMF");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&pmf);
rec_name_list.append("SPD");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&spd);
rec_name_list.append("RIR");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&rir);
rec_name_list.append("FTCN");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&ftcn);
rec_name_list.append("FSHMX");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&fshmx);
rec_name_list.append("VXYZMX");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&vxyzmx);
rec_name_list.append("delxyz");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&delxyz);
rec_name_list.append("VKNMX");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&vknmx);
rec_name_list.append("DELK");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&delk);
rec_name_list.append("DELKS");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&delks);
rec_name_list.append("NLMMT");
rec_type_list.append("INT32_(1)");
rec_ptr_list.append(&nlmmt);
rec_name_list.append("NRVC");
rec_type_list.append("INT32_(1)");
rec_ptr_list.append(&nrvc);
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++)
status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
status = hdf_file->close();
delete[] rec_names;
delete[] rec_types; delete[] rec_pointers;
delete hdf_file;
}
void Swap2::write_legacy(const std::string& file_name) {
fstream output;
double value;
output.open(file_name.c_str(), ios::out | ios::binary);
if (output.is_open()) {
output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
for (int j = 0; j < nkv; j++){
value = vkv[j];
output.write(reinterpret_cast<const char*>(&value), sizeof(double));
}
for (int mi = 0; mi < nkv; mi++) {
for (int mj = 0; mj < nkv; mj++) {
value = vkzm[mi][mj];
output.write(reinterpret_cast<const char*>(&value), sizeof(double));
}
}
output.write(reinterpret_cast<const char*>(&apfafa), sizeof(double));
output.write(reinterpret_cast<const char*>(&pmf), sizeof(double));
output.write(reinterpret_cast<const char*>(&spd), sizeof(double));
output.write(reinterpret_cast<const char*>(&rir), sizeof(double));
output.write(reinterpret_cast<const char*>(&ftcn), sizeof(double));
output.write(reinterpret_cast<const char*>(&fshmx), sizeof(double));
output.write(reinterpret_cast<const char*>(&vxyzmx), sizeof(double));
output.write(reinterpret_cast<const char*>(&delxyz), sizeof(double));
output.write(reinterpret_cast<const char*>(&vknmx), sizeof(double));
output.write(reinterpret_cast<const char*>(&delk), sizeof(double));
output.write(reinterpret_cast<const char*>(&delks), sizeof(double));
output.write(reinterpret_cast<const char*>(&nlmmt), sizeof(int));
output.write(reinterpret_cast<const char*>(&nrvc), sizeof(int));
output.close();
} else { // Should never occur
printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
}
}
bool Swap2::operator ==(Swap2 &other) {
if (nlmmt != other.nlmmt) {
return false;
}
if (nrvc != other.nrvc) {
return false;
}
if (nkv != other.nkv) {
return false;
}
if (apfafa != other.apfafa) {
return false;
}
if (pmf != other.pmf) {
return false;
}
if (spd != other.spd) {
return false;
}
if (rir != other.rir) {
return false;
}
if (ftcn != other.ftcn) {
return false;
}
if (fshmx != other.fshmx) {
return false;
}
if (vxyzmx != other.vxyzmx) {
return false;
} if (delxyz != other.delxyz) {
return false;
}
if (vknmx != other.vknmx) {
return false;
}
if (delk != other.delk) {
return false;
}
if (delks != other.delks) {
return false;
}
for (int vi = 0; vi < nkv; vi++) {
if (vkv[vi] != other.vkv[vi]) {
return false;
}
}
for (int mi = 0; mi < nkv; mi++) {
for (int mj = 0; mj < nkv; mj++) {
if (vkzm[mi][mj] != other.vkzm[mi][mj]) {
return false;
}
}
}
return true;
}
// >>> END OF Swap2 CLASS IMPLEMENTATION <<<
// >>> START OF TFRFME CLASS IMPLEMENTATION <<<
TFRFME::TFRFME(int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv) {
lmode = _lmode;
lm = _lm;
nkv = _nkv;
nxv = _nxv;
nyv = _nyv;
nzv = _nzv;
vk = 0.0;
exri = 0.0;
an = 0.0;
ff = 0.0;
tra = 0.0;
spd = 0.0;
frsh = 0.0;
exril = 0.0;
// Array initialization
xv = new double[nxv]();
yv = new double[nyv]();
zv = new double[nzv]();
nlmmt = lm * (lm + 2) * 2;
nrvc = nxv * nyv * nzv;
wsum = new dcomplex*[nlmmt];
for (int wi = 0; wi < nlmmt; wi++) wsum[wi] = new dcomplex[nrvc]();
}
TFRFME::~TFRFME() {
delete[] xv;
delete[] yv;
delete[] zv;
for (int wi = nlmmt - 1; wi > -1; wi--) delete[] wsum[wi];
delete[] wsum;
}
TFRFME* TFRFME::from_binary(const std::string& file_name, const std::string& mode) {
TFRFME *instance = NULL;
if (mode.compare("LEGACY") == 0) {
instance = from_legacy(file_name);
} else if (mode.compare("HDF5") == 0) {
instance = from_hdf5(file_name);
} else { string message = "Unknown format mode: \"" + mode + "\"";
throw UnrecognizedFormatException(message);
}
return instance;
}
TFRFME* TFRFME::from_hdf5(const std::string& file_name) {
TFRFME *instance = NULL;
unsigned int flags = H5F_ACC_RDONLY;
HDFFile *hdf_file = new HDFFile(file_name, flags);
herr_t status = hdf_file->get_status();
double *elements;
string str_type;
int _nlmmt, _nrvc, num_elements;
dcomplex value;
dcomplex **_wsum = NULL;
if (status == 0) {
int lmode, lm, nkv, nxv, nyv, nzv;
double vk, exri, an, ff, tra, spd, frsh, exril;
status = hdf_file->read("LMODE", "INT32", &lmode);
status = hdf_file->read("LM", "INT32", &lm);
status = hdf_file->read("NKV", "INT32", &nkv);
status = hdf_file->read("NXV", "INT32", &nxv);
status = hdf_file->read("NYV", "INT32", &nyv);
status = hdf_file->read("NZV", "INT32", &nzv);
status = hdf_file->read("VK", "FLOAT64", &vk);
status = hdf_file->read("EXRI", "FLOAT64", &exri);
status = hdf_file->read("AN", "FLOAT64", &an);
status = hdf_file->read("FF", "FLOAT64", &ff);
status = hdf_file->read("TRA", "FLOAT64", &tra);
status = hdf_file->read("SPD", "FLOAT64", &spd);
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);
instance->set_param("ff", ff);
instance->set_param("tra", tra);
instance->set_param("spd", spd);
instance->set_param("frsh", frsh);
instance->set_param("exril", exril);
str_type = "FLOAT64_(" + to_string(nxv) + ")";
status = hdf_file->read("XVEC", str_type, instance->xv);
str_type = "FLOAT64_(" + to_string(nyv) + ")";
status = hdf_file->read("YVEC", str_type, instance->yv);
str_type = "FLOAT64_(" + to_string(nzv) + ")";
status = hdf_file->read("ZVEC", str_type, instance->zv);
_nlmmt = 2 * lm * (lm + 2);
_nrvc = nxv * nyv * nzv;
num_elements = 2 * _nlmmt * _nrvc;
elements = new double[num_elements]();
str_type = "FLOAT64_(" + to_string(num_elements) + ")";
status = hdf_file->read("WSUM", str_type, elements);
int index = 0;
for (int wj = 0; wj < _nrvc; wj++) {
for (int wi = 0; wi < _nlmmt; wi++) {
value = elements[index] + elements[index + 1] * I;
_wsum[wi][wj] = value;
index += 2;
} // wi loop
} // wj loop
delete[] elements;
status = hdf_file->close();
delete hdf_file;
}
return instance;
}
TFRFME* TFRFME::from_legacy(const std::string& file_name) {
fstream input;
TFRFME *instance = NULL;
dcomplex **_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;
double vk, exri, an, ff, tra, spd, frsh, exril;
double dval, rval, ival;
input.read(reinterpret_cast<char *>(&lmode), sizeof(int));
input.read(reinterpret_cast<char *>(&lm), sizeof(int));
input.read(reinterpret_cast<char *>(&nkv), sizeof(int));
input.read(reinterpret_cast<char *>(&nxv), sizeof(int));
input.read(reinterpret_cast<char *>(&nyv), sizeof(int));
input.read(reinterpret_cast<char *>(&nzv), sizeof(int));
input.read(reinterpret_cast<char *>(&vk), sizeof(double));
input.read(reinterpret_cast<char *>(&exri), sizeof(double));
input.read(reinterpret_cast<char *>(&an), sizeof(double));
input.read(reinterpret_cast<char *>(&ff), sizeof(double));
input.read(reinterpret_cast<char *>(&tra), sizeof(double));
input.read(reinterpret_cast<char *>(&spd), sizeof(double));
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);
instance->set_param("ff", ff);
instance->set_param("tra", tra);
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));
coord_vec[xi] = dval;
}
coord_vec = instance->get_y();
for (int yi = 0; yi < nyv; yi++) {
input.read(reinterpret_cast<char *>(&dval), sizeof(double));
coord_vec[yi] = dval;
}
coord_vec = instance->get_z();
for (int zi = 0; zi < nzv; zi++) {
input.read(reinterpret_cast<char *>(&dval), sizeof(double));
coord_vec[zi] = dval;
}
int _nlmmt = 2 * lm * (lm + 2);
int _nrvc = nxv * nyv * nzv;
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));
dcomplex value = rval + ival * I;
_wsum[wi][wj] = value;
} // wi loop
} // wj loop
input.close();
} else {
printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
}
return instance;
}
long TFRFME::get_memory_requirement(
int _lmode, int _lm, int _nkv, int _nxv,
int _nyv, int _nzv
) { int _nlmmt = 2 * _lm * (_lm + 2);
int _nrvc = _nxv * _nyv * _nzv;
long size = (long)(8 * sizeof(int) + 8 * sizeof(double));
size += (long)((_nxv + _nyv + _nzv) * sizeof(double));
size += (long)(_nlmmt * _nrvc * sizeof(dcomplex));
return size;
}
double TFRFME::get_param(const std::string& param_name) {
double value;
if (param_name.compare("vk") == 0) value = vk;
else if (param_name.compare("exri") == 0) value = exri;
else if (param_name.compare("an") == 0) value = an;
else if (param_name.compare("ff") == 0) value = ff;
else if (param_name.compare("tra") == 0) value = tra;
else if (param_name.compare("spd") == 0) value = spd;
else if (param_name.compare("frsh") == 0) value = frsh;
else if (param_name.compare("exril") == 0) value = exril;
else if (param_name.compare("lmode") == 0) value = 1.0 * lmode;
else if (param_name.compare("lm") == 0) value = 1.0 * lm;
else if (param_name.compare("nkv") == 0) value = 1.0 * nkv;
else if (param_name.compare("nxv") == 0) value = 1.0 * nxv;
else if (param_name.compare("nyv") == 0) value = 1.0 * nyv;
else if (param_name.compare("nzv") == 0) value = 1.0 * nzv;
else {
string message = "Unrecognized parameter name \"" + param_name + "\"";
throw UnrecognizedParameterException(message);
}
return value;
}
void TFRFME::set_param(const std::string& param_name, double value) {
if (param_name.compare("vk") == 0) vk = value;
else if (param_name.compare("exri") == 0) exri = value;
else if (param_name.compare("an") == 0) an = value;
else if (param_name.compare("ff") == 0) ff = value;
else if (param_name.compare("tra") == 0) tra = value;
else if (param_name.compare("spd") == 0) spd = value;
else if (param_name.compare("frsh") == 0) frsh = value;
else if (param_name.compare("exril") == 0) exril = value;
else {
string message = "Unrecognized parameter name \"" + param_name + "\"";
throw UnrecognizedParameterException(message);
}
}
void TFRFME::write_binary(const std::string& file_name, const std::string& mode) {
if (mode.compare("LEGACY") == 0) {
write_legacy(file_name);
} else if (mode.compare("HDF5") == 0) {
write_hdf5(file_name);
} else {
string message = "Unknown format mode: \"" + mode + "\"";
throw UnrecognizedFormatException(message);
}
}
void TFRFME::write_hdf5(const std::string& file_name) {
List<string> rec_name_list(1);
List<string> rec_type_list(1);
List<void *> rec_ptr_list(1);
herr_t status;
string str_type;
rec_name_list.set(0, "LMODE");
rec_type_list.set(0, "INT32_(1)");
rec_ptr_list.set(0, &lmode);
rec_name_list.append("LM");
rec_type_list.append("INT32_(1)");
rec_ptr_list.append(&lm);
rec_name_list.append("NKV"); rec_type_list.append("INT32_(1)");
rec_ptr_list.append(&nkv);
rec_name_list.append("NXV");
rec_type_list.append("INT32_(1)");
rec_ptr_list.append(&nxv);
rec_name_list.append("NYV");
rec_type_list.append("INT32_(1)");
rec_ptr_list.append(&nyv);
rec_name_list.append("NZV");
rec_type_list.append("INT32_(1)");
rec_ptr_list.append(&nzv);
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("AN");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&an);
rec_name_list.append("FF");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&ff);
rec_name_list.append("TRA");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&tra);
rec_name_list.append("SPD");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&spd);
rec_name_list.append("FRSH");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&frsh);
rec_name_list.append("EXRIL");
rec_type_list.append("FLOAT64_(1)");
rec_ptr_list.append(&exril);
str_type = "FLOAT64_(" + to_string(nxv) + ")";
rec_name_list.append("XVEC");
rec_type_list.append(str_type);
rec_ptr_list.append(xv);
str_type = "FLOAT64_(" + to_string(nyv) + ")";
rec_name_list.append("YVEC");
rec_type_list.append(str_type);
rec_ptr_list.append(yv);
str_type = "FLOAT64_(" + to_string(nzv) + ")";
rec_name_list.append("ZVEC");
rec_type_list.append(str_type);
rec_ptr_list.append(zv);
rec_name_list.append("WSUM");
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]();
int index = 0;
for (int wj = 0; wj < nrvc; wj++) {
for (int wi = 0; wi < nlmmt; wi++) {
ptr_elements[index++] = real(wsum[wi][wj]);
ptr_elements[index++] = imag(wsum[wi][wj]);
} // wi loop
} // wj loop
rec_ptr_list.append(ptr_elements);
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++) status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
status = hdf_file->close();
delete[] ptr_elements;
delete[] rec_names;
delete[] rec_types;
delete[] rec_pointers;
delete hdf_file;
}
void TFRFME::write_legacy(const std::string& file_name) {
fstream output;
output.open(file_name.c_str(), ios::out | ios::binary);
if (output.is_open()) {
output.write(reinterpret_cast<char *>(&lmode), sizeof(int));
output.write(reinterpret_cast<char *>(&lm), sizeof(int));
output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
output.write(reinterpret_cast<char *>(&nxv), sizeof(int));
output.write(reinterpret_cast<char *>(&nyv), sizeof(int));
output.write(reinterpret_cast<char *>(&nzv), sizeof(int));
output.write(reinterpret_cast<char *>(&vk), sizeof(double));
output.write(reinterpret_cast<char *>(&exri), sizeof(double));
output.write(reinterpret_cast<char *>(&an), sizeof(double));
output.write(reinterpret_cast<char *>(&ff), sizeof(double));
output.write(reinterpret_cast<char *>(&tra), sizeof(double));
output.write(reinterpret_cast<char *>(&spd), sizeof(double));
output.write(reinterpret_cast<char *>(&frsh), sizeof(double));
output.write(reinterpret_cast<char *>(&exril), sizeof(double));
for (int xi = 0; xi < nxv; xi++)
output.write(reinterpret_cast<char *>(&(xv[xi])), sizeof(double));
for (int yi = 0; yi < nyv; yi++)
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 wj = 0; wj < nrvc; wj++) {
for (int wi = 0; wi < nlmmt; wi++) {
double rval = real(wsum[wi][wj]);
double ival = imag(wsum[wi][wj]);
output.write(reinterpret_cast<char *>(&rval), sizeof(double));
output.write(reinterpret_cast<char *>(&ival), sizeof(double));
} // wi loop
} // wj loop
output.close();
} else { // Should never occur
printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
}
}
bool TFRFME::operator ==(const TFRFME& other) {
if (lmode != other.lmode) {
return false;
}
if (lm != other.lm) {
return false;
}
if (nkv != other.nkv) {
return false;
}
if (nxv != other.nxv) {
return false;
}
if (nyv != other.nyv) {
return false;
}
if (nzv != other.nzv) {
return false;
}
if (vk != other.vk) {
return false;
} if (exri != other.exri) {
return false;
}
if (an != other.an) {
return false;
}
if (ff != other.ff) {
return false;
}
if (tra != other.tra) {
return false;
}
if (spd != other.spd) {
return false;
}
if (frsh != other.frsh) {
return false;
}
if (exril != other.exril) {
return false;
}
for (int xi = 0; xi < nxv; xi++) {
if (xv[xi] != other.xv[xi]) {
return false;
}
}
for (int yi = 0; yi < nyv; yi++) {
if (yv[yi] != other.yv[yi]) {
return false;
}
}
for (int zi = 0; zi < nzv; zi++) {
if (zv[zi] != other.zv[zi]) {
return false;
}
}
for (int wi = 0; wi < nlmmt; wi++) {
for (int wj = 0; wj < nrvc; wj++) {
if (wsum[wi][wj] != other.wsum[wi][wj]) {
return false;
}
} // wj loop
} // wi loop
return true;
}
// >>> END OF TFRFME CLASS IMPLEMENTATION <<<