Kahan 2x2 determinant computation reduces numerical imprecision.

src/fixedSizeSquareMatrixOpsImpl.hpp

@@ -88,7 +88,10 @@ struct SquareMatrixOps< 2 >
   static auto determinant( MATRIX const & matrix )
   {
     checkSizes< 2, 2 >( matrix );
-    return matrix[ 0 ][ 0 ] * matrix[ 1 ][ 1 ] - matrix[ 0 ][ 1 ] * matrix[ 1 ][ 0 ];
+
+    // Using the more precise Kahan method to compute the 2x2 determinant.
+    auto const tmp = matrix[0][1] * matrix[1][0];
+    return math::fma( -matrix[0][1], matrix[1][0], tmp ) + math::fma( matrix[0][0], matrix[1][1], -tmp );
   }
 
   /**

src/math.hpp

@@ -28,7 +28,7 @@ namespace LvArray
 {
 
 /**
- * @brief Contains protable wrappers around cmath functions and some cuda specific functions.
+ * @brief Contains portable wrappers around cmath functions and some cuda specific functions.
  */
 namespace math
 {
@@ -317,6 +317,46 @@ max( T const a, T const b )
 #endif
 }
 
+/**
+ * @return Returns @p x * @p y + @p z.
+ * @tparam T A floating point type.
+ * @param x The first number.
+ * @param y The second number.
+ * @param z The third number.
+ * @note fma stands for fused multiply add.
+ *
+ * The function computes the result without losing precision in any intermediate result.
+ * Integer arguments cast to double.
+ * We want to avoid that by providing a version dedicated to integers.
+ */
+template< typename T >
+LVARRAY_HOST_DEVICE inline constexpr
+std::enable_if_t< std::is_floating_point< T >::value, T >
+fma( T const x, T const y, T const z )
+{
+#if defined(__CUDA_ARCH__)
+  return ::fma( x, y, z );
+#else
+  return std::fma( x, y, z );
+#endif
+}
+
+/**
+ * @return Returns @p x * @p y + @p z.
+ * @tparam T A floating point type.
+ * @param x The first number.
+ * @param y The second number.
+ * @param z The third number.
+ * @note This is a dummy implementation for integers (in order to not cast to doubles).
+ */
+template< typename T >
+LVARRAY_HOST_DEVICE inline constexpr
+std::enable_if_t< std::is_integral< T >::value, T >
+fma( T const x, T const y, T const z )
+{
+  return x * y + z;
+}
+
 #if defined( LVARRAY_USE_CUDA )
 
 /// @copydoc max( T, T )