From da24d16c669e11b81bfcdbeb9a781170c3db6221 Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Tue, 9 Jan 2024 13:35:22 +0100 Subject: [PATCH 1/2] [MathLib] Avoid the division by zero in componentwiseDivide --- MathLib/LinAlg/LinAlg.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/MathLib/LinAlg/LinAlg.cpp b/MathLib/LinAlg/LinAlg.cpp index 0e1ae63c312..a06ab49642f 100644 --- a/MathLib/LinAlg/LinAlg.cpp +++ b/MathLib/LinAlg/LinAlg.cpp @@ -230,8 +230,9 @@ template <> void componentwiseDivide(EigenVector& w, EigenVector const& x, EigenVector const& y) { - w.getRawVector().noalias() = - x.getRawVector().cwiseQuotient(y.getRawVector()); + w.getRawVector().noalias() = x.getRawVector().binaryExpr( + y.getRawVector(), + [](auto const x, auto const y) { return y == 0 ? 0.0 : x / y; }); } // Explicit specialization From 2d6b4db7acc297eb6cd36525f57a3aed6e8445f3 Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Tue, 9 Jan 2024 13:42:41 +0100 Subject: [PATCH 2/2] [MathLib] Added a comment to PETSC::VecPointwiseDivide, which is called by componentwiseDivide --- MathLib/LinAlg/LinAlg.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/MathLib/LinAlg/LinAlg.cpp b/MathLib/LinAlg/LinAlg.cpp index a06ab49642f..5e688a3fd1e 100644 --- a/MathLib/LinAlg/LinAlg.cpp +++ b/MathLib/LinAlg/LinAlg.cpp @@ -70,6 +70,13 @@ void axpby(PETScVector& y, PetscScalar const a, PetscScalar const b, // Explicit specialization // Computes w = x/y componentwise. +// \note that VecPointwiseDivide avoids to divide by values that are +// identically zero such as +// for (int i=0; i void componentwiseDivide(PETScVector& w, PETScVector const& x, PETScVector const& y)