Merge pull request #12905 from brian-kelley/Fix12898

Tpetra: Fix #12898, Teuchos: absolute float array compare

packages/teuchos/core/src/Teuchos_LocalTestingHelpers.hpp

@@ -164,7 +164,7 @@
 
 /** \brief Assert that a1.size()==a2.size() and a[i]==b[i], i=0....
  *
- * Works for any object types that support a1[i], a1.size(), a2[j], and
+ * Works for all object types that support a1[i], a1.size(), a2[j], and
  * a2.size() and types a1 and a2 can be different types!
  *
  * \ingroup Teuchos_UnitTestAssertMacros_grp
@@ -178,7 +178,7 @@
 
 /** \brief Assert that a1.size()==a2.size() and rel_error(a[i],b[i]) <= tol, i=0....
  *
- * Works for any object types that support a1[i], a1.size(), a2[j], and
+ * Works for all object types that support a1[i], a1.size(), a2[j], and
  * a2.size() and types a1 and a2 can be different types!
  *
  * \ingroup Teuchos_UnitTestAssertMacros_grp
@@ -189,6 +189,19 @@
     if (!result) success = false; \
   }
 
+/** \brief Assert that a1.size()==a2.size() and |a[i]-b[i]| <= tol, i=0....
+ *
+ * Works for all object types that support a1[i], a1.size(), a2[j], and
+ * a2.size(), but the element types of a1 and a2 must be the same and Teuchos::ScalarTraits
+ * must have a specialization for this element type.
+ *
+ * \ingroup Teuchos_UnitTestAssertMacros_grp
+ */
+#define TEST_ABSOLUTE_COMPARE_FLOATING_ARRAYS( a1, a2, tol ) \
+  { \
+    const bool result = compareFloatingArraysAbsolute(a1,#a1,a2,#a2,tol,out); \
+    if (!result) success = false; \
+  }
 
 /** \brief Assert that the statement 'code' throws the exception 'ExceptType'
  * (otherwise the test fails).

packages/teuchos/core/src/Teuchos_TestingHelpers.hpp

@@ -221,12 +221,14 @@ bool testRelErr(
 
 /** \brief Compare if two array objects are the same or not.
  *
- * This function works with any two array objects are the same size and have
- * the same element value types.  The funtion is templated on the container
- * types and therefore can compare any two objects that have size() and
- * operator[](i) defined.
+ * This function works with any two array-like objects:
+ * <tt>Array1</tt> and <tt>Array2</tt> must have <tt>size()</tt>
+ * and <tt>operator[](i)</tt> methods defined. Their element types
+ * may be different, but it must be possible to compare them with
+ * <tt>==</tt>.
  *
- * \returns Returns <tt>true</tt> if the compare and <tt>false</tt> otherwise.
+ * \returns Returns <tt>true</tt> if <tt>a1.size() == a2.size()</tt> and
+ *   their elements are equal. Otherwise returns <tt>false</tt>.
  *
  * \ingroup teuchos_testing_grp
  */
@@ -241,12 +243,15 @@ bool compareArrays(
 /** \brief Compare if two array objects are the same or not up to a relative
  * floating point precision.
  *
- * This function works with any two array objects are the same size and have
- * the same element value types.  The funtion is templated on the container
- * types and therefore can compare any two objects that have size() and
- * operator[](i) defined.
+ * This function works with any two array-like objects:
+ * <tt>Array1</tt> and <tt>Array2</tt> must have <tt>size()</tt>
+ * and <tt>operator[](i)</tt> methods defined. Their element
+ * types may be different, as long as <tt>Teuchos::relErr(a1[i], a2[i])</tt>
+ * can be called to determine the relative difference between them.
  *
- * \returns Returns <tt>true</tt> if the compare and <tt>false</tt> otherwise.
+ * \returns Returns <tt>true</tt> if <tt>a1.size() == a2.size()</tt> and
+ *   <tt>relErr(a1[i], a2[i]) <= tol</tt> for all <tt>0 <= i < a1.size()</tt>.
+ *   Otherwise returns <tt>false</tt>.
  *
  * \ingroup teuchos_testing_grp
  */
@@ -258,6 +263,30 @@ bool compareFloatingArrays(
   Teuchos::FancyOStream &out
   );
 
+/** \brief Compare if two array objects are the same up to an absolute
+ * tolerance: elements <tt>a1[i]</tt> and <tt>a2[i]</tt> are considered
+ * the same if <tt>|a1[i]-a2[i]| <= tol</tt>.
+ *
+ * This function works with any two array-like objects
+ * with the same element types. <tt>Array1</tt> and <tt>Array2</tt> must have
+ * <tt>size()</tt> and <tt>operator[](i)</tt> methods defined.
+ * <tt>Teuchos::ScalarTraits</tt> must also have a specialization
+ * for the element type.
+ *
+ * \returns Returns <tt>true</tt> if <tt>a1.size() == a2.size()</tt> and
+ *   <tt>|a1[i]-a2[i]| <= tol</tt> for all <tt>0 <= i < a1.size()</tt>.
+ *   Otherwise returns <tt>false</tt>.
+ *
+ * \ingroup teuchos_testing_grp
+ */
+template<class Array1, class Array2, class ScalarMag>
+bool compareFloatingArraysAbsolute(
+  const Array1 &a1, const std::string &a1_name,
+  const Array2 &a2, const std::string &a2_name,
+  const ScalarMag &tol,
+  Teuchos::FancyOStream &out
+  );
+
 
 } // namespace Teuchos
 
@@ -654,6 +683,7 @@ bool Teuchos::compareArrays(
   )
 {
   using Teuchos::as;
+
   bool success = true;
 
   out << "Comparing " << a1_name << " == " << a2_name << " ... ";
@@ -694,6 +724,14 @@ bool Teuchos::compareFloatingArrays(
   )
 {
   using Teuchos::as;
+
+  // Determine the element types of Array1 and Array2
+  using Elem1 = std::decay_t<decltype(std::declval<Array1>()[0])>;
+  using Elem2 = std::decay_t<decltype(std::declval<Array2>()[0])>;
+
+  static_assert(std::is_same_v<Elem1, Elem2>,
+      "Teuchos::compareFloatingArrays: element types of Array1 and Array2 must be the same.");
+
   bool success = true;
 
   out << "Comparing " << a1_name << " == " << a2_name << " ... ";
@@ -726,5 +764,52 @@ bool Teuchos::compareFloatingArrays(
 
 }
 
+template<class Array1, class Array2, class ScalarMag>
+bool Teuchos::compareFloatingArraysAbsolute(
+  const Array1 &a1, const std::string &a1_name,
+  const Array2 &a2, const std::string &a2_name,
+  const ScalarMag &tol,
+  Teuchos::FancyOStream &out
+  )
+{
+  using Teuchos::as;
+  using Teuchos::ScalarTraits;
+
+  // Determine the element types of Array1 and Array2
+  using Elem1 = std::decay_t<decltype(std::declval<Array1>()[0])>;
+  using Elem2 = std::decay_t<decltype(std::declval<Array2>()[0])>;
+
+  static_assert(std::is_same_v<Elem1, Elem2>,
+      "Teuchos::compareFloatingArraysAbsolute: element types of Array1 and Array2 must be the same.");
+
+  bool success = true;
+
+  out << "Comparing " << a1_name << " == " << a2_name << " ... ";
+
+  const int n = a1.size();
+
+  // Compare sizes
+  if (as<int>(a2.size()) != n) {
+    out << "\nError, "<<a1_name<<".size() = "<<a1.size()<<" == "
+        << a2_name<<".size() = "<<a2.size()<<" : failed!\n";
+    return false;
+  }
+
+  // Compare elements
+  for( int i = 0; i < n; ++i ) {
+    const ScalarMag err = ScalarTraits<Elem1>::magnitude( a1[i] - a2[i] );
+    if ( !(err <= tol) ) {
+      out
+        <<"\nError, ||"<<a1_name<<"["<<i<<"] - " << a2_name<<"["<<i<<"]|| = "
+        <<err<<" <= tol = "<<tol<<": failed!\n";
+      success = false;
+    }
+  }
+  if (success) {
+    out << "passed\n";
+  }
+  return success;
+}
+
 
 #endif  // TEUCHOS_TESTING_HELPERS_HPP

packages/tpetra/core/test/CrsMatrix/CrsMatrix_UnitTests4.cpp

@@ -50,6 +50,7 @@
 #include "Tpetra_Details_scaleBlockDiagonal.hpp"
 #include "Tpetra_Apply_Helpers.hpp"
 #include "TpetraUtils_MatrixGenerator.hpp"
+#include "KokkosBlas1_nrm2.hpp"
 #include <type_traits> // std::is_same
 
 
@@ -1060,18 +1061,28 @@ inline void tupleToArray(Array<T> &arr, const tuple &tup)
     auto lclMtx1 = A1->getLocalMatrixDevice();
     values_type A1_values = lclMtx1.values;
     values_type A2_values("values",A1_values.extent(0));
+    values_type A2_abs_values("absolute values",A1_values.extent(0));
     Kokkos::parallel_for( range_policy(0,A1_values.extent(0)),KOKKOS_LAMBDA(const int i){
         A2_values[i] = A1_values[i] + A1_values[i];
+        A2_abs_values[i] = Kokkos::abs(A2_values[i]);
       });
 
     RCP<MAT> A2 = rcp(new MAT(map,A1->getColMap(),lclMtx1.graph.row_map,lclMtx1.graph.entries,A2_values));
     A2->expertStaticFillComplete(A1->getDomainMap(),A1->getRangeMap(),A1->getGraph()->getImporter(),A1->getGraph()->getExporter());
+    // Entrywise absolute value of A2, used for choosing a tolerance
+    RCP<MAT> A2_abs = rcp(new MAT(map,A1->getColMap(),lclMtx1.graph.row_map,lclMtx1.graph.entries,A2_abs_values));
+    A2_abs->expertStaticFillComplete(A1->getDomainMap(),A1->getRangeMap(),A1->getGraph()->getImporter(),A1->getGraph()->getExporter());
 
     /* Allocate multivectors */
     MV X(map,2), Y1a(map,2), Y1b(map,2), Y2a(map,2), Y2b(map,2), compare(map,2);
+    MV Y_abs(map,2);
     Array<Mag> norm(2), exact(2,MT_ZERO);
     X.putScalar(ONE);
 
+    A2_abs->apply(X, Y_abs);
+    // both columns of X are filled with 1, so just use the first column here
+    Mag absNorm = Y_abs.getVector(0)->norm1();
+
     // Do a std::vector version 
     {
       Y1a.putScalar(ZERO);Y1b.putScalar(ZERO);
@@ -1090,13 +1101,14 @@ inline void tupleToArray(Array<T> &arr, const tuple &tup)
       // Compare
       compare.update(ONE,Y1a,-ONE,Y1b,ZERO);
       compare.norm1(norm());
+
       if (ST::isOrdinal) {TEST_COMPARE_ARRAYS(exact,norm);}
-      else{TEST_COMPARE_FLOATING_ARRAYS(exact,norm,2.0*testingTol<Mag>());}
+      else{TEST_ABSOLUTE_COMPARE_FLOATING_ARRAYS(exact,norm,absNorm*testingTol<Mag>());}
 
       compare.update(ONE,Y2a,-ONE,Y2b,ZERO);
       compare.norm1(norm());
       if (ST::isOrdinal) {TEST_COMPARE_ARRAYS(exact,norm);}
-      else {TEST_COMPARE_FLOATING_ARRAYS(exact,norm,2.0*testingTol<Mag>());}
+      else {TEST_ABSOLUTE_COMPARE_FLOATING_ARRAYS(exact,norm,absNorm*testingTol<Mag>());}
     }
 
     // Do a std::vector version w/ a nice combination of options
@@ -1127,12 +1139,12 @@ inline void tupleToArray(Array<T> &arr, const tuple &tup)
         compare.update(ONE,Y1a,-ONE,Y1b,ZERO);
         compare.norm1(norm());
         if (ST::isOrdinal) {TEST_COMPARE_ARRAYS(exact,norm);}
-        else{TEST_COMPARE_FLOATING_ARRAYS(exact,norm,2.0*testingTol<Mag>());}
+        else{TEST_ABSOLUTE_COMPARE_FLOATING_ARRAYS(exact,norm,absNorm*testingTol<Mag>());}
 
         compare.update(ONE,Y2a,-ONE,Y2b,ZERO);
         compare.norm1(norm());
         if (ST::isOrdinal) {TEST_COMPARE_ARRAYS(exact,norm);}
-        else {TEST_COMPARE_FLOATING_ARRAYS(exact,norm,2.0*testingTol<Mag>());}
+        else {TEST_ABSOLUTE_COMPARE_FLOATING_ARRAYS(exact,norm,absNorm*testingTol<Mag>());}
       }
     }
 
@@ -1156,12 +1168,12 @@ inline void tupleToArray(Array<T> &arr, const tuple &tup)
       compare.update(ONE,Y1a,-ONE,Y1b,ZERO);
       compare.norm1(norm());
       if (ST::isOrdinal) {TEST_COMPARE_ARRAYS(exact,norm);}
-      else{TEST_COMPARE_FLOATING_ARRAYS(exact,norm,2.0*testingTol<Mag>());}
+      else{TEST_ABSOLUTE_COMPARE_FLOATING_ARRAYS(exact,norm,absNorm*testingTol<Mag>());}
 
       compare.update(ONE,Y2a,-ONE,Y2b,ZERO);
       compare.norm1(norm());
       if (ST::isOrdinal) {TEST_COMPARE_ARRAYS(exact,norm);}
-      else {TEST_COMPARE_FLOATING_ARRAYS(exact,norm,2.0*testingTol<Mag>());}
+      else {TEST_ABSOLUTE_COMPARE_FLOATING_ARRAYS(exact,norm,absNorm*testingTol<Mag>());}
     }
 
   }