Fix DC0M matrix to adapt to multi-layer case (0236f4c2) · Commits · Giacomo Mulas / NP_TMcode

src/include/Configuration.h

+4 −0

Original line number	Original line	Diff line number	Diff line
	@@ -312,6 +312,8 @@ protected:
	double _wp;		double _wp;
	//! \brief Peak XI. QUESTION: correct?		//! \brief Peak XI. QUESTION: correct?
	double _xip;		double _xip;
			//! \brief Maximum number of layers for the particle components.
			int _max_layers;
	//! \brief Flag to control whether to add an external layer.		//! \brief Flag to control whether to add an external layer.
	bool _use_external_sphere;		bool _use_external_sphere;

	@@ -377,6 +379,8 @@ public:
	const double& wp = _wp;		const double& wp = _wp;
	//! \brief Read only view on peak XI.		//! \brief Read only view on peak XI.
	const double& xip = _xip;		const double& xip = _xip;
			//! \brief Read only view on the maximum number of layers for the particle components.
			const int& max_layers = _max_layers;
	//! \brief Read only view on flag to control whether to add an external layer.		//! \brief Read only view on flag to control whether to add an external layer.
	const bool& use_external_sphere = _use_external_sphere;		const bool& use_external_sphere = _use_external_sphere;

src/libnptm/Configuration.cpp

+52 −51

Original line number	Original line	Diff line number	Diff line
	@@ -356,6 +356,10 @@ ScattererConfiguration::ScattererConfiguration(
	_iog_vec = iog_vector;		_iog_vec = iog_vector;
	_radii_of_spheres = ros_vector;		_radii_of_spheres = ros_vector;
	_nshl_vec = nshl_vector;		_nshl_vec = nshl_vector;
			_max_layers = 0;
			for (int li = 0; li < _configurations; li++) {
			if (_max_layers < _nshl_vec[li]) _max_layers = _nshl_vec[li];
			}
	_rcf = rcf_vector;		_rcf = rcf_vector;
	_idfc = dielectric_func_type;		_idfc = dielectric_func_type;
	_dc0_matrix = dc_matrix;		_dc0_matrix = dc_matrix;
	@@ -391,12 +395,13 @@ ScattererConfiguration::ScattererConfiguration(const ScattererConfiguration& rhs
	_exdc = rhs._exdc;		_exdc = rhs._exdc;
	_wp = rhs._wp;		_wp = rhs._wp;
	_xip = rhs._xip;		_xip = rhs._xip;
			_max_layers = rhs._max_layers;
	_iog_vec = new int[_number_of_spheres]();		_iog_vec = new int[_number_of_spheres]();
	_radii_of_spheres = new double[_number_of_spheres]();		_radii_of_spheres = new double[_number_of_spheres]();
	_nshl_vec = new int[_configurations]();		_nshl_vec = new int[_configurations]();
	_rcf = new double*[_configurations];		_rcf = new double*[_configurations];
	_scale_vec = new double[_number_of_scales]();		_scale_vec = new double[_number_of_scales]();
	_dc0_matrix = new dcomplex**[_configurations];		_dc0_matrix = new dcomplex**[_max_layers];
	for (int si = 0; si < _number_of_scales; si++) _scale_vec[si] = rhs._scale_vec[si];		for (int si = 0; si < _number_of_scales; si++) _scale_vec[si] = rhs._scale_vec[si];
	for (int si = 0; si < _number_of_spheres; si++) {		for (int si = 0; si < _number_of_spheres; si++) {
	_iog_vec[si] = rhs._iog_vec[si];		_iog_vec[si] = rhs._iog_vec[si];
	@@ -406,11 +411,13 @@ ScattererConfiguration::ScattererConfiguration(const ScattererConfiguration& rhs
	_radii_of_spheres[si] = rhs._radii_of_spheres[si];		_radii_of_spheres[si] = rhs._radii_of_spheres[si];
	_nshl_vec[si] = rhs._nshl_vec[si];		_nshl_vec[si] = rhs._nshl_vec[si];
	_rcf[si] = new double[_nshl_vec[si]]();		_rcf[si] = new double[_nshl_vec[si]]();
	_dc0_matrix[si] = new dcomplex*[_number_of_spheres];
	for (int sj = 0; sj < _nshl_vec[si]; sj++) _rcf[si][sj] = rhs._rcf[si][sj];		for (int sj = 0; sj < _nshl_vec[si]; sj++) _rcf[si][sj] = rhs._rcf[si][sj];
	for (int sj = 0; sj < _number_of_spheres; sj++) {		}
	_dc0_matrix[si][sj] = new dcomplex[dim3]();		for (int li = 0; li < _max_layers; li++) {
	for (int sk = 0; sk < dim3; sk++) _dc0_matrix[si][sj][sk] = rhs._dc0_matrix[si][sj][sk];		_dc0_matrix[li] = new dcomplex*[_number_of_spheres];
			for (int lj = 0; lj < _number_of_spheres; lj++) {
			_dc0_matrix[li][lj] = new dcomplex[dim3]();
			for (int lk = 0; lk < dim3; lk++) _dc0_matrix[li][lj][lk] = rhs._dc0_matrix[li][lj][lk];
	}		}
	}		}
	}		}
	@@ -431,6 +438,7 @@ ScattererConfiguration::ScattererConfiguration(const mixMPI *mpidata)
	itemp = sizeof(bool);		itemp = sizeof(bool);
	char ptemp = (char ) &_use_external_sphere;		char ptemp = (char ) &_use_external_sphere;
	MPI_Bcast(ptemp, itemp, MPI_CHAR, 0, MPI_COMM_WORLD);		MPI_Bcast(ptemp, itemp, MPI_CHAR, 0, MPI_COMM_WORLD);
			MPI_Bcast(&_max_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_exdc, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_exdc, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_wp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_wp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	@@ -442,16 +450,18 @@ ScattererConfiguration::ScattererConfiguration(const mixMPI *mpidata)
	MPI_Bcast(_nshl_vec, _configurations, MPI_INT, 0, MPI_COMM_WORLD);		MPI_Bcast(_nshl_vec, _configurations, MPI_INT, 0, MPI_COMM_WORLD);
	_rcf = new double*[_configurations];		_rcf = new double*[_configurations];
	_scale_vec = new double[_number_of_scales]();		_scale_vec = new double[_number_of_scales]();
	_dc0_matrix = new dcomplex**[_configurations];		_dc0_matrix = new dcomplex**[_max_layers];
	MPI_Bcast(_scale_vec, _number_of_scales, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(_scale_vec, _number_of_scales, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	int dim3 = (_idfc == 0) ? _number_of_scales : 1;		int dim3 = (_idfc == 0) ? _number_of_scales : 1;
	for (int si = 0; si < _configurations; si++) {		for (int si = 0; si < _configurations; si++) {
	_rcf[si] = new double[_nshl_vec[si]]();		_rcf[si] = new double[_nshl_vec[si]]();
	MPI_Bcast(_rcf[si], _nshl_vec[si], MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(_rcf[si], _nshl_vec[si], MPI_DOUBLE, 0, MPI_COMM_WORLD);
	_dc0_matrix[si] = new dcomplex*[_number_of_spheres];		}
	for (int sj = 0; sj < _number_of_spheres; sj++) {		for (int li = 0; li < _max_layers; li++) {
	_dc0_matrix[si][sj] = new dcomplex[dim3]();		_dc0_matrix[li] = new dcomplex*[_number_of_spheres];
	MPI_Bcast(_dc0_matrix[si][sj], dim3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);		for (int lj = 0; lj < _number_of_spheres; lj++) {
			_dc0_matrix[li][lj] = new dcomplex[dim3]();
			MPI_Bcast(_dc0_matrix[li][lj], dim3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
	}		}
	}		}
	}		}
	@@ -469,6 +479,7 @@ void ScattererConfiguration::mpibcast(const mixMPI *mpidata) {
	itemp = sizeof(bool);		itemp = sizeof(bool);
	char ptemp = (char ) &_use_external_sphere;		char ptemp = (char ) &_use_external_sphere;
	MPI_Bcast(ptemp, itemp, MPI_CHAR, 0, MPI_COMM_WORLD);		MPI_Bcast(ptemp, itemp, MPI_CHAR, 0, MPI_COMM_WORLD);
			MPI_Bcast(&_max_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_exdc, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_exdc, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_wp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_wp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	@@ -479,8 +490,10 @@ void ScattererConfiguration::mpibcast(const mixMPI *mpidata) {
	int dim3 = (_idfc == 0) ? _number_of_scales : 1;		int dim3 = (_idfc == 0) ? _number_of_scales : 1;
	for (int si = 0; si < _configurations; si++) {		for (int si = 0; si < _configurations; si++) {
	MPI_Bcast(_rcf[si], _nshl_vec[si], MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(_rcf[si], _nshl_vec[si], MPI_DOUBLE, 0, MPI_COMM_WORLD);
	for (int sj = 0; sj < _number_of_spheres; sj++) {		}
	MPI_Bcast(_dc0_matrix[si][sj], dim3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);		for (int li = 0; li < _max_layers; li++) {
			for (int lj = 0; lj < _number_of_spheres; lj++) {
			MPI_Bcast(_dc0_matrix[li][lj], dim3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
	}		}
	}		}
	}		}
	@@ -697,6 +710,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de
	double *ros_vector = new double[configurations]();		double *ros_vector = new double[configurations]();
	int *nshl_vector = new int[configurations]();		int *nshl_vector = new int[configurations]();
	double *rcf_vector = new double[configurations];		double *rcf_vector = new double[configurations];
			int max_layers = 0;
	last_configuration = 0;		last_configuration = 0;
	for (int i113 = 1; i113 <= fnsph; i113++) {		for (int i113 = 1; i113 <= fnsph; i113++) {
	if (iog_vector[i113 - 1] < i113) continue;		if (iog_vector[i113 - 1] < i113) continue;
	@@ -705,6 +719,8 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de
	re = regex("[0-9]+");		re = regex("[0-9]+");
	regex_search(str_target, m, re);		regex_search(str_target, m, re);
	nshl_vector[last_configuration - 1] = stoi(m.str());		nshl_vector[last_configuration - 1] = stoi(m.str());
			if (max_layers < nshl_vector[last_configuration - 1] + 1)
			max_layers = nshl_vector[last_configuration - 1] + 1;
	str_target = m.suffix().str();		str_target = m.suffix().str();
	re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");		re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
	regex_search(str_target, m, re);		regex_search(str_target, m, re);
	@@ -725,10 +741,10 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de
	} // ns112 loop		} // ns112 loop
	} // i113 loop		} // i113 loop

	dcomplex *dc0m = new dcomplex[configurations];		dcomplex *dc0m = new dcomplex[max_layers];
	dcomplex *dc0 = new dcomplex[configurations]();		dcomplex *dc0 = new dcomplex[max_layers]();
	int dim3 = (fidfc == 0) ? fnxi : 1;		int dim3 = (fidfc == 0) ? fnxi : 1;
	for (int di = 0; di < configurations; di++) {		for (int di = 0; di < max_layers; di++) {
	dc0m[di] = new dcomplex*[fnsph];		dc0m[di] = new dcomplex*[fnsph];
	for (int dj = 0; dj < fnsph; dj++) dc0m[di][dj] = new dcomplex[dim3]();		for (int dj = 0; dj < fnsph; dj++) dc0m[di][dj] = new dcomplex[dim3]();
	} // di loop		} // di loop
	@@ -763,21 +779,9 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de
	delete[] dc0;		delete[] dc0;

	ScattererConfiguration *config = new ScattererConfiguration(		ScattererConfiguration *config = new ScattererConfiguration(
	fnsph,		fnsph, configurations, variable_vector, fnxi, variable_name,
	configurations,		iog_vector, ros_vector, nshl_vector, rcf_vector, fidfc,
	variable_vector,		dc0m, (fies > 0), fexdc, fwp, fxip
	fnxi,
	variable_name,
	iog_vector,
	ros_vector,
	nshl_vector,
	rcf_vector,
	fidfc,
	dc0m,
	(fies > 0),
	fexdc,
	fwp,
	fxip
	);		);
	delete[] file_lines;		delete[] file_lines;
	delete[] variable_vector;		delete[] variable_vector;
	@@ -864,21 +868,9 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(const std::string& fil
	status = hdf_file->close();		status = hdf_file->close();
	delete hdf_file;		delete hdf_file;
	conf = new ScattererConfiguration(		conf = new ScattererConfiguration(
	nsph,		nsph, configuration_count, xi_vec, nxi, "XIV",
	configuration_count,		iog, ros_vector, nshl_vector, rcf_vector, _idfc,
	xi_vec,		dc0m, (ies == 1), _exdc, _wp, _xip
	nxi,
	"XIV",
	iog,
	ros_vector,
	nshl_vector,
	rcf_vector,
	_idfc,
	dc0m,
	(ies == 1),
	_exdc,
	_wp,
	_xip
	);		);
	delete[] xi_vec;		delete[] xi_vec;
	}		}
	@@ -932,9 +924,13 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f
	for (int nsi = 0; nsi < nsh; nsi++)		for (int nsi = 0; nsi < nsh; nsi++)
	input.read(reinterpret_cast<char *>(&(_rcf_vec[last_configuration - 1][nsi])), sizeof(double));		input.read(reinterpret_cast<char *>(&(_rcf_vec[last_configuration - 1][nsi])), sizeof(double));
	}		}
			int max_layers = 0;
			for (int li = 0; li < configurations; li++) {
			if (max_layers < _nshl_vec[li]) max_layers = _nshl_vec[li];
			}
	_dc0m = new dcomplex**[configurations];		_dc0m = new dcomplex**[configurations];
	int dim3 = (_idfc == 0) ? _nxi : 1;		int dim3 = (_idfc == 0) ? _nxi : 1;
	for (int di = 0; di < configurations; di++) {		for (int di = 0; di < max_layers; di++) {
	_dc0m[di] = new dcomplex*[_nsph];		_dc0m[di] = new dcomplex*[_nsph];
	for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new dcomplex[dim3]();		for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new dcomplex[dim3]();
	} // di loop		} // di loop
	@@ -995,7 +991,7 @@ void ScattererConfiguration::print() {
	for (int i = 0; i < _number_of_scales; i++) printf("\t%lg", _scale_vec[i]);		for (int i = 0; i < _number_of_scales; i++) printf("\t%lg", _scale_vec[i]);
	printf("\t]\n");		printf("\t]\n");
	printf("DC0M = [\n");		printf("DC0M = [\n");
	for (int i = 0; i < _configurations; i++) {		for (int i = 0; i < _max_layers; i++) {
	printf(" [\n");		printf(" [\n");
	for (int j = 0; j < _number_of_spheres; j++) {		for (int j = 0; j < _number_of_spheres; j++) {
	printf(" [");		printf(" [");
	@@ -1079,7 +1075,7 @@ void ScattererConfiguration::write_hdf5(const std::string& file_name) {
	}		}

	int dim3 = (idfc == 0) ? _number_of_scales : 1;		int dim3 = (idfc == 0) ? _number_of_scales : 1;
	int dc0m_size = 2 * dim3 * _number_of_spheres * _configurations;		int dc0m_size = 2 * dim3 * _number_of_spheres * _max_layers;
	double *dc0m = new double[dc0m_size]();		double *dc0m = new double[dc0m_size]();
	int dc0_index = 0;		int dc0_index = 0;
	for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) {		for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) {
	@@ -1376,6 +1372,9 @@ bool ScattererConfiguration::operator ==(const ScattererConfiguration &other) {
	if (_idfc != other._idfc) {		if (_idfc != other._idfc) {
	return false;		return false;
	}		}
			if (_max_layers != other._max_layers) {
			return false;
			}
	if (_exdc != other._exdc) {		if (_exdc != other._exdc) {
	return false;		return false;
	}		}
	@@ -1411,13 +1410,15 @@ bool ScattererConfiguration::operator ==(const ScattererConfiguration &other) {
	return false;		return false;
	}		}
	} // rj loop		} // rj loop
			} // ri loop
			for (int li = 0; li < _max_layers; li++) {
	for (int dj = 0; dj < _number_of_spheres; dj++) {		for (int dj = 0; dj < _number_of_spheres; dj++) {
	for (int dk = 0; dk < dim3; dk++) {		for (int dk = 0; dk < dim3; dk++) {
	if (_dc0_matrix[ri][dj][dk] != other._dc0_matrix[ri][dj][dk]) {		if (_dc0_matrix[li][dj][dk] != other._dc0_matrix[li][dj][dk]) {
	return false;		return false;
	}		}
	} // dk loop		} // dk loop
	} // dj loop		} // dj loop
	} // ri loop		}
	return true;		return true;
	}		}