Merge pull request #460 from toxa81/develop

debug floating-point exceptions

Author: toxa81

apps/dft_loop/sirius.scf.cpp

@@ -1,6 +1,8 @@
 #include "utils/profiler.hpp"
 #include <sirius.h>
 #include <utils/json.hpp>
+#include <cfenv>
+#include <fenv.h>
 
 using namespace sirius;
 using json = nlohmann::json;
@@ -334,13 +336,15 @@ void run_tasks(cmd_args const& args)
 
 int main(int argn, char** argv)
 {
+    std::feclearexcept(FE_ALL_EXCEPT);
     cmd_args args;
     args.register_key("--input=", "{string} input file name");
     args.register_key("--output=", "{string} output file name");
     args.register_key("--task=", "{int} task id");
     args.register_key("--aiida_output", "write output for AiiDA");
     args.register_key("--test_against=", "{string} json file with reference values");
     args.register_key("--repeat_update=", "{int} number of times to repeat update()");
+    args.register_key("--fpe", "enable check of floating-point exceptions using GNUC library");
     args.register_key("--control.processing_unit=", "");
     args.register_key("--control.verbosity=", "");
     args.register_key("--control.verification=", "");
@@ -361,6 +365,12 @@ int main(int argn, char** argv)
         return 0;
     }
 
+#if defined(_GNU_SOURCE)
+    if (args.exist("fpe")) {
+        feenableexcept(FE_DIVBYZERO|FE_INVALID|FE_OVERFLOW);
+    }
+#endif
+
     sirius::initialize(1);
 
     run_tasks(args);
@@ -375,6 +385,18 @@ int main(int argn, char** argv)
         std::ofstream ofs("timers.json", std::ofstream::out | std::ofstream::trunc);
         ofs << timing_result.json();
     }
+    if (std::fetestexcept(FE_DIVBYZERO)) {
+        std::cout << "FE_DIVBYZERO exception\n";
+    }
+    if (std::fetestexcept(FE_INVALID)) {
+        std::cout << "FE_INVALID exception\n";
+    }
+    if (std::fetestexcept(FE_UNDERFLOW)) {
+        std::cout << "FE_UNDERFLOW exception\n";
+    }
+    if (std::fetestexcept(FE_OVERFLOW)) {
+        std::cout << "FE_OVERFLOW exception\n";
+    }
 
     return 0;
 }

src/Potential/poisson.cpp

@@ -103,34 +103,38 @@ void Potential::poisson_add_pseudo_pw(mdarray<double_complex, 2>& qmt__,
             }
         }
 
+        double fourpi_omega = fourpi / unit_cell_.omega();
+
         /* add pseudo_density to interstitial charge density so that rho(G) has the correct
          * multipole moments in the muffin-tins */
         #pragma omp parallel for schedule(static)
-        for (int igloc = 0; igloc < ctx_.gvec().count(); igloc++) {
+        for (int igloc = ctx_.gvec().skip_g0(); igloc < ctx_.gvec().count(); igloc++) {
             int ig = ctx_.gvec().offset() + igloc;
 
             double gR = ctx_.gvec().gvec_len(ig) * R;
             double gRn = std::pow(2.0 / gR, pseudo_density_order_ + 1);
 
             double_complex rho_G(0, 0);
-            if (ig) { // G!=0
-                /* loop over atoms of the same type */
-                for (int l = 0, lm = 0; l <= ctx_.lmax_rho(); l++) {
-                    double_complex zt1(0, 0);
-                    for (int m = -l; m <= l; m++, lm++) {
-                        zt1 += gvec_ylm_(lm, igloc) * qapf(lm, igloc);
-                    }
-                    rho_G += (fourpi / unit_cell_.omega()) * std::conj(zil_[l]) * zt1 * gamma_factors_R_(l, iat) *
-                             sbessel_mt_(l + pseudo_density_order_ + 1, igloc, iat) * gRn;
-                } // l
-            } else { // G=0
-                for (int i = 0; i < unit_cell_.atom_type(iat).num_atoms(); i++) {
-                    int ia = unit_cell_.atom_type(iat).atom_id(i);
-                    rho_G += (fourpi / unit_cell_.omega()) * y00 * (qmt__(0, ia) - qit__(0, ia));
+            /* loop over atoms of the same type */
+            for (int l = 0, lm = 0; l <= ctx_.lmax_rho(); l++) {
+                double_complex zt1(0, 0);
+                for (int m = -l; m <= l; m++, lm++) {
+                    zt1 += gvec_ylm_(lm, igloc) * qapf(lm, igloc);
                 }
+                rho_G += fourpi_omega * std::conj(zil_[l]) * zt1 * gamma_factors_R_(l, iat) *
+                         sbessel_mt_(l + pseudo_density_order_ + 1, igloc, iat) * gRn;
             }
             rho_pw__[igloc] += rho_G;
         }
+        /* for G=0 case */
+        if (ctx_.comm().rank() == 0) {
+            double_complex rho_G(0, 0);
+            for (int i = 0; i < unit_cell_.atom_type(iat).num_atoms(); i++) {
+                int ia = unit_cell_.atom_type(iat).atom_id(i);
+                rho_G += fourpi_omega * y00 * (qmt__(0, ia) - qit__(0, ia));
+            }
+            rho_pw__[0] += rho_G;
+        }
     }
 }