diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 4e0d1e9a77b0dded1cda3a5316153db3e6909458..f87c318c7eb10f050f9b383cc27a8c92f12cf7cb 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -143,7 +143,7 @@ building_stage: - cat /etc/os-release - cd build - echo "Configuring with default compilers (MAGMA disabled)..." - - ./configure --without-magma + - ./configure --without-magma --disable-offload - make clean - echo "Building the default configuration..." - make -j diff --git a/build/Makefile.in b/build/Makefile.in index 01a51f5d476101c0e49dd0f87a08ae56e7d04e58..21198afa6b69ced519eecafa68e1fc1dd6beaf08 100644 --- a/build/Makefile.in +++ b/build/Makefile.in @@ -502,6 +502,7 @@ OBJDUMP = @OBJDUMP@ OBJEXT = @OBJEXT@ OFFLOADFLAGS = @OFFLOADFLAGS@ OMPFLAGS = @OMPFLAGS@ +OPTFLAGS = @OPTFLAGS@ OTOOL = @OTOOL@ OTOOL64 = @OTOOL64@ PACKAGE = @PACKAGE@ diff --git a/build/configure b/build/configure index fbc162d9ee626804521110f09db57e80d901857f..44ef1e6879e8804001c91c1fc199c2ab88307817 100755 --- a/build/configure +++ b/build/configure @@ -663,6 +663,7 @@ MAGMALDFLAGS MAGMAFLAGS LAPACKLDFLAGS LAPACKFLAGS +OPTFLAGS OMPFLAGS OFFLOADFLAGS BUILDFORTRAN_FALSE @@ -816,6 +817,7 @@ enable_libtool_lock enable_fortran enable_offload enable_openmp +enable_optimize with_lapack with_magma enable_nvtx @@ -1481,6 +1483,7 @@ Optional Features: --enable-offload enable target offloading (requires g++ version >= 13) [default=auto] --enable-openmp enable OpneMP multi-threading [default=yes] + --enable-optimize=LEVEL use optimization level LEVEL [default=3] --enable-nvtx use NVTX profiling [default=no] Optional Packages: @@ -25062,6 +25065,34 @@ esac fi +# Check whether --enable-optimize was given. +if test ${enable_optimize+y} +then : + enableval=$enable_optimize; + if test "x$enableval" = "x0"; then + OPTFLAGS="-O0" + + elif test "x$enableval" = "x1"; then + OPTFLAGS="-O1" + + elif test "x$enableval" = "x2"; then + OPTFLAGS="-O2" + + elif test "x$enableval" = "x3"; then + OPTFLAGS="-O3" + + else + as_fn_error $? "Invalid optimization flag!" "$LINENO" 5 + fi + +else case e in #( + e) OPTFLAGS="-O3" + + ;; +esac +fi + + # Check whether --with-lapack was given. if test ${with_lapack+y} @@ -25332,7 +25363,7 @@ else case e in #( ;; esac fi -CXXFLAGS="$CLANGFLAGS -O3 -ggdb $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS" +CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS" SUBDIRS="cluster libnptm sphere testing trapping" # Generate the output diff --git a/build/configure.ac b/build/configure.ac index 28f20c9ed4b3af85488b4e9448999a1525fcba84..0238f91edac77c5f9a6a61c1f3db6b3a86ec5fac 100644 --- a/build/configure.ac +++ b/build/configure.ac @@ -322,6 +322,25 @@ AC_ARG_ENABLE( ] ) +AC_ARG_ENABLE( + [optimize], + [AS_HELP_STRING([--enable-optimize=LEVEL], [use optimization level LEVEL [default=3]])], + [ + if test "x$enableval" = "x0"; then + AC_SUBST([OPTFLAGS], ["-O0"]) + elif test "x$enableval" = "x1"; then + AC_SUBST([OPTFLAGS], ["-O1"]) + elif test "x$enableval" = "x2"; then + AC_SUBST([OPTFLAGS], ["-O2"]) + elif test "x$enableval" = "x3"; then + AC_SUBST([OPTFLAGS], ["-O3"]) + else + AC_MSG_ERROR([Invalid optimization flag!]) + fi + ], + [AC_SUBST([OPTFLAGS], ["-O3"])] +) + AC_ARG_WITH( [lapack], [AS_HELP_STRING([--with-lapack], [use LAPACK @<:@default=auto@:>@])], @@ -444,7 +463,7 @@ AS_IF( [AC_SUBST([OMPFLAGS], [""])], [AC_SUBST([OMPFLAGS], [$OMPFLAGS])] ) -CXXFLAGS="$CLANGFLAGS -O3 -ggdb $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS" +CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS" SUBDIRS="cluster libnptm sphere testing trapping" # Generate the output diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp index ad48fe8eeede3fcf1841896cbd6720a4aa04bc01..73c24667a7ce36a478f20795aab0a93849d514b9 100644 --- a/src/cluster/cluster.cpp +++ b/src/cluster/cluster.cpp @@ -681,6 +681,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf int jwtm = gconf->jwtm; np_int ndit = 2 * nsph * cid->c4->nlim; int isq, ibf; + int last_configuration; #ifdef USE_NVTX nvtxRangePush("Prepare matrix calculation"); @@ -710,6 +711,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf return jer; // break; // rewrite this to go to the end of the function, to free locally allocated variables and return jer } + last_configuration = 0; for (int i132 = 1; i132 <= nsph; i132++) { int iogi = cid->c1->iog[i132 - 1]; if (iogi != i132) { @@ -718,7 +720,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf cid->c1->rei[l123 - 1][i132 - 1] = cid->c1->rei[l123 - 1][iogi - 1]; } // l123 loop } else { - int nsh = cid->c1->nshl[i132 - 1]; + last_configuration++; + int nsh = cid->c1->nshl[last_configuration - 1]; int ici = (nsh + 1) / 2; if (idfc == 0) { for (int ic = 0; ic < ici; ic++) @@ -732,7 +735,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc; dme( cid->c4->li, i132, npnt, npntts, vkarg, exdc, exri, - cid->c1, cid->c2, jer, lcalc, cid->arg + cid->c1, cid->c2, jer, lcalc, cid->arg, last_configuration ); if (jer != 0) { sprintf(virtual_line, " STOP IN DME\n"); @@ -813,16 +816,18 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf sprintf(virtual_line, " SPHERES; LMX=LI\n"); output->append_line(virtual_line); } + last_configuration = 0; for (int i170 = 1; i170 <= nsph; i170++) { if (cid->c1->iog[i170 - 1] >= i170) { int i = i170 - 1; + last_configuration++; double albeds = cid->c1->sscs[i] / cid->c1->sexs[i]; cid->c1->sqscs[i] *= cid->sqsfi; cid->c1->sqabs[i] *= cid->sqsfi; cid->c1->sqexs[i] *= cid->sqsfi; sprintf(virtual_line, " SPHERE %2d\n", i170); output->append_line(virtual_line); - if (cid->c1->nshl[i] != 1) { + if (cid->c1->nshl[last_configuration - 1] != 1) { sprintf(virtual_line, " SIZE=%15.7lE\n", cid->c2->vsz[i]); output->append_line(virtual_line); } else { // label 162 diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h index 66a3f6307c8465a3b49d1b5f8881bfa85a8cb466..bfcf4a720f44ce5e595d05581c60c25f02936f05 100644 --- a/src/include/sph_subs.h +++ b/src/include/sph_subs.h @@ -112,10 +112,11 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2); * \param jer: `int &` Reference to integer error code variable. * \param lcalc: `int &` Reference to integer variable recording the maximum expansion order accounted for. * \param arg: `complex double &` + * \param last_conf: `int` Last sphere configuration (used by CLUSTER) */ void dme( int li, int i, int npnt, int npntts, double vk, double exdc, double exri, - C1 *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg + C1 *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg, int last_conf=0 ); /*! \brief Bessel function calculation control parameters. diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp index 0ed4ed745f7df55d13005138e4b9caff35b46f82..6f36e84ff12defdd3cb078f333e31c23ea2be0fa 100644 --- a/src/libnptm/Configuration.cpp +++ b/src/libnptm/Configuration.cpp @@ -429,7 +429,7 @@ ScattererConfiguration::ScattererConfiguration(const mixMPI *mpidata) MPI_Bcast(&_xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); _iog_vec = new int[_number_of_spheres](); MPI_Bcast(_iog_vec, _number_of_spheres, MPI_INT, 0, MPI_COMM_WORLD); - _radii_of_spheres = new double[_number_of_spheres](); + _radii_of_spheres = new double[_configurations](); MPI_Bcast(_radii_of_spheres, _configurations, MPI_DOUBLE, 0, MPI_COMM_WORLD); _nshl_vec = new int[_configurations](); MPI_Bcast(_nshl_vec, _configurations, MPI_INT, 0, MPI_COMM_WORLD); @@ -513,6 +513,7 @@ ScattererConfiguration* ScattererConfiguration::from_binary(const std::string& f ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& dedfb_file_name) { int num_lines = 0; int last_read_line = 0; + int last_configuration; string *file_lines; regex re; smatch m; @@ -687,32 +688,34 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de } } int index = 0; - double *ros_vector = new double[fnsph](); - int *nshl_vector = new int[fnsph](); - double **rcf_vector = new double*[fnsph]; + double *ros_vector = new double[configurations](); + int *nshl_vector = new int[configurations](); + double **rcf_vector = new double*[configurations]; + last_configuration = 0; for (int i113 = 1; i113 <= fnsph; i113++) { if (iog_vector[i113 - 1] < i113) continue; + else last_configuration++; str_target = file_lines[++last_read_line]; re = regex("[0-9]+"); regex_search(str_target, m, re); - nshl_vector[i113 - 1] = stoi(m.str()); + nshl_vector[last_configuration - 1] = stoi(m.str()); str_target = m.suffix().str(); re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?"); regex_search(str_target, m, re); string str_number = m.str(); str_number = regex_replace(str_number, regex("D"), "e"); str_number = regex_replace(str_number, regex("d"), "e"); - ros_vector[i113 - 1] = stod(str_number); - int nsh = nshl_vector[i113 - 1]; + ros_vector[last_configuration - 1] = stod(str_number); + int nsh = nshl_vector[last_configuration - 1]; if (i113 == 1) nsh += fies; - rcf_vector[i113 - 1] = new double[nsh](); + rcf_vector[last_configuration - 1] = new double[nsh](); for (int ns112 = 0; ns112 < nsh; ns112++) { str_target = file_lines[++last_read_line]; regex_search(str_target, m, re); str_number = m.str(); str_number = regex_replace(str_number, regex("D"), "e"); str_number = regex_replace(str_number, regex("d"), "e"); - rcf_vector[i113 - 1][ns112] = stod(str_number); + rcf_vector[last_configuration - 1][ns112] = stod(str_number); } // ns112 loop } // i113 loop @@ -725,9 +728,11 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de } // di loop for (int jxi468 = 1; jxi468 <= fnxi; jxi468++) { if (fidfc != 0 && jxi468 > 1) continue; // jxi468 loop + last_configuration = 0; for (int i162 = 1; i162 <= fnsph; i162++) { if (iog_vector[i162 - 1] >= i162) { - int nsh = nshl_vector[i162 - 1]; + last_configuration++; + int nsh = nshl_vector[last_configuration - 1]; int ici = (nsh + 1) / 2; if (i162 == 1) ici += fies; for (int ic157 = 0; ic157 < ici; ic157++) { @@ -781,6 +786,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(const std::string& fil string str_name, str_type; if (status == 0) { int nsph, ies; + int last_configuration; int *iog; double _exdc, _wp, _xip; int _idfc, nxi; @@ -804,22 +810,24 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(const std::string& fil if (iog[ci] < ci + 1) continue; configuration_count++; } - nshl_vector = new int[nsph](); - ros_vector = new double[nsph](); - rcf_vector = new double*[nsph]; - for (int i115 = 1; i115 <= nsph; i115++) { + nshl_vector = new int[configuration_count](); + ros_vector = new double[configuration_count](); + rcf_vector = new double*[configuration_count]; + last_configuration = 0; + for (int i115 = 1; i115 <= configuration_count; i115++) { if (iog[i115 - 1] < i115) continue; - str_name = "NSHL_" + to_string(iog[i115 - 1]); + else last_configuration++; + str_name = "NSHL_" + to_string(last_configuration); str_type = "INT32_(1)"; - status = hdf_file->read(str_name, str_type, (nshl_vector + i115 - 1)); - str_name = "ROS_" + to_string(iog[i115 - 1]); + status = hdf_file->read(str_name, str_type, (nshl_vector + last_configuration - 1)); + str_name = "ROS_" + to_string(last_configuration); str_type = "FLOAT64_(1)"; - status = hdf_file->read(str_name, str_type, (ros_vector + i115 - 1)); - int nsh = nshl_vector[i115 - 1]; - rcf_vector[i115 - 1] = new double[nsh](); - str_name = "RCF_" + to_string(iog[i115 - 1]); + status = hdf_file->read(str_name, str_type, (ros_vector + last_configuration - 1)); + int nsh = nshl_vector[last_configuration - 1]; + rcf_vector[last_configuration - 1] = new double[nsh](); + str_name = "RCF_" + to_string(last_configuration); str_type = "FLOAT64_(" + to_string(nsh) + ")"; - status = hdf_file->read(str_name, str_type, (rcf_vector[i115 - 1])); + status = hdf_file->read(str_name, str_type, (rcf_vector[last_configuration - 1])); } xi_vec = new double[nxi](); str_name = "XIVEC"; @@ -874,6 +882,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f int _nsph; int *_iog_vec; int _ies; + int last_configuration; double _exdc, _wp, _xip; int _idfc, _nxi; int *_nshl_vec; @@ -901,18 +910,20 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f for (int i = 1; i <= _nsph; i++) { if (_iog_vec[i - 1] >= i) configurations++; } - _nshl_vec = new int[_nsph](); - _ros_vec = new double[_nsph](); - _rcf_vec = new double*[_nsph]; + _nshl_vec = new int[configurations](); + _ros_vec = new double[configurations](); + _rcf_vec = new double*[configurations]; + last_configuration = 0; for (int i115 = 1; i115 <= _nsph; i115++) { if (_iog_vec[i115 - 1] < i115) continue; - input.read(reinterpret_cast(&(_nshl_vec[i115 - 1])), sizeof(int)); - input.read(reinterpret_cast(&(_ros_vec[i115 - 1])), sizeof(double)); - int nsh = _nshl_vec[i115 - 1]; + else last_configuration++; + input.read(reinterpret_cast(&(_nshl_vec[last_configuration - 1])), sizeof(int)); + input.read(reinterpret_cast(&(_ros_vec[last_configuration - 1])), sizeof(double)); + int nsh = _nshl_vec[last_configuration - 1]; if (i115 == 1) nsh += _ies; - _rcf_vec[i115 - 1] = new double[nsh](); + _rcf_vec[last_configuration - 1] = new double[nsh](); for (int nsi = 0; nsi < nsh; nsi++) - input.read(reinterpret_cast(&(_rcf_vec[i115 - 1][nsi])), sizeof(double)); + input.read(reinterpret_cast(&(_rcf_vec[last_configuration - 1][nsi])), sizeof(double)); } _dc0m = new dcomplex**[configurations]; int dim3 = (_idfc == 0) ? _nxi : 1; @@ -922,9 +933,11 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f } // di loop for (int jxi468 = 1; jxi468 <= _nxi; jxi468++) { if (_idfc != 0 && jxi468 > 1) continue; + last_configuration = 0; for (int i162 = 1; i162 <= _nsph; i162++) { if (_iog_vec[i162 - 1] < i162) continue; - int nsh = _nshl_vec[i162 - 1]; + else last_configuration++; + int nsh = _nshl_vec[last_configuration - 1]; int ici = (nsh + 1) / 2; if (i162 == 1) ici = ici + _ies; for (int i157 = 0; i157 < ici; i157++) { @@ -957,26 +970,6 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f return conf; } -/* - double ScattererConfiguration::get_param(const std::string& param_name) { - double value; - if (param_name.compare("number_of_spheres") == 0) value = (double)number_of_spheres; - else if (param_name.compare("nsph") == 0) value = (double)number_of_spheres; - else if (param_name.compare("configurations") == 0) value = (double)configurations; - else if (param_name.compare("number_of_scales") == 0) value = (double)number_of_scales; - else if (param_name.compare("nxi") == 0) value = (double)number_of_scales; - else if (param_name.compare("idfc") == 0) value = (double)idfc; - else if (param_name.compare("exdc") == 0) value = exdc; - else if (param_name.compare("wp") == 0) value = wp; - else if (param_name.compare("xip") == 0) value = xip; - else { - string message = "unrecognized parameter \"" + param_name + "\""; - throw UnrecognizedParameterException(message); - } - return value; - } -*/ - void ScattererConfiguration::print() { int ies = (_use_external_sphere)? 1 : 0; printf("### CONFIGURATION DATA ###\n"); @@ -1035,6 +1028,7 @@ void ScattererConfiguration::write_binary(const std::string& file_name, const st void ScattererConfiguration::write_hdf5(const std::string& file_name) { int ies = (use_external_sphere)? 1 : 0; + int last_configuration; List *rec_name_list = new List(1); List *rec_type_list = new List(1); List *rec_ptr_list = new List(1); @@ -1068,23 +1062,25 @@ void ScattererConfiguration::write_hdf5(const std::string& file_name) { str_type = "FLOAT64_(" + to_string(number_of_scales) + ")"; rec_type_list->append(str_type); rec_ptr_list->append(_scale_vec); + last_configuration = 0; for (int i115 = 1; i115 <= _number_of_spheres; i115++) { if (_iog_vec[i115 - 1] < i115) continue; - str_name = "NSHL_" + to_string(i115); + else last_configuration++; + str_name = "NSHL_" + to_string(last_configuration); rec_name_list->append(str_name); rec_type_list->append("INT32_(1)"); - rec_ptr_list->append(&(_nshl_vec[i115 - 1])); // was not from IOG - str_name = "ROS_" + to_string(i115); + rec_ptr_list->append(&(_nshl_vec[last_configuration - 1])); // was not from IOG + str_name = "ROS_" + to_string(last_configuration); rec_name_list->append(str_name); rec_type_list->append("FLOAT64_(1)"); - rec_ptr_list->append(&(_radii_of_spheres[i115 - 1])); // was not from IOG - int nsh = _nshl_vec[i115 - 1]; // was not from IOG + rec_ptr_list->append(&(_radii_of_spheres[last_configuration - 1])); // was not from IOG + int nsh = _nshl_vec[last_configuration - 1]; // was not from IOG if (i115 == 1) nsh += ies; - str_name = "RCF_" + to_string(i115); // was not from IOG + str_name = "RCF_" + to_string(last_configuration); // was not from IOG str_type = "FLOAT64_(" + to_string(nsh) + ")"; rec_name_list->append(str_name); rec_type_list->append(str_type); - rec_ptr_list->append(&(_rcf[i115 - 1][0])); // was not from IOG + rec_ptr_list->append(&(_rcf[last_configuration - 1][0])); // was not from IOG } int dim3 = (idfc == 0) ? _number_of_scales : 1; @@ -1093,9 +1089,11 @@ void ScattererConfiguration::write_hdf5(const std::string& file_name) { int dc0_index = 0; for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) { if (_idfc != 0 && jxi468 > 1) continue; + last_configuration = 0; for (int i162 = 1; i162 <= _number_of_spheres; i162++) { if (_iog_vec[i162 - 1] < i162) continue; - int nsh = _nshl_vec[i162 - 1]; // was not from IOG + else last_configuration++; + int nsh = _nshl_vec[last_configuration - 1]; // was not from IOG int ici = (nsh + 1) / 2; if (i162 == 1) ici = ici + ies; for (int i157 = 0; i157 < ici; i157++) { @@ -1138,6 +1136,7 @@ void ScattererConfiguration::write_hdf5(const std::string& file_name) { void ScattererConfiguration::write_legacy(const std::string& file_name) { fstream output; int ies = (_use_external_sphere)? 1 : 0; + int last_configuration; output.open(file_name.c_str(), ios::out | ios::binary); output.write(reinterpret_cast(&_number_of_spheres), sizeof(int)); output.write(reinterpret_cast(&ies), sizeof(int)); @@ -1150,20 +1149,24 @@ void ScattererConfiguration::write_legacy(const std::string& file_name) { output.write(reinterpret_cast(&_number_of_scales), sizeof(int)); for (int i = 0; i < _number_of_scales; i++) output.write(reinterpret_cast(&(_scale_vec[i])), sizeof(double)); + last_configuration = 0; for (int i115 = 1; i115 <= _number_of_spheres; i115++) { if (_iog_vec[i115 - 1] < i115) continue; - output.write(reinterpret_cast(&(_nshl_vec[i115 - 1])), sizeof(int)); // was not from IOG - output.write(reinterpret_cast(&(_radii_of_spheres[i115 - 1])), sizeof(double)); // was not from IOG - int nsh = _nshl_vec[i115 - 1]; // was not from IOG + else last_configuration++; + output.write(reinterpret_cast(&(_nshl_vec[last_configuration - 1])), sizeof(int)); // was not from IOG + output.write(reinterpret_cast(&(_radii_of_spheres[last_configuration - 1])), sizeof(double)); // was not from IOG + int nsh = _nshl_vec[last_configuration - 1]; // was not from IOG if (i115 == 1) nsh += ies; for (int nsi = 0; nsi < nsh; nsi++) - output.write(reinterpret_cast(&(_rcf[i115 - 1][nsi])), sizeof(double)); // was not from IOG + output.write(reinterpret_cast(&(_rcf[last_configuration - 1][nsi])), sizeof(double)); // was not from IOG } for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) { if (idfc != 0 && jxi468 > 1) continue; + last_configuration = 0; for (int i162 = 1; i162 <= _number_of_spheres; i162++) { if (_iog_vec[i162 - 1] < i162) continue; - int nsh = _nshl_vec[i162 - 1]; // was not from IOG + else last_configuration++; + int nsh = _nshl_vec[last_configuration - 1]; // was not from IOG int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here? if (i162 == 1) ici = ici + ies; for (int i157 = 0; i157 < ici; i157++) { @@ -1186,6 +1189,7 @@ void ScattererConfiguration::write_formatted(const std::string& file_name) { const double two_pi = acos(0.0) * 4.0; double *xi_vec; int ici; + int last_configuration; int ies = (use_external_sphere)? 1: 0; FILE *output = fopen(file_name.c_str(), "w"); int scale_type = -1; @@ -1332,9 +1336,11 @@ void ScattererConfiguration::write_formatted(const std::string& file_name) { } if (_idfc != 0) { fprintf(output, " DIELECTRIC CONSTANTS\n"); + last_configuration = 0; for (int i473 = 1; i473 <= _number_of_spheres; i473++) { if (_iog_vec[i473 - 1] != i473) continue; - ici = (_nshl_vec[i473 - 1] + 1) / 2; + else last_configuration++; + ici = (_nshl_vec[last_configuration - 1] + 1) / 2; if (i473 == 1) ici += ies; fprintf(output, " SPHERE N.%4d\n", i473); for (int ic472 = 0; ic472 < ici; ic472++) { @@ -1345,9 +1351,11 @@ void ScattererConfiguration::write_formatted(const std::string& file_name) { } } else { fprintf(output, " DIELECTRIC FUNCTIONS\n"); + last_configuration = 0; for (int i478 = 1; i478 <= _number_of_spheres; i478++) { if (_iog_vec[i478 - 1] != i478) continue; - ici = (_nshl_vec[i478 - 1] + 1) / 2; + last_configuration++; + ici = (_nshl_vec[last_configuration - 1] + 1) / 2; if (i478 == 1) ici += ies; fprintf(output, " SPHERE N.%4d\n", i478); for (int ic477 = 1; ic477 <= ici; ic477++) { @@ -1418,4 +1426,3 @@ bool ScattererConfiguration::operator ==(const ScattererConfiguration &other) { } // ri loop return true; } - diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp index f10d3f788ac77b67adfc917f11d15895ace741a8..1712036ade1d3cb770ce91d1f03d29b0e5565925 100644 --- a/src/libnptm/clu_subs.cpp +++ b/src/libnptm/clu_subs.cpp @@ -2439,16 +2439,19 @@ void scr2( void str(ScattererConfiguration *sconf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) { dcomplex *ylm; const double pi = acos(-1.0); + int last_configuration; c3->gcs = 0.0; double gcss = 0.0; + last_configuration = 0; for (int i18 = 1; i18 <= c4->nsph; i18++) { int iogi = c1->iog[i18 - 1]; if (iogi >= i18) { - gcss = pi * c1->ros[i18 - 1] * c1->ros[i18 - 1]; + last_configuration++; + gcss = pi * c1->ros[last_configuration - 1] * c1->ros[last_configuration - 1]; c1->gcsv[i18 - 1] = gcss; - int nsh = c1->nshl[i18 - 1]; + int nsh = c1->nshl[last_configuration - 1]; for (int j16 = 1; j16 <= nsh; j16++) { - c1->rc[i18 - 1][j16 - 1] = sconf->get_rcf(i18 - 1, j16 - 1) * c1->ros[i18 - 1]; + c1->rc[last_configuration - 1][j16 - 1] = sconf->get_rcf(last_configuration - 1, j16 - 1) * c1->ros[last_configuration - 1]; } // j16 loop } c3->gcs += gcss; diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp index 1b5a8e6aa2ea61ce52c8b893d7f3666cfecc0726..4a63948008c71a3362bc6b25d44818a25e7efa0f 100644 --- a/src/libnptm/sph_subs.cpp +++ b/src/libnptm/sph_subs.cpp @@ -197,20 +197,6 @@ double cg1(int lmpml, int mu, int l, int m) { return result; } -/* -#ifdef USE_TARGET_OFFLOAD -#pragma omp begin declare target device_type(any) -#endif -dcomplex dconjg(dcomplex z) { - double zreal = real(z); - double zimag = imag(z); - return (zreal - zimag * I); -} -#ifdef USE_TARGET_OFFLOAD -#pragma omp end declare target -#endif -*/ - void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) { const double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1]; const double half_step = 0.5 * dif / npntmo; @@ -232,7 +218,7 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) { void dme( int li, int i, int npnt, int npntts, double vk, double exdc, double exri, - C1 *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg + C1 *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg, int last_conf ) { const int lipo = li + 1; const int lipt = li + 2; @@ -243,13 +229,14 @@ void dme( dcomplex dfbi, dfb, dfn, ccna, ccnb, ccnc, ccnd; dcomplex y1, dy1, y2, dy2, arin, cri; const dcomplex uim = 1.0 * I; + int sph_index = (last_conf == 0) ? i : last_conf; jer = 0; int nstp = npnt - 1; int nstpts = npntts - 1; - double sz = vk * c1->ros[i - 1]; + double sz = vk * c1->ros[sph_index - 1]; c2->vsz[i - 1] = sz; - double vkr1 = vk * c1->rc[i - 1][0]; - int nsh = c1->nshl[i - 1]; + double vkr1 = vk * c1->rc[sph_index - 1][0]; + int nsh = c1->nshl[sph_index - 1]; c2->vkt[i - 1] = csqrt(c2->dc0[0]); arg = vkr1 * c2->vkt[i - 1]; arin = arg; diff --git a/src/testing/test_TEDF.cpp b/src/testing/test_TEDF.cpp index 1b3632e886fc93240f2a5a0068d79fe6f4aa4b95..0c50ee74eff9468e0095a1d698d2fb0759df63b4 100644 --- a/src/testing/test_TEDF.cpp +++ b/src/testing/test_TEDF.cpp @@ -74,9 +74,21 @@ int main(int argc, char **argv) { hdf5_file = string(argv[3]); } ScattererConfiguration *a, *b, *c; - a = ScattererConfiguration::from_dedfb(dedfb_file); - b = ScattererConfiguration::from_binary(legacy_file); - c = ScattererConfiguration::from_binary(hdf5_file, "HDF5"); + try { a = ScattererConfiguration::from_dedfb(dedfb_file); } + catch (...) { + printf("Failed to open legacy configuration file.\n"); + return 2; + } + try { b = ScattererConfiguration::from_binary(legacy_file); } + catch (...) { + printf("Failed to open legacy binary file.\n"); + return 20; + } + try { c = ScattererConfiguration::from_binary(hdf5_file, "HDF5"); } + catch (...) { + printf("Failed to open HDF5 configuration file.\n"); + return 200; + } if (*a == *b) printf("Configuration objects a and b are equal.\n"); else { printf("Configuration objects a and b are different.\n");