Cholesky: Keep blocked complex LLT inverse diagonals exactly real

libeigen/eigen!3031

Closes #3149

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

Eigen/src/Cholesky/LLT.h

diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 65545b3..64340d4 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -636,8 +636,8 @@
       llt.solveInPlace(dst);
     }
     // Mirror; (i, j) reads (j, i), which lies in the computed triangle and is never written here, so
-    // the aliasing is benign. The computed diagonal is exactly real (a squaredNorm above the
-    // threshold, a real scalar below it), so the result is exactly self-adjoint.
+    // the aliasing is benign. LAUUM makes the computed diagonal exactly real, as does the real-scalar
+    // solve below the threshold, so the result is exactly self-adjoint.
     dst.template triangularView<kMirrorMode>() = dst.adjoint();
   }
 };

Eigen/src/Core/TriangularInPlace.h

diff --git a/Eigen/src/Core/TriangularInPlace.h b/Eigen/src/Core/TriangularInPlace.h
index a5e28d7..69499ba 100644
--- a/Eigen/src/Core/TriangularInPlace.h
+++ b/Eigen/src/Core/TriangularInPlace.h
@@ -176,7 +176,13 @@
       if (rs > 0) B.noalias() += mat.block(k + bs, k, rs, bs).adjoint() * mat.block(k + bs, 0, rs, k);
     }
     triangular_adjoint_square_unblocked(L11);
-    if (rs > 0) L11.template selfadjointView<Lower>().rankUpdate(mat.block(k + bs, k, rs, bs).adjoint());
+    if (rs > 0) {
+      L11.template selfadjointView<Lower>().rankUpdate(mat.block(k + bs, k, rs, bs).adjoint());
+      // Fused complex products can leave a rounding residual in the imaginary diagonal.
+      EIGEN_IF_CONSTEXPR (NumTraits<typename MatrixType::Scalar>::IsComplex) {
+        L11.diagonal() = L11.diagonal().real().template cast<typename MatrixType::Scalar>();
+      }
+    }
   }
 }
 

test/triangular_in_place.cpp

diff --git a/test/triangular_in_place.cpp b/test/triangular_in_place.cpp
index 1b06ea5..e1e510b 100644
--- a/test/triangular_in_place.cpp
+++ b/test/triangular_in_place.cpp
@@ -145,6 +145,7 @@
   const MatrixType t = ill_conditioned_triangular<MatrixType>(n);
   MatrixType a = t;
   internal::triangular_adjoint_square_in_place<UpLo>(a);
+  VERIFY((a.diagonal().imag().array() == RealScalar(0)).all());
 
   const DenseType td = t.template triangularView<UpLo>().toDenseMatrix();
   const DenseType expected = UpLo == Lower ? DenseType(td.adjoint() * td) : DenseType(td * td.adjoint());
@@ -180,6 +181,7 @@
       CALL_SUBTEST_5(triangular_inverse_on_strided_map<MatrixXd>(n));
       CALL_SUBTEST_7(adjoint_square_both_triangles<MatrixXd>(n));
       CALL_SUBTEST_8(adjoint_square_both_triangles<MatrixXcd>(n));
+      CALL_SUBTEST_8((adjoint_square_both_triangles<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>(n)));
       CALL_SUBTEST_9(adjoint_square_both_triangles<MatrixXf>(n));
       CALL_SUBTEST_10((adjoint_square_both_triangles<Matrix<double, Dynamic, Dynamic, RowMajor>>(n)));
     }