diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 39e36dad06329deff757d046ef360ef211c575e1..a77e21b6db9cb73ffab0a0709391b0f9849ef7b6 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -78,13 +78,13 @@ compatibility_stage: - CXX=g++-14 FC=gfortran-14 ./configure - make wipe - make -j - - echo "Running make with refinement with gnu compilers version 14..." + - echo "Running make with gnu compilers version 14..." - cd .. - rm -rf build_gnu14 - mkdir build_gnu14_refine - cd build_gnu14_refine - cp -r ../build/* . - - CXX=g++-14 FC=gfortran-14 ./configure --enable-refinement + - CXX=g++-14 FC=gfortran-14 ./configure - make wipe - make -j #- echo "Running make with flang version 16 and clang version 16..." @@ -173,7 +173,7 @@ building_stage: - cat /etc/os-release - cd build - echo "Configuring with default compilers (MAGMA disabled)..." - - ./configure --without-magma --without-cublas --disable-offload --enable-refinement --enable-shared + - ./configure --without-magma --without-cublas --disable-offload --enable-shared - make wipe - echo "Building the default configuration..." - make -j diff --git a/COPYING b/COPYING index 281d399f13dfd23087acc56303dd38d68162587c..92995c440ab669f6b4d1df7ee4cac2f809b03925 100644 --- a/COPYING +++ b/COPYING @@ -617,3 +617,5 @@ reviewing courts shall apply local law that most closely approximates an absolute waiver of all civil liability in connection with the Program, unless a warranty or assumption of liability accompanies a copy of the Program in return for a fee. + + diff --git a/build/configure.sh b/build/configure.sh index 161ddb933460e7f32952856381a652c2a5c86c25..3ab2082371d505b2acbff59954eba026f356e5bc 100755 --- a/build/configure.sh +++ b/build/configure.sh @@ -21,21 +21,20 @@ NVTXFLAGS="" OMPMODE="auto" OFFLOAD="auto" OFFLOADFLAGS="" -REFINEFLAGS="" # End of default configuration settings # Function declarations function guess_cxx { # Guess the name of the C++ compiler - result=$(which mpicxx) + result=$(which mpicxx 2>/dev/null) if [ "x$result" = "x" ]; then - result=$(which g++) + result=$(which g++ 2>/dev/null) fi if [ "x$result" = "x" ]; then - result=$(which clang) + result=$(which clang 2>/dev/null) fi if [ "x$result" = "x" ]; then - result=$(which icpx) + result=$(which icpx 2>/dev/null) fi if [ "x$result" = "x" ]; then result="none" @@ -45,21 +44,21 @@ function guess_cxx { function guess_fc { # Guess the name of the FORTRAN compiler - result=$(which mpif90) + result=$(which mpif90 2>/dev/null) if [ "x$result" = "x" ]; then - result=$(which gfortran) + result=$(which gfortran 2>/dev/null) fi if [ "x$result" = "x" ]; then - result=$(which f77) + result=$(which f77 2>/dev/null) fi if [ "x$result" = "x" ]; then - result=$(which flang) + result=$(which flang 2>/dev/null) fi if [ "x$result" = "x" ]; then - result=$(which flang-new) + result=$(which flang-new 2>/dev/null) fi if [ "x$result" = "x" ]; then - result=$(which ifx) + result=$(which ifx 2>/dev/null) fi if [ "x$result" = "x" ]; then result="none" @@ -87,9 +86,6 @@ function print_help { echo "--enable-nvtx Enable NVTX profiling tools. " echo "--disable-nvtx Disable NVTX profiling tools (DEFAULT). " echo "--enable-optimize=OPT Use OPT-level compiler optimization. " - echo "--enable-refinement Use iterative refinement of matrix inversion. " - echo "--disable-refinement Do not iterate refinement of matrix inversion " - echo " (DEFAULT). " echo "--enable-shared Use shared libraries (default is static). " echo "--disable-shared Use static libraries (DEFAULT). " echo "--help Print this help and exit. " @@ -198,10 +194,6 @@ do fi FC_OPT=$opt_level CXX_OPT=$opt_level - elif [ "x$cut_arg" = "x--enable-refinement" ]; then - REFINEFLAGS=" -DUSE_REFINEMENT" - elif [ "x$cut_arg" = "x--disable-refinement" ]; then - REFINEFLAGS="" elif [ "x$cut_arg" = "x--enable-shared" ]; then LIBMODE="shared" elif [ "x$cut_arg" = "x--disable-shared" ]; then @@ -634,6 +626,11 @@ if [ "x$CUBLAS" != "xno" ]; then cuda_pkg=$(for i in "${pkg_array[@]}"; do echo "$i" | cut --delimiter=" " -f1; done | grep cublas) CUDAFLAGS=$(pkg-config --cflags ${cuda_pkg}) CUDALDFLAGS=$(pkg-config --libs ${cuda_pkg}) + else + # CUBLAS not detected + CUBLAS="no" + CUDAFLAGS="" + CUDALDFLAGS="" fi # end of CUBLAS runtime decision tree declare -a pkg_array=$(pkg-config --list-all | grep cudart) for i in "${pkg_array[@]}"; do echo "$i" | cut --delimiter=" " -f1; done | grep cudart > /dev/null @@ -647,12 +644,16 @@ if [ "x$CUBLAS" != "xno" ]; then echo $CUDALDFLAGS | grep cublas > /dev/null cudart_check=$? if [ "x${cudart_check}" != "x0" ]; then - CUDALDFLAGS="$CUDALDFLAGS -lcublas" + if [ "x${CUBLAS}" != "xno" ]; then + CUDALDFLAGS="$CUDALDFLAGS -lcublas" + fi fi echo $CUDALDFLAGS | grep cudart > /dev/null cudart_check=$? if [ "x${cudart_check}" != "x0" ]; then - CUDALDFLAGS="$CUDALDFLAGS -lcudart" + if [ "x${CUBLAS}" != "xno" ]; then + CUDALDFLAGS="$CUDALDFLAGS -lcudart" + fi fi else # pkg-config is not available @@ -813,7 +814,7 @@ EOF fi if [ "x$result" = "x0" ]; then echo "yes" - OFFLOADFLAGS=" -fcf-protection=check -foffload=nvptx-none=\"-O${CXX_OPT}${CXX_DBG} -fcf-protection=check -fopt-info -lm -latomic -gomp\"" + OFFLOADFLAGS=" -DUSE_TARGET_OFFLOAD -fcf-protection=check -foffload=nvptx-none=\"-O${CXX_OPT}${CXX_DBG} -fcf-protection=check -fopt-info -lm -latomic -lgomp\"" if [ "x${OMPFLAGS}" = "x" ]; then OFFLOADFLAGS="${OFFLOADFLAGS} -fopenmp" fi @@ -826,7 +827,7 @@ else fi # End of offload checks if [ "x$CXXFLAGS" = "x" ]; then - CXXFLAGS="-O${CXX_OPT}${CXX_DBG}${CLANGFLAGS}${INCLUDEFLAGS}${HDF5FLAGS}${OMPFLAGS}${MPIFLAGS}${LAPACKFLAGS}${CUBLASFLAGS}${MAGMAFLAGS}${REFINEFLAGS}${DEBUGFLAGS}${OFFLOADFLAGS}${NVTXFLAGS}" + CXXFLAGS="-O${CXX_OPT}${CXX_DBG}${CLANGFLAGS}${INCLUDEFLAGS}${HDF5FLAGS}${OMPFLAGS}${MPIFLAGS}${LAPACKFLAGS}${CUBLASFLAGS}${MAGMAFLAGS}${DEBUGFLAGS}${OFFLOADFLAGS}${NVTXFLAGS}" fi if [ "x$CXXLDFLAGS" = "x" ]; then if [ "x$LIBMODE" = "xstatic" ]; then @@ -885,11 +886,6 @@ if [ "x${NVTXFLAGS}" = "x" ]; then else echo "INFO: NVTX profiling is enabled." fi -if [ "x${REFINEFLAGS}" = "x" ]; then - echo "INFO: iterative matrix inversion refinement is disabled." -else - echo "INFO: iterative matrix inversion refinement is enabled." -fi if [ "x${OFFLOADFLAGS}" = "x" ]; then echo "INFO: GPU offload through OpenMP is disabled." else diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp index f9e43625f7c346bf32dedc9661d734002346f65b..68892d7a5db66d97d94fd0a01f6958cd6e19ade7 100644 --- a/src/cluster/cluster.cpp +++ b/src/cluster/cluster.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -757,9 +757,15 @@ int cluster_jxi488_cycle( + ").\n"; logger->log(message, LOG_WARN); } else if (recommended_li < cid->c1->li) { - message = "INFO: lowering internal order from " + to_string(cid->c1->li) + " to " - + to_string(recommended_li) + " for scale iteration " + to_string(jxi488) + ".\n"; - logger->log(message, LOG_INFO); + if (gconf->dyn_order_flag > 0) { + message = "INFO: lowering internal order from " + to_string(cid->c1->li) + " to " + + to_string(recommended_li) + " for scale iteration " + to_string(jxi488) + ".\n"; + logger->log(message, LOG_INFO); + } else { + message = "WARNING: internal order " + to_string(cid->c1->li) + " too high for scale iteration " + + to_string(jxi488) + ".\n"; + logger->log(message, LOG_WARN); + } } if (recommended_le > cid->c1->le) { message = "WARNING: external order " + to_string(cid->c1->le) + " for scale iteration " @@ -767,17 +773,25 @@ int cluster_jxi488_cycle( + ").\n"; logger->log(message, LOG_WARN); } else if (recommended_le < cid->c1->le) { - message = "INFO: lowering external order from " + to_string(cid->c1->le) + " to " - + to_string(recommended_le) + " for scale iteration " + to_string(jxi488) + ".\n"; - logger->log(message, LOG_INFO); + if (gconf->dyn_order_flag > 0) { + message = "INFO: lowering external order from " + to_string(cid->c1->le) + " to " + + to_string(recommended_le) + " for scale iteration " + to_string(jxi488) + ".\n"; + logger->log(message, LOG_INFO); + } else { + message = "INFO: external order " + to_string(cid->c1->le) + " too high for scale iteration " + + to_string(jxi488) + ".\n"; + logger->log(message, LOG_WARN); + } } if (recommended_li < max_li || recommended_le < max_le) { - int new_li = (recommended_li < max_li) ? recommended_li : max_li; - int new_le = (recommended_le < max_le) ? recommended_le : max_le; - cid->update_orders(sconf->_rcf, new_li, new_le); - is_first_scale = true; - jaw = 1; - cid->refinemode = 2; + if (gconf->dyn_order_flag > 0) { + int new_li = (recommended_li < max_li) ? recommended_li : max_li; + int new_le = (recommended_le < max_le) ? recommended_le : max_le; + cid->update_orders(sconf->_rcf, new_li, new_le); + is_first_scale = true; + jaw = 1; + cid->refinemode = 2; + } } } int li = cid->c1->li; @@ -883,16 +897,18 @@ int cluster_jxi488_cycle( double actualaccuracy = cid->accuracygoal; invert_matrix(cid->am, ndit, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, output_path, jxi488, mxndm, cid->proc_device); // in principle, we should check whether the returned actualaccuracy is indeed lower than the accuracygoal, and do something about it if not -#ifdef USE_REFINEMENT - if (cid->refinemode==2) { - message = "INFO: calibration obtained accuracy " + to_string(actualaccuracy) + " (" + to_string(cid->accuracygoal) + " requested) in " + to_string(cid->maxrefiters) + " refinement iterations\n"; - logger->log(message); - if (actualaccuracy > 1e-1) { - printf("Accuracy worse than 0.1, stopping"); - exit(1); + if (gconf->refine_flag > 0) { + if (cid->refinemode==2) { + message = "DEBUG: iterative refinement enabled at run-time.\n"; + logger->log(message, LOG_DEBG); + message = "INFO: calibration obtained accuracy " + to_string(actualaccuracy) + " (" + to_string(cid->accuracygoal) + " requested) in " + to_string(cid->maxrefiters) + " refinement iterations\n"; + logger->log(message); + if (actualaccuracy > 1e-1) { + printf("Accuracy worse than 0.1, stopping"); + exit(1); + } } } -#endif cid->refinemode = 0; #ifdef DEBUG_AM VirtualAsciiFile *outam2 = new VirtualAsciiFile(); diff --git a/src/cluster/np_cluster.cpp b/src/cluster/np_cluster.cpp index 7beb785cf6003945889b2d9263bbf2fa23934942..4da183e7cedee60e01efe75fcd808c3af1b845a0 100644 --- a/src/cluster/np_cluster.cpp +++ b/src/cluster/np_cluster.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/Commons.h b/src/include/Commons.h index 54f062e5d35adc5a808193ed6b1441d23ecac927..b7245d03a2c58eb5f6c3bc459879cfd4d8f4ca2b 100644 --- a/src/include/Commons.h +++ b/src/include/Commons.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/Configuration.h b/src/include/Configuration.h index 4235f56d9d42cb4849380282f66b90871b5bf1cb..0e1113ac7be97b6287af6ed3a7455de7eadca711 100644 --- a/src/include/Configuration.h +++ b/src/include/Configuration.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -109,54 +109,62 @@ protected: double *_sph_y; //! \brief Vector of spherical components Z coordinates. double *_sph_z; + //! \brief Flag for matrix inversion refinement. + short _refine_flag; + //! \brief Flag for dynamic order management. + short _dyn_order_flag; public: - //! \brief Read only view on number of spherical components. + //! \brief Read-only view on number of spherical components. const int& number_of_spheres = _number_of_spheres; - //! \brief Read only view on maximum field expansion order. + //! \brief Read-only view on maximum field expansion order. const int& l_max = _l_max; - //! \brief Read only view on maximum internal field expansion order. + //! \brief Read-only view on maximum internal field expansion order. const int& li = _li; - //! \brief Read only view on maximum external field expansion order. + //! \brief Read-only view on maximum external field expansion order. const int& le = _le; - //! \brief Read only view on maximum dimension of allocated matrix allowance (deprecated). + //! \brief Read-only view on maximum dimension of allocated matrix allowance (deprecated). const np_int& mxndm = _mxndm; - //! \brief Read only view on the intensity mode flag. + //! \brief Read-only view on the intensity mode flag. const int& iavm = _iavm; - //! \brief Read only view on incident field polarization status (0 - linear, 1 - circular). + //! \brief Read-only view on incident field polarization status (0 - linear, 1 - circular). const int& in_pol = _in_pol; - //! \brief Read only view on number of points for transition layer integration. + //! \brief Read-only view on number of points for transition layer integration. const int& npnt = _npnt; - //! \brief Read only view on number of points for non-transition layer integration. + //! \brief Read-only view on number of points for non-transition layer integration. const int& npntts = _npntts; - //! \brief Read only view on type of meridional plane definition. + //! \brief Read-only view on type of meridional plane definition. const int& isam = _isam; - //! \brief Read only view on scale index for T-matrix output. + //! \brief Read-only view on scale index for T-matrix output. const int& jwtm = _jwtm; - //! \brief Read only view on incident field initial azimuth. + //! \brief Read-only view on incident field initial azimuth. const double& in_theta_start = _in_theta_start; - //! \brief Read only view on incident field azimuth step. + //! \brief Read-only view on incident field azimuth step. const double& in_theta_step = _in_theta_step; - //! \brief Read only view on incident field final azimuth. + //! \brief Read-only view on incident field final azimuth. const double& in_theta_end = _in_theta_end; - //! \brief Read only view on scattered field initial azimuth. + //! \brief Read-only view on scattered field initial azimuth. const double& sc_theta_start = _sc_theta_start; - //! \brief Read only view on scattered field azimuth step. + //! \brief Read-only view on scattered field azimuth step. const double& sc_theta_step = _sc_theta_step; - //! \brief Read only view on scattered field final azimuth. + //! \brief Read-only view on scattered field final azimuth. const double& sc_theta_end = _sc_theta_end; - //! \brief Read only view on incident field initial elevation. + //! \brief Read-only view on incident field initial elevation. const double& in_phi_start = _in_phi_start; - //! \brief Read only view on incident field elevation step. + //! \brief Read-only view on incident field elevation step. const double& in_phi_step = _in_phi_step; - //! \brief Read only view on incident field final elevation. + //! \brief Read-only view on incident field final elevation. const double& in_phi_end = _in_phi_end; - //! \brief Read only view on scattered field initial elevation. + //! \brief Read-only view on scattered field initial elevation. const double& sc_phi_start = _sc_phi_start; - //! \brief Read only view on scattered field elevation step. + //! \brief Read-only view on scattered field elevation step. const double& sc_phi_step = _sc_phi_step; - //! \brief Read only view on scattered field final elevation. + //! \brief Read-only view on scattered field final elevation. const double& sc_phi_end = _sc_phi_end; + //! \brief Read-only view on flag for matrix inversion refinement. + const short& refine_flag = _refine_flag; + //! \brief Read-only view on flag for dynamic order management. + const short& dyn_order_flag = _dyn_order_flag; /*! \brief Build a scattering geometry configuration structure. * @@ -361,25 +369,25 @@ protected: */ void write_legacy(const std::string& file_name); public: - //! \brief Read only view on name of the reference variable type. + //! \brief Read-only view on name of the reference variable type. const std::string& reference_variable_name = _reference_variable_name; - //! \brief Read only view on number of spherical components. + //! \brief Read-only view on number of spherical components. const int& number_of_spheres = _number_of_spheres; - //! \brief Read only view on number of configurations. + //! \brief Read-only view on number of configurations. const int& configurations = _configurations; - //! \brief Read only view on number of scales to use in calculation. + //! \brief Read-only view on number of scales to use in calculation. const int& number_of_scales = _number_of_scales; - //! \brief Read only view on type of dielectric functions. + //! \brief Read-only view on type of dielectric functions. const int& idfc = _idfc; - //! \brief Read only view on external medium dielectric constant. + //! \brief Read-only view on external medium dielectric constant. const double& exdc = _exdc; - //! \brief Read only view on WP. + //! \brief Read-only view on WP. const double& wp = _wp; - //! \brief Read only view on peak XI. + //! \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. + //! \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; //! \brief Matrix of fractional transition radii with size [CONFIGURATIONS x LAYERS]. double **_rcf; diff --git a/src/include/IterationData.h b/src/include/IterationData.h index a8d82dbf9e2d5a02d583fc94cf3f5f80442437dd..effb94faf480d503675babebbf27c7f489507380 100644 --- a/src/include/IterationData.h +++ b/src/include/IterationData.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/TransitionMatrix.h b/src/include/TransitionMatrix.h index a0be4c77475a685a8efbcf676acf76c787a50cb9..1ff69a1c7b5653a6a66e5de2b5a9a64048066bad 100644 --- a/src/include/TransitionMatrix.h +++ b/src/include/TransitionMatrix.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/algebraic.h b/src/include/algebraic.h index c59dbae2c43a452313970d1e3cca215b407cc55a..f7b8bdacfb5ba50b108a6e0437b402d32e4d6c6f 100644 --- a/src/include/algebraic.h +++ b/src/include/algebraic.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h index 40c5f0e322702b1a3e3c7801cdaf400fc792e852..a02d51b70bd51aecd373c08183ee929b54b8ab5c 100644 --- a/src/include/clu_subs.h +++ b/src/include/clu_subs.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/cublas_calls.h b/src/include/cublas_calls.h index 3d379fdef70ae6bbef2e27d28c1c4496296e8552..8c711175826d79c4858a330da795ce96c5bb6095 100644 --- a/src/include/cublas_calls.h +++ b/src/include/cublas_calls.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/inclu_subs.h b/src/include/inclu_subs.h index c7af5e26c70594e4a0d9184acf33fffc1de9113e..aaa89a5460a5131949499e03c41ee186b2b85032 100644 --- a/src/include/inclu_subs.h +++ b/src/include/inclu_subs.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/lapack_calls.h b/src/include/lapack_calls.h index 45a18733b3432697dd46bc3f763df74a0129a3a1..e0dc6e8222fe6b291ebe17bbb54ed973fab7e6e1 100644 --- a/src/include/lapack_calls.h +++ b/src/include/lapack_calls.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/logging.h b/src/include/logging.h index f5fc061535419a63e2498744b376177e2cbbcada..689c9ca6d6c5740d98853f50edade993d0372391 100644 --- a/src/include/logging.h +++ b/src/include/logging.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/magma_calls.h b/src/include/magma_calls.h index 679b07f37c70b0cc88036f0af4e4cca9163ada78..4767a3c488c0b3914bea263233f170672a6de9ba 100644 --- a/src/include/magma_calls.h +++ b/src/include/magma_calls.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/outputs.h b/src/include/outputs.h index 82ee40ae0c2d19dbaa18135e680c1f001a94ecd6..33d2a2c67eb9483d919bee2b5553759e59756fb5 100644 --- a/src/include/outputs.h +++ b/src/include/outputs.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h index 1bac56f53536c89c113255e661ab1733c4396d20..41c000a6ea4fc893376cb3c3eff1258369174e91 100644 --- a/src/include/sph_subs.h +++ b/src/include/sph_subs.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/tfrfme.h b/src/include/tfrfme.h index ded5377a19f0a4879e57e3733601648b4db12913..25666aa9d8086fdf1c2afc925c2b3545de56908a 100644 --- a/src/include/tfrfme.h +++ b/src/include/tfrfme.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -27,14 +27,14 @@ class Swap1 { protected: //! Index of the last element to be filled. - int last_index; + int _last_index; //! Number of vector coordinates. QUESTION: correct? - int nkv; + int _nkv; //! NLMMT = 2 * LM * (LM + 2) - int nlmmt; + int _nlmmt; //! QUESTION: definition? - dcomplex *wk; + dcomplex *_wk; /*! \brief Load a Swap1 instance from a HDF5 binary file. * @@ -64,24 +64,24 @@ protected: public: //! \brief Read only view on WK. - const dcomplex *vec_wk; + const dcomplex *wk; /*! \brief Swap1 instance constructor. * * \param lm: `int` Maximum field expansion order. * \param _nkv: `int` Number of vector coordinates. QUESTION: correct? */ - Swap1(int lm, int _nkv); + Swap1(int lm, int nkv); /*! \brief Swap1 instance destroyer. */ - ~Swap1() { delete[] wk; } + ~Swap1() { delete[] _wk; } /*! \brief Append an element at the end of the vector. * * \param value: `complex double` The value to be added to the vector. */ - void append(dcomplex value) { wk[last_index++] = value; } + void append(dcomplex value) { _wk[_last_index++] = value; } /*! \brief Load a Swap1 instance from binary file. * @@ -95,14 +95,14 @@ public: /*! \brief Calculate the necessary amount of memory to create a new instance. * * \param lm: `int` Maximum field expansion order. - * \param _nkv: `int` Number of vector coordinates. QUESTION: correct? + * \param nkv: `int` Number of vector coordinates. QUESTION: correct? * \return size: `long` The necessary memory size in bytes. */ - static long get_memory_requirement(int lm, int _nkv); + static long get_size(int lm, int nkv); /*! \brief Bring the pointer to the next element at the start of vector. */ - void reset() { last_index = 0; } + void reset() { _last_index = 0; } /*! \brief Write a Swap1 instance to binary file. * @@ -127,41 +127,39 @@ public: class Swap2 { protected: //! Index of the last vector element to be filled. - int last_vector; + int _last_vector; //! Index of the last matrix element to be filled. - int last_matrix; + int _last_matrix; //! Number of vector coordinates. QUESTION: correct? - int nkv; - //! QUESTION: definition? - double *vkv; + int _nkv; + //! Contiguous vector of VKZM matrix. + double *vec_vkzm; //! QUESTION: definition? - double **vkzm; + double _apfafa; //! QUESTION: definition? - double apfafa; + double _pmf; //! QUESTION: definition? - double pmf; + double _spd; //! QUESTION: definition? - double spd; + double _rir; //! QUESTION: definition? - double rir; + double _ftcn; //! QUESTION: definition? - double ftcn; + double _fshmx; //! QUESTION: definition? - double fshmx; - //! QUESTION: definition? - double vxyzmx; + double _vxyzmx; //! Cartesian displacement. QUESTION: correct? - double delxyz; + double _delxyz; //! QUESTION: definition? - double vknmx; + double _vknmx; //! QUESTION: definition? - double delk; + double _delk; //! QUESTION: definition? - double delks; + double _delks; //! NLMMT = LM * (LM + 2) * 2 - int nlmmt; + int _nlmmt; //! Number of radial vector coordinates. QUESTION: correct? - int nrvc; + int _nrvc; /*! \brief Load a Swap2 instance from a HDF5 binary file. * @@ -190,11 +188,48 @@ protected: void write_legacy(const std::string& file_name); public: + //! Read-only view onthe index of the last vector element to be filled. + const int &last_vector = _last_vector; + //! Read-only view on the index of the last matrix element to be filled. + const int &last_matrix = _last_matrix; + //! Read-only view on the number of vector coordinates. QUESTION: correct? + const int &nkv = _nkv; + //! QUESTION: definition? + double *vkv; + //! QUESTION: definition? + double **vkzm; + //! QUESTION: definition? + const double &apfafa = _apfafa; + //! QUESTION: definition? + const double &pmf = _pmf; + //! QUESTION: definition? + const double &spd = _spd; + //! QUESTION: definition? + const double &rir = _rir; + //! QUESTION: definition? + const double &ftcn = _ftcn; + //! QUESTION: definition? + const double &fshmx = _fshmx; + //! QUESTION: definition? + const double &vxyzmx = _vxyzmx; + //! Cartesian displacement. QUESTION: correct? + const double &delxyz = _delxyz; + //! QUESTION: definition? + const double &vknmx = _vknmx; + //! QUESTION: definition? + const double &delk = _delk; + //! QUESTION: definition? + const double &delks = _delks; + //! NLMMT = LM * (LM + 2) * 2 + const int &nlmmt = _nlmmt; + //! Number of radial vector coordinates. QUESTION: correct? + const int &nrvc = _nrvc; + /*! \brief Swap2 instance constructor. * - * \param _nkv: `int` Number of vector coordinates. QUESTION: correct? + * \param nkv: `int` Number of vector coordinates. QUESTION: correct? */ - Swap2(int _nkv); + Swap2(int nkv); /*! \brief Swap2 instance destroyer. */ @@ -217,17 +252,10 @@ public: /*! \brief Calculate the necessary amount of memory to create a new instance. * - * \param _nkv: `int` Number of radial vector coordinates. QUESTION: correct? + * \param nkv: `int` Number of radial vector coordinates. QUESTION: correct? * \return size: `long` The necessary memory size in bytes. */ - static long get_memory_requirement(int _nkv); - - /*! \brief Get a parameter by its name. - * - * \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 nkv); /*! \brief Get the pointer to the VKV vector. * @@ -245,15 +273,15 @@ public: * * \param value: `double` The value to be pushed in the vector. */ - void push_vector(double value) { vkv[last_vector++] = value; } + void push_vector(double value) { vkv[_last_vector++] = value; } /*! \brief Bring the matrix pointer to the start of the array. */ - void reset_matrix() { last_matrix = 0; } + void reset_matrix() { _last_matrix = 0; } /*! \brief Bring the vector pointer to the start of the array. */ - void reset_vector() { last_vector = 0; } + void reset_vector() { _last_vector = 0; } /*! \brief Set a parameter by its name and value. * @@ -285,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 @@ -324,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. * @@ -355,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. */ @@ -380,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. * diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h index 5ac6c1f5e2c98c8a5b3e017c45405fca8c2c86e2..01f3e0b2ec4691ed30ad09c8582ca4205acfecb7 100644 --- a/src/include/tra_subs.h +++ b/src/include/tra_subs.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/types.h b/src/include/types.h index ba73e2a1edee19c6a8b66090a55adef6fb566e47..2da04679bda9d6475c8e6ea08a8f18a23ab82c01 100644 --- a/src/include/types.h +++ b/src/include/types.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/include/utils.h b/src/include/utils.h index 82779beb33fa14026b8a20c3cbf9f47a01a0c65c..6678c1ad8adff28cb2d8a1c206a41a8e26c6101d 100644 --- a/src/include/utils.h +++ b/src/include/utils.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/inclusion/inclusion.cpp b/src/inclusion/inclusion.cpp index b70306c044b424f3d719d8df501b98f605964d44..1d4433a0b7d49fc7a1d342cc26975927df496e26 100644 --- a/src/inclusion/inclusion.cpp +++ b/src/inclusion/inclusion.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -743,9 +743,15 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo + ").\n"; logger->log(message, LOG_WARN); } else if (recommended_li < cid->c1->li) { - message = "INFO: lowering internal order from " + to_string(cid->c1->li) + " to " - + to_string(recommended_li) + " for scale iteration " + to_string(jxi488) + ".\n"; - logger->log(message, LOG_INFO); + if (gconf->dyn_order_flag > 0) { + message = "INFO: lowering internal order from " + to_string(cid->c1->li) + " to " + + to_string(recommended_li) + " for scale iteration " + to_string(jxi488) + ".\n"; + logger->log(message, LOG_INFO); + } else { + message = "WARNING: internal order " + to_string(cid->c1->li) + " too high for scale iteration " + + to_string(jxi488) + ".\n"; + logger->log(message, LOG_WARN); + } } if (recommended_le > cid->c1->le) { message = "WARNING: external order " + to_string(cid->c1->le) + " for scale iteration " @@ -753,17 +759,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo + ").\n"; logger->log(message, LOG_WARN); } else if (recommended_le < cid->c1->le) { - message = "INFO: lowering external order from " + to_string(cid->c1->le) + " to " - + to_string(recommended_le) + " for scale iteration " + to_string(jxi488) + ".\n"; - logger->log(message, LOG_INFO); + if (gconf->dyn_order_flag > 0) { + message = "INFO: lowering external order from " + to_string(cid->c1->le) + " to " + + to_string(recommended_le) + " for scale iteration " + to_string(jxi488) + ".\n"; + logger->log(message, LOG_INFO); + } else { + message = "WARNING: external order " + to_string(cid->c1->le) + " too high for scale iteration " + + to_string(jxi488) + ".\n"; + logger->log(message, LOG_WARN); + } } if (recommended_li < max_li || recommended_le < max_le) { - int new_li = (recommended_li < max_li) ? recommended_li : max_li; - int new_le = (recommended_le < max_le) ? recommended_le : max_le; - cid->update_orders(sconf->_rcf, new_li, new_le); - is_first_scale = true; - jaw = 1; - cid->refinemode = 2; + if (gconf->dyn_order_flag > 0) { + int new_li = (recommended_li < max_li) ? recommended_li : max_li; + int new_le = (recommended_le < max_le) ? recommended_le : max_le; + cid->update_orders(sconf->_rcf, new_li, new_le); + is_first_scale = true; + jaw = 1; + cid->refinemode = 2; + } } } int li = cid->c1->li; @@ -873,16 +887,18 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo double actualaccuracy = cid->accuracygoal; invert_matrix(cid->am, cid->c1->ndm, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, output_path, jxi488, mxndm, cid->proc_device); // in principle, we should check whether the returned actualaccuracy is indeed lower than the accuracygoal, and do something about it if not -#ifdef USE_REFINEMENT - if (cid->refinemode==2) { - message = "INFO: calibration obtained accuracy " + to_string(actualaccuracy) + " (" + to_string(cid->accuracygoal) + " requested) in " + to_string(cid->maxrefiters) + " refinement iterations\n"; - logger->log(message); - if (actualaccuracy > 1e-2) { - printf("Accuracy worse than 0.01, stopping"); - exit(1); + if (gconf->refine_flag > 0) { + if (cid->refinemode==2) { + message = "DEBUG: iterative refinement enabled at run-time.\n"; + logger->log(message, LOG_DEBG); + message = "INFO: calibration obtained accuracy " + to_string(actualaccuracy) + " (" + to_string(cid->accuracygoal) + " requested) in " + to_string(cid->maxrefiters) + " refinement iterations\n"; + logger->log(message); + if (actualaccuracy > 1e-2) { + printf("Accuracy worse than 0.01, stopping"); + exit(1); + } } } -#endif // USE_REFINEMENT #ifdef USE_NVTX nvtxRangePop(); #endif diff --git a/src/inclusion/np_inclusion.cpp b/src/inclusion/np_inclusion.cpp index fde9371b53744de005bfdd6302226ab84a28a948..875fe094287dbfd7741f7d7ba9835c0b6d9c6c55 100644 --- a/src/inclusion/np_inclusion.cpp +++ b/src/inclusion/np_inclusion.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp index b93abab2ce2fa19c341256ed4a9847901eb99f4c..03dd266352cf63495d425be6e25905d87cced0f9 100644 --- a/src/libnptm/Commons.cpp +++ b/src/libnptm/Commons.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp index 2aa61343aaf91e99294d2d786c72ddae2f290c5b..d348ef167e2c74404fa1feef883fa7ff74aa4c9d 100644 --- a/src/libnptm/Configuration.cpp +++ b/src/libnptm/Configuration.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -64,15 +64,13 @@ using namespace std; GeometryConfiguration::GeometryConfiguration( - int nsph, int lm, int in_pol, int npnt, int npntts, int isam, - int li, int le, np_int mxndm, int iavm, - double *x, double *y, double *z, - double in_th_start, double in_th_step, double in_th_end, - double sc_th_start, double sc_th_step, double sc_th_end, - double in_ph_start, double in_ph_step, double in_ph_end, - double sc_ph_start, double sc_ph_step, double sc_ph_end, - int jwtm - ) { + int nsph, int lm, int in_pol, int npnt, int npntts, int isam, + int li, int le, np_int mxndm, int iavm, double *x, double *y, + double *z, double in_th_start, double in_th_step, double in_th_end, + double sc_th_start, double sc_th_step, double sc_th_end, + double in_ph_start, double in_ph_step, double in_ph_end, + double sc_ph_start, double sc_ph_step, double sc_ph_end, int jwtm +) { _number_of_spheres = nsph; _l_max = lm; _in_pol = in_pol; @@ -99,6 +97,8 @@ GeometryConfiguration::GeometryConfiguration( _sph_x = x; _sph_y = y; _sph_z = z; + _refine_flag = 0; + _dyn_order_flag = 1; } GeometryConfiguration::GeometryConfiguration(const GeometryConfiguration& rhs) @@ -134,6 +134,8 @@ GeometryConfiguration::GeometryConfiguration(const GeometryConfiguration& rhs) _sph_y[ni] = rhs._sph_y[ni]; _sph_z[ni] = rhs._sph_z[ni]; } + _refine_flag = rhs._refine_flag; + _dyn_order_flag = rhs._dyn_order_flag; } #ifdef MPI_VERSION @@ -169,6 +171,8 @@ GeometryConfiguration::GeometryConfiguration(const mixMPI *mpidata) { MPI_Bcast(_sph_x, _number_of_spheres, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(_sph_y, _number_of_spheres, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(_sph_z, _number_of_spheres, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&_refine_flag, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast(&_dyn_order_flag, 1, MPI_SHORT, 0, MPI_COMM_WORLD); } void GeometryConfiguration::mpibcast(const mixMPI *mpidata) { @@ -200,6 +204,8 @@ void GeometryConfiguration::mpibcast(const mixMPI *mpidata) { MPI_Bcast(_sph_x, _number_of_spheres, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(_sph_y, _number_of_spheres, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(_sph_z, _number_of_spheres, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&_refine_flag, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast(&_dyn_order_flag, 1, MPI_SHORT, 0, MPI_COMM_WORLD); } #endif @@ -221,6 +227,7 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(const std::string& fil OpenConfigurationFileException ex(file_name); throw ex; } + // Read the legacy FORTRAN mandatory configuration data int _nsph = 0, _lm = 0, _in_pol = 0, _npnt = 0, _npntts = 0, _isam = 0; int _li = 0, _le = 0, _iavm = 0, num_params = 0; np_int _mxndm = 0; @@ -319,6 +326,7 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(const std::string& fil str_target = file_lines[last_read_line++]; regex_search(str_target, m, re); fjwtm = stoi(m.str()); + // Mandatory configuration data were read. Create the configuration object. GeometryConfiguration *conf = new GeometryConfiguration( _nsph, _lm, _in_pol, _npnt, _npntts, _isam, _li, _le, _mxndm, _iavm, x, y, z, @@ -328,6 +336,25 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(const std::string& fil sc_ph_start, sc_ph_step, sc_ph_end, fjwtm ); + + // Read optional configuration data used only by the C++ code. + while (num_lines > last_read_line) { + str_target = file_lines[last_read_line++]; + if (str_target.size() > 0) { + if (str_target.substr(0, 15).compare("USE_REFINEMENT=") == 0) { + regex_search(str_target, m, re); + short refine_flag = (short)stoi(m.str()); + conf->_refine_flag = refine_flag; + } + else if (str_target.substr(0, 14).compare("USE_DYN_ORDER=") == 0) { + regex_search(str_target, m, re); + short dyn_order_flag = (short)stoi(m.str()); + conf->_dyn_order_flag = dyn_order_flag; + } + } + } + + // Clean up memory and return configuration object. delete[] file_lines; return conf; } diff --git a/src/libnptm/TransitionMatrix.cpp b/src/libnptm/TransitionMatrix.cpp index 8c72f899998f3860167977661587f776ab365db2..754951ad636e39a5979cb7c7e974db24e10e0c77 100644 --- a/src/libnptm/TransitionMatrix.cpp +++ b/src/libnptm/TransitionMatrix.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/algebraic.cpp b/src/libnptm/algebraic.cpp index 33158f0c0dec340389d8fb61ca34823b13ddb73f..e2fe94f849429f35c8f7f810112884821525f5e4 100644 --- a/src/libnptm/algebraic.cpp +++ b/src/libnptm/algebraic.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp index 3548af3d656cff48a18fddafafb2e61a04542852..7bd73b106914ffd80ef5ab02a57f5cde17b6c938 100644 --- a/src/libnptm/clu_subs.cpp +++ b/src/libnptm/clu_subs.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -1341,6 +1341,8 @@ void pcros(double vk, double exri, ParticleDescriptor *c1) { #endif #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd reduction(+:sum, sump, sum1, sum2, sum3, sum4) +#else +#pragma omp parallel for simd reduction(+:sum, sump, sum1, sum2, sum3, sum4) #endif for (int i12 = 0; i12 < nlemt; i12++) { // int i = i12 - 1; @@ -1408,6 +1410,8 @@ void pcrsm0(double vk, double exri, int inpol, ParticleDescriptor *c1) { sum3 = cc0; #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd reduction(+:sum2,sum3) +#else +#pragma omp parallel for simd reduction(+:sum2,sum3) #endif for (int i14 = 0; i14 < c1->nlem; i14++) { int ie = i14 + c1->nlem; @@ -1418,6 +1422,8 @@ void pcrsm0(double vk, double exri, int inpol, ParticleDescriptor *c1) { dcomplex sumpd = cc0; #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd collapse(2) reduction(+:sumpi,sumpd) +#else +#pragma omp parallel for simd collapse(2) reduction(+:sumpi,sumpd) #endif for (int i16 = 0; i16 < nlemt; i16++) { for (int j16 = 0; j16 < c1->nlem; j16++) { @@ -2001,6 +2007,8 @@ void raba( #endif #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd reduction(+:c1, c2) +#else +#pragma omp parallel for simd reduction(+:c1, c2) #endif for (int j10 = 1; j10 <= nlemt; j10++) { int j = j10 - 1; @@ -2021,6 +2029,8 @@ void raba( #endif #ifdef USE_TARGET_OFFLOAD #pragma omp teams distribute parallel for +#else +#pragma omp parallel for #endif for (int ipo = 0; ipo < 2; ipo++) { int jpo = 1 - ipo; @@ -2051,12 +2061,14 @@ void raba( int kmax = le*(le+2); // for efficiency I should also linearise array w, but I cannot easily since I do not know for sure its major dimension (changes to containing class needed) #ifdef USE_NVTX - nvtxRangePush("raba inner loop 2"); + nvtxRangePush("raba inner loop 2"); #endif #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd reduction(+:ctqce0, ctqce1, ctqce2, ctqcs0, ctqcs1, ctqcs2, tqcpe0, tqcpe1, tqcpe2, tqcps0, tqcps1, tqcps2) +#else +#pragma omp parallel for simd reduction(+:ctqce0, ctqce1, ctqce2, ctqcs0, ctqcs1, ctqcs2, tqcpe0, tqcpe1, tqcpe2, tqcps0, tqcps1, tqcps2) #endif - for (int k = 1; k<=kmax; k++) { + for (int k = 1; k<=kmax; k++) { int l60 = (int) sqrt(k+1); int im60 = k - (l60*l60) + 1; if (im60 == 0) { @@ -2079,44 +2091,44 @@ void raba( dcomplex acwp; dcomplex aca; dcomplex acap; - if (mmmmu <= l60) { - int immu = mmmu + il - 1; - int immue = immu + nlem; - rmu = -sqrt(1.0 * (l60 + mmmu) * (l60 - m)) * sq2i; - acw = dconjg(vec_a[2*i+ipo]) * vec_w[4*immu+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*immue+ipo]; - acwp = dconjg(vec_a[2*i+ipo]) * vec_w[4*immu+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*immue+jpo]; - aca = dconjg(vec_a[2*i+ipo]) * vec_a[2*immu+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*immue+ipo]; - acap = dconjg(vec_a[2*i+ipo]) * vec_a[2*immu+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*immue+jpo]; - ctqce0 += (acw * rmu); - tqcpe0 += (acwp * rmu); - ctqcs0 += (aca * rmu); - tqcps0 += (acap * rmu); - } - // label 30 - rmu = -1.0 * m; - acw = dconjg(vec_a[2*i+ipo]) * vec_w[4*i+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*ie+ipo]; - acwp = dconjg(vec_a[2*i+ipo]) * vec_w[4*i+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*ie+jpo]; - aca = dconjg(vec_a[2*i+ipo]) * vec_a[2*i+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*ie+ipo]; - acap = dconjg(vec_a[2*i+ipo]) * vec_a[2*i+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*ie+jpo]; - ctqce1 += (acw * rmu); - tqcpe1 += (acwp * rmu); - ctqcs1 += (aca * rmu); - tqcps1 += (acap * rmu); - mmmu = m - 1; - mmmmu = (mmmu > 0) ? mmmu : -mmmu; - if (mmmmu <= l60) { - int immu = mmmu + il - 1; - int immue = immu + nlem; - rmu = sqrt(1.0 * (l60 - mmmu) * (l60 + m)) * sq2i; - acw = dconjg(vec_a[2*i+ipo]) * vec_w[4*immu+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*immue+ipo]; - acwp = dconjg(vec_a[2*i+ipo]) * vec_w[4*immu+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*immue+jpo]; - aca = dconjg(vec_a[2*i+ipo]) * vec_a[2*immu+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*immue+ipo]; - acap = dconjg(vec_a[2*i+ipo]) * vec_a[2*immu+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*immue+jpo]; - ctqce2 += (acw * rmu); - tqcpe2 += (acwp * rmu); - ctqcs2 += (aca * rmu); - tqcps2 += (acap * rmu); - } // ends if clause + if (mmmmu <= l60) { + int immu = mmmu + il - 1; + int immue = immu + nlem; + rmu = -sqrt(1.0 * (l60 + mmmu) * (l60 - m)) * sq2i; + acw = dconjg(vec_a[2*i+ipo]) * vec_w[4*immu+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*immue+ipo]; + acwp = dconjg(vec_a[2*i+ipo]) * vec_w[4*immu+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*immue+jpo]; + aca = dconjg(vec_a[2*i+ipo]) * vec_a[2*immu+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*immue+ipo]; + acap = dconjg(vec_a[2*i+ipo]) * vec_a[2*immu+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*immue+jpo]; + ctqce0 += (acw * rmu); + tqcpe0 += (acwp * rmu); + ctqcs0 += (aca * rmu); + tqcps0 += (acap * rmu); + } + // label 30 + rmu = -1.0 * m; + acw = dconjg(vec_a[2*i+ipo]) * vec_w[4*i+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*ie+ipo]; + acwp = dconjg(vec_a[2*i+ipo]) * vec_w[4*i+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*ie+jpo]; + aca = dconjg(vec_a[2*i+ipo]) * vec_a[2*i+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*ie+ipo]; + acap = dconjg(vec_a[2*i+ipo]) * vec_a[2*i+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*ie+jpo]; + ctqce1 += (acw * rmu); + tqcpe1 += (acwp * rmu); + ctqcs1 += (aca * rmu); + tqcps1 += (acap * rmu); + mmmu = m - 1; + mmmmu = (mmmu > 0) ? mmmu : -mmmu; + if (mmmmu <= l60) { + int immu = mmmu + il - 1; + int immue = immu + nlem; + rmu = sqrt(1.0 * (l60 - mmmu) * (l60 + m)) * sq2i; + acw = dconjg(vec_a[2*i+ipo]) * vec_w[4*immu+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*immue+ipo]; + acwp = dconjg(vec_a[2*i+ipo]) * vec_w[4*immu+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_w[4*immue+jpo]; + aca = dconjg(vec_a[2*i+ipo]) * vec_a[2*immu+ipo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*immue+ipo]; + acap = dconjg(vec_a[2*i+ipo]) * vec_a[2*immu+jpo] + dconjg(vec_a[2*ie+ipo]) * vec_a[2*immue+jpo]; + ctqce2 += (acw * rmu); + tqcpe2 += (acwp * rmu); + ctqcs2 += (aca * rmu); + tqcps2 += (acap * rmu); + } // ends if clause } // k loop (previously the l60 and im60 loops #ifdef USE_NVTX nvtxRangePop(); @@ -2129,7 +2141,9 @@ void raba( nvtxRangePush("raba loop 3"); #endif #ifdef USE_TARGET_OFFLOAD -#pragma omp teams distribute parallel for simd +#pragma omp target teams distribute parallel for simd +#else +#pragma omp parallel for simd #endif for (int ipo78 = 1; ipo78 <= 2; ipo78++) { int ipo = ipo78 - 1; @@ -2202,6 +2216,8 @@ void scr0(double vk, double exri, ParticleDescriptor *c1) { #endif #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd reduction(+:sums, sum21) +#else +#pragma omp parallel for simd reduction(+:sums, sum21) #endif for (int l10 = 1; l10 <= c1->li; l10++) { double fl = 1.0 * (l10 + l10 + 1); @@ -2248,6 +2264,8 @@ void scr0(double vk, double exri, ParticleDescriptor *c1) { #endif #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd reduction(+:scs, ecs, acs, tfsas) +#else +#pragma omp parallel for simd reduction(+:scs, ecs, acs, tfsas) #endif for (int i14 = 1; i14 <= c1->nsph; i14++) { int iogi = c1->iog[i14 - 1]; @@ -2312,6 +2330,8 @@ void scr2( #endif #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd reduction(-:s11, s21, s12, s22) +#else +#pragma omp parallel for simd reduction(-:s11, s21, s12, s22) #endif for (int k = 1; k<=kmax; k++) { int l10 = (int) sqrt(k+1); @@ -2366,6 +2386,8 @@ void scr2( #endif #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd reduction(+:tsas00, tsas10, tsas01, tsas11) +#else +#pragma omp parallel for simd reduction(+:tsas00, tsas10, tsas01, tsas11) #endif for (int i14 = 1; i14 <= c1->nsph; i14++) { int i = i14 - 1; @@ -2398,6 +2420,8 @@ void scr2( #endif #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd collapse(4) +#else +#pragma omp parallel for simd collapse(4) #endif for (int ipo1 = 1; ipo1 <=2; ipo1++) { for (int jpo1 = 1; jpo1 <= 2; jpo1++) { @@ -2422,7 +2446,9 @@ void scr2( nvtxRangePush("scr2 loop 4"); #endif #ifdef USE_TARGET_OFFLOAD -#pragma omp target parallel for collapse(4) +#pragma omp target teams distribute parallel for collapse(4) +#else +#pragma omp parallel for collapse(4) #endif for (int ipo1 = 1; ipo1 <=2; ipo1++) { for (int jpo1 = 1; jpo1 <= 2; jpo1++) { @@ -2534,7 +2560,7 @@ void tqr( } void ztm(dcomplex **am, ParticleDescriptor *c1) { - dcomplex gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4; + // dcomplex gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4; const dcomplex cc0 = 0.0 + 0.0 * I; // int i2 = 0; // old implementation #ifdef USE_NVTX @@ -2546,7 +2572,6 @@ void ztm(dcomplex **am, ParticleDescriptor *c1) { // C9 *c9_para = new C9(*c9); dcomplex *gis_v = c1->gis[0]; dcomplex *gls_v = c1->gls[0]; - double *rac3j_local = (double *) malloc(c1->lmtpo*sizeof(double)); int k2max = c1->li*(c1->li+2); int k3max = c1->le*(c1->le+2); // To parallelise, I run a linearised loop directly over k @@ -2559,6 +2584,8 @@ void ztm(dcomplex **am, ParticleDescriptor *c1) { // im = im -(2*l+1) and l = l+1 (there was a rounding error in a nearly exact root) #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd collapse(3) +#else +#pragma omp parallel for simd collapse(3) #endif for (int n2 = 1; n2 <= c1->nsph; n2++) { // GPU portable? for (int k2 = 1; k2<=k2max; k2++) { @@ -2590,12 +2617,13 @@ void ztm(dcomplex **am, ParticleDescriptor *c1) { int i3 = l3 * l3 + im3 - 1; int m3 = -l3 - 1 + im3; int vecindex = (i2 - 1) * c1->nlem + i3 - 1; + double *rac3j_local = (double *) malloc(c1->lmtpo*sizeof(double)); gis_v[vecindex] = ghit_d(2, 0, n2, l2, m2, l3, m3, c1, rac3j_local); gls_v[vecindex] = ghit_d(2, 1, n2, l2, m2, l3, m3, c1, rac3j_local); + free(rac3j_local); } // close k3 loop, former l3 + im3 loops } // close k2 loop, former l2 + im2 loops } // close n2 loop - free(rac3j_local); #ifdef USE_NVTX nvtxRangePop(); #endif @@ -2606,6 +2634,8 @@ void ztm(dcomplex **am, ParticleDescriptor *c1) { dcomplex *sam_v = c1->sam[0]; #ifdef USE_TARGET_OFFLOAD #pragma omp target teams distribute parallel for simd collapse(2) +#else +#pragma omp parallel for simd collapse(2) #endif for (int i1 = 1; i1 <= c1->ndi; i1++) { // GPU portable? for (int i3 = 1; i3 <= c1->nlem; i3++) { @@ -2615,7 +2645,6 @@ void ztm(dcomplex **am, ParticleDescriptor *c1) { dcomplex sum4 = cc0; int i1e = i1 + c1->ndi; int i3e = i3 + c1->nlem; -#pragma parallel for simd reduction(+:sum1,sum2,sum3,sum4) for (int i2 = 1; i2 <= c1->ndi; i2++) { int i2e = i2 + c1->ndi; int vecindg_23 = (i2 - 1) * c1->nlem + i3 - 1; @@ -2640,7 +2669,11 @@ void ztm(dcomplex **am, ParticleDescriptor *c1) { sam_v[vecind1e + i3e - 1] = sum4; } // i3 loop } // i1 loop -#pragma omp parallel for collapse(2) +#ifdef USE_TARGET_OFFLOAD +#pragma omp target teams distribute parallel for simd collapse(2) +#else +#pragma omp parallel for simd collapse(2) +#endif for (int i1 = 1; i1 <= c1->ndi; i1++) { for (int i0 = 1; i0 <= c1->nlem; i0++) { int vecindex = (i1 - 1) * c1->nlem + i0 - 1; @@ -2650,7 +2683,9 @@ void ztm(dcomplex **am, ParticleDescriptor *c1) { } // i1 loop dcomplex *vec_am0m = c1->am0m[0]; #ifdef USE_TARGET_OFFLOAD -#pragma omp target parallel for collapse(2) +#pragma omp target teams distribute parallel for simd collapse(2) +#else +#pragma omp parallel for simd collapse(2) #endif for (int i0 = 1; i0 <= c1->nlem; i0++) { for (int i3 = 1; i3 <= c1->nlemt; i3++) { @@ -2661,11 +2696,11 @@ void ztm(dcomplex **am, ParticleDescriptor *c1) { int i1e = i1 + c1->ndi; int vecind1 = (i1 - 1) * c1->nlemt; int vecind1e = (i1e - 1) * c1->nlemt; - a1 = sam_v[vecind1 + i3 - 1]; - a2 = sam_v[vecind1e + i3 - 1]; + dcomplex a1 = sam_v[vecind1 + i3 - 1]; + dcomplex a2 = sam_v[vecind1e + i3 - 1]; int vecindex = (i1 - 1) * c1->nlem + i0 - 1; - gie = gis_v[vecindex]; - gle = gls_v[vecindex]; + dcomplex gie = gis_v[vecindex]; + dcomplex gle = gls_v[vecindex]; sum1 += (a1 * gie + a2 * gle); sum2 += (a1 * gle + a2 * gie); } // i1 loop diff --git a/src/libnptm/cublas_calls.cpp b/src/libnptm/cublas_calls.cpp index caf483a61251d24e9991be7752f181c782c96dc0..7e070762917a6d13a9aff917164812523cff54da 100644 --- a/src/libnptm/cublas_calls.cpp +++ b/src/libnptm/cublas_calls.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/inclu_subs.cpp b/src/libnptm/inclu_subs.cpp index 12859ab2a1953eb2387fa5d3c9eb804fbb8c4f0b..9b1125e515c99068e08475197be88a1012733b86 100644 --- a/src/libnptm/inclu_subs.cpp +++ b/src/libnptm/inclu_subs.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/lapack_calls.cpp b/src/libnptm/lapack_calls.cpp index 0b9af94ed13be9f3d9cf408f887bae334738b49e..407a28d3c3774c87744897fc39e97621a7e7b567 100644 --- a/src/libnptm/lapack_calls.cpp +++ b/src/libnptm/lapack_calls.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/logging.cpp b/src/libnptm/logging.cpp index 514eeaeed947f5012d2ac2d0948f67884522c7a5..a208e03d944032a9816c011bc12355da00b0470b 100644 --- a/src/libnptm/logging.cpp +++ b/src/libnptm/logging.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/magma_calls.cpp b/src/libnptm/magma_calls.cpp index 88a5b811088659bfd213b4d56dcb8406d48a322a..11441645ed86f64a515614fcb8edbfea11377757 100644 --- a/src/libnptm/magma_calls.cpp +++ b/src/libnptm/magma_calls.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/outputs.cpp b/src/libnptm/outputs.cpp index 8c4ecf7df21e3ed3d04320eeac8ec58666d2c510..4b6e7265c4cbc62aea148ac2a799db00cb4353c0 100644 --- a/src/libnptm/outputs.cpp +++ b/src/libnptm/outputs.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp index a67d2adda340d6ccb6ff78b36297784fbe0d6d7f..10cc50a5eb70a6f225884fee7a2d3e1348f5509c 100644 --- a/src/libnptm/sph_subs.cpp +++ b/src/libnptm/sph_subs.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/libnptm/tfrfme.cpp b/src/libnptm/tfrfme.cpp index 2b81215780af3610fd3c535a49326a5663780823..fc3651b637e7227028767ad192cd029765e8a2ec 100644 --- a/src/libnptm/tfrfme.cpp +++ b/src/libnptm/tfrfme.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -47,13 +47,13 @@ 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](); - vec_wk = wk; - last_index = 0; +Swap1::Swap1(int lm, int nkv) { + _nkv = nkv; + _nlmmt = 2 * lm * (lm + 2); + const int size = nkv * nkv * _nlmmt; + _wk = new dcomplex[size](); + wk = _wk; + _last_index = 0; } Swap1* Swap1::from_binary(const std::string& file_name, const std::string& mode) { @@ -76,21 +76,21 @@ Swap1* Swap1::from_hdf5(const std::string& file_name) { herr_t status = hdf_file->get_status(); double *elements; string str_type; - int _nlmmt, _nkv, lm, num_elements, index; + int nlmmt, nkv, lm, num_elements, index; dcomplex value; if (status == 0) { - status = hdf_file->read("NLMMT", "INT32", &_nlmmt); - status = hdf_file->read("NKV", "INT32", &_nkv); - lm = (int)(sqrt(4.0 + 2.0 * _nlmmt) / 2.0) - 1; - num_elements = 2 * _nlmmt * _nkv * _nkv; - instance = new Swap1(lm, _nkv); + status = hdf_file->read("NLMMT", "INT32", &nlmmt); + status = hdf_file->read("NKV", "INT32", &nkv); + lm = (int)(sqrt(4.0 + 2.0 * nlmmt) / 2.0) - 1; + num_elements = 2 * nlmmt * nkv * nkv; + instance = new Swap1(lm, nkv); 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; - instance->wk[wi] = value; + instance->_wk[wi] = value; } // wi loop delete[] elements; status = hdf_file->close(); @@ -102,19 +102,19 @@ Swap1* Swap1::from_hdf5(const std::string& file_name) { Swap1* Swap1::from_legacy(const std::string& file_name) { fstream input; Swap1 *instance = NULL; - int _nlmmt, _nkv, lm; + 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(&_nlmmt), sizeof(int)); - lm = (int)(sqrt(4.0 + 2.0 * _nlmmt) / 2.0) - 1; - input.read(reinterpret_cast(&_nkv), sizeof(int)); - instance = new Swap1(lm, _nkv); - int num_elements = _nlmmt * _nkv * _nkv; + input.read(reinterpret_cast(&nlmmt), sizeof(int)); + lm = (int)(sqrt(4.0 + 2.0 * nlmmt) / 2.0) - 1; + input.read(reinterpret_cast(&nkv), sizeof(int)); + instance = new Swap1(lm, nkv); + int num_elements = nlmmt * nkv * nkv; for (int j = 0; j < num_elements; j++) { input.read(reinterpret_cast(&rval), sizeof(double)); input.read(reinterpret_cast(&ival), sizeof(double)); - instance->wk[j] = rval + ival * I; + instance->_wk[j] = rval + ival * I; } input.close(); } else { @@ -123,9 +123,9 @@ Swap1* Swap1::from_legacy(const std::string& file_name) { return instance; } -long Swap1::get_memory_requirement(int lm, int _nkv) { +long Swap1::get_size(int lm, int nkv) { long size = (long)(3 * sizeof(int)); - size += (long)(sizeof(dcomplex) * 2 * lm * (lm + 2) * _nkv * _nkv); + size += (long)(sizeof(dcomplex) * 2 * lm * (lm + 2) * nkv * nkv); return size; } @@ -146,20 +146,20 @@ void Swap1::write_hdf5(const std::string& file_name) { List rec_ptr_list(1); herr_t status; string str_type; - int num_elements = 2 * nlmmt * nkv * nkv; + 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_ptr_list.set(0, &_nlmmt); 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("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]); + ptr_elements[2 * wi] = real(_wk[wi]); + ptr_elements[2 * wi + 1] = imag(_wk[wi]); } rec_ptr_list.append(ptr_elements); @@ -185,12 +185,12 @@ void Swap1::write_legacy(const std::string& file_name) { 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(&nlmmt), sizeof(int)); - output.write(reinterpret_cast(&nkv), sizeof(int)); + int num_elements = _nlmmt * _nkv * _nkv; + output.write(reinterpret_cast(&_nlmmt), sizeof(int)); + output.write(reinterpret_cast(&_nkv), sizeof(int)); for (int j = 0; j < num_elements; j++) { - rval = real(wk[j]); - ival = imag(wk[j]); + rval = real(_wk[j]); + ival = imag(_wk[j]); output.write(reinterpret_cast(&rval), sizeof(double)); output.write(reinterpret_cast(&ival), sizeof(double)); } @@ -201,15 +201,15 @@ void Swap1::write_legacy(const std::string& file_name) { } bool Swap1::operator ==(Swap1 &other) { - if (nlmmt != other.nlmmt) { + if (_nlmmt != other._nlmmt) { return false; } - if (nkv != other.nkv) { + if (_nkv != other._nkv) { return false; } - const int num_elements = nlmmt * nkv * nkv; + const int num_elements = _nlmmt * _nkv * _nkv; for (int i = 0; i < num_elements; i++) { - if (wk[i] != other.wk[i]) { + if (_wk[i] != other._wk[i]) { return false; } } @@ -218,18 +218,19 @@ bool Swap1::operator ==(Swap1 &other) { // >>> END OF Swap1 CLASS IMPLEMENTATION <<< // >>> START OF Swap2 CLASS IMPLEMENTATION <<< -Swap2::Swap2(int _nkv) { - nkv = _nkv; - vkv = new double[nkv](); +Swap2::Swap2(int nkv) { + _nkv = nkv; + vkv = new double[_nkv](); + vec_vkzm = new double[_nkv * _nkv](); vkzm = new double*[nkv]; - for (int vi = 0; vi < nkv; vi++) vkzm[vi] = new double[nkv](); - last_vector = 0; - last_matrix = 0; + for (int vi = 0; vi < _nkv; vi++) vkzm[vi] = vec_vkzm + vi * _nkv; + _last_vector = 0; + _last_matrix = 0; } Swap2::~Swap2() { delete[] vkv; - for (int vi = nkv - 1; vi > -1; vi--) delete[] vkzm[vi]; + delete[] vec_vkzm; delete[] vkzm; } @@ -252,15 +253,15 @@ Swap2* Swap2::from_hdf5(const std::string& file_name) { HDFFile *hdf_file = new HDFFile(file_name, flags); herr_t status = hdf_file->get_status(); string str_type; - int _nkv, _nlmmt, _nrvc; + int fnkv, fnlmmt, fnrvc; 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("NKV", "INT32", &fnkv); + instance = new Swap2(fnkv); + str_type = "FLOAT64_(" + to_string(fnkv) + ")"; 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); + str_type = "FLOAT64_(" + to_string(fnkv * fnkv) + ")"; + status = hdf_file->read("VKZM", str_type, instance->vec_vkzm); status = hdf_file->read("APFAFA", "FLOAT64", &value); instance->set_param("apfafa", value); status = hdf_file->read("PMF", "FLOAT64", &value); @@ -283,10 +284,10 @@ Swap2* Swap2::from_hdf5(const std::string& file_name) { 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->read("NLMMT", "INT32", &fnlmmt); + instance->set_param("nlmmt", 1.0 * fnlmmt); + status = hdf_file->read("NRVC", "INT32", &fnrvc); + instance->set_param("nrvc", 1.0 * fnrvc); status = hdf_file->close(); delete hdf_file; } @@ -296,24 +297,24 @@ Swap2* Swap2::from_hdf5(const std::string& file_name) { Swap2* Swap2::from_legacy(const std::string& file_name) { fstream input; Swap2 *instance = NULL; - int _nkv, _nlmmt, _nrvc; - double **_vkzm = NULL; - double *_vkv = NULL; + int fnkv, fnlmmt, fnrvc; + double **fvkzm = NULL; + double *fvkv = NULL; double value; input.open(file_name.c_str(), ios::in | ios::binary); if (input.is_open()) { - input.read(reinterpret_cast(&_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(&fnkv), sizeof(int)); + instance = new Swap2(fnkv); + fvkzm = instance->get_matrix(); + fvkv = instance->get_vector(); + for (int vj = 0; vj < fnkv; vj++) { input.read(reinterpret_cast(&value), sizeof(double)); - _vkv[vj] = value; + fvkv[vj] = value; } - for (int mi = 0; mi < _nkv; mi++) { - for (int mj = 0; mj < _nkv; mj++) { + for (int mi = 0; mi < fnkv; mi++) { + for (int mj = 0; mj < fnkv; mj++) { input.read(reinterpret_cast(&value), sizeof(double)); - _vkzm[mi][mj] = value; + fvkzm[mi][mj] = value; } } input.read(reinterpret_cast(&value), sizeof(double)); @@ -338,10 +339,10 @@ Swap2* Swap2::from_legacy(const std::string& file_name) { instance->set_param("delk", value); input.read(reinterpret_cast(&value), sizeof(double)); instance->set_param("delks", value); - input.read(reinterpret_cast(&_nlmmt), sizeof(int)); - instance->set_param("nlmmt", 1.0 * _nlmmt); - input.read(reinterpret_cast(&_nrvc), sizeof(int)); - instance->set_param("nrvc", 1.0 * _nrvc); + input.read(reinterpret_cast(&fnlmmt), sizeof(int)); + instance->set_param("nlmmt", 1.0 * fnlmmt); + input.read(reinterpret_cast(&fnrvc), sizeof(int)); + instance->set_param("nrvc", 1.0 * fnrvc); input.close(); } else { printf("ERROR: could not open input file \"%s\"\n", file_name.c_str()); @@ -349,57 +350,34 @@ Swap2* Swap2::from_legacy(const std::string& file_name) { 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)); +long Swap2::get_size(int nkv) { + long size = (long)(10 * sizeof(int) + 22 * 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); + int col = _last_matrix % (_nkv - 1); + int row = _last_matrix - (_nkv * row); vkzm[row][col] = value; - last_matrix++; + _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; + 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); @@ -425,54 +403,54 @@ void Swap2::write_hdf5(const std::string& file_name) { string str_type; rec_name_list.set(0, "NKV"); rec_type_list.set(0, "INT32_(1)"); - rec_ptr_list.set(0, &nkv); + rec_ptr_list.set(0, &_nkv); rec_name_list.append("VKV"); - str_type = "FLOAT64_(" + to_string(nkv) + ")"; + 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) + ")"; + str_type = "FLOAT64_(" + to_string(_nkv * _nkv) + ")"; rec_type_list.append(str_type); - rec_ptr_list.append(vkzm); + rec_ptr_list.append(vec_vkzm); rec_name_list.append("APFAFA"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&apfafa); + rec_ptr_list.append(&_apfafa); rec_name_list.append("PMF"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&pmf); + rec_ptr_list.append(&_pmf); 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("RIR"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&rir); + rec_ptr_list.append(&_rir); rec_name_list.append("FTCN"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&ftcn); + rec_ptr_list.append(&_ftcn); rec_name_list.append("FSHMX"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&fshmx); + rec_ptr_list.append(&_fshmx); rec_name_list.append("VXYZMX"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&vxyzmx); + rec_ptr_list.append(&_vxyzmx); rec_name_list.append("delxyz"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&delxyz); + rec_ptr_list.append(&_delxyz); rec_name_list.append("VKNMX"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&vknmx); + rec_ptr_list.append(&_vknmx); rec_name_list.append("DELK"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&delk); + rec_ptr_list.append(&_delk); rec_name_list.append("DELKS"); rec_type_list.append("FLOAT64_(1)"); - rec_ptr_list.append(&delks); + rec_ptr_list.append(&_delks); rec_name_list.append("NLMMT"); rec_type_list.append("INT32_(1)"); - rec_ptr_list.append(&nlmmt); + rec_ptr_list.append(&_nlmmt); rec_name_list.append("NRVC"); rec_type_list.append("INT32_(1)"); - rec_ptr_list.append(&nrvc); + rec_ptr_list.append(&_nrvc); string *rec_names = rec_name_list.to_array(); string *rec_types = rec_type_list.to_array(); @@ -495,30 +473,30 @@ void Swap2::write_legacy(const std::string& file_name) { double value; output.open(file_name.c_str(), ios::out | ios::binary); if (output.is_open()) { - output.write(reinterpret_cast(&nkv), sizeof(int)); - for (int j = 0; j < nkv; j++){ + output.write(reinterpret_cast(&_nkv), sizeof(int)); + for (int j = 0; j < _nkv; j++){ value = vkv[j]; output.write(reinterpret_cast(&value), sizeof(double)); } - for (int mi = 0; mi < nkv; mi++) { - for (int mj = 0; mj < nkv; mj++) { + for (int mi = 0; mi < _nkv; mi++) { + for (int mj = 0; mj < _nkv; mj++) { value = vkzm[mi][mj]; output.write(reinterpret_cast(&value), sizeof(double)); } } - output.write(reinterpret_cast(&apfafa), sizeof(double)); - output.write(reinterpret_cast(&pmf), sizeof(double)); - output.write(reinterpret_cast(&spd), sizeof(double)); - output.write(reinterpret_cast(&rir), sizeof(double)); - output.write(reinterpret_cast(&ftcn), sizeof(double)); - output.write(reinterpret_cast(&fshmx), sizeof(double)); - output.write(reinterpret_cast(&vxyzmx), sizeof(double)); - output.write(reinterpret_cast(&delxyz), sizeof(double)); - output.write(reinterpret_cast(&vknmx), sizeof(double)); - output.write(reinterpret_cast(&delk), sizeof(double)); - output.write(reinterpret_cast(&delks), sizeof(double)); - output.write(reinterpret_cast(&nlmmt), sizeof(int)); - output.write(reinterpret_cast(&nrvc), sizeof(int)); + output.write(reinterpret_cast(&_apfafa), sizeof(double)); + output.write(reinterpret_cast(&_pmf), sizeof(double)); + output.write(reinterpret_cast(&_spd), sizeof(double)); + output.write(reinterpret_cast(&_rir), sizeof(double)); + output.write(reinterpret_cast(&_ftcn), sizeof(double)); + output.write(reinterpret_cast(&_fshmx), sizeof(double)); + output.write(reinterpret_cast(&_vxyzmx), sizeof(double)); + output.write(reinterpret_cast(&_delxyz), sizeof(double)); + output.write(reinterpret_cast(&_vknmx), sizeof(double)); + output.write(reinterpret_cast(&_delk), sizeof(double)); + output.write(reinterpret_cast(&_delks), sizeof(double)); + output.write(reinterpret_cast(&_nlmmt), sizeof(int)); + output.write(reinterpret_cast(&_nrvc), sizeof(int)); output.close(); } else { // Should never occur printf("ERROR: could not open output file \"%s\"\n", file_name.c_str()); @@ -526,55 +504,55 @@ void Swap2::write_legacy(const std::string& file_name) { } bool Swap2::operator ==(Swap2 &other) { - if (nlmmt != other.nlmmt) { + if (_nlmmt != other._nlmmt) { return false; } - if (nrvc != other.nrvc) { + if (_nrvc != other._nrvc) { return false; } - if (nkv != other.nkv) { + if (_nkv != other._nkv) { return false; } - if (apfafa != other.apfafa) { + if (_apfafa != other._apfafa) { return false; } - if (pmf != other.pmf) { + if (_pmf != other._pmf) { return false; } - if (spd != other.spd) { + if (_spd != other._spd) { return false; } - if (rir != other.rir) { + if (_rir != other._rir) { return false; } - if (ftcn != other.ftcn) { + if (_ftcn != other._ftcn) { return false; } - if (fshmx != other.fshmx) { + if (_fshmx != other._fshmx) { return false; } - if (vxyzmx != other.vxyzmx) { + if (_vxyzmx != other._vxyzmx) { return false; } - if (delxyz != other.delxyz) { + if (_delxyz != other._delxyz) { return false; } - if (vknmx != other.vknmx) { + if (_vknmx != other._vknmx) { return false; } - if (delk != other.delk) { + if (_delk != other._delk) { return false; } - if (delks != other.delks) { + if (_delks != other._delks) { return false; } - for (int vi = 0; vi < nkv; vi++) { + 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++) { + for (int mi = 0; mi < _nkv; mi++) { + for (int mj = 0; mj < _nkv; mj++) { if (vkzm[mi][mj] != other.vkzm[mi][mj]) { return false; } @@ -585,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; } @@ -639,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; @@ -660,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); @@ -675,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 @@ -699,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()) { @@ -721,7 +697,6 @@ TFRFME* TFRFME::from_legacy(const std::string& file_name) { input.read(reinterpret_cast(&frsh), sizeof(double)); input.read(reinterpret_cast(&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); @@ -745,14 +720,14 @@ TFRFME* TFRFME::from_legacy(const std::string& file_name) { input.read(reinterpret_cast(&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(&rval), sizeof(double)); input.read(reinterpret_cast(&ival), sizeof(double)); dcomplex value = rval + ival * I; - _wsum[wi][wj] = value; + instance->wsum[wi][wj] = value; } // wi loop } // wj loop input.close(); @@ -762,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); @@ -831,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); @@ -920,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(&lmode), sizeof(int)); - output.write(reinterpret_cast(&lm), sizeof(int)); - output.write(reinterpret_cast(&nkv), sizeof(int)); - output.write(reinterpret_cast(&nxv), sizeof(int)); - output.write(reinterpret_cast(&nyv), sizeof(int)); - output.write(reinterpret_cast(&nzv), sizeof(int)); - output.write(reinterpret_cast(&vk), sizeof(double)); - output.write(reinterpret_cast(&exri), sizeof(double)); - output.write(reinterpret_cast(&an), sizeof(double)); - output.write(reinterpret_cast(&ff), sizeof(double)); - output.write(reinterpret_cast(&tra), sizeof(double)); - output.write(reinterpret_cast(&spd), sizeof(double)); - output.write(reinterpret_cast(&frsh), sizeof(double)); - output.write(reinterpret_cast(&exril), sizeof(double)); - for (int xi = 0; xi < nxv; xi++) + output.write(reinterpret_cast(&_lmode), sizeof(int)); + output.write(reinterpret_cast(&_lm), sizeof(int)); + output.write(reinterpret_cast(&_nkv), sizeof(int)); + output.write(reinterpret_cast(&_nxv), sizeof(int)); + output.write(reinterpret_cast(&_nyv), sizeof(int)); + output.write(reinterpret_cast(&_nzv), sizeof(int)); + output.write(reinterpret_cast(&_vk), sizeof(double)); + output.write(reinterpret_cast(&_exri), sizeof(double)); + output.write(reinterpret_cast(&_an), sizeof(double)); + output.write(reinterpret_cast(&_ff), sizeof(double)); + output.write(reinterpret_cast(&_tra), sizeof(double)); + output.write(reinterpret_cast(&_spd), sizeof(double)); + output.write(reinterpret_cast(&_frsh), sizeof(double)); + output.write(reinterpret_cast(&_exril), sizeof(double)); + for (int xi = 0; xi < _nxv; xi++) output.write(reinterpret_cast(&(xv[xi])), sizeof(double)); - for (int yi = 0; yi < nyv; yi++) + for (int yi = 0; yi < _nyv; yi++) output.write(reinterpret_cast(&(yv[yi])), sizeof(double)); - for (int zi = 0; zi < nzv; zi++) + for (int zi = 0; zi < _nzv; zi++) output.write(reinterpret_cast(&(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(&rval), sizeof(double)); @@ -955,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; } diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp index 59496de9b3d642b013c36d394db50e5989ce9017..4332b85fcc467f82166ce8e9837146e0eaf142bd 100644 --- a/src/libnptm/tra_subs.cpp +++ b/src/libnptm/tra_subs.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -291,12 +291,9 @@ dcomplex *frfmer( double vkz = sqrt(sq - sqn); wamff(wk, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf); for (int j = 0; j < nlemt; j++) tt1->append(wk[j]); - //tt2.write(reinterpret_cast(&vkz), sizeof(double)); _vkzm[jx80][jy90] = vkz; } else { // label 50 for (int j = 0; j < nlemt; j++) tt1->append(cc0); - //double vkz = 0.0; - //tt2.write(reinterpret_cast(&vkz), sizeof(double)); _vkzm[jx80][jy90] = 0.0; } } // jx80 loop diff --git a/src/libnptm/utils.cpp b/src/libnptm/utils.cpp index f989f07565b33e18d93ed6390675468a63b41fb3..679c69622b64026aaeebe2cea8ac1ec6863c86ce 100644 --- a/src/libnptm/utils.cpp +++ b/src/libnptm/utils.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -19,27 +19,6 @@ * \brief Implementation of auxiliary code utilities. */ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - A copy of the GNU General Public License is distributed along with - this program in the COPYING file. If not, see: . - */ - -/*! \file utils.h - * - * \brief Definition of auxiliary code utilities. - */ - #include #ifndef INCLUDE_TYPES_H_ #include "../include/types.h" diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py index 4475a826aafa9a2cf053ac9c73ce94c59a7c7996..f5cda8f15652e181d1ca935d7fc0fdd7a48c920d 100755 --- a/src/scripts/model_maker.py +++ b/src/scripts/model_maker.py @@ -328,6 +328,16 @@ def load_model(model_file): for j in range(expected_radii): sconf['rcf'][i][j] = float(model['particle_settings']['rad_frac'][i][j]) # Create the gconf dict + use_refinement = False + dyn_orders = True + try: + use_refinement = False if model['system_settings']['refinement'] == "0" else True + except KeyError: + use_refinement = False + try: + dyn_orders = False if model['radiation_settings']['dyn_orders'] == "0" else True + except KeyError: + dyn_orders = True str_polar = model['radiation_settings']['polarization'] if (str_polar not in ["LINEAR", "CIRCULAR"]): print("ERROR: %s is not a recognized polarization state."%str_polar) @@ -338,6 +348,8 @@ def load_model(model_file): model['input_settings']['geometry_file'] ) } + gconf['use_refinement'] = use_refinement + gconf['dyn_orders'] = dyn_orders gconf['nsph'] = sconf['nsph'] gconf['application'] = model['particle_settings']['application'] gconf['li'] = int(model['geometry_settings']['li']) @@ -865,6 +877,12 @@ def write_legacy_gconf(conf): output.write(str_line) str_line = " {0:d}\n0\n".format(conf['jwtm']) output.write(str_line) + if (conf['use_refinement']): + str_line = "USE_REFINEMENT=1\n" + output.write(str_line) + if (not conf['dyn_orders']): + str_line = "USE_DYN_ORDER=0\n" + output.write(str_line) output.close() return result diff --git a/src/sphere/np_sphere.cpp b/src/sphere/np_sphere.cpp index e83f1edeedefdc422b30f89540fd1d0b8642ee93..ff9da9a59b23003e35f29e305a988814f18afa8d 100644 --- a/src/sphere/np_sphere.cpp +++ b/src/sphere/np_sphere.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp index 7d26091f21762d4f1edf166e1a65cf6dd7c92163..07dbd7b28f67418764a92b7c411c0c0d9296df60 100644 --- a/src/sphere/sphere.cpp +++ b/src/sphere/sphere.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -593,14 +593,21 @@ int sphere_jxi488_cycle( + ").\n"; logger->log(message, LOG_WARN); } else if (recommended_lm < max_lm) { - int new_lm = recommended_lm; - message = "INFO: lowering internal order from " + to_string(max_lm) + " to " - + to_string(recommended_lm) + " for scale iteration " + to_string(jxi488) + ".\n"; - logger->log(message, LOG_INFO); - sid->update_order(new_lm); - is_first_scale = true; - // jw = 1; - l_max = new_lm; + if (gconf->dyn_order_flag > 0) { + int new_lm = recommended_lm; + message = "INFO: lowering internal order from " + to_string(max_lm) + " to " + + to_string(recommended_lm) + " for scale iteration " + to_string(jxi488) + ".\n"; + logger->log(message, LOG_INFO); + sid->update_order(new_lm); + is_first_scale = true; + // jw = 1; + l_max = new_lm; + } else { + message = "WARNING: internal order " + to_string(max_lm) + " for scale iteration " + + to_string(jxi488) + " too high (recommended order is " + to_string(recommended_lm) + + ").\n"; + logger->log(message, LOG_WARN); + } } } // End of dynamic order check diff --git a/src/testing/test_TEDF.cpp b/src/testing/test_TEDF.cpp index eeb975221b36c2c2eb03cfce301b324a3ae0c221..804d55d79bc22d8252ae6bc477bb583a3876b065 100644 --- a/src/testing/test_TEDF.cpp +++ b/src/testing/test_TEDF.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/testing/test_TTMS.cpp b/src/testing/test_TTMS.cpp index e662c2c708101481f1e93fb9e839ab9a16479f7e..7687e63466ec9fd2bb0870ddb90f9e5352b4a643 100644 --- a/src/testing/test_TTMS.cpp +++ b/src/testing/test_TTMS.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/testing/test_file_io.cpp b/src/testing/test_file_io.cpp index 549866af0286ffc2919224754a0cdf7e48999d2e..2156c1d5b1d34d6371a87a514168bf9846a89801 100644 --- a/src/testing/test_file_io.cpp +++ b/src/testing/test_file_io.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp index 5901be17f7a617910bc6b3e4413a5bc8eaa879cb..cc5ca373e0e7c488cb161a054832c060bb767524 100644 --- a/src/trapping/cfrfme.cpp +++ b/src/trapping/cfrfme.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -18,6 +18,7 @@ * * C++ implementation of FRFME functions. */ +#include #include #include #include @@ -60,6 +61,14 @@ #include #endif +#ifdef _OPENMP +#include +#endif + +#ifdef USE_TARGET_OFFLOAD +#pragma omp requires unified_shared_memory +#endif + using namespace std; /*! \brief C++ implementation of FRFME @@ -71,6 +80,11 @@ void frfme(string data_file, string output_path) { #ifdef USE_NVTX nvtxRangePush("Running frfme()"); #endif + chrono::time_point t_start = chrono::high_resolution_clock::now(); + chrono::duration elapsed; + char buffer[256]; + string message = "INIT"; + Logger logger(LOG_INFO); string tfrfme_name = output_path + "/c_TFRFME.hd5"; TFRFME *tfrfme = NULL; Swap1 *tt1 = NULL; @@ -110,43 +124,45 @@ void frfme(string data_file, string output_path) { int nxv, nyv, nzv; if (tfrfme == NULL) tfrfme = TFRFME::from_binary(tfrfme_name, "HDF5"); if (tfrfme != NULL) { - lmode = (int)tfrfme->get_param("lmode"); - lm = (int)tfrfme->get_param("lm"); - nkv = (int)tfrfme->get_param("nkv"); - nxv = (int)tfrfme->get_param("nxv"); - nyv = (int)tfrfme->get_param("nyv"); - nzv = (int)tfrfme->get_param("nzv"); - vk = tfrfme->get_param("vk"); - exri = tfrfme->get_param("exri"); - an = tfrfme->get_param("an"); - ff = tfrfme->get_param("ff"); - tra = tfrfme->get_param("tra"); - spd = tfrfme->get_param("spd"); - frsh = tfrfme->get_param("frsh"); - exril = tfrfme->get_param("exril"); + lmode = tfrfme->lmode; + lm = tfrfme->lm; + nkv = tfrfme->nkv; + nxv = tfrfme->nxv; + nyv = tfrfme->nyv; + nzv = tfrfme->nzv; + vk = tfrfme->vk; + exri = tfrfme->exri; + an = tfrfme->an; + ff = tfrfme->ff; + tra = tfrfme->tra; + spd = tfrfme->spd; + frsh = tfrfme->frsh; + exril = tfrfme->exril; string tempname2 = output_path + "/c_TEMPTAPE2.hd5"; if (tt2 == NULL) tt2 = Swap2::from_binary(tempname2, "HDF5"); if (tt2 != NULL) { vkv = tt2->get_vector(); vkzm = tt2->get_matrix(); - apfafa = tt2->get_param("apfafa"); - pmf = tt2->get_param("pmf"); - spd = tt2->get_param("spd"); - rir = tt2->get_param("rir"); - ftcn = tt2->get_param("ftcn"); - fshmx = tt2->get_param("fshmx"); - vxyzmx = tt2->get_param("vxyzmx"); - delxyz = tt2->get_param("delxyz"); - vknmx = tt2->get_param("vknmx"); - delk = tt2->get_param("delk"); - delks = tt2->get_param("delks"); - nlmmt = (int)tt2->get_param("nlmmt"); - nrvc = (int)tt2->get_param("nrvc"); + apfafa = tt2->apfafa; + pmf = tt2->pmf; + spd = tt2->spd; + rir = tt2->rir; + ftcn = tt2->ftcn; + fshmx = tt2->fshmx; + vxyzmx = tt2->vxyzmx; + delxyz = tt2->delxyz; + vknmx = tt2->vknmx; + delk = tt2->delk; + delks = tt2->delks; + nlmmt = tt2->nlmmt; + nrvc = tt2->nrvc; } else { - printf("ERROR: could not open TEMPTAPE2 file.\n"); + message = "ERROR: could not open TEMPTAPE2 file.\n"; + logger.err(message); } } else { - printf("ERROR: could not open TFRFME file.\n"); + message = "ERROR: could not open TFRFME file.\n"; + logger.err(message); } nks = nkv - 1; #ifdef USE_NVTX @@ -272,13 +288,16 @@ void frfme(string data_file, string output_path) { // Array initialization long swap1_size, swap2_size, tfrfme_size; double size_mb; - printf("INFO: calculating memory requirements\n"); - swap1_size = Swap1::get_memory_requirement(lm, nkv); + message = "INFO: calculating memory requirements...\n"; + logger.log(message); + swap1_size = Swap1::get_size(lm, nkv); size_mb = 1.0 * swap1_size / 1024.0 / 1024.0; printf("Swap 1: %.2lg MB\n", size_mb); - swap2_size = Swap2::get_memory_requirement(nkv); + swap2_size = Swap2::get_size(nkv); size_mb = 1.0 * swap2_size / 1024.0 / 1024.0; - printf("Swap 2: %.2lg MB\n", size_mb); + sprintf(buffer, "Swap 2: %.2lg MB\n", size_mb); + message = string(buffer); + logger.log(message); tt2 = new Swap2(nkv); vkv = tt2->get_vector(); vkzm = tt2->get_matrix(); @@ -299,16 +318,19 @@ void frfme(string data_file, string output_path) { int nzs = nzsh * 2; int nzv = nzs + 1; int nzshpo = nzsh + 1; - tfrfme_size = TFRFME::get_memory_requirement(lmode, lm, nkv, nxv, nyv, nzv); + tfrfme_size = TFRFME::get_size(lm, nkv, nxv, nyv, nzv); size_mb = 1.0 * tfrfme_size / 1024.0 / 1024.0; - printf("TFRFME: %.2lg MB\n", size_mb); + sprintf(buffer, "TFRFME: %.2lg MB\n", size_mb); + message = string(buffer); + logger.log(message); size_mb = 1.0 * (swap1_size + swap2_size + tfrfme_size) / 1024.0 / 1024.0; - printf("TOTAL: %.2lg MB\n", size_mb); + sprintf(buffer, "TOTAL: %.2lg MB\n", size_mb); + message = string(buffer); + logger.log(message); tfrfme = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv); double *_xv = tfrfme->get_x(); double *_yv = tfrfme->get_y(); double *_zv = tfrfme->get_z(); - dcomplex **_wsum = tfrfme->get_matrix(); #ifdef USE_NVTX nvtxRangePop(); #endif @@ -374,54 +396,95 @@ void frfme(string data_file, string output_path) { #ifdef USE_NVTX nvtxRangePop(); #endif - dcomplex *vec_w = new dcomplex[nkv * nkv](); - dcomplex **w = new dcomplex*[nkv]; - for (int wi = 0; wi < nkv; wi++) w[wi] = vec_w + wi * nkv; #ifdef USE_NVTX nvtxRangePush("j80 loop"); #endif - for (int j80 = jlmf; j80 <= jlml; j80++) { + dcomplex *vec_wsum = tfrfme->wsum[0]; + double *vec_vkzm = vkzm[0]; +#ifdef USE_TARGET_OFFLOAD +#pragma omp target teams distribute parallel for simd +#else +#pragma omp parallel for simd +#endif +// #pragma omp parallel for + for (int j80 = jlmf-1; j80 < jlml; j80++) { + int nkvs = nkv * nkv; + dcomplex *vec_w = (dcomplex *) calloc(nkvs, sizeof(dcomplex)); + dcomplex **w = (dcomplex **) calloc(nkv, sizeof(dcomplex *)); + // dcomplex *wk_local = new dcomplex[nlmmt](); + for (int wi = 0; wi < nkv; wi++) w[wi] = vec_w + wi * nkv; + dcomplex wk_value; int wk_index = 0; - for (int jy50 = 0; jy50 < nkv; jy50++) { - for (int jx50 = 0; jx50 < nkv; jx50++) { - for (int wi = 0; wi < nlmmt; wi++) wk[wi] = tt1->vec_wk[wk_index++]; - w[jx50][jy50] = wk[j80 - 1]; - } // jx50 - } // jy50 loop + // for (int jy50 = 0; jy50 < nkv; jy50++) { + // for (int jx50 = 0; jx50 < nkv; jx50++) { +// #ifdef USE_TARGET_OFFLOAD +// #pragma omp target teams distribute parallel for simd +// #else +// #pragma omp parallel for simd +// #endif + for (int jxy50 = 0; jxy50 < nkvs; jxy50++) { + // for (int wi = 0; wi < nlmmt; wi++) wk_local[wi] = tt1->wk[wk_index++]; + // w[jx50][jy50] = wk_local[j80]; + wk_index = nlmmt*jxy50; + wk_value = tt1->wk[wk_index + j80]; + // wk_index += nlmmt; + int jy50 = jxy50 / nkv; + int jx50 = jxy50 % nkv; + vec_w[(nkv*jx50) + jy50] = wk_value; + // w[jx50][jy50] = wk_value; + } // jxy50 loop + // } // jx50 + // } // jy50 loop int ixyz = 0; - for (int wj = 0; wj < nrvc; wj++) _wsum[j80 - 1][wj] = cc0; - for (int iz75 = 0; iz75 < nzv; iz75++) { - double z = _zv[iz75] + frsh; - for (int iy70 = 0; iy70 < nyv; iy70++) { - double y = _yv[iy70]; - for (int ix65 = 0; ix65 < nxv; ix65++) { - double x = _xv[ix65]; - ixyz++; - dcomplex sumy = cc0; - for (int jy60 = 0; jy60 < nkv; jy60++) { - double vky = vkv[jy60]; - double vkx = vkv[nkv - 1]; - double vkzf = vkzm[0][jy60]; - dcomplex phasf = cexp(uim * (-vkx * x + vky * y + vkzf * z)); - double vkzl = vkzm[nkv - 1][jy60]; - dcomplex phasl = cexp(uim * (vkx * x + vky * y + vkzl * z)); - dcomplex sumx = 0.5 * (w[0][jy60] * phasf + w[nkv - 1][jy60] * phasl); - for (int jx55 = 2; jx55 <= nks; jx55++) { - vkx = vkv[jx55 - 1]; - double vkz = vkzm[jx55 - 1][jy60]; - dcomplex phas = cexp(uim * (vkx * x + vky * y + vkz * z)); - sumx += (w[jx55 - 1][jy60] * phas); - } // jx55 loop - if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5; - sumy += sumx; - } // jy60 loop - _wsum[j80 - 1][ixyz - 1] = sumy * delks; - } // ix65 loop - } // iy70 loop - } // iz75 loop + for (int wj = 0; wj < nrvc; wj++) vec_wsum[(j80*nrvc)+wj] = cc0; + int nvtot = nxv*nyv*nzv; + int nvxy = nxv *nyv; +// #ifdef USE_TARGET_OFFLOAD +// #pragma omp target teams distribute parallel for +// #else +// #pragma omp parallel for +// #endif + for (int ixyz = 0; ixyz < nvtot; ixyz++) { + int iz75 = ixyz / nvxy; + int iy70 = (ixyz % nvxy) / nxv; + int ix65 = ixyz % nxv; + // for (int iz75 = 0; iz75 < nzv; iz75++) { + double z = _zv[iz75] + frsh; + // for (int iy70 = 0; iy70 < nyv; iy70++) { + double y = _yv[iy70]; + // for (int ix65 = 0; ix65 < nxv; ix65++) { + double x = _xv[ix65]; + // ixyz++; + dcomplex sumy = cc0; +// #ifdef USE_TARGET_OFFLOAD +// #pragma omp target teams distribute parallel for simd reduction(+:sumy) +// #else +// #pragma omp parallel for simd reduction(+:sumy) +// #endif + for (int jy60x55 = 0; jy60x55 < nkvs ; jy60x55++) { + int jy60 = jy60x55 / nkv; + int jx55 = jy60x55 % nkv; + double vky = vkv[jy60]; + double vkx = (jx55==0) ? vkv[nkv - 1] : vkv[jx55]; + double vkz = vec_vkzm[(jx55*nkv)+jy60]; + dcomplex phas = (jx55==0) ? + cexp(uim * (-vkx * x + vky * y + vkz * z)): + cexp(uim * (vkx * x + vky * y + vkz * z)); + dcomplex sumx = vec_w[(jx55*nkv)+jy60] * phas; + double factor1 = ((jx55==0)||(jx55==(nkv-1))) ? 0.5 : 1.0; + double factor2 = ((jy60==0)||(jy60==(nkv-1))) ? 0.5 : 1.0; + sumx *= factor1*factor2; + sumy += sumx; + } // jy60x55 loop + vec_wsum[((j80)*nrvc)+ixyz] = sumy * delks; + // } // ix65 loop + // } // iy70 loop + // } // iz75 loop + } // ixyz loop + free(vec_w); + free(w); + // delete[] wk_local; } // j80 loop - delete[] vec_w; - delete[] w; #ifdef USE_NVTX nvtxRangePop(); #endif @@ -441,10 +504,12 @@ void frfme(string data_file, string output_path) { nvtxRangePop(); #endif } else { // Should never happen. - printf("ERROR: could not open TFRFME file for output.\n"); + message = "ERROR: could not open TFRFME file for output.\n"; + logger.err(message); } } else { - printf("ERROR: could not open TEDF file.\n"); + message = "ERROR: could not open TEDF file.\n"; + logger.err(message); } } else { // label 98 string output_name = output_path + "/c_OFRFME"; @@ -468,8 +533,11 @@ void frfme(string data_file, string output_path) { #ifdef USE_NVTX nvtxRangePop(); #endif - printf("FRFME: Done.\n"); + elapsed = chrono::high_resolution_clock::now() - t_start; + message = "INFO: FRFME took " + to_string(elapsed.count()) + "s.\n"; + logger.log(message); #ifdef USE_NVTX nvtxRangePop(); #endif } + diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp index 4d66d8bad883686c89c169c7afb90641735f1729..39d4a3a3760c4d43b1d5ece1f0d4032312e82730 100644 --- a/src/trapping/clffft.cpp +++ b/src/trapping/clffft.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -18,6 +18,7 @@ * * \brief C++ implementation of LFFFT functions. */ +#include #include #include #include @@ -71,10 +72,15 @@ void lffft(string data_file, string output_path) { #ifdef USE_NVTX nvtxRangePush("Running lffft()"); #endif + chrono::time_point t_start = chrono::high_resolution_clock::now(); + chrono::duration elapsed; const dcomplex uim = 0.0 + 1.0 * I; const double sq2i = 1.0 / sqrt(2.0); const dcomplex sq2iti = sq2i * uim; + char buffer[256]; + string message = "INIT"; + Logger logger(LOG_INFO); fstream tlfff, tlfft; double ****zpv = NULL; dcomplex *ac = NULL, *ws = NULL, *wsl = NULL; @@ -82,7 +88,6 @@ void lffft(string data_file, string output_path) { dcomplex **amd = NULL; int **indam = NULL; dcomplex *tmsm = NULL, *tmse = NULL, **tms = NULL; - dcomplex **_wsum = NULL; int jft, jss, jtw; int is, le, nvam = 0; double vks, exris; @@ -188,26 +193,25 @@ void lffft(string data_file, string output_path) { int lmode, lm, nkv, nxv, nyv, nzv; double vk, exri, an, ff, tra; double spd, frsh, exril; - lmode = (int)tfrfme->get_param("lmode"); - lm = (int)tfrfme->get_param("lm"); - nkv = (int)tfrfme->get_param("nkv"); - nxv = (int)tfrfme->get_param("nxv"); - nyv = (int)tfrfme->get_param("nyv"); - nzv = (int)tfrfme->get_param("nzv"); - _wsum = tfrfme->get_matrix(); + lmode = tfrfme->lmode; + lm = tfrfme->lm; + nkv = tfrfme->nkv; + nxv = tfrfme->nxv; + nyv = tfrfme->nyv; + nzv = tfrfme->nzv; double *_xv = tfrfme->get_x(); double *_yv = tfrfme->get_y(); double *_zv = tfrfme->get_z(); if (lm >= le) { - vk = tfrfme->get_param("vk"); - exri = tfrfme->get_param("exri"); - an = tfrfme->get_param("an"); - ff = tfrfme->get_param("ff"); - tra = tfrfme->get_param("tra"); + vk = tfrfme->vk; + exri = tfrfme->exri; + an = tfrfme->an; + ff = tfrfme->ff; + tra = tfrfme->tra; if (vk == vks && exri == exris) { - spd = tfrfme->get_param("spd"); - frsh = tfrfme->get_param("frsh"); - exril = tfrfme->get_param("exril"); + spd = tfrfme->spd; + frsh = tfrfme->frsh; + exril = tfrfme->exril; bool goto160 = false; if (jft <= 0) { zpv = new double***[le]; @@ -316,7 +320,7 @@ void lffft(string data_file, string output_path) { //binary_input.read(reinterpret_cast(&vimag), sizeof(double)); int row = i; int col = (nyv * nxv * iz475) + (nxv * iy475) + ix475; - dcomplex value = _wsum[row][col]; + dcomplex value = tfrfme->wsum[row][col]; if (lm <= le) { ws[i] = value; } else { // label 170 @@ -482,7 +486,9 @@ void lffft(string data_file, string output_path) { delete cil; delete ccr; delete[] file_lines; - printf("LFFT: Done.\n"); + elapsed = chrono::high_resolution_clock::now() - t_start; + message = "INFO: LFFT took " + to_string(elapsed.count()) + "s.\n"; + logger.log(message); #ifdef USE_NVTX nvtxRangePop(); #endif diff --git a/src/trapping/np_trapping.cpp b/src/trapping/np_trapping.cpp index 99b1f7898e900cd852113f427b926152cce08c91..053fa6abf5e7476a3385609987804372282ccfd3 100644 --- a/src/trapping/np_trapping.cpp +++ b/src/trapping/np_trapping.cpp @@ -1,10 +1,31 @@ +/* Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + A copy of the GNU General Public License is distributed along with + this program in the COPYING file. If not, see: . + */ + /*! \file np_trapping.cpp * * \brief Trapping problem handler. */ +#include #include #include +#ifndef INCLUDE_LOGGING_H_ +#include "../include/logging.h" +#endif + using namespace std; extern void frfme(string data_file, string output_path); @@ -17,9 +38,13 @@ extern void lffft(string data_file, string output_path); * \return result: `int` */ int main(int argc, char **argv) { + chrono::time_point t_start = chrono::high_resolution_clock::now(); + chrono::duration elapsed; string frfme_data_file = "../../test_data/trapping/DFRFME"; string lffft_data_file = "../../test_data/trapping/DLFFFT"; string output_path = "."; + string message; + Logger logger(LOG_DEBG); if (argc == 4) { frfme_data_file = string(argv[1]); lffft_data_file = string(argv[2]); @@ -27,5 +52,8 @@ int main(int argc, char **argv) { } frfme(frfme_data_file, output_path); lffft(lffft_data_file, output_path); + elapsed = chrono::high_resolution_clock::now() - t_start; + message = "INFO: calculation lasted " + to_string(elapsed.count()) + "s.\n"; + logger.log(message); return 0; } diff --git a/test_data/cluster/DCLU b/test_data/cluster/DCLU index dd6d963fb87a48906e1958603fcd6d45a18c4471..f819ca7377a0dbba0071e92e430d794adb71510d 100644 --- a/test_data/cluster/DCLU +++ b/test_data/cluster/DCLU @@ -7,3 +7,4 @@ 3.40D01 1.50D01 4.90D01 5.00D00 5.00D00 1.00D01 0 0 +# This is a comment line diff --git a/test_data/cluster/DCLU_48 b/test_data/cluster/DCLU_48 index 6a4dba14577f0ec92d13346bf036cf1557860a41..0a6975bb9153a9fb84774ad3d43bbb880eb38dfe 100644 --- a/test_data/cluster/DCLU_48 +++ b/test_data/cluster/DCLU_48 @@ -51,3 +51,5 @@ 0.00D01 0.00D00 0.00D01 0.00D00 0.00D01 0.00D01 1 0 +# COMMENT: next line enables iterative refinement at run-time +USE_REFINEMENT=1