Make TFRFME matrices and vectors contiguous (64b58343) · Commits · Giacomo Mulas / NP_TMcode

src/include/tfrfme.h

+70 −53

Original line number	Diff line number	Diff line
		@@ -313,38 +313,37 @@ public:
		class TFRFME {
		protected:
		//! NLMMT = 2 * LM * (LM + 2)
		int nlmmt;
		int _nlmmt;
		//! NRVC = NXV * NYV * NZV
		int nrvc;

		int _nrvc;
		//! Field expansion mode identifier.
		int lmode;
		//! Maximim field expansion order;
		int lm;
		int _lmode;
		//! Maximum field expansion order.
		int _lm;
		//! QUESTION: definition?
		int nkv;
		int _nkv;
		//! Number of computed X coordinates.
		int nxv;
		int _nxv;
		//! Number of computed Y coordinates.
		int nyv;
		int _nyv;
		//! Number of computed Z coordinates.
		int nzv;
		//! Wave number in scale units
		double vk;
		int _nzv;
		//! Vacuum wave number.
		double _vk;
		//! External medium refractive index
		double exri;
		double _exri;
		//! QUESTION: definition?
		double an;
		double _an;
		//! QUESTION: definition?
		double ff;
		double _ff;
		//! QUESTION: definition?
		double tra;
		double _tra;
		//! QUESTION: definition?
		double spd;
		double _spd;
		//! QUESTION: definition?
		double frsh;
		double _frsh;
		//! QUESTION: definition?
		double exril;
		double _exril;
		//! Vector of computed x positions
		double *xv;
		//! Vector of computed y positions
		@@ -352,7 +351,7 @@ protected:
		//! Vector of computed z positions
		double *zv;
		//! QUESTION: definition?
		dcomplex **wsum;
		dcomplex *vec_wsum;

		/*! \brief Load a configuration instance from a HDF5 binary file.
		*
		@@ -383,18 +382,51 @@ protected:
		void write_legacy(const std::string& file_name);

		public:
		//! Read-only view on NLMMT.
		const int& nlmmt = _nlmmt;
		//! Read-only view on NRVC.
		const int& nrvc = _nrvc;
		//! Read-only view on field expansion mode identifier.
		const int& lmode = _lmode;
		//! Read-only view on maximum field expansion order.
		const int& lm = _lm;
		//! QUESTION: definition?
		const int& nkv = _nkv;
		//! Read-only view on number of computed X coordinates.
		const int& nxv = _nxv;
		//! Read-only view on number of computed Y coordinates.
		const int& nyv = _nyv;
		//! Read-only view on number of computed Z coordinates.
		const int& nzv = _nzv;
		//! Read-only view on vacuum wave number.
		const double& vk = _vk;
		//! Read-only view on external medium refractive index
		const double& exri = _exri;
		//! QUESTION: definition?
		const double& an = _an;
		//! QUESTION: definition?
		const double& ff = _ff;
		//! QUESTION: definition?
		const double& tra = _tra;
		//! QUESTION: definition?
		const double& spd = _spd;
		//! QUESTION: definition?
		const double& frsh = _frsh;
		//! QUESTION: definition?
		const double& exril = _exril;
		//! QUESTION: definition?
		dcomplex **wsum;

		/*! \brief Trapping configuration instance constructor.
		*
		* \param _lmode: `int` Order expansion mode flag.
		* \param _lm: `int` Maximum field expansion order.
		* \param _nkv: `int` Number of wave vector coordinates. QUESTION: correct?
		* \param _nxv: `int` Number of computed X coordinates.
		* \param _nyv: `int` Number of computed Y coordinates.
		* \param _nzv: `int` Number of computed Z coordinates.
		* \param lmode: `int` Order expansion mode flag.
		* \param lm: `int` Maximum field expansion order.
		* \param nkv: `int` Number of wave vector coordinates. QUESTION: correct?
		* \param nxv: `int` Number of computed X coordinates.
		* \param nyv: `int` Number of computed Y coordinates.
		* \param nzv: `int` Number of computed Z coordinates.
		*/
		TFRFME(
		int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv
		);
		TFRFME(int lmode, int lm, int nkv, int nxv, int nyv, int nzv);

		/*! \brief Trapping configuration instance destroyer.
		*/
		@@ -408,35 +440,20 @@ public:
		* \return instance: `TFRFME *` Pointer to a newly created configuration
		* instance.
		*/
		static TFRFME* from_binary(const std::string& file_name, const std::string& mode="LEGACY");

		/*! \brief Get the pointer to the WSUM matrix.
		*
		* \return value: `complex double **` Pointer to the WSUM matrix.
		*/
		dcomplex **get_matrix() { return wsum; }
		static TFRFME* from_binary(
		const std::string& file_name, const std::string& mode="LEGACY"
		);

		/*! \brief Calculate the necessary amount of memory to create a new instance.
		*
		* \param _lmode: `int` Order expansion mode flag.
		* \param _lm: `int` Maximum field expansion order.
		* \param _nkv: `int` Number of radial vector coordinates. QUESTION: correct?
		* \param _nxv: `int` Number of computed X coordinates.
		* \param _nyv: `int` Number of computed Y coordinates.
		* \param _nzv: `int` Number of computed Z coordinates.
		* \param lm: `int` Maximum field expansion order.
		* \param nkv: `int` Number of radial vector coordinates. QUESTION: correct?
		* \param nxv: `int` Number of computed X coordinates.
		* \param nyv: `int` Number of computed Y coordinates.
		* \param nzv: `int` Number of computed Z coordinates.
		* \return size: `long` The necessary memory size in bytes.
		*/
		static long get_memory_requirement(
		int _lmode, int _lm, int _nkv, int _nxv,
		int _nyv, int _nzv
		);

		/*! \brief Get a configuration parameter.
		*
		* \param param_name: `string` Name of the parameter.
		* \return value: `double` The value of the requested parameter.
		*/
		double get_param(const std::string& param_name);
		static long get_size(int lm, int nkv, int nxv, int nyv, int nzv);

		/*! \brief Get the pointer to the X coordinates vector.
		*

src/libnptm/tfrfme.cpp

+98 −127

Original line number	Diff line number	Diff line
		@@ -563,37 +563,38 @@ bool Swap2::operator ==(Swap2 &other) {
		// >>> 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;
		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;
		_nlmmt = _lm * (_lm + 2) * 2;
		_nrvc = _nxv * _nyv * _nzv;
		vec_wsum = new dcomplex[nrvc * nlmmt]();
		wsum = new dcomplex*[nlmmt];
		for (int wi = 0; wi < nlmmt; wi++) wsum[wi] = new dcomplex[nrvc]();
		for (int wi = 0; wi < nlmmt; wi++) wsum[wi] = vec_wsum + wi * nrvc;
		}

		TFRFME::~TFRFME() {
		delete[] xv;
		delete[] yv;
		delete[] zv;
		for (int wi = nlmmt - 1; wi > -1; wi--) delete[] wsum[wi];
		delete[] vec_wsum;
		delete[] wsum;
		}

		@@ -617,9 +618,8 @@ TFRFME* TFRFME::from_hdf5(const std::string& file_name) {
		herr_t status = hdf_file->get_status();
		double *elements;
		string str_type;
		int _nlmmt, _nrvc, num_elements;
		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;
		@@ -638,7 +638,6 @@ TFRFME* TFRFME::from_hdf5(const std::string& file_name) {
		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);
		@@ -653,17 +652,17 @@ TFRFME* TFRFME::from_hdf5(const std::string& file_name) {
		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;
		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++) {
		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;
		instance->wsum[wi][wj] = value;
		index += 2;
		} // wi loop
		} // wj loop
		@@ -677,7 +676,6 @@ TFRFME* TFRFME::from_hdf5(const std::string& file_name) {
		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()) {
		@@ -699,7 +697,6 @@ TFRFME* TFRFME::from_legacy(const std::string& file_name) {
		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);
		@@ -723,14 +720,14 @@ TFRFME* TFRFME::from_legacy(const std::string& file_name) {
		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++) {
		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;
		instance->wsum[wi][wj] = value;
		} // wi loop
		} // wj loop
		input.close();
		@@ -740,50 +737,24 @@ TFRFME* TFRFME::from_legacy(const std::string& file_name) {
		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));
		long TFRFME::get_size(int lm, int nkv, int nxv, int nyv, int nzv) {
		int nlmmt = 2 * lm * (lm + 2);
		int nrvc = nxv * nyv * nzv;
		long size = (long)(16 * sizeof(int) + 16 * 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;
		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);
		@@ -809,46 +780,46 @@ void TFRFME::write_hdf5(const std::string& file_name) {
		string str_type;
		rec_name_list.set(0, "LMODE");
		rec_type_list.set(0, "INT32_(1)");
		rec_ptr_list.set(0, &lmode);
		rec_ptr_list.set(0, &_lmode);
		rec_name_list.append("LM");
		rec_type_list.append("INT32_(1)");
		rec_ptr_list.append(&lm);
		rec_ptr_list.append(&_lm);
		rec_name_list.append("NKV");
		rec_type_list.append("INT32_(1)");
		rec_ptr_list.append(&nkv);
		rec_ptr_list.append(&_nkv);
		rec_name_list.append("NXV");
		rec_type_list.append("INT32_(1)");
		rec_ptr_list.append(&nxv);
		rec_ptr_list.append(&_nxv);
		rec_name_list.append("NYV");
		rec_type_list.append("INT32_(1)");
		rec_ptr_list.append(&nyv);
		rec_ptr_list.append(&_nyv);
		rec_name_list.append("NZV");
		rec_type_list.append("INT32_(1)");
		rec_ptr_list.append(&nzv);
		rec_ptr_list.append(&_nzv);
		rec_name_list.append("VK");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&vk);
		rec_ptr_list.append(&_vk);
		rec_name_list.append("EXRI");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&exri);
		rec_ptr_list.append(&_exri);
		rec_name_list.append("AN");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&an);
		rec_ptr_list.append(&_an);
		rec_name_list.append("FF");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&ff);
		rec_ptr_list.append(&_ff);
		rec_name_list.append("TRA");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&tra);
		rec_ptr_list.append(&_tra);
		rec_name_list.append("SPD");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&spd);
		rec_ptr_list.append(&_spd);
		rec_name_list.append("FRSH");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&frsh);
		rec_ptr_list.append(&_frsh);
		rec_name_list.append("EXRIL");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&exril);
		rec_ptr_list.append(&_exril);
		str_type = "FLOAT64_(" + to_string(nxv) + ")";
		rec_name_list.append("XVEC");
		rec_type_list.append(str_type);
		@@ -898,28 +869,28 @@ 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 *>(&_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++)
		for (int yi = 0; yi < _nyv; yi++)
		output.write(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));
		for (int zi = 0; zi < nzv; zi++)
		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++) {
		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));
		@@ -933,65 +904,65 @@ void TFRFME::write_legacy(const std::string& file_name) {
		}

		bool TFRFME::operator ==(const TFRFME& other) {
		if (lmode != other.lmode) {
		if (_lmode != other._lmode) {
		return false;
		}
		if (lm != other.lm) {
		if (_lm != other._lm) {
		return false;
		}
		if (nkv != other.nkv) {
		if (_nkv != other._nkv) {
		return false;
		}
		if (nxv != other.nxv) {
		if (_nxv != other._nxv) {
		return false;
		}
		if (nyv != other.nyv) {
		if (_nyv != other._nyv) {
		return false;
		}
		if (nzv != other.nzv) {
		if (_nzv != other._nzv) {
		return false;
		}
		if (vk != other.vk) {
		if (_vk != other._vk) {
		return false;
		}
		if (exri != other.exri) {
		if (_exri != other._exri) {
		return false;
		}
		if (an != other.an) {
		if (_an != other._an) {
		return false;
		}
		if (ff != other.ff) {
		if (_ff != other._ff) {
		return false;
		}
		if (tra != other.tra) {
		if (_tra != other._tra) {
		return false;
		}
		if (spd != other.spd) {
		if (_spd != other._spd) {
		return false;
		}
		if (frsh != other.frsh) {
		if (_frsh != other._frsh) {
		return false;
		}
		if (exril != other.exril) {
		if (_exril != other._exril) {
		return false;
		}
		for (int xi = 0; xi < nxv; xi++) {
		for (int xi = 0; xi < _nxv; xi++) {
		if (xv[xi] != other.xv[xi]) {
		return false;
		}
		}
		for (int yi = 0; yi < nyv; yi++) {
		for (int yi = 0; yi < _nyv; yi++) {
		if (yv[yi] != other.yv[yi]) {
		return false;
		}
		}
		for (int zi = 0; zi < nzv; zi++) {
		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++) {
		for (int wi = 0; wi < _nlmmt; wi++) {
		for (int wj = 0; wj < _nrvc; wj++) {
		if (wsum[wi][wj] != other.wsum[wi][wj]) {
		return false;
		}