Cleaning pass on rcond estimator.

1 files changed, 60 insertions, 101 deletions

diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h
index f53c2a837..68c5e918e 100644
--- a/Eigen/src/Core/ConditionEstimator.h
+++ b/Eigen/src/Core/ConditionEstimator.h
@@ -13,139 +13,97 @@
 namespace Eigen {
 
 namespace internal {
-template <typename MatrixType>
-inline typename MatrixType::RealScalar MatrixL1Norm(const MatrixType& matrix) {
-  return matrix.cwiseAbs().colwise().sum().maxCoeff();
-}
-
-template <typename Vector>
-inline typename Vector::RealScalar VectorL1Norm(const Vector& v) {
-  return v.template lpNorm<1>();
-}
 
 template <typename Vector, typename RealVector, bool IsComplex>
-struct SignOrUnity {
+struct rcond_compute_sign {
   static inline Vector run(const Vector& v) {
     const RealVector v_abs = v.cwiseAbs();
     return (v_abs.array() == static_cast<typename Vector::RealScalar>(0))
-        .select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs));
+            .select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs));
   }
 };
 
 // Partial specialization to avoid elementwise division for real vectors.
 template <typename Vector>
-struct SignOrUnity<Vector, Vector, false> {
+struct rcond_compute_sign<Vector, Vector, false> {
   static inline Vector run(const Vector& v) {
     return (v.array() < static_cast<typename Vector::RealScalar>(0))
-        .select(-Vector::Ones(v.size()), Vector::Ones(v.size()));
+           .select(-Vector::Ones(v.size()), Vector::Ones(v.size()));
   }
 };
 
-}  // namespace internal
-
-/** \class ConditionEstimator
-    * \ingroup Core_Module
-    *
-    * \brief Condition number estimator.
-    *
-    * Computing a decomposition of a dense matrix takes O(n^3) operations, while
-    * this method estimates the condition number quickly and reliably in O(n^2)
-    * operations.
-    *
-    * \returns an estimate of the reciprocal condition number
-    * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given the matrix and
-    * its decomposition. Supports the following decompositions: FullPivLU,
-    * PartialPivLU, LDLT, and LLT.
-    *
-    * \sa FullPivLU, PartialPivLU, LDLT, LLT.
-    */
+/** \brief Reciprocal condition number estimator.
+  *
+  * Computing a decomposition of a dense matrix takes O(n^3) operations, while
+  * this method estimates the condition number quickly and reliably in O(n^2)
+  * operations.
+  *
+  * \returns an estimate of the reciprocal condition number
+  * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given ||matrix||_1 and
+  * its decomposition. Supports the following decompositions: FullPivLU,
+  * PartialPivLU, LDLT, and LLT.
+  *
+  * \sa FullPivLU, PartialPivLU, LDLT, LLT.
+  */
 template <typename Decomposition>
-typename Decomposition::RealScalar ReciprocalConditionNumberEstimate(
-    const typename Decomposition::MatrixType& matrix,
-    const Decomposition& dec) {
-  eigen_assert(matrix.rows() == dec.rows());
-  eigen_assert(matrix.cols() == dec.cols());
-  if (dec.rows() == 0) return typename Decomposition::RealScalar(1);
-  return ReciprocalConditionNumberEstimate(MatrixL1Norm(matrix), dec);
-}
-
-/** \class ConditionEstimator
- * \ingroup Core_Module
- *
- * \brief Condition number estimator.
- *
- * Computing a decomposition of a dense matrix takes O(n^3) operations, while
- * this method estimates the condition number quickly and reliably in O(n^2)
- * operations.
- *
- * \returns an estimate of the reciprocal condition number
- * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given ||matrix||_1 and
- * its decomposition. Supports the following decompositions: FullPivLU,
- * PartialPivLU, LDLT, and LLT.
- *
- * \sa FullPivLU, PartialPivLU, LDLT, LLT.
- */
-template <typename Decomposition>
-typename Decomposition::RealScalar ReciprocalConditionNumberEstimate(
-    typename Decomposition::RealScalar matrix_norm, const Decomposition& dec) {
+typename Decomposition::RealScalar
+rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Decomposition& dec)
+{
   typedef typename Decomposition::RealScalar RealScalar;
   eigen_assert(dec.rows() == dec.cols());
-  if (dec.rows() == 0) return RealScalar(1);
+  if (dec.rows() == 0)              return RealScalar(1);
   if (matrix_norm == RealScalar(0)) return RealScalar(0);
-  if (dec.rows() == 1) return RealScalar(1);
-  const typename Decomposition::RealScalar inverse_matrix_norm =
-      InverseMatrixL1NormEstimate(dec);
-  return (inverse_matrix_norm == RealScalar(0)
-              ? RealScalar(0)
-              : (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
+  if (dec.rows() == 1)              return RealScalar(1);
+  const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec);
+  return (inverse_matrix_norm == RealScalar(0) ? RealScalar(0)
+                                               : (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
 }
 
 /**
- * \returns an estimate of ||inv(matrix)||_1 given a decomposition of
- * matrix that implements .solve() and .adjoint().solve() methods.
- *
- * The method implements Algorithms 4.1 and 5.1 from
- *   http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf
- * which also forms the basis for the condition number estimators in
- * LAPACK. Since at most 10 calls to the solve method of dec are
- * performed, the total cost is O(dims^2), as opposed to O(dims^3)
- * needed to compute the inverse matrix explicitly.
- *
- * The most common usage is in estimating the condition number
- * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be
- * computed directly in O(n^2) operations.
- *
- * Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, and
- * LLT.
- *
- * \sa FullPivLU, PartialPivLU, LDLT, LLT.
- */
+  * \returns an estimate of ||inv(matrix)||_1 given a decomposition of
+  * \a matrix that implements .solve() and .adjoint().solve() methods.
+  *
+  * This function implements Algorithms 4.1 and 5.1 from
+  *   http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf
+  * which also forms the basis for the condition number estimators in
+  * LAPACK. Since at most 10 calls to the solve method of dec are
+  * performed, the total cost is O(dims^2), as opposed to O(dims^3)
+  * needed to compute the inverse matrix explicitly.
+  *
+  * The most common usage is in estimating the condition number
+  * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be
+  * computed directly in O(n^2) operations.
+  *
+  * Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, and
+  * LLT.
+  *
+  * \sa FullPivLU, PartialPivLU, LDLT, LLT.
+  */
 template <typename Decomposition>
-typename Decomposition::RealScalar InverseMatrixL1NormEstimate(
-    const Decomposition& dec) {
+typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomposition& dec)
+{
   typedef typename Decomposition::MatrixType MatrixType;
   typedef typename Decomposition::Scalar Scalar;
   typedef typename Decomposition::RealScalar RealScalar;
   typedef typename internal::plain_col_type<MatrixType>::type Vector;
-  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type
-      RealVector;
+  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVector;
   const bool is_complex = (NumTraits<Scalar>::IsComplex != 0);
 
   eigen_assert(dec.rows() == dec.cols());
   const Index n = dec.rows();
-  if (n == 0) {
+  if (n == 0)
     return 0;
-  }
+
   Vector v = dec.solve(Vector::Ones(n) / Scalar(n));
 
   // lower_bound is a lower bound on
   //   ||inv(matrix)||_1  = sup_v ||inv(matrix) v||_1 / ||v||_1
   // and is the objective maximized by the ("super-") gradient ascent
   // algorithm below.
-  RealScalar lower_bound = internal::VectorL1Norm(v);
-  if (n == 1) {
+  RealScalar lower_bound = v.template lpNorm<1>();
+  if (n == 1)
     return lower_bound;
-  }
+
   // Gradient ascent algorithm follows: We know that the optimum is achieved at
   // one of the simplices v = e_i, so in each iteration we follow a
   // super-gradient to move towards the optimal one.
@@ -154,8 +112,9 @@ typename Decomposition::RealScalar InverseMatrixL1NormEstimate(
   Vector old_sign_vector;
   Index v_max_abs_index = -1;
   Index old_v_max_abs_index = v_max_abs_index;
-  for (int k = 0; k < 4; ++k) {
-    sign_vector = internal::SignOrUnity<Vector, RealVector, is_complex>::run(v);
+  for (int k = 0; k < 4; ++k)
+  {
+    sign_vector = internal::rcond_compute_sign<Vector, RealVector, is_complex>::run(v);
     if (k > 0 && !is_complex && sign_vector == old_sign_vector) {
       // Break if the solution stagnated.
       break;
@@ -169,7 +128,7 @@ typename Decomposition::RealScalar InverseMatrixL1NormEstimate(
     }
     // Move to the new simplex e_j, where j = v_max_abs_index.
     v = dec.solve(Vector::Unit(n, v_max_abs_index));  // v = inv(matrix) * e_j.
-    lower_bound = internal::VectorL1Norm(v);
+    lower_bound = v.template lpNorm<1>();
     if (lower_bound <= old_lower_bound) {
       // Break if the gradient step did not increase the lower_bound.
       break;
@@ -192,16 +151,16 @@ typename Decomposition::RealScalar InverseMatrixL1NormEstimate(
   // Hager's algorithm to vastly underestimate ||matrix||_1.
   Scalar alternating_sign(RealScalar(1));
   for (Index i = 0; i < n; ++i) {
-    v[i] = alternating_sign *
-           (RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1))));
+    v[i] = alternating_sign * (RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1))));
     alternating_sign = -alternating_sign;
   }
   v = dec.solve(v);
-  const RealScalar alternate_lower_bound =
-      (2 * internal::VectorL1Norm(v)) / (3 * RealScalar(n));
+  const RealScalar alternate_lower_bound = (2 * v.template lpNorm<1>()) / (3 * RealScalar(n));
   return numext::maxi(lower_bound, alternate_lower_bound);
 }
 
+}  // namespace internal
+
 }  // namespace Eigen
 
 #endif