Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
FullPivLU.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_LU_H
11 #define EIGEN_LU_H
12 
13 #include "./InternalHeaderCheck.h"
14 
15 namespace Eigen {
16 
17 namespace internal {
18 template<typename MatrixType_> struct traits<FullPivLU<MatrixType_> >
19  : traits<MatrixType_>
20 {
21  typedef MatrixXpr XprKind;
22  typedef SolverStorage StorageKind;
23  typedef int StorageIndex;
24  enum { Flags = 0 };
25 };
26 
27 } // end namespace internal
28 
62 template<typename MatrixType_> class FullPivLU
63  : public SolverBase<FullPivLU<MatrixType_> >
64 {
65  public:
66  typedef MatrixType_ MatrixType;
68  friend class SolverBase<FullPivLU>;
69 
70  EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivLU)
71  enum {
72  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
73  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
74  };
75  typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType;
76  typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType;
79  typedef typename MatrixType::PlainObject PlainObject;
80 
87  FullPivLU();
88 
95  FullPivLU(Index rows, Index cols);
96 
102  template<typename InputType>
103  explicit FullPivLU(const EigenBase<InputType>& matrix);
104 
111  template<typename InputType>
112  explicit FullPivLU(EigenBase<InputType>& matrix);
113 
121  template<typename InputType>
123  m_lu = matrix.derived();
124  computeInPlace();
125  return *this;
126  }
127 
134  inline const MatrixType& matrixLU() const
135  {
136  eigen_assert(m_isInitialized && "LU is not initialized.");
137  return m_lu;
138  }
139 
147  inline Index nonzeroPivots() const
148  {
149  eigen_assert(m_isInitialized && "LU is not initialized.");
150  return m_nonzero_pivots;
151  }
152 
156  RealScalar maxPivot() const { return m_maxpivot; }
157 
162  EIGEN_DEVICE_FUNC inline const PermutationPType& permutationP() const
163  {
164  eigen_assert(m_isInitialized && "LU is not initialized.");
165  return m_p;
166  }
167 
172  inline const PermutationQType& permutationQ() const
173  {
174  eigen_assert(m_isInitialized && "LU is not initialized.");
175  return m_q;
176  }
177 
192  inline const internal::kernel_retval<FullPivLU> kernel() const
193  {
194  eigen_assert(m_isInitialized && "LU is not initialized.");
195  return internal::kernel_retval<FullPivLU>(*this);
196  }
197 
217  inline const internal::image_retval<FullPivLU>
218  image(const MatrixType& originalMatrix) const
219  {
220  eigen_assert(m_isInitialized && "LU is not initialized.");
221  return internal::image_retval<FullPivLU>(*this, originalMatrix);
222  }
223 
224  #ifdef EIGEN_PARSED_BY_DOXYGEN
244  template<typename Rhs>
245  inline const Solve<FullPivLU, Rhs>
246  solve(const MatrixBase<Rhs>& b) const;
247  #endif
248 
252  inline RealScalar rcond() const
253  {
254  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
255  return internal::rcond_estimate_helper(m_l1_norm, *this);
256  }
257 
273  typename internal::traits<MatrixType>::Scalar determinant() const;
274 
292  FullPivLU& setThreshold(const RealScalar& threshold)
293  {
294  m_usePrescribedThreshold = true;
295  m_prescribedThreshold = threshold;
296  return *this;
297  }
298 
308  {
309  m_usePrescribedThreshold = false;
310  return *this;
311  }
312 
317  RealScalar threshold() const
318  {
319  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
320  return m_usePrescribedThreshold ? m_prescribedThreshold
321  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
322  // and turns out to be identical to Higham's formula used already in LDLt.
323  : NumTraits<Scalar>::epsilon() * RealScalar(m_lu.diagonalSize());
324  }
325 
332  inline Index rank() const
333  {
334  using std::abs;
335  eigen_assert(m_isInitialized && "LU is not initialized.");
336  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
337  Index result = 0;
338  for(Index i = 0; i < m_nonzero_pivots; ++i)
339  result += (abs(m_lu.coeff(i,i)) > premultiplied_threshold);
340  return result;
341  }
342 
349  inline Index dimensionOfKernel() const
350  {
351  eigen_assert(m_isInitialized && "LU is not initialized.");
352  return cols() - rank();
353  }
354 
362  inline bool isInjective() const
363  {
364  eigen_assert(m_isInitialized && "LU is not initialized.");
365  return rank() == cols();
366  }
367 
375  inline bool isSurjective() const
376  {
377  eigen_assert(m_isInitialized && "LU is not initialized.");
378  return rank() == rows();
379  }
380 
387  inline bool isInvertible() const
388  {
389  eigen_assert(m_isInitialized && "LU is not initialized.");
390  return isInjective() && (m_lu.rows() == m_lu.cols());
391  }
392 
400  inline const Inverse<FullPivLU> inverse() const
401  {
402  eigen_assert(m_isInitialized && "LU is not initialized.");
403  eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
404  return Inverse<FullPivLU>(*this);
405  }
406 
407  MatrixType reconstructedMatrix() const;
408 
409  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
410  inline Index rows() const EIGEN_NOEXCEPT { return m_lu.rows(); }
411  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
412  inline Index cols() const EIGEN_NOEXCEPT { return m_lu.cols(); }
413 
414  #ifndef EIGEN_PARSED_BY_DOXYGEN
415  template<typename RhsType, typename DstType>
416  void _solve_impl(const RhsType &rhs, DstType &dst) const;
417 
418  template<bool Conjugate, typename RhsType, typename DstType>
419  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
420  #endif
421 
422  protected:
423 
424  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
425 
426  void computeInPlace();
427 
428  MatrixType m_lu;
429  PermutationPType m_p;
430  PermutationQType m_q;
431  IntColVectorType m_rowsTranspositions;
432  IntRowVectorType m_colsTranspositions;
433  Index m_nonzero_pivots;
434  RealScalar m_l1_norm;
435  RealScalar m_maxpivot, m_prescribedThreshold;
436  signed char m_det_pq;
437  bool m_isInitialized, m_usePrescribedThreshold;
438 };
439 
440 template<typename MatrixType>
442  : m_isInitialized(false), m_usePrescribedThreshold(false)
443 {
444 }
445 
446 template<typename MatrixType>
448  : m_lu(rows, cols),
449  m_p(rows),
450  m_q(cols),
451  m_rowsTranspositions(rows),
452  m_colsTranspositions(cols),
453  m_isInitialized(false),
454  m_usePrescribedThreshold(false)
455 {
456 }
457 
458 template<typename MatrixType>
459 template<typename InputType>
461  : m_lu(matrix.rows(), matrix.cols()),
462  m_p(matrix.rows()),
463  m_q(matrix.cols()),
464  m_rowsTranspositions(matrix.rows()),
465  m_colsTranspositions(matrix.cols()),
466  m_isInitialized(false),
467  m_usePrescribedThreshold(false)
468 {
469  compute(matrix.derived());
470 }
471 
472 template<typename MatrixType>
473 template<typename InputType>
475  : m_lu(matrix.derived()),
476  m_p(matrix.rows()),
477  m_q(matrix.cols()),
478  m_rowsTranspositions(matrix.rows()),
479  m_colsTranspositions(matrix.cols()),
480  m_isInitialized(false),
481  m_usePrescribedThreshold(false)
482 {
483  computeInPlace();
484 }
485 
486 template<typename MatrixType>
488 {
489  // the permutations are stored as int indices, so just to be sure:
490  eigen_assert(m_lu.rows()<=NumTraits<int>::highest() && m_lu.cols()<=NumTraits<int>::highest());
491 
492  m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
493 
494  const Index size = m_lu.diagonalSize();
495  const Index rows = m_lu.rows();
496  const Index cols = m_lu.cols();
497 
498  // will store the transpositions, before we accumulate them at the end.
499  // can't accumulate on-the-fly because that will be done in reverse order for the rows.
500  m_rowsTranspositions.resize(m_lu.rows());
501  m_colsTranspositions.resize(m_lu.cols());
502  Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
503 
504  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
505  m_maxpivot = RealScalar(0);
506 
507  for(Index k = 0; k < size; ++k)
508  {
509  // First, we need to find the pivot.
510 
511  // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
512  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
513  typedef internal::scalar_score_coeff_op<Scalar> Scoring;
514  typedef typename Scoring::result_type Score;
515  Score biggest_in_corner;
516  biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
517  .unaryExpr(Scoring())
518  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
519  row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
520  col_of_biggest_in_corner += k; // need to add k to them.
521 
522  if(numext::is_exactly_zero(biggest_in_corner))
523  {
524  // before exiting, make sure to initialize the still uninitialized transpositions
525  // in a sane state without destroying what we already have.
526  m_nonzero_pivots = k;
527  for(Index i = k; i < size; ++i)
528  {
529  m_rowsTranspositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
530  m_colsTranspositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
531  }
532  break;
533  }
534 
535  RealScalar abs_pivot = internal::abs_knowing_score<Scalar>()(m_lu(row_of_biggest_in_corner, col_of_biggest_in_corner), biggest_in_corner);
536  if(abs_pivot > m_maxpivot) m_maxpivot = abs_pivot;
537 
538  // Now that we've found the pivot, we need to apply the row/col swaps to
539  // bring it to the location (k,k).
540 
541  m_rowsTranspositions.coeffRef(k) = internal::convert_index<StorageIndex>(row_of_biggest_in_corner);
542  m_colsTranspositions.coeffRef(k) = internal::convert_index<StorageIndex>(col_of_biggest_in_corner);
543  if(k != row_of_biggest_in_corner) {
544  m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
545  ++number_of_transpositions;
546  }
547  if(k != col_of_biggest_in_corner) {
548  m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
549  ++number_of_transpositions;
550  }
551 
552  // Now that the pivot is at the right location, we update the remaining
553  // bottom-right corner by Gaussian elimination.
554 
555  if(k<rows-1)
556  m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
557  if(k<size-1)
558  m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).tail(rows-k-1) * m_lu.row(k).tail(cols-k-1);
559  }
560 
561  // the main loop is over, we still have to accumulate the transpositions to find the
562  // permutations P and Q
563 
564  m_p.setIdentity(rows);
565  for(Index k = size-1; k >= 0; --k)
566  m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
567 
568  m_q.setIdentity(cols);
569  for(Index k = 0; k < size; ++k)
570  m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
571 
572  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
573 
574  m_isInitialized = true;
575 }
576 
577 template<typename MatrixType>
578 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
579 {
580  eigen_assert(m_isInitialized && "LU is not initialized.");
581  eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
582  return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
583 }
584 
588 template<typename MatrixType>
590 {
591  eigen_assert(m_isInitialized && "LU is not initialized.");
592  const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
593  // LU
594  MatrixType res(m_lu.rows(),m_lu.cols());
595  // FIXME the .toDenseMatrix() should not be needed...
596  res = m_lu.leftCols(smalldim)
597  .template triangularView<UnitLower>().toDenseMatrix()
598  * m_lu.topRows(smalldim)
599  .template triangularView<Upper>().toDenseMatrix();
600 
601  // P^{-1}(LU)
602  res = m_p.inverse() * res;
603 
604  // (P^{-1}LU)Q^{-1}
605  res = res * m_q.inverse();
606 
607  return res;
608 }
609 
610 /********* Implementation of kernel() **************************************************/
611 
612 namespace internal {
613 template<typename MatrixType_>
614 struct kernel_retval<FullPivLU<MatrixType_> >
615  : kernel_retval_base<FullPivLU<MatrixType_> >
616 {
617  EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<MatrixType_>)
618 
619  enum { MaxSmallDimAtCompileTime = min_size_prefer_fixed(
620  MatrixType::MaxColsAtCompileTime,
621  MatrixType::MaxRowsAtCompileTime)
622  };
623 
624  template<typename Dest> void evalTo(Dest& dst) const
625  {
626  using std::abs;
627  const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
628  if(dimker == 0)
629  {
630  // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
631  // avoid crashing/asserting as that depends on floating point calculations. Let's
632  // just return a single column vector filled with zeros.
633  dst.setZero();
634  return;
635  }
636 
637  /* Let us use the following lemma:
638  *
639  * Lemma: If the matrix A has the LU decomposition PAQ = LU,
640  * then Ker A = Q(Ker U).
641  *
642  * Proof: trivial: just keep in mind that P, Q, L are invertible.
643  */
644 
645  /* Thus, all we need to do is to compute Ker U, and then apply Q.
646  *
647  * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
648  * Thus, the diagonal of U ends with exactly
649  * dimKer zero's. Let us use that to construct dimKer linearly
650  * independent vectors in Ker U.
651  */
652 
653  Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
654  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
655  Index p = 0;
656  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
657  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
658  pivots.coeffRef(p++) = i;
659  eigen_internal_assert(p == rank());
660 
661  // we construct a temporaty trapezoid matrix m, by taking the U matrix and
662  // permuting the rows and cols to bring the nonnegligible pivots to the top of
663  // the main diagonal. We need that to be able to apply our triangular solvers.
664  // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
665  Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
666  MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
667  m(dec().matrixLU().block(0, 0, rank(), cols));
668  for(Index i = 0; i < rank(); ++i)
669  {
670  if(i) m.row(i).head(i).setZero();
671  m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
672  }
673  m.block(0, 0, rank(), rank());
674  m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
675  for(Index i = 0; i < rank(); ++i)
676  m.col(i).swap(m.col(pivots.coeff(i)));
677 
678  // ok, we have our trapezoid matrix, we can apply the triangular solver.
679  // notice that the math behind this suggests that we should apply this to the
680  // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
681  m.topLeftCorner(rank(), rank())
682  .template triangularView<Upper>().solveInPlace(
683  m.topRightCorner(rank(), dimker)
684  );
685 
686  // now we must undo the column permutation that we had applied!
687  for(Index i = rank()-1; i >= 0; --i)
688  m.col(i).swap(m.col(pivots.coeff(i)));
689 
690  // see the negative sign in the next line, that's what we were talking about above.
691  for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
692  for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
693  for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
694  }
695 };
696 
697 /***** Implementation of image() *****************************************************/
698 
699 template<typename MatrixType_>
700 struct image_retval<FullPivLU<MatrixType_> >
701  : image_retval_base<FullPivLU<MatrixType_> >
702 {
703  EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<MatrixType_>)
704 
705  enum { MaxSmallDimAtCompileTime = min_size_prefer_fixed(
706  MatrixType::MaxColsAtCompileTime,
707  MatrixType::MaxRowsAtCompileTime)
708  };
709 
710  template<typename Dest> void evalTo(Dest& dst) const
711  {
712  using std::abs;
713  if(rank() == 0)
714  {
715  // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
716  // avoid crashing/asserting as that depends on floating point calculations. Let's
717  // just return a single column vector filled with zeros.
718  dst.setZero();
719  return;
720  }
721 
722  Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
723  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
724  Index p = 0;
725  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
726  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
727  pivots.coeffRef(p++) = i;
728  eigen_internal_assert(p == rank());
729 
730  for(Index i = 0; i < rank(); ++i)
731  dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
732  }
733 };
734 
735 /***** Implementation of solve() *****************************************************/
736 
737 } // end namespace internal
738 
739 #ifndef EIGEN_PARSED_BY_DOXYGEN
740 template<typename MatrixType_>
741 template<typename RhsType, typename DstType>
742 void FullPivLU<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
743 {
744  /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
745  * So we proceed as follows:
746  * Step 1: compute c = P * rhs.
747  * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
748  * Step 3: replace c by the solution x to Ux = c. May or may not exist.
749  * Step 4: result = Q * c;
750  */
751 
752  const Index rows = this->rows(),
753  cols = this->cols(),
754  nonzero_pivots = this->rank();
755  const Index smalldim = (std::min)(rows, cols);
756 
757  if(nonzero_pivots == 0)
758  {
759  dst.setZero();
760  return;
761  }
762 
763  typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
764 
765  // Step 1
766  c = permutationP() * rhs;
767 
768  // Step 2
769  m_lu.topLeftCorner(smalldim,smalldim)
770  .template triangularView<UnitLower>()
771  .solveInPlace(c.topRows(smalldim));
772  if(rows>cols)
773  c.bottomRows(rows-cols) -= m_lu.bottomRows(rows-cols) * c.topRows(cols);
774 
775  // Step 3
776  m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
777  .template triangularView<Upper>()
778  .solveInPlace(c.topRows(nonzero_pivots));
779 
780  // Step 4
781  for(Index i = 0; i < nonzero_pivots; ++i)
782  dst.row(permutationQ().indices().coeff(i)) = c.row(i);
783  for(Index i = nonzero_pivots; i < m_lu.cols(); ++i)
784  dst.row(permutationQ().indices().coeff(i)).setZero();
785 }
786 
787 template<typename MatrixType_>
788 template<bool Conjugate, typename RhsType, typename DstType>
789 void FullPivLU<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
790 {
791  /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1},
792  * and since permutations are real and unitary, we can write this
793  * as A^T = Q U^T L^T P,
794  * So we proceed as follows:
795  * Step 1: compute c = Q^T rhs.
796  * Step 2: replace c by the solution x to U^T x = c. May or may not exist.
797  * Step 3: replace c by the solution x to L^T x = c.
798  * Step 4: result = P^T c.
799  * If Conjugate is true, replace "^T" by "^*" above.
800  */
801 
802  const Index rows = this->rows(), cols = this->cols(),
803  nonzero_pivots = this->rank();
804  const Index smalldim = (std::min)(rows, cols);
805 
806  if(nonzero_pivots == 0)
807  {
808  dst.setZero();
809  return;
810  }
811 
812  typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
813 
814  // Step 1
815  c = permutationQ().inverse() * rhs;
816 
817  // Step 2
818  m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
819  .template triangularView<Upper>()
820  .transpose()
821  .template conjugateIf<Conjugate>()
822  .solveInPlace(c.topRows(nonzero_pivots));
823 
824  // Step 3
825  m_lu.topLeftCorner(smalldim, smalldim)
826  .template triangularView<UnitLower>()
827  .transpose()
828  .template conjugateIf<Conjugate>()
829  .solveInPlace(c.topRows(smalldim));
830 
831  // Step 4
832  PermutationPType invp = permutationP().inverse().eval();
833  for(Index i = 0; i < smalldim; ++i)
834  dst.row(invp.indices().coeff(i)) = c.row(i);
835  for(Index i = smalldim; i < rows; ++i)
836  dst.row(invp.indices().coeff(i)).setZero();
837 }
838 
839 #endif
840 
841 namespace internal {
842 
843 
844 /***** Implementation of inverse() *****************************************************/
845 template<typename DstXprType, typename MatrixType>
846 struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivLU<MatrixType>::Scalar>, Dense2Dense>
847 {
848  typedef FullPivLU<MatrixType> LuType;
849  typedef Inverse<LuType> SrcXprType;
850  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename MatrixType::Scalar> &)
851  {
852  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
853  }
854 };
855 } // end namespace internal
856 
857 /******* MatrixBase methods *****************************************************************/
858 
865 template<typename Derived>
866 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
868 {
869  return FullPivLU<PlainObject>(eval());
870 }
871 
872 } // end namespace Eigen
873 
874 #endif // EIGEN_LU_H
LU decomposition of a matrix with complete pivoting, and related features.
Definition: FullPivLU.h:64
RealScalar maxPivot() const
Definition: FullPivLU.h:156
bool isInvertible() const
Definition: FullPivLU.h:387
MatrixType reconstructedMatrix() const
Definition: FullPivLU.h:589
const PermutationQType & permutationQ() const
Definition: FullPivLU.h:172
const Inverse< FullPivLU > inverse() const
Definition: FullPivLU.h:400
Index dimensionOfKernel() const
Definition: FullPivLU.h:349
FullPivLU & compute(const EigenBase< InputType > &matrix)
Definition: FullPivLU.h:122
const MatrixType & matrixLU() const
Definition: FullPivLU.h:134
const PermutationPType & permutationP() const
Definition: FullPivLU.h:162
const Solve< FullPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
internal::traits< MatrixType >::Scalar determinant() const
Definition: FullPivLU.h:578
Index rank() const
Definition: FullPivLU.h:332
bool isInjective() const
Definition: FullPivLU.h:362
const internal::image_retval< FullPivLU > image(const MatrixType &originalMatrix) const
Definition: FullPivLU.h:218
Index nonzeroPivots() const
Definition: FullPivLU.h:147
FullPivLU & setThreshold(const RealScalar &threshold)
Definition: FullPivLU.h:292
const internal::kernel_retval< FullPivLU > kernel() const
Definition: FullPivLU.h:192
RealScalar threshold() const
Definition: FullPivLU.h:317
RealScalar rcond() const
Definition: FullPivLU.h:252
FullPivLU & setThreshold(Default_t)
Definition: FullPivLU.h:307
FullPivLU()
Default Constructor.
Definition: FullPivLU.h:441
bool isSurjective() const
Definition: FullPivLU.h:375
Expression of the inverse of another expression.
Definition: Inverse.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const FullPivLU< PlainObject > fullPivLu() const
Definition: FullPivLU.h:867
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SolverBase.h:107
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231