diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 7cc4d636f..390225f55 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -26,7 +26,7 @@ void selfadjointeigensolver_essential_check(const MatrixType& m) {
   VERIFY_IS_EQUAL(eiSymm.info(), Success);
 
   RealScalar scaling = m.cwiseAbs().maxCoeff();
-  RealScalar unitary_error_factor = RealScalar(16);
+  RealScalar unitary_error_factor = RealScalar(32);
 
   if (scaling < (std::numeric_limits<RealScalar>::min)()) {
     VERIFY(eiSymm.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
@@ -80,7 +80,7 @@ void selfadjointeigensolver(const MatrixType& m) {
   MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
   MatrixType symmC = symmA;
 
-  svd_fill_random(symmA, Symmetric);
+  svd_fill_random(symmA, SelfAdjoint);
 
   symmA.template triangularView<StrictlyUpper>().setZero();
   symmC.template triangularView<StrictlyUpper>().setZero();
diff --git a/test/product.h b/test/product.h
index 645d40aa9..00e67efe3 100644
--- a/test/product.h
+++ b/test/product.h
@@ -111,8 +111,12 @@ void product(const MatrixType& m) {
   VERIFY_IS_APPROX(m3, m1 * (m1.transpose() * m2));
 
   // continue testing Product.h: distributivity
-  VERIFY_IS_APPROX(square * (m1 + m2), square * m1 + square * m2);
-  VERIFY_IS_APPROX(square * (m1 - m2), square * m1 - square * m2);
+  {
+    // Increase tolerance, since coefficients here can get relatively large.
+    RealScalar tol = RealScalar(2) * get_test_precision(m1);
+    VERIFY(verifyIsApprox(square * (m1 + m2), square * m1 + square * m2, tol));
+    VERIFY(verifyIsApprox(square * (m1 - m2), square * m1 - square * m2, tol));
+  }
 
   // continue testing Product.h: compatibility with ScalarMultiple.h
   VERIFY_IS_APPROX(s1 * (square * m1), (s1 * square) * m1);
@@ -130,7 +134,10 @@ void product(const MatrixType& m) {
 
   // test the previous tests were not screwed up because operator* returns 0
   // (we use the more accurate default epsilon)
-  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
+  // Skip non-commutativity check for very low-precision types (e.g. bfloat16) where
+  // random small matrices can produce products that are approximately equal.
+  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1 &&
+      NumTraits<RealScalar>::dummy_precision() < RealScalar(0.04)) {
     VERIFY(areNotApprox(m1.transpose() * m2, m2.transpose() * m1, not_approx_epsilon));
   }
 
@@ -138,7 +145,8 @@ void product(const MatrixType& m) {
   res = square;
   res.noalias() += m1 * m2.transpose();
   VERIFY_IS_APPROX(res, square + m1 * m2.transpose());
-  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
+  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1 &&
+      NumTraits<RealScalar>::dummy_precision() < RealScalar(0.04)) {
     VERIFY(areNotApprox(res, square + m2 * m1.transpose(), not_approx_epsilon));
   }
   vcres = vc2;
@@ -149,7 +157,8 @@ void product(const MatrixType& m) {
   res = square;
   res.noalias() -= m1 * m2.transpose();
   VERIFY_IS_APPROX(res, square - (m1 * m2.transpose()));
-  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
+  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1 &&
+      NumTraits<RealScalar>::dummy_precision() < RealScalar(0.04)) {
     VERIFY(areNotApprox(res, square - m2 * m1.transpose(), not_approx_epsilon));
   }
   vcres = vc2;
@@ -183,21 +192,35 @@ void product(const MatrixType& m) {
   res.noalias() -= square - m1 * m2.transpose();
   VERIFY_IS_APPROX(res, square - m1 * m2.transpose());
 
-  tm1 = m1;
-  VERIFY_IS_APPROX(tm1.transpose() * v1, m1.transpose() * v1);
-  VERIFY_IS_APPROX(v1.transpose() * tm1, v1.transpose() * m1);
+  // Row-major vs col-major GEMV: different accumulation orders produce
+  // rounding differences bounded by O(sqrt(n)) * epsilon per inner product.
+  // tm1.transpose() * v1 has inner dimension rows; v1^T * tm1 has inner dim cols.
+  {
+    RealScalar gemv_tol = (std::max)(get_test_precision(m1),
+                                     numext::sqrt(RealScalar((std::max)(rows, cols))) * NumTraits<Scalar>::epsilon());
+    tm1 = m1;
+    VERIFY(verifyIsApprox(tm1.transpose() * v1, m1.transpose() * v1, gemv_tol));
+    VERIFY(verifyIsApprox(v1.transpose() * tm1, v1.transpose() * m1, gemv_tol));
+  }
 
   // test submatrix and matrix/vector product
-  for (int i = 0; i < rows; ++i) res.row(i) = m1.row(i) * m2.transpose();
-  VERIFY_IS_APPROX(res, m1 * m2.transpose());
-  // the other way round:
-  for (int i = 0; i < rows; ++i) res.col(i) = m1 * m2.transpose().col(i);
-  VERIFY_IS_APPROX(res, m1 * m2.transpose());
+  // Row-by-row vs full GEMM: different evaluation strategies can differ
+  // by O(sqrt(n)) * epsilon for low-precision types.
+  {
+    RealScalar prod_tol =
+        (std::max)(get_test_precision(m1), numext::sqrt(RealScalar(cols)) * NumTraits<Scalar>::epsilon());
+    for (int i = 0; i < rows; ++i) res.row(i) = m1.row(i) * m2.transpose();
+    VERIFY(verifyIsApprox(res, m1 * m2.transpose(), prod_tol));
+    // the other way round:
+    for (int i = 0; i < rows; ++i) res.col(i) = m1 * m2.transpose().col(i);
+    VERIFY(verifyIsApprox(res, m1 * m2.transpose(), prod_tol));
+  }
 
   res2 = square2;
   res2.noalias() += m1.transpose() * m2;
   VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2);
-  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
+  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1 &&
+      NumTraits<RealScalar>::dummy_precision() < RealScalar(0.04)) {
     VERIFY(areNotApprox(res2, square2 + m2.transpose() * m1, not_approx_epsilon));
   }