diff --git a/src/include/Configuration.h b/src/include/Configuration.h index 81477ab6b223ef61f158137702a3a4c5753ed3e4..d23aa6ff178ee60c8f251ab9b88a86349b6e1762 100644 --- a/src/include/Configuration.h +++ b/src/include/Configuration.h @@ -312,6 +312,8 @@ protected: double _wp; //! \brief Peak XI. QUESTION: correct? double _xip; + //! \brief Maximum number of layers for the particle components. + int _max_layers; //! \brief Flag to control whether to add an external layer. bool _use_external_sphere; @@ -377,6 +379,8 @@ public: const double& wp = _wp; //! \brief Read only view on peak XI. 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. const bool& use_external_sphere = _use_external_sphere; diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp index fe78967c8e6a18705b076d266cd29fe4686c8f69..9d52636290029330e08e830e1bfa6fd2e4a577a9 100644 --- a/src/libnptm/Configuration.cpp +++ b/src/libnptm/Configuration.cpp @@ -356,6 +356,10 @@ ScattererConfiguration::ScattererConfiguration( _iog_vec = iog_vector; _radii_of_spheres = ros_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; _idfc = dielectric_func_type; _dc0_matrix = dc_matrix; @@ -391,12 +395,13 @@ ScattererConfiguration::ScattererConfiguration(const ScattererConfiguration& rhs _exdc = rhs._exdc; _wp = rhs._wp; _xip = rhs._xip; + _max_layers = rhs._max_layers; _iog_vec = new int[_number_of_spheres](); _radii_of_spheres = new double[_number_of_spheres](); _nshl_vec = new int[_configurations](); _rcf = new double*[_configurations]; _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_spheres; 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]; _nshl_vec[si] = rhs._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 < _number_of_spheres; sj++) { - _dc0_matrix[si][sj] = new dcomplex[dim3](); - for (int sk = 0; sk < dim3; sk++) _dc0_matrix[si][sj][sk] = rhs._dc0_matrix[si][sj][sk]; + } + for (int li = 0; li < _max_layers; li++) { + _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); char *ptemp = (char *) &_use_external_sphere; 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(&_wp, 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); _rcf = new double*[_configurations]; _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); int dim3 = (_idfc == 0) ? _number_of_scales : 1; for (int si = 0; si < _configurations; si++) { _rcf[si] = new double[_nshl_vec[si]](); 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++) { - _dc0_matrix[si][sj] = new dcomplex[dim3](); - MPI_Bcast(_dc0_matrix[si][sj], dim3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD); + } + for (int li = 0; li < _max_layers; li++) { + _dc0_matrix[li] = new dcomplex*[_number_of_spheres]; + 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); char *ptemp = (char *) &_use_external_sphere; 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(&_wp, 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; for (int si = 0; si < _configurations; si++) { 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](); int *nshl_vector = new int[configurations](); double **rcf_vector = new double*[configurations]; + int max_layers = 0; last_configuration = 0; for (int i113 = 1; i113 <= fnsph; i113++) { if (iog_vector[i113 - 1] < i113) continue; @@ -705,6 +719,8 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de re = regex("[0-9]+"); regex_search(str_target, m, re); 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(); re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?"); regex_search(str_target, m, re); @@ -725,10 +741,10 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de } // ns112 loop } // i113 loop - dcomplex ***dc0m = new dcomplex**[configurations]; - dcomplex *dc0 = new dcomplex[configurations](); + dcomplex ***dc0m = new dcomplex**[max_layers]; + dcomplex *dc0 = new dcomplex[max_layers](); 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]; for (int dj = 0; dj < fnsph; dj++) dc0m[di][dj] = new dcomplex[dim3](); } // di loop @@ -763,22 +779,10 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de delete[] dc0; ScattererConfiguration *config = new ScattererConfiguration( - fnsph, - configurations, - variable_vector, - fnxi, - variable_name, - iog_vector, - ros_vector, - nshl_vector, - rcf_vector, - fidfc, - dc0m, - (fies > 0), - fexdc, - fwp, - fxip - ); + fnsph, configurations, variable_vector, fnxi, variable_name, + iog_vector, ros_vector, nshl_vector, rcf_vector, fidfc, + dc0m, (fies > 0), fexdc, fwp, fxip + ); delete[] file_lines; delete[] variable_vector; return config; @@ -864,21 +868,9 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(const std::string& fil status = hdf_file->close(); delete hdf_file; conf = new ScattererConfiguration( - nsph, - configuration_count, - xi_vec, - nxi, - "XIV", - iog, - ros_vector, - nshl_vector, - rcf_vector, - _idfc, - dc0m, - (ies == 1), - _exdc, - _wp, - _xip + nsph, configuration_count, xi_vec, nxi, "XIV", + iog, ros_vector, nshl_vector, rcf_vector, _idfc, + dc0m, (ies == 1), _exdc, _wp, _xip ); delete[] xi_vec; } @@ -932,9 +924,13 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f for (int nsi = 0; nsi < nsh; nsi++) 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]; 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]; for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new dcomplex[dim3](); } // di loop @@ -995,7 +991,7 @@ void ScattererConfiguration::print() { for (int i = 0; i < _number_of_scales; i++) printf("\t%lg", _scale_vec[i]); printf("\t]\n"); printf("DC0M = [\n"); - for (int i = 0; i < _configurations; i++) { + for (int i = 0; i < _max_layers; i++) { printf(" [\n"); for (int j = 0; j < _number_of_spheres; j++) { printf(" ["); @@ -1079,7 +1075,7 @@ void ScattererConfiguration::write_hdf5(const std::string& file_name) { } 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](); int dc0_index = 0; for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) { @@ -1376,6 +1372,9 @@ bool ScattererConfiguration::operator ==(const ScattererConfiguration &other) { if (_idfc != other._idfc) { return false; } + if (_max_layers != other._max_layers) { + return false; + } if (_exdc != other._exdc) { return false; } @@ -1411,13 +1410,15 @@ bool ScattererConfiguration::operator ==(const ScattererConfiguration &other) { return false; } } // rj loop + } // ri loop + for (int li = 0; li < _max_layers; li++) { for (int dj = 0; dj < _number_of_spheres; dj++) { 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; } } // dk loop } // dj loop - } // ri loop + } return true; }