Update Eigen to commit:954e21152e204b1960aca802eb9f16d054d70fd9

CHANGELOG
=========
954e21152 - Include <limits> in test main.h
e15cd620a - Remove select class
1c0048a08 - Fix inconsistency between ptrue and pcmp_* in HVX
ddce1d7d1 - Fixes #2952
8b9dbcdaa - Fix numext::bit_cast() compilation failure in C++20
975a5aba4 - Fix TODO: Use std::bit_cast or __builtin_bit_cast if available.

PiperOrigin-RevId: 794397153
Change-Id: I6acf902aeeafe4f69dd9fafefa9d823615d01970
diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h
index ec731ac..c9b2d2d 100644
--- a/Eigen/src/Core/CoreEvaluators.h
+++ b/Eigen/src/Core/CoreEvaluators.h
@@ -1565,50 +1565,6 @@
   }
 };
 
-// -------------------- Select --------------------
-// NOTE shall we introduce a ternary_evaluator?
-
-// TODO enable vectorization for Select
-template <typename ConditionMatrixType, typename ThenMatrixType, typename ElseMatrixType>
-struct evaluator<Select<ConditionMatrixType, ThenMatrixType, ElseMatrixType>>
-    : evaluator_base<Select<ConditionMatrixType, ThenMatrixType, ElseMatrixType>> {
-  typedef Select<ConditionMatrixType, ThenMatrixType, ElseMatrixType> XprType;
-  enum {
-    CoeffReadCost = evaluator<ConditionMatrixType>::CoeffReadCost +
-                    plain_enum_max(evaluator<ThenMatrixType>::CoeffReadCost, evaluator<ElseMatrixType>::CoeffReadCost),
-
-    Flags = (unsigned int)evaluator<ThenMatrixType>::Flags & evaluator<ElseMatrixType>::Flags & HereditaryBits,
-
-    Alignment = plain_enum_min(evaluator<ThenMatrixType>::Alignment, evaluator<ElseMatrixType>::Alignment)
-  };
-
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit evaluator(const XprType& select)
-      : m_conditionImpl(select.conditionMatrix()), m_thenImpl(select.thenMatrix()), m_elseImpl(select.elseMatrix()) {
-    EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
-  }
-
-  typedef typename XprType::CoeffReturnType CoeffReturnType;
-
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index row, Index col) const {
-    if (m_conditionImpl.coeff(row, col))
-      return m_thenImpl.coeff(row, col);
-    else
-      return m_elseImpl.coeff(row, col);
-  }
-
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const {
-    if (m_conditionImpl.coeff(index))
-      return m_thenImpl.coeff(index);
-    else
-      return m_elseImpl.coeff(index);
-  }
-
- protected:
-  evaluator<ConditionMatrixType> m_conditionImpl;
-  evaluator<ThenMatrixType> m_thenImpl;
-  evaluator<ElseMatrixType> m_elseImpl;
-};
-
 // -------------------- Replicate --------------------
 
 template <typename ArgType, int RowFactor, int ColFactor>
diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h
index 5e4e5c2..bf41c3b 100644
--- a/Eigen/src/Core/NumTraits.h
+++ b/Eigen/src/Core/NumTraits.h
@@ -95,9 +95,22 @@
 }  // end namespace internal
 
 namespace numext {
-/** \internal bit-wise cast without changing the underlying bit representation. */
 
-// TODO: Replace by std::bit_cast (available in C++20)
+/** \internal bit-wise cast without changing the underlying bit representation. */
+#if defined(__cpp_lib_bit_cast) && __cpp_lib_bit_cast >= 201806L
+template <typename Tgt, typename Src>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC constexpr Tgt bit_cast(const Src& src) {
+  return std::bit_cast<Tgt>(src);
+}
+#elif EIGEN_HAS_BUILTIN(__builtin_bit_cast)
+template <typename Tgt, typename Src>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC constexpr Tgt bit_cast(const Src& src) {
+  EIGEN_STATIC_ASSERT(std::is_trivially_copyable<Src>::value, THIS_TYPE_IS_NOT_SUPPORTED)
+  EIGEN_STATIC_ASSERT(std::is_trivially_copyable<Tgt>::value, THIS_TYPE_IS_NOT_SUPPORTED)
+  EIGEN_STATIC_ASSERT(sizeof(Src) == sizeof(Tgt), THIS_TYPE_IS_NOT_SUPPORTED)
+  return __builtin_bit_cast(Tgt, src);
+}
+#else
 template <typename Tgt, typename Src>
 EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Tgt bit_cast(const Src& src) {
   // The behaviour of memcpy is not specified for non-trivially copyable types
@@ -113,6 +126,7 @@
   memcpy(static_cast<void*>(&tgt), static_cast<const void*>(&staged), sizeof(Tgt));
   return tgt;
 }
+#endif
 }  // namespace numext
 
 // clang-format off
diff --git a/Eigen/src/Core/Select.h b/Eigen/src/Core/Select.h
index 0fa5f1e..61a67c2 100644
--- a/Eigen/src/Core/Select.h
+++ b/Eigen/src/Core/Select.h
@@ -15,7 +15,7 @@
 
 namespace Eigen {
 
-/** \class Select
+/** \typedef Select
  * \ingroup Core_Module
  *
  * \brief Expression of a coefficient wise version of the C++ ternary operator ?:
@@ -24,73 +24,16 @@
  * \tparam ThenMatrixType the type of the \em then expression
  * \tparam ElseMatrixType the type of the \em else expression
  *
- * This class represents an expression of a coefficient wise version of the C++ ternary operator ?:.
+ * This type represents an expression of a coefficient wise version of the C++ ternary operator ?:.
  * It is the return type of DenseBase::select() and most of the time this is the only way it is used.
  *
  * \sa DenseBase::select(const DenseBase<ThenDerived>&, const DenseBase<ElseDerived>&) const
  */
-
-namespace internal {
 template <typename ConditionMatrixType, typename ThenMatrixType, typename ElseMatrixType>
-struct traits<Select<ConditionMatrixType, ThenMatrixType, ElseMatrixType> > : traits<ThenMatrixType> {
-  typedef typename traits<ThenMatrixType>::Scalar Scalar;
-  typedef Dense StorageKind;
-  typedef typename traits<ThenMatrixType>::XprKind XprKind;
-  typedef typename ConditionMatrixType::Nested ConditionMatrixNested;
-  typedef typename ThenMatrixType::Nested ThenMatrixNested;
-  typedef typename ElseMatrixType::Nested ElseMatrixNested;
-  enum {
-    RowsAtCompileTime = ConditionMatrixType::RowsAtCompileTime,
-    ColsAtCompileTime = ConditionMatrixType::ColsAtCompileTime,
-    MaxRowsAtCompileTime = ConditionMatrixType::MaxRowsAtCompileTime,
-    MaxColsAtCompileTime = ConditionMatrixType::MaxColsAtCompileTime,
-    Flags = (unsigned int)ThenMatrixType::Flags & ElseMatrixType::Flags & RowMajorBit
-  };
-};
-}  // namespace internal
-
-template <typename ConditionMatrixType, typename ThenMatrixType, typename ElseMatrixType>
-class Select : public internal::dense_xpr_base<Select<ConditionMatrixType, ThenMatrixType, ElseMatrixType> >::type,
-               internal::no_assignment_operator {
- public:
-  typedef typename internal::dense_xpr_base<Select>::type Base;
-  EIGEN_DENSE_PUBLIC_INTERFACE(Select)
-
-  inline EIGEN_DEVICE_FUNC Select(const ConditionMatrixType& a_conditionMatrix, const ThenMatrixType& a_thenMatrix,
-                                  const ElseMatrixType& a_elseMatrix)
-      : m_condition(a_conditionMatrix), m_then(a_thenMatrix), m_else(a_elseMatrix) {
-    eigen_assert(m_condition.rows() == m_then.rows() && m_condition.rows() == m_else.rows());
-    eigen_assert(m_condition.cols() == m_then.cols() && m_condition.cols() == m_else.cols());
-  }
-
-  EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept { return m_condition.rows(); }
-  EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return m_condition.cols(); }
-
-  inline EIGEN_DEVICE_FUNC const Scalar coeff(Index i, Index j) const {
-    if (m_condition.coeff(i, j))
-      return m_then.coeff(i, j);
-    else
-      return m_else.coeff(i, j);
-  }
-
-  inline EIGEN_DEVICE_FUNC const Scalar coeff(Index i) const {
-    if (m_condition.coeff(i))
-      return m_then.coeff(i);
-    else
-      return m_else.coeff(i);
-  }
-
-  inline EIGEN_DEVICE_FUNC const ConditionMatrixType& conditionMatrix() const { return m_condition; }
-
-  inline EIGEN_DEVICE_FUNC const ThenMatrixType& thenMatrix() const { return m_then; }
-
-  inline EIGEN_DEVICE_FUNC const ElseMatrixType& elseMatrix() const { return m_else; }
-
- protected:
-  typename ConditionMatrixType::Nested m_condition;
-  typename ThenMatrixType::Nested m_then;
-  typename ElseMatrixType::Nested m_else;
-};
+using Select = CwiseTernaryOp<internal::scalar_boolean_select_op<typename DenseBase<ThenMatrixType>::Scalar,
+                                                                 typename DenseBase<ElseMatrixType>::Scalar,
+                                                                 typename DenseBase<ConditionMatrixType>::Scalar>,
+                              ThenMatrixType, ElseMatrixType, ConditionMatrixType>;
 
 /** \returns a matrix where each coefficient (i,j) is equal to \a thenMatrix(i,j)
  * if \c *this(i,j) != Scalar(0), and \a elseMatrix(i,j) otherwise.
@@ -98,7 +41,7 @@
  * Example: \include MatrixBase_select.cpp
  * Output: \verbinclude MatrixBase_select.out
  *
- * \sa DenseBase::bitwiseSelect(const DenseBase<ThenDerived>&, const DenseBase<ElseDerived>&)
+ * \sa typedef Select
  */
 template <typename Derived>
 template <typename ThenDerived, typename ElseDerived>
@@ -107,15 +50,12 @@
                                        typename DenseBase<Derived>::Scalar>,
     ThenDerived, ElseDerived, Derived>
 DenseBase<Derived>::select(const DenseBase<ThenDerived>& thenMatrix, const DenseBase<ElseDerived>& elseMatrix) const {
-  using Op = internal::scalar_boolean_select_op<typename DenseBase<ThenDerived>::Scalar,
-                                                typename DenseBase<ElseDerived>::Scalar, Scalar>;
-  return CwiseTernaryOp<Op, ThenDerived, ElseDerived, Derived>(thenMatrix.derived(), elseMatrix.derived(), derived(),
-                                                               Op());
+  return Select<Derived, ThenDerived, ElseDerived>(thenMatrix.derived(), elseMatrix.derived(), derived());
 }
 /** Version of DenseBase::select(const DenseBase&, const DenseBase&) with
  * the \em else expression being a scalar value.
  *
- * \sa DenseBase::booleanSelect(const DenseBase<ThenDerived>&, const DenseBase<ElseDerived>&) const, class Select
+ * \sa typedef Select
  */
 template <typename Derived>
 template <typename ThenDerived>
@@ -126,15 +66,13 @@
 DenseBase<Derived>::select(const DenseBase<ThenDerived>& thenMatrix,
                            const typename DenseBase<ThenDerived>::Scalar& elseScalar) const {
   using ElseConstantType = typename DenseBase<ThenDerived>::ConstantReturnType;
-  using Op = internal::scalar_boolean_select_op<typename DenseBase<ThenDerived>::Scalar,
-                                                typename DenseBase<ThenDerived>::Scalar, Scalar>;
-  return CwiseTernaryOp<Op, ThenDerived, ElseConstantType, Derived>(
-      thenMatrix.derived(), ElseConstantType(rows(), cols(), elseScalar), derived(), Op());
+  return Select<Derived, ThenDerived, ElseConstantType>(thenMatrix.derived(),
+                                                        ElseConstantType(rows(), cols(), elseScalar), derived());
 }
 /** Version of DenseBase::select(const DenseBase&, const DenseBase&) with
  * the \em then expression being a scalar value.
  *
- * \sa DenseBase::booleanSelect(const DenseBase<ThenDerived>&, const DenseBase<ElseDerived>&) const, class Select
+ * \sa typedef Select
  */
 template <typename Derived>
 template <typename ElseDerived>
@@ -145,10 +83,8 @@
 DenseBase<Derived>::select(const typename DenseBase<ElseDerived>::Scalar& thenScalar,
                            const DenseBase<ElseDerived>& elseMatrix) const {
   using ThenConstantType = typename DenseBase<ElseDerived>::ConstantReturnType;
-  using Op = internal::scalar_boolean_select_op<typename DenseBase<ElseDerived>::Scalar,
-                                                typename DenseBase<ElseDerived>::Scalar, Scalar>;
-  return CwiseTernaryOp<Op, ThenConstantType, ElseDerived, Derived>(ThenConstantType(rows(), cols(), thenScalar),
-                                                                    elseMatrix.derived(), derived(), Op());
+  return Select<Derived, ThenConstantType, ElseDerived>(ThenConstantType(rows(), cols(), thenScalar),
+                                                        elseMatrix.derived(), derived());
 }
 
 }  // end namespace Eigen
diff --git a/Eigen/src/Core/VectorwiseOp.h b/Eigen/src/Core/VectorwiseOp.h
index 9ccbf7d..9e34d8c 100644
--- a/Eigen/src/Core/VectorwiseOp.h
+++ b/Eigen/src/Core/VectorwiseOp.h
@@ -146,6 +146,22 @@
   const BinaryOp& binaryFunc() const { return m_functor; }
   const BinaryOp m_functor;
 };
+
+template <typename Scalar>
+struct scalar_replace_zero_with_one_op {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar operator()(const Scalar& x) const {
+    return numext::is_exactly_zero(x) ? Scalar(1) : x;
+  }
+  template <typename Packet>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+    return pselect(pcmp_eq(x, pzero(x)), pset1<Packet>(Scalar(1)), x);
+  }
+};
+template <typename Scalar>
+struct functor_traits<scalar_replace_zero_with_one_op<Scalar>> {
+  enum { Cost = 1, PacketAccess = packet_traits<Scalar>::HasCmp };
+};
+
 }  // namespace internal
 
 /** \class VectorwiseOp
@@ -624,18 +640,28 @@
     return m_matrix / extendedTo(other.derived());
   }
 
+  using Normalized_NonzeroNormType =
+      CwiseUnaryOp<internal::scalar_replace_zero_with_one_op<Scalar>, const NormReturnType>;
+  using NormalizedReturnType = CwiseBinaryOp<internal::scalar_quotient_op<Scalar>, const ExpressionTypeNestedCleaned,
+                                             const typename OppositeExtendedType<Normalized_NonzeroNormType>::Type>;
+
   /** \returns an expression where each column (or row) of the referenced matrix are normalized.
    * The referenced matrix is \b not modified.
+   *
+   * \warning If the input columns (or rows) are too small (i.e., their norm equals to 0), they remain unchanged in the
+   *          resulting expression.
+   *
    * \sa MatrixBase::normalized(), normalize()
    */
-  EIGEN_DEVICE_FUNC CwiseBinaryOp<internal::scalar_quotient_op<Scalar>, const ExpressionTypeNestedCleaned,
-                                  const typename OppositeExtendedType<NormReturnType>::Type>
-  normalized() const {
-    return m_matrix.cwiseQuotient(extendedToOpposite(this->norm()));
+  EIGEN_DEVICE_FUNC NormalizedReturnType normalized() const {
+    return m_matrix.cwiseQuotient(extendedToOpposite(Normalized_NonzeroNormType(this->norm())));
   }
 
   /** Normalize in-place each row or columns of the referenced matrix.
-   * \sa MatrixBase::normalize(), normalized()
+   *
+   * \warning If the input columns (or rows) are too small (i.e., their norm equals to 0), they are left unchanged.
+   *
+   * \sa MatrixBase::normalized(), normalize()
    */
   EIGEN_DEVICE_FUNC void normalize() { m_matrix = this->normalized(); }
 
diff --git a/Eigen/src/Core/arch/HVX/PacketMath.h b/Eigen/src/Core/arch/HVX/PacketMath.h
index b9080d9..9b6ceb3 100644
--- a/Eigen/src/Core/arch/HVX/PacketMath.h
+++ b/Eigen/src/Core/arch/HVX/PacketMath.h
@@ -400,8 +400,25 @@
 }
 
 template <HVXPacketSize T>
+EIGEN_STRONG_INLINE HVXPacket<T> ptrue_hvx(const HVXPacket<T>& a) {
+  return HVXPacket<T>::Create(Q6_V_vsplat_R(0x3f800000));
+}
+template <>
+EIGEN_STRONG_INLINE Packet32f ptrue(const Packet32f& a) {
+  return ptrue_hvx(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet16f ptrue(const Packet16f& a) {
+  return ptrue_hvx(a);
+}
+template <>
+EIGEN_STRONG_INLINE Packet8f ptrue(const Packet8f& a) {
+  return ptrue_hvx(a);
+}
+
+template <HVXPacketSize T>
 EIGEN_STRONG_INLINE HVXPacket<T> pcmp_le_hvx(const HVXPacket<T>& a, const HVXPacket<T>& b) {
-  HVX_Vector v_true = Q6_V_vsplat_R(0x3f800000);
+  HVX_Vector v_true = ptrue(a).Get();
   HVX_VectorPred pred = Q6_Q_vcmp_gt_VsfVsf(a.Get(), b.Get());
   return HVXPacket<T>::Create(Q6_V_vmux_QVV(pred, Q6_V_vzero(), v_true));
 }
@@ -420,7 +437,7 @@
 
 template <HVXPacketSize T>
 EIGEN_STRONG_INLINE HVXPacket<T> pcmp_eq_hvx(const HVXPacket<T>& a, const HVXPacket<T>& b) {
-  HVX_Vector v_true = Q6_V_vsplat_R(0x3f800000);
+  HVX_Vector v_true = ptrue(a).Get();
   HVX_VectorPred pred = Q6_Q_vcmp_eq_VwVw(a.Get(), b.Get());
   return HVXPacket<T>::Create(Q6_V_vmux_QVV(pred, v_true, Q6_V_vzero()));
 }
@@ -439,7 +456,7 @@
 
 template <HVXPacketSize T>
 EIGEN_STRONG_INLINE HVXPacket<T> pcmp_lt_hvx(const HVXPacket<T>& a, const HVXPacket<T>& b) {
-  HVX_Vector v_true = Q6_V_vsplat_R(0x3f800000);
+  HVX_Vector v_true = ptrue(a).Get();
   HVX_VectorPred pred = Q6_Q_vcmp_gt_VsfVsf(b.Get(), a.Get());
   return HVXPacket<T>::Create(Q6_V_vmux_QVV(pred, v_true, Q6_V_vzero()));
 }
@@ -458,7 +475,7 @@
 
 template <HVXPacketSize T>
 EIGEN_STRONG_INLINE HVXPacket<T> pcmp_lt_or_nan_hvx(const HVXPacket<T>& a, const HVXPacket<T>& b) {
-  HVX_Vector v_true = Q6_V_vsplat_R(0x3f800000);
+  HVX_Vector v_true = ptrue(a).Get();
   HVX_VectorPred pred = Q6_Q_vcmp_gt_VsfVsf(b.Get(), a.Get());
   return HVXPacket<T>::Create(Q6_V_vmux_QVV(pred, v_true, Q6_V_vzero()));
 }
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 3c0bc46..4eb134f 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -397,8 +397,6 @@
                                                                    : EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION),
           int MaxRows_ = Rows_, int MaxCols_ = Cols_>
 class Array;
-template <typename ConditionMatrixType, typename ThenMatrixType, typename ElseMatrixType>
-class Select;
 template <typename MatrixType, typename BinaryOp, int Direction>
 class PartialReduxExpr;
 template <typename ExpressionType, int Direction>
diff --git a/test/array_cwise.cpp b/test/array_cwise.cpp
index 6ff8d67..d37e47a 100644
--- a/test/array_cwise.cpp
+++ b/test/array_cwise.cpp
@@ -734,7 +734,7 @@
   VERIFY_IS_CWISE_EQUAL(m1.abs().cwiseLessOrEqual(NumTraits<Scalar>::highest()), bool_true);
   VERIFY_IS_CWISE_EQUAL(m1.abs().cwiseGreaterOrEqual(Scalar(0)), bool_true);
 
-  // test Select
+  // test select
   VERIFY_IS_APPROX((m1 < m2).select(m1, m2), m1.cwiseMin(m2));
   VERIFY_IS_APPROX((m1 > m2).select(m1, m2), m1.cwiseMax(m2));
   Scalar mid = (m1.cwiseAbs().minCoeff() + m1.cwiseAbs().maxCoeff()) / Scalar(2);
diff --git a/test/array_for_matrix.cpp b/test/array_for_matrix.cpp
index e94398a..e8d2c7d 100644
--- a/test/array_for_matrix.cpp
+++ b/test/array_for_matrix.cpp
@@ -120,7 +120,7 @@
   VERIFY((m1.array() == m1(r, c)).any());
   VERIFY(m1.cwiseEqual(m1(r, c)).any());
 
-  // test Select
+  // test select
   VERIFY_IS_APPROX((m1.array() < m2.array()).select(m1, m2), m1.cwiseMin(m2));
   VERIFY_IS_APPROX((m1.array() > m2.array()).select(m1, m2), m1.cwiseMax(m2));
   Scalar mid = m1.cwiseAbs().minCoeff() / Scalar(2) + m1.cwiseAbs().maxCoeff() / Scalar(2);
diff --git a/test/evaluators.cpp b/test/evaluators.cpp
index d96a0a1..5a4ab37 100644
--- a/test/evaluators.cpp
+++ b/test/evaluators.cpp
@@ -340,7 +340,7 @@
     matXcd_ref.imag() = mat2;
     VERIFY_IS_APPROX(matXcd, matXcd_ref);
 
-    // test Select
+    // test select
     VERIFY_IS_APPROX_EVALUATOR(aX, (aXsrc > 0).select(aXsrc, -aXsrc));
 
     // test Replicate
diff --git a/test/main.h b/test/main.h
index 2288778..db4d484 100644
--- a/test/main.h
+++ b/test/main.h
@@ -14,6 +14,7 @@
 #include <iostream>
 #include <iomanip>
 #include <fstream>
+#include <limits>
 #include <string>
 #include <sstream>
 #include <vector>
diff --git a/test/vectorwiseop.cpp b/test/vectorwiseop.cpp
index 6d0e5cb..d037bb4 100644
--- a/test/vectorwiseop.cpp
+++ b/test/vectorwiseop.cpp
@@ -114,6 +114,8 @@
   RealColVectorType rcres;
   RealRowVectorType rrres;
 
+  Scalar small_scalar = (std::numeric_limits<RealScalar>::min)();
+
   // test broadcast assignment
   m2 = m1;
   m2.colwise() = colvec;
@@ -171,18 +173,26 @@
   VERIFY_IS_APPROX(m1.cwiseAbs().colwise().sum().x(), m1.col(0).cwiseAbs().sum());
 
   // test normalized
-  m2 = m1.colwise().normalized();
-  VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized());
-  m2 = m1.rowwise().normalized();
-  VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized());
+  m2 = m1;
+  m2.col(c).fill(small_scalar);
+  m3 = m2.colwise().normalized();
+  for (Index k = 0; k < cols; ++k) VERIFY_IS_APPROX(m3.col(k), m2.col(k).normalized());
+  m2 = m1;
+  m2.row(r).setZero();
+  m3 = m2.rowwise().normalized();
+  for (Index k = 0; k < rows; ++k) VERIFY_IS_APPROX(m3.row(k), m2.row(k).normalized());
 
   // test normalize
   m2 = m1;
-  m2.colwise().normalize();
-  VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized());
+  m2.col(c).setZero();
+  m3 = m2;
+  m3.colwise().normalize();
+  for (Index k = 0; k < cols; ++k) VERIFY_IS_APPROX(m3.col(k), m2.col(k).normalized());
   m2 = m1;
-  m2.rowwise().normalize();
-  VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized());
+  m2.row(r).fill(small_scalar);
+  m3 = m2;
+  m3.rowwise().normalize();
+  for (Index k = 0; k < rows; ++k) VERIFY_IS_APPROX(m3.row(k), m2.row(k).normalized());
 
   // test with partial reduction of products
   Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> m1m1 = m1 * m1.transpose();