Eigen::internal::ldlt_inplace< Lower > Struct Reference

#include <LDLT.h>

Static Public Member Functions
template<typename MatrixType , typename TranspositionType , typename Workspace >
static bool	unblocked (MatrixType &mat, TranspositionType &transpositions, Workspace &temp, SignMatrix &sign)
template<typename MatrixType , typename WDerived >
static bool	updateInPlace (MatrixType &mat, MatrixBase< WDerived > &w, const typename MatrixType::RealScalar &sigma=1)
template<typename MatrixType , typename TranspositionType , typename Workspace , typename WType >
static bool	update (MatrixType &mat, const TranspositionType &transpositions, Workspace &tmp, const WType &w, const typename MatrixType::RealScalar &sigma=1)

Member Function Documentation

unblocked()

template<typename MatrixType , typename TranspositionType , typename Workspace >

static bool Eigen::internal::ldlt_inplace< Lower >::unblocked ( MatrixType &  mat, TranspositionType &  transpositions, Workspace &  temp, SignMatrix &  sign  )

inlinestatic

                                                                                                               {
    using std::abs;
    typedef typename MatrixType::Scalar Scalar;
    typedef typename MatrixType::RealScalar RealScalar;
    typedef typename TranspositionType::StorageIndex IndexType;
    eigen_assert(mat.rows() == mat.cols());
    const Index size = mat.rows();
    bool found_zero_pivot = false;
    bool ret = true;
 
    if (size <= 1) {
      transpositions.setIdentity();
      if (size == 0)
        sign = ZeroSign;
      else if (numext::real(mat.coeff(0, 0)) > static_cast<RealScalar>(0))
        sign = PositiveSemiDef;
      else if (numext::real(mat.coeff(0, 0)) < static_cast<RealScalar>(0))
        sign = NegativeSemiDef;
      else
        sign = ZeroSign;
      return true;
    }
 
    for (Index k = 0; k < size; ++k) {
      // Find largest diagonal element
      Index index_of_biggest_in_corner;
      mat.diagonal().tail(size - k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
      index_of_biggest_in_corner += k;
 
      transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner);
      if (k != index_of_biggest_in_corner) {
        // apply the transposition while taking care to consider only
        // the lower triangular part
        Index s = size - index_of_biggest_in_corner - 1;  // trailing size after the biggest element
        mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
        mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
        std::swap(mat.coeffRef(k, k), mat.coeffRef(index_of_biggest_in_corner, index_of_biggest_in_corner));
        for (Index i = k + 1; i < index_of_biggest_in_corner; ++i) {
          Scalar tmp = mat.coeffRef(i, k);
          mat.coeffRef(i, k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner, i));
          mat.coeffRef(index_of_biggest_in_corner, i) = numext::conj(tmp);
        }
        if (NumTraits<Scalar>::IsComplex)
          mat.coeffRef(index_of_biggest_in_corner, k) = numext::conj(mat.coeff(index_of_biggest_in_corner, k));
      }
 
      // partition the matrix:
      //       A00 |  -  |  -
      // lu  = A10 | A11 |  -
      //       A20 | A21 | A22
      Index rs = size - k - 1;
      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);
 
      if (k > 0) {
        temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint();
        mat.coeffRef(k, k) -= (A10 * temp.head(k)).value();
        if (rs > 0) A21.noalias() -= A20 * temp.head(k);
      }
 
      // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot
      // was smaller than the cutoff value. However, since LDLT is not rank-revealing
      // we should only make sure that we do not introduce INF or NaN values.
      // Remark that LAPACK also uses 0 as the cutoff value.
      RealScalar realAkk = numext::real(mat.coeffRef(k, k));
      bool pivot_is_valid = (abs(realAkk) > RealScalar(0));
 
      if (k == 0 && !pivot_is_valid) {
        // The entire diagonal is zero, there is nothing more to do
        // except filling the transpositions, and checking whether the matrix is zero.
        sign = ZeroSign;
        for (Index j = 0; j < size; ++j) {
          transpositions.coeffRef(j) = IndexType(j);
          ret = ret && (mat.col(j).tail(size - j - 1).array() == Scalar(0)).all();
        }
        return ret;
      }
 
      if ((rs > 0) && pivot_is_valid)
        A21 /= realAkk;
      else if (rs > 0)
        ret = ret && (A21.array() == Scalar(0)).all();
 
      if (found_zero_pivot && pivot_is_valid)
        ret = false;  // factorization failed
      else if (!pivot_is_valid)
        found_zero_pivot = true;
 
      if (sign == PositiveSemiDef) {
        if (realAkk < static_cast<RealScalar>(0)) sign = Indefinite;
      } else if (sign == NegativeSemiDef) {
        if (realAkk > static_cast<RealScalar>(0)) sign = Indefinite;
      } else if (sign == ZeroSign) {
        if (realAkk > static_cast<RealScalar>(0))
          sign = PositiveSemiDef;
        else if (realAkk < static_cast<RealScalar>(0))
          sign = NegativeSemiDef;
      }
    }
 
    return ret;
  }

References abs(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::coeff(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::coeffRef(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::cols(), conj(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::diagonal(), eigen_assert, i, j, k, Eigen::internal::NegativeSemiDef, Eigen::internal::PositiveSemiDef, ret, Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::rows(), s, SYCL::sign(), Eigen::EigenBase< Derived >::size(), swap(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::swap(), tmp, Eigen::value, and Eigen::internal::ZeroSign.

update()

template<typename MatrixType , typename TranspositionType , typename Workspace , typename WType >

static bool Eigen::internal::ldlt_inplace< Lower >::update ( MatrixType &  mat, const TranspositionType &  transpositions, Workspace &  tmp, const WType &  w, const typename MatrixType::RealScalar &  sigma = 1  )

inlinestatic

                                                                     {
    // Apply the permutation to the input w
    tmp = transpositions * w;
 
    return ldlt_inplace<Lower>::updateInPlace(mat, tmp, sigma);
  }

References calibrate::sigma, tmp, and w.

Referenced by smc.smc::recursiveBayesian().

updateInPlace()

template<typename MatrixType , typename WDerived >

static bool Eigen::internal::ldlt_inplace< Lower >::updateInPlace ( MatrixType &  mat, MatrixBase< WDerived > &  w, const typename MatrixType::RealScalar &  sigma = 1  )

inlinestatic

                                                                            {
    using numext::isfinite;
    typedef typename MatrixType::Scalar Scalar;
    typedef typename MatrixType::RealScalar RealScalar;
 
    const Index size = mat.rows();
    eigen_assert(mat.cols() == size && w.size() == size);
 
    RealScalar alpha = 1;
 
    // Apply the update
    for (Index j = 0; j < size; j++) {
      // Check for termination due to an original decomposition of low-rank
      if (!(isfinite)(alpha)) break;
 
      // Update the diagonal terms
      RealScalar dj = numext::real(mat.coeff(j, j));
      Scalar wj = w.coeff(j);
      RealScalar swj2 = sigma * numext::abs2(wj);
      RealScalar gamma = dj * alpha + swj2;
 
      mat.coeffRef(j, j) += swj2 / alpha;
      alpha += swj2 / dj;
 
      // Update the terms of L
      Index rs = size - j - 1;
      w.tail(rs) -= wj * mat.col(j).tail(rs);
      if (!numext::is_exactly_zero(gamma)) mat.col(j).tail(rs) += (sigma * numext::conj(wj) / gamma) * w.tail(rs);
    }
    return true;
  }

References Eigen::numext::abs2(), alpha, Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::coeff(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::coeffRef(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::cols(), conj(), eigen_assert, mathsFunc::gamma(), Eigen::numext::is_exactly_zero(), Eigen::numext::isfinite(), isfinite, j, Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::rows(), calibrate::sigma, Eigen::EigenBase< Derived >::size(), and w.

---

The documentation for this struct was generated from the following file:

• LDLT.h