Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime > Struct Template Reference

#include <PartialPivLU.h>

Public Types
typedef Matrix< Scalar, ActualSizeAtCompileTime, ActualSizeAtCompileTime, StorageOrder >	MatrixType
typedef Ref< MatrixType >	MatrixTypeRef
typedef Ref< Matrix< Scalar, Dynamic, Dynamic, StorageOrder > >	BlockType
typedef MatrixType::RealScalar	RealScalar

Static Public Member Functions
static Index	unblocked_lu (MatrixTypeRef &lu, PivIndex *row_transpositions, PivIndex &nb_transpositions)
static Index	blocked_lu (Index rows, Index cols, Scalar *lu_data, Index luStride, PivIndex *row_transpositions, PivIndex &nb_transpositions, Index maxBlockSize=256)

Static Public Attributes
static constexpr int	UnBlockedBound = 16
static constexpr bool	UnBlockedAtCompileTime = SizeAtCompileTime != Dynamic && SizeAtCompileTime <= UnBlockedBound
static constexpr int	ActualSizeAtCompileTime = UnBlockedAtCompileTime ? SizeAtCompileTime : Dynamic
static constexpr int	RRows = SizeAtCompileTime == 2 ? 1 : Dynamic
static constexpr int	RCols = SizeAtCompileTime == 2 ? 1 : Dynamic

Detailed Description

template<typename Scalar, int StorageOrder, typename PivIndex, int SizeAtCompileTime = Dynamic>
struct Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >

This is the blocked version of fullpivlu_unblocked()

Member Typedef Documentation

BlockType

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

typedef Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::BlockType

MatrixType

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

typedef Matrix<Scalar, ActualSizeAtCompileTime, ActualSizeAtCompileTime, StorageOrder> Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::MatrixType

MatrixTypeRef

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

typedef Ref<MatrixType> Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::MatrixTypeRef

RealScalar

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

typedef MatrixType::RealScalar Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::RealScalar

Member Function Documentation

blocked_lu()

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

static Index Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::blocked_lu ( Index  rows, Index  cols, Scalar *  lu_data, Index  luStride, PivIndex *  row_transpositions, PivIndex &  nb_transpositions, Index  maxBlockSize = 256  )

inlinestatic

performs the LU decomposition in-place of the matrix represented by the variables rows, cols, lu_data, and lu_stride using a recursive, blocked algorithm.

In addition, this function returns the row transpositions in the vector row_transpositions which must have a size equal to the number of columns of the matrix lu, and an integer nb_transpositions which returns the actual number of transpositions.

Returns

The index of the first pivot which is exactly zero if any, or a negative number otherwise.

Note

This very low level interface using pointers, etc. is to: 1 - reduce the number of instantiations to the strict minimum 2 - avoid infinite recursion of the instantiations with Block<Block<Block<...> > >

                                                                                 {
    MatrixTypeRef lu = MatrixType::Map(lu_data, rows, cols, OuterStride<>(luStride));
 
    const Index size = (std::min)(rows, cols);
 
    // if the matrix is too small, no blocking:
    if (UnBlockedAtCompileTime || size <= UnBlockedBound) {
      return unblocked_lu(lu, row_transpositions, nb_transpositions);
    }
 
    // automatically adjust the number of subdivisions to the size
    // of the matrix so that there is enough sub blocks:
    Index blockSize;
    {
      blockSize = size / 8;
      blockSize = (blockSize / 16) * 16;
      blockSize = (std::min)((std::max)(blockSize, Index(8)), maxBlockSize);
    }
 
    nb_transpositions = 0;
    Index first_zero_pivot = -1;
    for (Index k = 0; k < size; k += blockSize) {
      Index bs = (std::min)(size - k, blockSize);  // actual size of the block
      Index trows = rows - k - bs;                 // trailing rows
      Index tsize = size - k - bs;                 // trailing size
 
      // partition the matrix:
      //                          A00 | A01 | A02
      // lu  = A_0 | A_1 | A_2 =  A10 | A11 | A12
      //                          A20 | A21 | A22
      BlockType A_0 = lu.block(0, 0, rows, k);
      BlockType A_2 = lu.block(0, k + bs, rows, tsize);
      BlockType A11 = lu.block(k, k, bs, bs);
      BlockType A12 = lu.block(k, k + bs, bs, tsize);
      BlockType A21 = lu.block(k + bs, k, trows, bs);
      BlockType A22 = lu.block(k + bs, k + bs, trows, tsize);
 
      PivIndex nb_transpositions_in_panel;
      // recursively call the blocked LU algorithm on [A11^T A21^T]^T
      // with a very small blocking size:
      Index ret = blocked_lu(trows + bs, bs, &lu.coeffRef(k, k), luStride, row_transpositions + k,
                             nb_transpositions_in_panel, 16);
      if (ret >= 0 && first_zero_pivot == -1) first_zero_pivot = k + ret;
 
      nb_transpositions += nb_transpositions_in_panel;
      // update permutations and apply them to A_0
      for (Index i = k; i < k + bs; ++i) {
        Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k));
        A_0.row(i).swap(A_0.row(piv));
      }
 
      if (trows) {
        // apply permutations to A_2
        for (Index i = k; i < k + bs; ++i) A_2.row(i).swap(A_2.row(row_transpositions[i]));
 
        // A12 = A11^-1 A12
        A11.template triangularView<UnitLower>().solveInPlace(A12);
 
        A22.noalias() -= A21 * A12;
      }
    }
    return first_zero_pivot;
  }

References cols, i, k, lu(), Eigen::PlainObjectBase< Matrix< Scalar_, Rows_, Cols_, Options_, MaxRows_, MaxCols_ > >::Map(), max, min, ret, rows, size, Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::unblocked_lu(), Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::UnBlockedAtCompileTime, and Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::UnBlockedBound.

unblocked_lu()

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

static Index Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::unblocked_lu ( MatrixTypeRef &  lu, PivIndex *  row_transpositions, PivIndex &  nb_transpositions  )

inlinestatic

performs the LU decomposition in-place of the matrix lu using an unblocked algorithm.

In addition, this function returns the row transpositions in the vector row_transpositions which must have a size equal to the number of columns of the matrix lu, and an integer nb_transpositions which returns the actual number of transpositions.

Returns

The index of the first pivot which is exactly zero if any, or a negative number otherwise.

                                                                                                          {
    typedef scalar_score_coeff_op<Scalar> Scoring;
    typedef typename Scoring::result_type Score;
    const Index rows = lu.rows();
    const Index cols = lu.cols();
    const Index size = (std::min)(rows, cols);
    // For small compile-time matrices it is worth processing the last row separately:
    //  speedup: +100% for 2x2, +10% for others.
    const Index endk = UnBlockedAtCompileTime ? size - 1 : size;
    nb_transpositions = 0;
    Index first_zero_pivot = -1;
    for (Index k = 0; k < endk; ++k) {
      int rrows = internal::convert_index<int>(rows - k - 1);
      int rcols = internal::convert_index<int>(cols - k - 1);
 
      Index row_of_biggest_in_col;
      Score biggest_in_corner = lu.col(k).tail(rows - k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
      row_of_biggest_in_col += k;
 
      row_transpositions[k] = PivIndex(row_of_biggest_in_col);
 
      if (!numext::is_exactly_zero(biggest_in_corner)) {
        if (k != row_of_biggest_in_col) {
          lu.row(k).swap(lu.row(row_of_biggest_in_col));
          ++nb_transpositions;
        }
 
        lu.col(k).tail(fix<RRows>(rrows)) /= lu.coeff(k, k);
      } else if (first_zero_pivot == -1) {
        // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
        // and continue the factorization such we still have A = PLU
        first_zero_pivot = k;
      }
 
      if (k < rows - 1)
        lu.bottomRightCorner(fix<RRows>(rrows), fix<RCols>(rcols)).noalias() -=
            lu.col(k).tail(fix<RRows>(rrows)) * lu.row(k).tail(fix<RCols>(rcols));
    }
 
    // special handling of the last entry
    if (UnBlockedAtCompileTime) {
      Index k = endk;
      row_transpositions[k] = PivIndex(k);
      if (numext::is_exactly_zero(Scoring()(lu(k, k))) && first_zero_pivot == -1) first_zero_pivot = k;
    }
 
    return first_zero_pivot;
  }

References cols, Eigen::numext::is_exactly_zero(), k, lu(), min, rows, size, and Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::UnBlockedAtCompileTime.

Referenced by Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::blocked_lu().

Member Data Documentation

ActualSizeAtCompileTime

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

constexpr int Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::ActualSizeAtCompileTime = UnBlockedAtCompileTime ? SizeAtCompileTime : Dynamic

staticconstexpr

RCols

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

constexpr int Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::RCols = SizeAtCompileTime == 2 ? 1 : Dynamic

staticconstexpr

RRows

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

constexpr int Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::RRows = SizeAtCompileTime == 2 ? 1 : Dynamic

staticconstexpr

UnBlockedAtCompileTime

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

constexpr bool Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::UnBlockedAtCompileTime = SizeAtCompileTime != Dynamic && SizeAtCompileTime <= UnBlockedBound

staticconstexpr

Referenced by Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::blocked_lu(), and Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::unblocked_lu().

UnBlockedBound

template<typename Scalar , int StorageOrder, typename PivIndex , int SizeAtCompileTime = Dynamic>

constexpr int Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::UnBlockedBound = 16

staticconstexpr

Referenced by Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex, SizeAtCompileTime >::blocked_lu().

---

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

• PartialPivLU.h