SelfAdjointView: add l1Norm() and use it in LLT / LDLT

libeigen/eigen!2456

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

Eigen/src/Cholesky/LDLT.h

diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 45912a2..63aa5bd 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -489,19 +489,8 @@
 
   m_matrix = a.derived();
 
-  // Compute matrix L1 norm = max abs column sum.
-  m_l1_norm = RealScalar(0);
-  // TODO: move this code to SelfAdjointView
-  for (Index col = 0; col < size; ++col) {
-    RealScalar abs_col_sum;
-    if (UpLo_ == Lower)
-      abs_col_sum =
-          m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
-    else
-      abs_col_sum =
-          m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
-    if (abs_col_sum > m_l1_norm) m_l1_norm = abs_col_sum;
-  }
+  // Compute matrix L1 norm = max abs column sum over the implicit self-adjoint matrix.
+  m_l1_norm = m_matrix.template selfadjointView<UpLo_>().l1Norm();
 
   m_transpositions.resize(size);
   m_isInitialized = false;

Eigen/src/Cholesky/LLT.h

diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 7066cd0..9bffeae 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -404,19 +404,8 @@
   m_matrix.resize(size, size);
   if (!internal::is_same_dense(m_matrix, a.derived())) m_matrix = a.derived();
 
-  // Compute matrix L1 norm = max abs column sum.
-  m_l1_norm = RealScalar(0);
-  // TODO: move this code to SelfAdjointView
-  for (Index col = 0; col < size; ++col) {
-    RealScalar abs_col_sum;
-    if (UpLo_ == Lower)
-      abs_col_sum =
-          m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
-    else
-      abs_col_sum =
-          m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
-    if (abs_col_sum > m_l1_norm) m_l1_norm = abs_col_sum;
-  }
+  // Compute matrix L1 norm = max abs column sum over the implicit self-adjoint matrix.
+  m_l1_norm = m_matrix.template selfadjointView<UpLo_>().l1Norm();
 
   m_isInitialized = true;
   bool ok = Traits::inplace_decomposition(m_matrix);

Eigen/src/Core/SelfAdjointView.h

diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h
index ff450eb..59731b8 100644
--- a/Eigen/src/Core/SelfAdjointView.h
+++ b/Eigen/src/Core/SelfAdjointView.h
@@ -223,6 +223,28 @@
     return typename MatrixType::ConstDiagonalReturnType(m_matrix);
   }
 
+  /** \returns the matrix 1-norm (maximum absolute column sum) of the implicit
+   * full self-adjoint matrix, reading only the stored triangle. For Hermitian
+   * (complex) scalars the unstored entries are conjugates of stored ones, and
+   * since |conj(x)| = |x| the result matches the L1 norm of the full matrix.
+   */
+  EIGEN_DEVICE_FUNC typename NumTraits<Scalar>::Real l1Norm() const {
+    typedef typename NumTraits<Scalar>::Real RealScalar_;
+    RealScalar_ norm = RealScalar_(0);
+    const Index n = m_matrix.rows();
+    for (Index col = 0; col < n; ++col) {
+      RealScalar_ abs_col_sum;
+      if (UpLo == Lower)
+        abs_col_sum =
+            m_matrix.col(col).tail(n - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
+      else
+        abs_col_sum =
+            m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(n - col).template lpNorm<1>();
+      if (abs_col_sum > norm) norm = abs_col_sum;
+    }
+    return norm;
+  }
+
   /////////// Cholesky module ///////////
 
   LLT<PlainObject, UpLo> llt() const;

test/selfadjoint.cpp

diff --git a/test/selfadjoint.cpp b/test/selfadjoint.cpp
index 65d4e79..87d3152 100644
--- a/test/selfadjoint.cpp
+++ b/test/selfadjoint.cpp
@@ -54,6 +54,24 @@
   VERIFY_IS_APPROX(m4, MatrixType((s * m1).template selfadjointView<Lower>()));
   m4 = m1.template selfadjointView<Lower>() * s;
   VERIFY_IS_APPROX(m4, MatrixType((m1 * s).template selfadjointView<Lower>()));
+
+  // l1Norm: reads only the stored triangle but must agree with the L1 norm of
+  // the materialized full matrix. Upper and Lower views of the same (stored-
+  // full) self-adjoint matrix must return the same value; complex scalars
+  // behave identically since |conj(x)| = |x|.
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  m3 = m1.template selfadjointView<Upper>();  // m3 is now fully self-adjoint
+  RealScalar ref_l1 = m3.cwiseAbs().colwise().sum().maxCoeff();
+  VERIFY_IS_APPROX(m3.template selfadjointView<Upper>().l1Norm(), ref_l1);
+  VERIFY_IS_APPROX(m3.template selfadjointView<Lower>().l1Norm(), ref_l1);
+  // Either triangle alone still gives the correct L1 norm even if the other
+  // half is zero (the view conjure it back via symmetry).
+  MatrixType upperOnly = MatrixType::Zero(rows, cols);
+  upperOnly.template triangularView<Upper>() = m3;
+  VERIFY_IS_APPROX(upperOnly.template selfadjointView<Upper>().l1Norm(), ref_l1);
+  MatrixType lowerOnly = MatrixType::Zero(rows, cols);
+  lowerOnly.template triangularView<Lower>() = m3;
+  VERIFY_IS_APPROX(lowerOnly.template selfadjointView<Lower>().l1Norm(), ref_l1);
 }
 
 void bug_159() {