Eigen/src/Cholesky/LLT.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0

#ifndef EIGEN_LLT_H
#define EIGEN_LLT_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

// Smallest n at which LLT::inverse() runs the POTRI sequence on a real scalar rather than solving
// against an explicit identity. POTRI's n^3/3 runs in the unblocked TRTRI and LAUUM kernels at small
// n, the solve's 3x at the blocked TRSM rate, so a wider ISA helps the competitor more and moves the
// crossover right: on Zen 4 it is n = 32 under SSE2 and AVX2, 128 (double) to 256 (float) under
// AVX-512. Complex is not thresholded: POTRI leads from n = 8 in all twelve configurations measured.
#ifndef EIGEN_LLT_INVERSE_POTRI_THRESHOLD
#if defined(EIGEN_VECTORIZE_AVX512)
#define EIGEN_LLT_INVERSE_POTRI_THRESHOLD 256
#else
#define EIGEN_LLT_INVERSE_POTRI_THRESHOLD 32
#endif
#endif

namespace Eigen {

namespace internal {

template <typename MatrixType_, int UpLo_>
struct traits<LLT<MatrixType_, UpLo_> > : traits<MatrixType_> {
  using XprKind = MatrixXpr;
  using StorageKind = SolverStorage;
  using StorageIndex = int;
  enum { Flags = 0 };
};

template <typename MatrixType, int UpLo>
struct LLT_Traits;
}  // namespace internal

/** \ingroup Cholesky_Module
 *
 * \class LLT
 *
 * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
 *
 * \tparam MatrixType_ the type of the matrix of which we are computing the LL^T Cholesky decomposition
 * \tparam UpLo_ the triangular part that will be used for the decomposition: Lower (default) or Upper.
 *               The other triangular part won't be read.
 *
 * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
 * matrix A such that A = LL^* = U^*U, where L is lower triangular.
 *
 * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
 * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
 * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
 * situations like generalised eigen problems with hermitian matrices.
 *
 * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive
 * definite matrices, use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine
 * whether a system of equations has a solution.
 *
 * Example: \include LLT_example.cpp
 * Output: \verbinclude LLT_example.out
 *
 * \b Performance: for best performance, it is recommended to use a column-major storage format
 * with the Lower triangular part (the default), or, equivalently, a row-major storage format
 * with the Upper triangular part. Otherwise, you might get a 20% slowdown for the full factorization
 * step, and rank-updates can be up to 3 times slower.
 *
 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
 *
 * Note that during the decomposition, only the lower (or upper, as defined by UpLo_) triangular part of A is
 * considered. Therefore, the strict upper part (or the strict lower part when UpLo_ is Upper) does not have to
 * store correct values.
 *
 * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT
 */
template <typename MatrixType_, int UpLo_>
class LLT : public SolverBase<LLT<MatrixType_, UpLo_> > {
 public:
  using MatrixType = MatrixType_;
  using Base = SolverBase<LLT>;
  friend class SolverBase<LLT>;

  EIGEN_GENERIC_PUBLIC_INTERFACE(LLT)
  enum { MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime };

  enum { PacketSize = internal::packet_traits<Scalar>::size, AlignmentMask = int(PacketSize) - 1, UpLo = UpLo_ };

  using Traits = internal::LLT_Traits<MatrixType, UpLo>;
  using PlainObject = typename MatrixType::PlainObject;

  /**
   * \brief Default Constructor.
   *
   * The default constructor is useful in cases in which the user intends to
   * perform decompositions via LLT::compute(const MatrixType&).
   */
  LLT() : m_matrix(), m_l1_norm(0), m_isInitialized(false), m_info(InvalidInput) {}

  /** \brief Default Constructor with memory preallocation
   *
   * Like the default constructor but with preallocation of the internal data
   * according to the specified problem \a size.
   * \sa LLT()
   */
  explicit LLT(Index size) : m_matrix(size, size), m_l1_norm(0), m_isInitialized(false), m_info(InvalidInput) {}

  template <typename InputType>
  explicit LLT(const EigenBase<InputType>& matrix)
      : m_matrix(matrix.rows(), matrix.cols()), m_l1_norm(0), m_isInitialized(false), m_info(InvalidInput) {
    compute(matrix.derived());
  }

  /** \brief Constructs a LLT factorization from a given matrix
   *
   * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when
   * \c MatrixType is a Eigen::Ref.
   *
   * \sa LLT(const EigenBase&)
   */
  template <typename InputType>
  explicit LLT(EigenBase<InputType>& matrix)
      : m_matrix(matrix.derived()), m_l1_norm(0), m_isInitialized(false), m_info(InvalidInput) {
    compute(matrix.derived());
  }

  /** \returns a view of the upper triangular matrix U */
  inline typename Traits::MatrixU matrixU() const {
    eigen_assert(m_isInitialized && "LLT is not initialized.");
    return Traits::getU(m_matrix);
  }

  /** \returns a view of the lower triangular matrix L */
  inline typename Traits::MatrixL matrixL() const {
    eigen_assert(m_isInitialized && "LLT is not initialized.");
    return Traits::getL(m_matrix);
  }

#ifdef EIGEN_PARSED_BY_DOXYGEN
  /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
   *
   * Since this LLT class assumes anyway that the matrix A is invertible, the solution
   * theoretically exists and is unique regardless of b.
   *
   * Example: \include LLT_solve.cpp
   * Output: \verbinclude LLT_solve.out
   *
   * \sa solveInPlace(), MatrixBase::llt(), SelfAdjointView::llt()
   */
  template <typename Rhs>
  inline Solve<LLT, Rhs> solve(const MatrixBase<Rhs>& b) const;
#endif

  template <typename Derived>
  void solveInPlace(const MatrixBase<Derived>& bAndX) const;

  template <typename InputType>
  LLT& compute(const EigenBase<InputType>& matrix);

  /** \returns an estimate of the reciprocal condition number of the matrix of
   *  which \c *this is the Cholesky decomposition.
   */
  RealScalar rcond() const {
    eigen_assert(m_isInitialized && "LLT is not initialized.");
    eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
    return internal::rcond_estimate_helper(m_l1_norm, *this);
  }

  /** \returns the LLT decomposition matrix
   *
   * TODO: document the storage layout
   */
  inline const MatrixType& matrixLLT() const {
    eigen_assert(m_isInitialized && "LLT is not initialized.");
    return m_matrix;
  }

  /** \returns the inverse of the matrix of which \c *this is the Cholesky decomposition.
   *
   * The result is computed as \f$ A^{-1} = L^{-*} L^{-1} \f$ by inverting the stored factor in place and
   * squaring it, the LAPACK \c *POTRI sequence, for 2n^3/3 flops against the 2n^3 of solving with an
   * explicit identity right hand side. Beyond the destination, only a block-sized scratch panel of at most
   * 128x128 coefficients is used, whatever the matrix size.
   *
   * Fewer flops is not fewer seconds at every size: for a real scalar below
   * \c EIGEN_LLT_INVERSE_POTRI_THRESHOLD (32, or 256 where AVX-512 is enabled) the POTRI sequence runs in
   * its unblocked kernels while the solve against an identity runs its 3x at the blocked TRSM rate, so that
   * is what this method does there. Complex scalars always take the POTRI path. The result is exactly
   * self-adjoint either way: one triangle is computed and mirrored onto the other.
   *
   * An in-place decomposition (see the class documentation) may overwrite its own factor with the result,
   * as in <tt>storage = llt.inverse()</tt>; like any other write to the referenced matrix, that leaves the
   * decomposition unusable afterwards.
   *
   * The matrix must be positive definite, that is, info() must be \c Success.
   *
   * \sa solve(), MatrixBase::inverse()
   */
  inline Inverse<LLT> inverse() const {
    eigen_assert(m_isInitialized && "LLT is not initialized.");
    return Inverse<LLT>(*this);
  }

  MatrixType reconstructedMatrix() const;

  /** \brief Reports whether previous computation was successful.
   *
   * \returns \c Success if computation was successful,
   *          \c NumericalIssue if the matrix appears not to be positive definite.
   */
  ComputationInfo info() const {
    eigen_assert(m_isInitialized && "LLT is not initialized.");
    return m_info;
  }

  /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix
   * is self-adjoint.
   *
   * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as:
   * \code x = decomposition.adjoint().solve(b) \endcode
   */
  const LLT& adjoint() const noexcept { return *this; }

  constexpr Index rows() const noexcept { return m_matrix.rows(); }
  constexpr Index cols() const noexcept { return m_matrix.cols(); }

  template <typename VectorType>
  LLT& rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);

  /** \returns the determinant of the matrix of which *this is the Cholesky decomposition.
   *
   * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix)
   * as the Cholesky decomposition has already been computed.
   *
   * \warning a determinant can be very big or small, so for matrices
   * of large enough dimension, there is a risk of overflow/underflow.
   * One way to work around that is to use logAbsDeterminant() instead.
   *
   * \pre info() returns \c Success. A failed factorization does not represent the input matrix.
   *
   * \sa absDeterminant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant()
   */
  Scalar determinant() const;

  /** \returns the absolute value of the determinant of the matrix of which *this is the Cholesky decomposition.
   *
   * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix)
   * as the Cholesky decomposition has already been computed.
   *
   * \note The decomposed matrix is positive definite, so this is the determinant itself.
   *
   * \warning a determinant can be very big or small, so for matrices
   * of large enough dimension, there is a risk of overflow/underflow.
   * One way to work around that is to use logAbsDeterminant() instead.
   *
   * \pre info() returns \c Success. A failed factorization does not represent the input matrix.
   *
   * \sa determinant(), logAbsDeterminant(), signDeterminant(), MatrixBase::determinant()
   */
  RealScalar absDeterminant() const;

  /** \returns the natural log of the absolute value of the determinant of the matrix of which *this is the Cholesky
   * decomposition.
   *
   * It has only linear complexity (that is, O(n) where n is the dimension of the square matrix)
   * as the Cholesky decomposition has already been computed.
   *
   * \note This method is useful to work around the risk of overflow/underflow that's inherent
   * to determinant computation.
   *
   * \pre info() returns \c Success. A failed factorization does not represent the input matrix.
   *
   * \sa determinant(), absDeterminant(), signDeterminant(), MatrixBase::determinant()
   */
  RealScalar logAbsDeterminant() const;

  /** \returns the sign of the determinant of the matrix of which *this is the Cholesky decomposition,
   * which is \c 1 since that matrix is positive definite.
   *
   * This method is provided for compatibility with the other decompositions, thus enabling generic code.
   *
   * \pre info() returns \c Success. A failed factorization does not represent the input matrix.
   *
   * \sa determinant(), absDeterminant(), logAbsDeterminant(), MatrixBase::determinant()
   */
  Scalar signDeterminant() const;

#ifndef EIGEN_PARSED_BY_DOXYGEN
  template <typename RhsType, typename DstType>
  void _solve_impl(const RhsType& rhs, DstType& dst) const;

  template <bool Conjugate, typename RhsType, typename DstType>
  void _solve_impl_transposed(const RhsType& rhs, DstType& dst) const;
#endif

 protected:
  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)

  /** \internal
   * Used to compute and store L, or U when UpLo_ is Upper.
   * The strict part of the other triangle is not used and even not initialized.
   */
  MatrixType m_matrix;
  RealScalar m_l1_norm;
  bool m_isInitialized;
  ComputationInfo m_info;
};

namespace internal {

template <typename Scalar, int UpLo>
struct llt_inplace;

template <typename MatrixType, typename VectorType>
static Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec,
                                   const typename MatrixType::RealScalar& sigma) {
  using std::sqrt;
  using Scalar = typename MatrixType::Scalar;
  using RealScalar = typename MatrixType::RealScalar;
  using ColXpr = typename MatrixType::ColXpr;
  using ColXprCleaned = internal::remove_all_t<ColXpr>;
  using ColXprSegment = typename ColXprCleaned::SegmentReturnType;
  using TempVectorType = Matrix<Scalar, MatrixType::RowsAtCompileTime, 1, 0, MatrixType::MaxRowsAtCompileTime, 1>;
  using TempVecSegment = typename TempVectorType::SegmentReturnType;

  Index n = mat.cols();
  eigen_assert(mat.rows() == n && vec.size() == n);

  TempVectorType temp;

  if (sigma > 0) {
    // This version is based on Givens rotations.
    // It is faster than the other one below, but only works for updates,
    // i.e., for sigma > 0
    temp = sqrt(sigma) * vec;

    for (Index i = 0; i < n; ++i) {
      JacobiRotation<Scalar> g;
      g.makeGivens(mat(i, i), -temp(i), &mat(i, i));

      Index rs = n - i - 1;
      if (rs > 0) {
        ColXprSegment x(mat.col(i).tail(rs));
        TempVecSegment y(temp.tail(rs));
        apply_rotation_in_the_plane(x, y, g);
      }
    }
  } else {
    temp = vec;
    RealScalar beta = 1;
    for (Index j = 0; j < n; ++j) {
      RealScalar Ljj = numext::real(mat.coeff(j, j));
      RealScalar dj = numext::abs2(Ljj);
      Scalar wj = temp.coeff(j);
      RealScalar swj2 = sigma * numext::abs2(wj);
      RealScalar gamma = dj * beta + swj2;

      RealScalar x = dj + swj2 / beta;
      if (x <= RealScalar(0)) return j;
      RealScalar nLjj = sqrt(x);
      mat.coeffRef(j, j) = nLjj;
      beta += swj2 / dj;

      // Update the terms of L
      Index rs = n - j - 1;
      if (rs) {
        temp.tail(rs) -= (wj / Ljj) * mat.col(j).tail(rs);
        if (!numext::is_exactly_zero(gamma))
          mat.col(j).tail(rs) =
              (nLjj / Ljj) * mat.col(j).tail(rs) + (nLjj * sigma * numext::conj(wj) / gamma) * temp.tail(rs);
      }
    }
  }
  return -1;
}

template <typename Scalar>
struct llt_inplace<Scalar, Lower> {
  using RealScalar = typename NumTraits<Scalar>::Real;
  template <typename MatrixType>
  static Index unblocked(MatrixType& mat) {
    using std::sqrt;

    eigen_assert(mat.rows() == mat.cols());
    const Index size = mat.rows();
    for (Index k = 0; k < size; ++k) {
      Index rs = size - k - 1;  // remaining size

      Block<MatrixType, Dynamic, 1> A21(mat, k + 1, k, rs, 1);
      Block<MatrixType, 1, Dynamic> A10(mat, k, 0, 1, k);
      Block<MatrixType, Dynamic, Dynamic> A20(mat, k + 1, 0, rs, k);

      RealScalar x = numext::real(mat.coeff(k, k));
      if (k > 0) x -= A10.squaredNorm();
      if (x <= RealScalar(0)) return k;
      mat.coeffRef(k, k) = x = sqrt(x);
      if (k > 0 && rs > 0) A21.noalias() -= A20 * A10.adjoint();
      if (rs > 0) A21 /= x;
    }
    return -1;
  }

  template <typename MatrixType>
  static Index blocked(MatrixType& m) {
    eigen_assert(m.rows() == m.cols());
    Index size = m.rows();
    if (size < 32) return unblocked(m);

    Index blockSize = size / 8;
    blockSize = (blockSize / 16) * 16;
    blockSize = (std::min)((std::max)(blockSize, Index(8)), Index(128));

    for (Index k = 0; k < size; k += blockSize) {
      // partition the matrix:
      //       A00 |  -  |  -
      // lu  = A10 | A11 |  -
      //       A20 | A21 | A22
      Index bs = (std::min)(blockSize, size - k);
      Index rs = size - k - bs;
      Block<MatrixType, Dynamic, Dynamic> A11(m, k, k, bs, bs);
      Block<MatrixType, Dynamic, Dynamic> A21(m, k + bs, k, rs, bs);
      Block<MatrixType, Dynamic, Dynamic> A22(m, k + bs, k + bs, rs, rs);

      Index ret;
      if ((ret = unblocked(A11)) >= 0) return k + ret;
      if (rs > 0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
      if (rs > 0)
        A22.template selfadjointView<Lower>().rankUpdate(A21,
                                                         typename NumTraits<RealScalar>::Literal(-1));  // bottleneck
    }
    return -1;
  }

  template <typename MatrixType, typename VectorType>
  static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) {
    return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
  }
};

template <typename Scalar>
struct llt_inplace<Scalar, Upper> {
  using RealScalar = typename NumTraits<Scalar>::Real;

  template <typename MatrixType>
  static EIGEN_STRONG_INLINE Index unblocked(MatrixType& mat) {
    Transpose<MatrixType> matt(mat);
    return llt_inplace<Scalar, Lower>::unblocked(matt);
  }
  template <typename MatrixType>
  static EIGEN_STRONG_INLINE Index blocked(MatrixType& mat) {
    Transpose<MatrixType> matt(mat);
    return llt_inplace<Scalar, Lower>::blocked(matt);
  }
  template <typename MatrixType, typename VectorType>
  static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) {
    Transpose<MatrixType> matt(mat);
    return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
  }
};

template <typename MatrixType>
struct LLT_Traits<MatrixType, Lower> {
  using MatrixL = const TriangularView<const MatrixType, Lower>;
  using MatrixU = const TriangularView<const typename MatrixType::AdjointReturnType, Upper>;
  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
  static bool inplace_decomposition(MatrixType& m) {
    return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m) == -1;
  }
};

template <typename MatrixType>
struct LLT_Traits<MatrixType, Upper> {
  using MatrixL = const TriangularView<const typename MatrixType::AdjointReturnType, Lower>;
  using MatrixU = const TriangularView<const MatrixType, Upper>;
  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
  static bool inplace_decomposition(MatrixType& m) {
    return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m) == -1;
  }
};

}  // end namespace internal

/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
 *
 * \returns a reference to *this
 *
 * Example: \include TutorialLinAlgComputeTwice.cpp
 * Output: \verbinclude TutorialLinAlgComputeTwice.out
 */
template <typename MatrixType, int UpLo_>
template <typename InputType>
LLT<MatrixType, UpLo_>& LLT<MatrixType, UpLo_>::compute(const EigenBase<InputType>& a) {
  eigen_assert(a.rows() == a.cols());
  const Index size = a.rows();
  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 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);
  m_info = ok ? Success : NumericalIssue;

  return *this;
}

/** Performs a rank one update (or downdate) of the current decomposition.
 * If A = LL^* before the rank one update,
 * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector
 * of same dimension.
 */
template <typename MatrixType_, int UpLo_>
template <typename VectorType>
LLT<MatrixType_, UpLo_>& LLT<MatrixType_, UpLo_>::rankUpdate(const VectorType& v, const RealScalar& sigma) {
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
  eigen_assert(v.size() == m_matrix.cols());
  eigen_assert(m_isInitialized);
  if (internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix, v, sigma) >= 0)
    m_info = NumericalIssue;
  else
    m_info = Success;

  return *this;
}

// A = L L^*, with L real and positive on the diagonal, so det(A) = prod(L_ii)^2 > 0.

template <typename MatrixType_, int UpLo_>
typename LLT<MatrixType_, UpLo_>::Scalar LLT<MatrixType_, UpLo_>::determinant() const {
  return Scalar(absDeterminant());
}

template <typename MatrixType_, int UpLo_>
typename LLT<MatrixType_, UpLo_>::RealScalar LLT<MatrixType_, UpLo_>::absDeterminant() const {
  eigen_assert(m_isInitialized && "LLT is not initialized.");
  eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
  return numext::abs2(m_matrix.diagonal().real().prod());
}

template <typename MatrixType_, int UpLo_>
typename LLT<MatrixType_, UpLo_>::RealScalar LLT<MatrixType_, UpLo_>::logAbsDeterminant() const {
  eigen_assert(m_isInitialized && "LLT is not initialized.");
  eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
  return RealScalar(2) * m_matrix.diagonal().real().array().log().sum();
}

template <typename MatrixType_, int UpLo_>
typename LLT<MatrixType_, UpLo_>::Scalar LLT<MatrixType_, UpLo_>::signDeterminant() const {
  eigen_assert(m_isInitialized && "LLT is not initialized.");
  eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
  return Scalar(1);
}

#ifndef EIGEN_PARSED_BY_DOXYGEN
template <typename MatrixType_, int UpLo_>
template <typename RhsType, typename DstType>
void LLT<MatrixType_, UpLo_>::_solve_impl(const RhsType& rhs, DstType& dst) const {
  _solve_impl_transposed<true>(rhs, dst);
}

template <typename MatrixType_, int UpLo_>
template <bool Conjugate, typename RhsType, typename DstType>
void LLT<MatrixType_, UpLo_>::_solve_impl_transposed(const RhsType& rhs, DstType& dst) const {
  dst = rhs;

  matrixL().template conjugateIf<!Conjugate>().solveInPlace(dst);
  matrixU().template conjugateIf<!Conjugate>().solveInPlace(dst);
}
#endif

/** \internal use x = llt_object.solve(x);
 *
 * This is the \em in-place version of solve().
 *
 * \param bAndX represents both the right-hand side matrix b and result x.
 *
 * This version avoids a copy when the right hand side matrix b is not needed anymore.
 *
 * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here.
 * This function will const_cast it, so constness isn't honored here.
 *
 * \sa LLT::solve(), MatrixBase::llt()
 */
template <typename MatrixType, int UpLo_>
template <typename Derived>
void LLT<MatrixType, UpLo_>::solveInPlace(const MatrixBase<Derived>& bAndX) const {
  eigen_assert(m_isInitialized && "LLT is not initialized.");
  eigen_assert(m_matrix.rows() == bAndX.rows());
  matrixL().solveInPlace(bAndX);
  matrixU().solveInPlace(bAndX);
}

namespace internal {

/***** Implementation of inverse() *****************************************************/
template <typename DstXprType, typename MatrixType, int UpLo_>
struct Assignment<DstXprType, Inverse<LLT<MatrixType, UpLo_> >,
                  internal::assign_op<typename DstXprType::Scalar, typename LLT<MatrixType, UpLo_>::Scalar>,
                  Dense2Dense> {
  using LltType = LLT<MatrixType, UpLo_>;
  using SrcXprType = Inverse<LltType>;
  static constexpr unsigned int kMirrorMode = UpLo_ == Lower ? StrictlyUpper : StrictlyLower;
  static void run(DstXprType& dst, const SrcXprType& src,
                  const internal::assign_op<typename DstXprType::Scalar, typename LltType::Scalar>&) {
    const LltType& llt = src.nestedExpression();
    eigen_assert(llt.info() == Success && "LLT::inverse(): the factorization failed.");
    const Index size = llt.rows();
    if ((dst.rows() != size) || (dst.cols() != size)) dst.resize(size, size);

    // Complex is never thresholded; see EIGEN_LLT_INVERSE_POTRI_THRESHOLD.
    constexpr Index kPotriThreshold =
        NumTraits<typename LltType::Scalar>::IsComplex ? Index(0) : Index(EIGEN_LLT_INVERSE_POTRI_THRESHOLD);
    // An in-place LLT<Ref<...>> may be asked to overwrite its own factor. The POTRI sequence does exactly
    // that, so it takes the alias at every size; the solve fallback would read a factor that setIdentity()
    // had already destroyed. extract_data() is null for a destination whose inner stride is not known to
    // be 1 at compile time, which leaves the alias unknown rather than excluded, so that goes the same way.
    const typename DstXprType::Scalar* dst_data = extract_data(dst);
    const bool overwrites_factor = dst_data == nullptr || dst_data == llt.matrixLLT().data();
    if (overwrites_factor || size >= kPotriThreshold) {
      // A = L L^*, hence A^-1 = L^-* L^-1: invert the factor (xTRTRI), then square it (xLAUUM).
      dst.template triangularView<UpLo_>() = llt.matrixLLT().template triangularView<UpLo_>();
      dst.template triangularView<UpLo_>().inverseInPlace();
      internal::triangular_adjoint_square_in_place<UpLo_>(dst);
    } else {
      dst.setIdentity();
      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.
    dst.template triangularView<kMirrorMode>() = dst.adjoint();
  }
};

}  // end namespace internal

/** \returns the matrix represented by the decomposition,
 * i.e., it returns the product: L L^*.
 * This function is provided for debug purpose. */
template <typename MatrixType, int UpLo_>
MatrixType LLT<MatrixType, UpLo_>::reconstructedMatrix() const {
  eigen_assert(m_isInitialized && "LLT is not initialized.");
  return matrixL() * matrixL().adjoint().toDenseMatrix();
}

/** \cholesky_module
 * \returns the LLT decomposition of \c *this
 * \sa SelfAdjointView::llt()
 */
template <typename Derived>
inline LLT<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::llt() const {
  return LLT<PlainObject>(derived());
}

/** \cholesky_module
 * \returns the LLT decomposition of \c *this
 * \sa SelfAdjointView::llt()
 */
template <typename MatrixType, unsigned int UpLo>
inline LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo> SelfAdjointView<MatrixType, UpLo>::llt()
    const {
  return LLT<PlainObject, UpLo>(m_matrix);
}

}  // end namespace Eigen

#endif  // EIGEN_LLT_H