Core: Stream and vectorize the self-adjoint 1-norm

libeigen/eigen!3060

Eigen/src/Core/SelfAdjointView.h

diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h
index b9e5e20..439e19c 100644
--- a/Eigen/src/Core/SelfAdjointView.h
+++ b/Eigen/src/Core/SelfAdjointView.h
@@ -33,8 +33,224 @@
  */
 
 namespace internal {
+
+// Column step of the self-adjoint 1-norm, on two columns sharing a range of rows: sums[i] +=
+// |m(i, j0)| + |m(i, j1)|, and each column's own sum of those rows goes to sums[j0] and sums[j1].
+// Walking two columns at once halves the traffic on sums, and one packet pass does everything, so
+// short columns pay no per-expression setup.
+//
+// Real scalars use pabs. Complex ones have no packet abs, so |z| = sqrt(re^2 + im^2) is computed
+// on the real lanes of the complex packet. That leaves |z|^2 in both lanes of each slot, so one
+// square root serves both columns: even lanes from j0 and odd lanes from j1 give |u0| |v0| |u1|
+// |v1| ..., whose sum with its flip is the update of sums, and whose reduction as a complex packet
+// is (sum |u|, sum |v|). The accumulator has the matrix's scalar type and only its real parts are
+// read, so what lands in the imaginary lanes is harmless. Squaring overflows above sqrt(max) and
+// loses precision below sqrt(min), so the pass also records the largest component it has seen and
+// the caller recomputes the norm through numext::abs when that is out of range. The record is
+// exact and the range test finite, so neither depends on infinities surviving fast-math.
+template <typename Scalar_>
+struct selfadjoint_l1norm_real_lanes {
+  using Scalar = Scalar_;
+  using Real = Scalar_;
+  using Packet = typename packet_traits<Scalar>::type;
+  using RPacket = Packet;
+  static constexpr Index PacketSize = unpacket_traits<Packet>::size;
+  // Up to this size the per-column form (mirrored term read as a row) beats the column pass, whose
+  // accumulator costs more to set up than these columns cost to read.
+  static constexpr Index PerColumnUpTo = 4;
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE RPacket lanes(const Packet& p) { return p; }
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real abs(const Scalar& x) { return numext::abs(x); }
+  // Nothing to record: pabs is exact.
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real component(const Scalar&) { return Real(0); }
+  static EIGEN_DEVICE_FUNC bool inRange(Real, Index) { return true; }
+
+  // The running sums of a pass over two columns.
+  struct Pass {
+    RPacket acc0 = pzero(RPacket());
+    RPacket acc1 = pzero(RPacket());
+    template <typename SumsEvaluator>
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void step(const RPacket& u, const RPacket& v, SumsEvaluator& s, Index i) {
+      RPacket a = pabs(u);
+      RPacket b = pabs(v);
+      acc0 = padd(acc0, a);
+      acc1 = padd(acc1, b);
+      s.template writePacket<Unaligned>(i, padd(s.template packet<Unaligned, Packet>(i), padd(a, b)));
+    }
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real sum0() const { return predux(acc0); }
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real sum1() const { return predux(acc1); }
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real peak() const { return Real(0); }
+  };
+};
+
+template <typename T>
+struct selfadjoint_l1norm_complex_lanes {
+  using Scalar = std::complex<T>;
+  using Real = T;
+  using Packet = typename packet_traits<Scalar>::type;
+  using RPacket = typename unpacket_traits<Packet>::as_real;
+  static constexpr Index PacketSize = unpacket_traits<Packet>::size;
+  static constexpr Index PerColumnUpTo = 0;
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE RPacket lanes(const Packet& p) { return p.v; }
+  // Same formula as the packets, for the diagonal and the tails: hypot costs more than the packets
+  // spend on the rest of a short column.
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real abs(const Scalar& z) { return numext::sqrt(numext::abs2(z)); }
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real component(const Scalar& z) {
+    return numext::maxi(numext::abs(numext::real(z)), numext::abs(numext::imag(z)));
+  }
+  // Components this large overflow when squared, and below the lower bound the squares lose
+  // precision the sum of n of them cannot hide.
+  static EIGEN_DEVICE_FUNC bool inRange(Real peak, Index n) {
+    Real tiny = Real(n) * numext::sqrt((std::numeric_limits<Real>::min)()) / NumTraits<Real>::epsilon();
+    Real huge = numext::sqrt(NumTraits<Real>::highest()) / Real(2);
+    return peak > tiny && peak < huge;
+  }
+
+  struct Pass {
+    RPacket acc = pzero(RPacket());    // |u| in the even lanes, |v| in the odd ones
+    RPacket peak_ = pzero(RPacket());  // the largest component seen
+    template <typename SumsEvaluator>
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void step(const RPacket& u, const RPacket& v, SumsEvaluator& s, Index i) {
+      peak_ = pmax(peak_, pmax(pabs(u), pabs(v)));
+      RPacket r = psqrt(pselect(peven_mask(u), abs2(u), abs2(v)));  // |u0| |v0| |u1| |v1| ...
+      acc = padd(acc, r);
+      s.template writePacket<Unaligned>(i,
+                                        Packet(padd(lanes(s.template packet<Unaligned, Packet>(i)), padd(r, flip(r)))));
+    }
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real sum0() const { return numext::real(predux(Packet(acc))); }
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real sum1() const { return numext::imag(predux(Packet(acc))); }
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real peak() const { return predux_max(peak_); }
+  };
+
+ private:
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE RPacket flip(const RPacket& r) { return pcplxflip(Packet(r)).v; }
+  // |z|^2 in both lanes of its slot.
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE RPacket abs2(const RPacket& v) {
+    RPacket s = pmul(v, v);
+    return padd(s, flip(s));
+  }
+};
+
+// The pass over two columns: packets over the shared rows, coefficients for the tail. An instance
+// remembers the largest component it has seen, for inRange().
+template <typename Lanes>
+struct selfadjoint_l1norm_packet_impl : Lanes {
+  using Lanes::PacketSize;
+  using typename Lanes::Packet;
+  using typename Lanes::Real;
+  using typename Lanes::RPacket;
+  using typename Lanes::Scalar;
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real abs(const Scalar& x) {
+    m_peak = numext::maxi(m_peak, Lanes::component(x));
+    return Lanes::abs(x);
+  }
+  EIGEN_DEVICE_FUNC bool inRange(Index n) const { return Lanes::inRange(m_peak, n); }
+
+  template <typename SumsDerived, typename Derived>
+  EIGEN_DEVICE_FUNC void accumulate(DenseBase<SumsDerived>& sums, const DenseBase<Derived>& m, Index j0, Index j1,
+                                    Index begin, Index end) {
+    accumulateCast(sums, j0, j1, begin, m.col(j0).segment(begin, end - begin).template cast<Scalar>(),
+                   m.col(j1).segment(begin, end - begin).template cast<Scalar>());
+  }
+
+ private:
+  Real m_peak = Real(0);
+
+  template <typename SumsDerived, typename Derived0, typename Derived1>
+  EIGEN_DEVICE_FUNC void accumulateCast(DenseBase<SumsDerived>& sums, Index j0, Index j1, Index begin,
+                                        const DenseBase<Derived0>& x0, const DenseBase<Derived1>& x1) {
+    using SumsEvaluator = evaluator<SumsDerived>;
+    using Evaluator0 = evaluator<Derived0>;
+    using Evaluator1 = evaluator<Derived1>;
+    constexpr int Needed = PacketAccessBit | LinearAccessBit;
+    constexpr bool Vectorize = (SumsEvaluator::Flags & Needed) == Needed && (Evaluator0::Flags & Needed) == Needed &&
+                               (Evaluator1::Flags & Needed) == Needed;
+    SumsEvaluator s(sums.derived());
+    accumulate(s, j0, j1, begin, Evaluator0(x0.derived()), Evaluator1(x1.derived()), x0.size(),
+               bool_constant<Vectorize>());
+  }
+  template <typename Evaluator>
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE RPacket load(const Evaluator& x, Index i) {
+    return Lanes::lanes(x.template packet<Unaligned, Packet>(i));
+  }
+  template <typename SumsEvaluator, typename Evaluator0, typename Evaluator1>
+  EIGEN_DEVICE_FUNC void accumulate(SumsEvaluator& s, Index j0, Index j1, Index begin, const Evaluator0& x0,
+                                    const Evaluator1& x1, Index from, Index to) {
+    Real sum0 = Real(0);
+    Real sum1 = Real(0);
+    for (Index i = from; i < to; ++i) {
+      Real a = abs(x0.coeff(i));
+      Real b = abs(x1.coeff(i));
+      s.coeffRef(begin + i) += a + b;
+      sum0 += a;
+      sum1 += b;
+    }
+    s.coeffRef(j0) += sum0;
+    s.coeffRef(j1) += sum1;
+  }
+  template <typename SumsEvaluator, typename Evaluator0, typename Evaluator1>
+  EIGEN_DEVICE_FUNC void accumulate(SumsEvaluator& s, Index j0, Index j1, Index begin, const Evaluator0& x0,
+                                    const Evaluator1& x1, Index n, std::false_type) {
+    accumulate(s, j0, j1, begin, x0, x1, Index(0), n);
+  }
+  template <typename SumsEvaluator, typename Evaluator0, typename Evaluator1>
+  EIGEN_DEVICE_FUNC void accumulate(SumsEvaluator& s, Index j0, Index j1, Index begin, const Evaluator0& x0,
+                                    const Evaluator1& x1, Index n, std::true_type) {
+    if (n < PacketSize) return accumulate(s, j0, j1, begin, x0, x1, Index(0), n);
+    typename Lanes::Pass pass;
+    Index i = 0;
+    for (; i + PacketSize <= n; i += PacketSize) pass.step(load(x0, i), load(x1, i), s, begin + i);
+    accumulate(s, j0, j1, begin, x0, x1, i, n);
+    s.coeffRef(j0) += pass.sum0();
+    s.coeffRef(j1) += pass.sum1();
+    m_peak = numext::maxi(m_peak, pass.peak());
+  }
+};
+
+// Coefficient fallback: custom complex types, or complex packets without a plain real view.
+template <typename Scalar_, typename Enable = void>
+struct selfadjoint_l1norm_impl {
+  using Scalar = Scalar_;
+  using Real = typename NumTraits<Scalar>::Real;
+  static constexpr Index PerColumnUpTo = 16;
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real abs(const Scalar& x) const { return numext::abs(x); }
+  EIGEN_DEVICE_FUNC bool inRange(Index) const { return true; }
+  template <typename SumsDerived, typename Derived>
+  EIGEN_DEVICE_FUNC void accumulate(DenseBase<SumsDerived>& sums, const DenseBase<Derived>& m, Index j0, Index j1,
+                                    Index begin, Index end) const {
+    Real sum0 = Real(0);
+    Real sum1 = Real(0);
+    for (Index i = begin; i < end; ++i) {
+      Real a = numext::abs(m.coeff(i, j0));
+      Real b = numext::abs(m.coeff(i, j1));
+      sums.coeffRef(i) += Scalar(a + b);
+      sum0 += a;
+      sum1 += b;
+    }
+    sums.coeffRef(j0) += Scalar(sum0);
+    sums.coeffRef(j1) += Scalar(sum1);
+  }
+};
+// half and bfloat16 accumulate in float, as stableNorm does.
+template <typename Scalar>
+struct selfadjoint_l1norm_impl<Scalar, std::enable_if_t<!NumTraits<Scalar>::IsComplex>>
+    : selfadjoint_l1norm_packet_impl<selfadjoint_l1norm_real_lanes<typename stable_norm_accumulator<Scalar>::type>> {};
+// The real view must lay the components out one per lane: Z13 stores four floats in two double
+// packets, which the lane masks do not describe.
+template <typename Packet, typename Enable = void>
+struct selfadjoint_l1norm_plain_real_view : std::false_type {};
+template <typename Packet>
+struct selfadjoint_l1norm_plain_real_view<Packet, void_t<typename unpacket_traits<Packet>::as_real>>
+    : bool_constant<sizeof(typename unpacket_traits<Packet>::as_real) ==
+                    unpacket_traits<Packet>::size * sizeof(typename unpacket_traits<Packet>::type)> {};
+template <typename T>
+struct selfadjoint_l1norm_impl<
+    std::complex<T>,
+    std::enable_if_t<selfadjoint_l1norm_plain_real_view<typename packet_traits<std::complex<T>>::type>::value>>
+    : selfadjoint_l1norm_packet_impl<selfadjoint_l1norm_complex_lanes<T>> {};
+
 template <typename MatrixType, unsigned int UpLo>
-struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType> {
+struct traits<SelfAdjointView<MatrixType, UpLo>> : traits<MatrixType> {
   using MatrixTypeNested = typename ref_selector<MatrixType>::non_const_type;
   using MatrixTypeNestedCleaned = remove_all_t<MatrixTypeNested>;
   using ExpressionType = MatrixType;
@@ -217,66 +433,81 @@
    */
   EIGEN_DEVICE_FUNC RealScalar l1Norm() const {
 #ifdef EIGEN_GPU_COMPILE_PHASE
-    // The panel accumulator below is per-thread local storage on a device, so it would cost every
-    // kernel instantiating this kPanelSize scalars of stack and the registers to address them.
+    // No per-thread accumulator on a device.
     return l1NormPerColumn();
 #else
-    // For a self-adjoint matrix |a_ij| = |a_ji|, so the stored triangle of a row-major matrix is
-    // the transposed, column-major, complementary one and yields the same norm read the fast way.
+    if (m_matrix.rows() <= L1NormImpl::PerColumnUpTo) return l1NormPerColumn();
+    // The stored triangle of a row-major matrix is the complementary triangle of its column-major
+    // transpose, which has the same norm.
     EIGEN_IF_CONSTEXPR (bool(MatrixType::IsRowMajor)) {
-      return l1NormColumnwise<TransposeMode>(m_matrix.transpose());
+      return l1NormStreaming<TransposeMode>(m_matrix.transpose());
     } else {
-      return l1NormColumnwise<UpLo>(m_matrix);
+      return l1NormStreaming<UpLo>(m_matrix);
     }
 #endif
   }
 
  private:
-  // Reading the mirrored term of column j as a row of the stored triangle costs a stride-n
-  // traversal of a column-major matrix. Instead accumulate column sums a panel at a time:
-  // |a_ij| from a column left of the panel is added to the sum of column i, which walks that
-  // column. Only the panel's diagonal block keeps the row traversal, where it is cache resident,
-  // and a panel-sized accumulator stays a stack object.
+  using L1NormImpl = internal::selfadjoint_l1norm_impl<Scalar>;
+  // float for half and bfloat16, Scalar otherwise.
+  using L1NormScalar = typename L1NormImpl::Scalar;
+  using L1NormAccumulator = typename L1NormImpl::Real;
+
+  // Each column is read once, top to bottom, two at a time: |a_ij| goes to column j's sum and, as
+  // the mirrored a_ji, to sums[i]. Lower walks the columns forward and Upper backward so that
+  // sums[j] is complete when column j is reached. Of a pair (j0, j1) only j0's element in row j1
+  // lies outside the rows the two share.
   template <int Mode, typename Mat>
-  EIGEN_DEVICE_FUNC static RealScalar l1NormColumnwise(const Mat& m) {
-    static constexpr int kPanelSize = 64;
-    RealScalar norm = RealScalar(0);
+  RealScalar l1NormStreaming(const Mat& m) const {
     const Index n = m.rows();
-    Matrix<RealScalar, kPanelSize, 1> sums;
-    for (Index p = 0; p < n; p += kPanelSize) {
-      const Index len = numext::mini(Index(kPanelSize), n - p);
-      EIGEN_IF_CONSTEXPR (Mode == Lower) {
-        for (Index j = 0; j < len; ++j)
-          sums.coeffRef(j) =
-              m.col(p + j).tail(n - p - j).template lpNorm<1>() + m.row(p + j).segment(p, j).template lpNorm<1>();
-        for (Index j = 0; j < p; ++j) sums.head(len) += m.col(j).segment(p, len).cwiseAbs();
-      } else {
-        for (Index j = 0; j < len; ++j)
-          sums.coeffRef(j) = m.col(p + j).head(p + j + 1).template lpNorm<1>() +
-                             m.row(p + j).segment(p + j + 1, len - j - 1).template lpNorm<1>();
-        for (Index j = p + len; j < n; ++j) sums.head(len) += m.col(j).segment(p, len).cwiseAbs();
-      }
-      norm = numext::maxi(norm, sums.head(len).maxCoeff());
+    // The accumulator lives in the object for bounded sizes and on the stack otherwise, so that
+    // neither fixed-size nor preallocated dynamic-size decompositions allocate.
+    internal::gemv_static_vector_if<L1NormScalar, Mat::RowsAtCompileTime, Mat::MaxRowsAtCompileTime, true> static_sums;
+    ei_declare_aligned_stack_constructed_variable(L1NormScalar, sums_data, n, static_sums.data());
+    Map<Matrix<L1NormScalar, Dynamic, 1>> sums(sums_data, n);
+    sums.setZero();
+    L1NormImpl impl;
+    L1NormAccumulator norm = L1NormAccumulator(0);
+    Index k = 0;
+    for (; k + 1 < n; k += 2) {
+      Index j0 = Mode == Lower ? k : n - 1 - k;
+      Index j1 = Mode == Lower ? j0 + 1 : j0 - 1;
+      Index rowBegin = Mode == Lower ? j1 + 1 : 0;
+      Index rowEnd = Mode == Lower ? n : j1;
+      impl.accumulate(sums, m, j0, j1, rowBegin, rowEnd);
+      // The element of j0 in row j1 lies outside the shared rows: it counts for both columns.
+      L1NormAccumulator boundary = impl.abs(m.coeff(j1, j0));
+      // Totals are materialized so that maxi compares two accumulators (an integer sum promotes,
+      // an autodiff sum is an expression).
+      L1NormAccumulator col0 = numext::real(sums.coeff(j0)) + impl.abs(m.coeff(j0, j0)) + boundary;
+      L1NormAccumulator col1 = numext::real(sums.coeff(j1)) + impl.abs(m.coeff(j1, j1)) + boundary;
+      norm = numext::maxi(norm, col0);
+      norm = numext::maxi(norm, col1);
     }
-    return norm;
+    if (k < n) {
+      Index j = Mode == Lower ? k : 0;
+      L1NormAccumulator col = numext::real(sums.coeff(j)) + impl.abs(m.coeff(j, j));
+      norm = numext::maxi(norm, col);
+    }
+    return impl.inRange(n) ? RealScalar(norm) : l1NormPerColumn();
   }
 
-  // Workspace-free form, one column sum at a time; the mirrored term is read as a row.
+  // One column at a time, the mirrored term read as a row; no workspace.
   EIGEN_DEVICE_FUNC RealScalar l1NormPerColumn() const {
-    RealScalar norm = RealScalar(0);
+    L1NormAccumulator norm = L1NormAccumulator(0);
     const Index n = m_matrix.rows();
     for (Index col = 0; col < n; ++col) {
-      RealScalar abs_col_sum;
+      L1NormAccumulator abs_col_sum;
       EIGEN_IF_CONSTEXPR (UpLo == Lower) {
-        abs_col_sum =
-            m_matrix.col(col).tail(n - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
+        abs_col_sum = m_matrix.col(col).tail(n - col).template cast<L1NormScalar>().template lpNorm<1>() +
+                      m_matrix.row(col).head(col).template cast<L1NormScalar>().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>();
+        abs_col_sum = m_matrix.col(col).head(col).template cast<L1NormScalar>().template lpNorm<1>() +
+                      m_matrix.row(col).tail(n - col).template cast<L1NormScalar>().template lpNorm<1>();
       }
       norm = numext::maxi(norm, abs_col_sum);
     }
-    return norm;
+    return RealScalar(norm);
   }
 
  public:

contrib/test/autodiff.cpp

diff --git a/contrib/test/autodiff.cpp b/contrib/test/autodiff.cpp
index 7769154..582fcfd 100644
--- a/contrib/test/autodiff.cpp
+++ b/contrib/test/autodiff.cpp
@@ -332,6 +332,19 @@
   return (y1 + y2 + y3).value();
 }
 
+// LLT stores the 1-norm of its input for rcond(); the sum of AutoDiffScalars is an expression, so
+// the norm must materialize its column totals before comparing them.
+void test_autodiff_selfadjoint_l1norm() {
+  typedef AutoDiffScalar<Vector2d> AD;
+  Matrix<AD, 8, 8> m = Matrix<AD, 8, 8>::Identity();
+  m(3, 3) = AD(2, Vector2d(1, 0));
+  AD norm = m.selfadjointView<Lower>().l1Norm();
+  VERIFY_IS_EQUAL(norm.value(), 2.0);
+  VERIFY_IS_EQUAL(norm.derivatives()(0), 1.0);
+  LLT<Matrix<AD, 8, 8>> llt(m);
+  VERIFY(llt.info() == Success);
+}
+
 EIGEN_DECLARE_TEST(autodiff) {
   for (int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1(test_autodiff_scalar<1>());
@@ -345,4 +358,5 @@
   CALL_SUBTEST_5(bug_1260());
   CALL_SUBTEST_5(bug_1261());
   CALL_SUBTEST_5(bug_1281());
+  CALL_SUBTEST_5(test_autodiff_selfadjoint_l1norm());
 }

test/cholesky.cpp

diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 856ffa5..802ceeb 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -9,6 +9,7 @@
 // SPDX-License-Identifier: MPL-2.0
 
 #define TEST_ENABLE_TEMPORARY_TRACKING
+#define EIGEN_RUNTIME_NO_MALLOC
 
 #include "main.h"
 #include <Eigen/Cholesky>
@@ -709,7 +710,7 @@
 
   const RealScalar absdet = d.prod();
   const RealScalar logabsdet = d.array().log().sum();
-  const MatrixType spd = q * d.template cast<Scalar>().asDiagonal() * q.adjoint();
+  MatrixType spd = q * d.template cast<Scalar>().asDiagonal() * q.adjoint();
 
   LLT<MatrixType, Lower> lltlo(spd);
   VERIFY(lltlo.info() == Success);
@@ -911,10 +912,32 @@
   }
 }
 
+// Preallocated decompositions of dynamic-size matrices stay allocation-free, including the 1-norm
+// they take for rcond(), whose workspace comes from the stack.
+template <typename Scalar>
+void cholesky_dynamic_preallocated_no_malloc() {
+  typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
+  Index size = 8;
+  MatrixType A = MatrixType::Random(size, size);
+  MatrixType spd = A * A.adjoint() + MatrixType::Identity(size, size) * Scalar(size);
+  LLT<MatrixType> llt(size);
+  LDLT<MatrixType> ldlt(size);
+  internal::set_is_malloc_allowed(false);
+  llt.compute(spd);
+  VERIFY_IS_EQUAL(llt.info(), Success);
+  ldlt.compute(spd);
+  VERIFY_IS_EQUAL(ldlt.info(), Success);
+  VERIFY(spd.template selfadjointView<Lower>().l1Norm() > 0);
+  internal::set_is_malloc_allowed(true);
+  VERIFY(llt.rcond() > 0);
+}
+
 EIGEN_DECLARE_TEST(cholesky) {
   int s = 0;
   for (int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1(cholesky(Matrix<double, 1, 1>()));
+    CALL_SUBTEST_1(cholesky_dynamic_preallocated_no_malloc<double>());
+    CALL_SUBTEST_1(cholesky_dynamic_preallocated_no_malloc<std::complex<double> >());
     CALL_SUBTEST_3(cholesky(Matrix2d()));
     CALL_SUBTEST_3(cholesky_bug241(Matrix2d()));
     CALL_SUBTEST_3(cholesky_definiteness(Matrix2d()));

test/fastmath.cpp

diff --git a/test/fastmath.cpp b/test/fastmath.cpp
index 9056b72..98fcbea 100644
--- a/test/fastmath.cpp
+++ b/test/fastmath.cpp
@@ -301,12 +301,31 @@
   VERIFY_IS_APPROX(mat, q * r);
 }
 
+// The complex self-adjoint 1-norm takes sqrt(re^2 + im^2) in packets and falls back to the scalar
+// abs when a component is too large to square. The decision must not depend on an infinity or
+// NaN surviving fast-math: an approximate packet sqrt can turn the overflow into a NaN that a
+// maximum then discards. Large entries both on the diagonal and inside a packet of a column.
+template <typename RealScalar>
+void check_complex_selfadjoint_l1norm() {
+  typedef std::complex<RealScalar> Scalar;
+  RealScalar big = numext::sqrt(NumTraits<RealScalar>::highest()) * RealScalar(1e3);
+  Matrix<Scalar, 8, 8> m = Matrix<Scalar, 8, 8>::Identity() * big;
+  VERIFY_IS_APPROX(m.template selfadjointView<Lower>().l1Norm(), big);
+  VERIFY_IS_APPROX(m.template selfadjointView<Upper>().l1Norm(), big);
+  m.setIdentity();
+  m(2, 0) = m(0, 2) = Scalar(big, big);
+  RealScalar expected = numext::abs(Scalar(big, big)) + RealScalar(1);
+  VERIFY_IS_APPROX(m.template selfadjointView<Lower>().l1Norm(), expected);
+  VERIFY_IS_APPROX(m.template selfadjointView<Upper>().l1Norm(), expected);
+}
+
 template <typename RealScalar>
 void check_complex_fastmath() {
   check_complex_rowmajor_adjoint_product<RealScalar>();
   check_complex_packet_arithmetic<RealScalar>();
   check_complex_packet_math_functions<RealScalar>();
   check_complex_householder_qr<RealScalar>();
+  check_complex_selfadjoint_l1norm<RealScalar>();
 }
 
 // The packet implementations of these functions manipulate signs of non-zero values through a

test/selfadjoint.cpp

diff --git a/test/selfadjoint.cpp b/test/selfadjoint.cpp
index 581eb18..708b102 100644
--- a/test/selfadjoint.cpp
+++ b/test/selfadjoint.cpp
@@ -75,13 +75,15 @@
   VERIFY_IS_APPROX(lowerOnly.template selfadjointView<Lower>().l1Norm(), ref_l1);
 }
 
-// l1Norm accumulates the column sums a panel at a time, so sweep sizes across two panel
-// boundaries: a mis-sized segment would otherwise hide between the random sizes above.
+// l1Norm switches from a per-column to a streaming form at a small size and its column pass has
+// packet tails, so sweep sizes around the switch and across packet boundaries: a mis-sized
+// segment would otherwise hide between the random sizes above.
 template <typename Scalar>
 void selfadjoint_l1norm_sizes() {
   typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
   typedef typename NumTraits<Scalar>::Real RealScalar;
-  for (Index n : {Index(0), Index(1), Index(2), Index(63), Index(64), Index(65), Index(127), Index(128), Index(129)}) {
+  for (Index n : {Index(0), Index(1), Index(2), Index(15), Index(16), Index(17), Index(63), Index(64), Index(65),
+                  Index(127), Index(128), Index(129)}) {
     MatrixType m = MatrixType::Random(n, n);
     MatrixType full = m.template selfadjointView<Lower>();
     RealScalar ref = n == 0 ? RealScalar(0) : full.cwiseAbs().colwise().sum().maxCoeff();
@@ -90,6 +92,55 @@
   }
 }
 
+// The vectorized complex path squares the parts, so entries beyond the square-root range of the
+// scalar must come back through the scalar fallback: squares that overflow, and squares that
+// land in the denormals and lose precision.
+template <typename Scalar>
+void selfadjoint_l1norm_range() {
+  typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  Index n = 70;
+  MatrixType m = MatrixType::Random(n, n);
+  RealScalar big = numext::sqrt(NumTraits<RealScalar>::highest()) * RealScalar(1e3);
+  RealScalar small = numext::sqrt((std::numeric_limits<RealScalar>::min)()) * RealScalar(1e-3);
+  for (RealScalar scale : {big, small}) {
+    MatrixType ms = m * scale;
+    MatrixType full = ms.template selfadjointView<Lower>();
+    RealScalar ref = full.cwiseAbs().colwise().sum().maxCoeff();
+    VERIFY_IS_APPROX(ms.template selfadjointView<Lower>().l1Norm(), ref);
+    VERIFY_IS_APPROX(full.template selfadjointView<Upper>().l1Norm(), ref);
+  }
+}
+
+// half and bfloat16 accumulate the norm in float, so only the final rounding to the scalar separates
+// the result from a float reference, not the size of the matrix.
+template <typename Scalar>
+void selfadjoint_l1norm_lowprec() {
+  typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
+  for (Index n : {Index(8), Index(64), Index(300)}) {
+    MatrixType m = MatrixType::Random(n, n).template selfadjointView<Lower>();
+    float ref = m.template cast<float>().cwiseAbs().colwise().sum().maxCoeff();
+    float tol = 2 * float(NumTraits<Scalar>::epsilon()) * ref;
+    VERIFY(numext::abs(float(m.template selfadjointView<Lower>().l1Norm()) - ref) <= tol);
+    VERIFY(numext::abs(float(m.template selfadjointView<Upper>().l1Norm()) - ref) <= tol);
+  }
+}
+
+// Narrow integers promote when added: the column totals must be materialized in the scalar type
+// before they are compared (n=2 takes the per-column form, n=8 the column pass).
+template <typename Scalar>
+void selfadjoint_l1norm_integer() {
+  for (Index n : {Index(2), Index(8)}) {
+    typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
+    // Random() spans the whole range; keep the column sums representable.
+    MatrixType m = (MatrixType::Random(n, n) / Scalar(NumTraits<Scalar>::highest() / 16)).eval();
+    m = m.template selfadjointView<Lower>();
+    Scalar ref = m.template cast<int>().cwiseAbs().colwise().sum().maxCoeff();
+    VERIFY_IS_EQUAL(m.template selfadjointView<Lower>().l1Norm(), ref);
+    VERIFY_IS_EQUAL(m.template selfadjointView<Upper>().l1Norm(), ref);
+  }
+}
+
 void bug_159() {
   Matrix3d m = Matrix3d::Random().selfadjointView<Lower>();
   EIGEN_UNUSED_VARIABLE(m);
@@ -104,6 +155,7 @@
     CALL_SUBTEST_3(selfadjoint(Matrix3cf()));
     CALL_SUBTEST_4(selfadjoint(MatrixXcd(s, s)));
     CALL_SUBTEST_5(selfadjoint(Matrix<float, Dynamic, Dynamic, RowMajor>(s, s)));
+    CALL_SUBTEST_6(selfadjoint(Matrix<std::complex<float>, Dynamic, Dynamic, RowMajor>(s, s)));
 
     TEST_SET_BUT_UNUSED_VARIABLE(s);
   }
@@ -111,4 +163,11 @@
   CALL_SUBTEST_1(bug_159());
   CALL_SUBTEST_4(selfadjoint_l1norm_sizes<double>());
   CALL_SUBTEST_4(selfadjoint_l1norm_sizes<std::complex<double> >());
+  CALL_SUBTEST_6(selfadjoint_l1norm_sizes<std::complex<float> >());
+  CALL_SUBTEST_4(selfadjoint_l1norm_range<std::complex<double> >());
+  CALL_SUBTEST_6(selfadjoint_l1norm_range<std::complex<float> >());
+  CALL_SUBTEST_1(selfadjoint_l1norm_integer<short>());
+  CALL_SUBTEST_1(selfadjoint_l1norm_integer<int>());
+  CALL_SUBTEST_7(selfadjoint_l1norm_lowprec<half>());
+  CALL_SUBTEST_7(selfadjoint_l1norm_lowprec<bfloat16>());
 }