Eigen  3.3.4
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 namespace Eigen {
14 
15 namespace internal {
16 template<typename _MatrixType> struct traits<FullPivLU<_MatrixType> >
17  : traits<_MatrixType>
18 {
19  typedef MatrixXpr XprKind;
20  typedef SolverStorage StorageKind;
21  enum { Flags = 0 };
22 };
23 
24 } // end namespace internal
25 
59 template<typename _MatrixType> class FullPivLU
60  : public SolverBase<FullPivLU<_MatrixType> >
61 {
62  public:
63  typedef _MatrixType MatrixType;
64  typedef SolverBase<FullPivLU> Base;
65 
66  EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivLU)
67  // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int
68  enum {
69  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
71  };
72  typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType;
73  typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType;
74  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
75  typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
76  typedef typename MatrixType::PlainObject PlainObject;
77 
84  FullPivLU();
85 
92  FullPivLU(Index rows, Index cols);
93 
99  template<typename InputType>
100  explicit FullPivLU(const EigenBase<InputType>& matrix);
101 
108  template<typename InputType>
109  explicit FullPivLU(EigenBase<InputType>& matrix);
110 
118  template<typename InputType>
120  m_lu = matrix.derived();
121  computeInPlace();
122  return *this;
123  }
124 
131  inline const MatrixType& matrixLU() const
132  {
133  eigen_assert(m_isInitialized && "LU is not initialized.");
134  return m_lu;
135  }
136 
144  inline Index nonzeroPivots() const
145  {
146  eigen_assert(m_isInitialized && "LU is not initialized.");
147  return m_nonzero_pivots;
148  }
149 
153  RealScalar maxPivot() const { return m_maxpivot; }
154 
159  EIGEN_DEVICE_FUNC inline const PermutationPType& permutationP() const
160  {
161  eigen_assert(m_isInitialized && "LU is not initialized.");
162  return m_p;
163  }
164 
169  inline const PermutationQType& permutationQ() const
170  {
171  eigen_assert(m_isInitialized && "LU is not initialized.");
172  return m_q;
173  }
174 
189  inline const internal::kernel_retval<FullPivLU> kernel() const
190  {
191  eigen_assert(m_isInitialized && "LU is not initialized.");
192  return internal::kernel_retval<FullPivLU>(*this);
193  }
194 
214  inline const internal::image_retval<FullPivLU>
215  image(const MatrixType& originalMatrix) const
216  {
217  eigen_assert(m_isInitialized && "LU is not initialized.");
218  return internal::image_retval<FullPivLU>(*this, originalMatrix);
219  }
220 
240  // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion.
241  template<typename Rhs>
242  inline const Solve<FullPivLU, Rhs>
243  solve(const MatrixBase<Rhs>& b) const
244  {
245  eigen_assert(m_isInitialized && "LU is not initialized.");
246  return Solve<FullPivLU, Rhs>(*this, b.derived());
247  }
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() * 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 inline Index rows() const { return m_lu.rows(); }
410  EIGEN_DEVICE_FUNC inline Index cols() const { return m_lu.cols(); }
411 
412  #ifndef EIGEN_PARSED_BY_DOXYGEN
413  template<typename RhsType, typename DstType>
414  EIGEN_DEVICE_FUNC
415  void _solve_impl(const RhsType &rhs, DstType &dst) const;
416 
417  template<bool Conjugate, typename RhsType, typename DstType>
418  EIGEN_DEVICE_FUNC
419  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
420  #endif
421 
422  protected:
423 
424  static void check_template_parameters()
425  {
426  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
427  }
428 
429  void computeInPlace();
430 
431  MatrixType m_lu;
432  PermutationPType m_p;
433  PermutationQType m_q;
434  IntColVectorType m_rowsTranspositions;
435  IntRowVectorType m_colsTranspositions;
436  Index m_nonzero_pivots;
437  RealScalar m_l1_norm;
438  RealScalar m_maxpivot, m_prescribedThreshold;
439  signed char m_det_pq;
440  bool m_isInitialized, m_usePrescribedThreshold;
441 };
442 
443 template<typename MatrixType>
445  : m_isInitialized(false), m_usePrescribedThreshold(false)
446 {
447 }
448 
449 template<typename MatrixType>
451  : m_lu(rows, cols),
452  m_p(rows),
453  m_q(cols),
454  m_rowsTranspositions(rows),
455  m_colsTranspositions(cols),
456  m_isInitialized(false),
457  m_usePrescribedThreshold(false)
458 {
459 }
460 
461 template<typename MatrixType>
462 template<typename InputType>
464  : m_lu(matrix.rows(), matrix.cols()),
465  m_p(matrix.rows()),
466  m_q(matrix.cols()),
467  m_rowsTranspositions(matrix.rows()),
468  m_colsTranspositions(matrix.cols()),
469  m_isInitialized(false),
470  m_usePrescribedThreshold(false)
471 {
472  compute(matrix.derived());
473 }
474 
475 template<typename MatrixType>
476 template<typename InputType>
478  : m_lu(matrix.derived()),
479  m_p(matrix.rows()),
480  m_q(matrix.cols()),
481  m_rowsTranspositions(matrix.rows()),
482  m_colsTranspositions(matrix.cols()),
483  m_isInitialized(false),
484  m_usePrescribedThreshold(false)
485 {
486  computeInPlace();
487 }
488 
489 template<typename MatrixType>
491 {
492  check_template_parameters();
493 
494  // the permutations are stored as int indices, so just to be sure:
495  eigen_assert(m_lu.rows()<=NumTraits<int>::highest() && m_lu.cols()<=NumTraits<int>::highest());
496 
497  m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
498 
499  const Index size = m_lu.diagonalSize();
500  const Index rows = m_lu.rows();
501  const Index cols = m_lu.cols();
502 
503  // will store the transpositions, before we accumulate them at the end.
504  // can't accumulate on-the-fly because that will be done in reverse order for the rows.
505  m_rowsTranspositions.resize(m_lu.rows());
506  m_colsTranspositions.resize(m_lu.cols());
507  Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
508 
509  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
510  m_maxpivot = RealScalar(0);
511 
512  for(Index k = 0; k < size; ++k)
513  {
514  // First, we need to find the pivot.
515 
516  // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
517  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
518  typedef internal::scalar_score_coeff_op<Scalar> Scoring;
519  typedef typename Scoring::result_type Score;
520  Score biggest_in_corner;
521  biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
522  .unaryExpr(Scoring())
523  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
524  row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
525  col_of_biggest_in_corner += k; // need to add k to them.
526 
527  if(biggest_in_corner==Score(0))
528  {
529  // before exiting, make sure to initialize the still uninitialized transpositions
530  // in a sane state without destroying what we already have.
531  m_nonzero_pivots = k;
532  for(Index i = k; i < size; ++i)
533  {
534  m_rowsTranspositions.coeffRef(i) = i;
535  m_colsTranspositions.coeffRef(i) = i;
536  }
537  break;
538  }
539 
540  RealScalar abs_pivot = internal::abs_knowing_score<Scalar>()(m_lu(row_of_biggest_in_corner, col_of_biggest_in_corner), biggest_in_corner);
541  if(abs_pivot > m_maxpivot) m_maxpivot = abs_pivot;
542 
543  // Now that we've found the pivot, we need to apply the row/col swaps to
544  // bring it to the location (k,k).
545 
546  m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
547  m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
548  if(k != row_of_biggest_in_corner) {
549  m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
550  ++number_of_transpositions;
551  }
552  if(k != col_of_biggest_in_corner) {
553  m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
554  ++number_of_transpositions;
555  }
556 
557  // Now that the pivot is at the right location, we update the remaining
558  // bottom-right corner by Gaussian elimination.
559 
560  if(k<rows-1)
561  m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
562  if(k<size-1)
563  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);
564  }
565 
566  // the main loop is over, we still have to accumulate the transpositions to find the
567  // permutations P and Q
568 
569  m_p.setIdentity(rows);
570  for(Index k = size-1; k >= 0; --k)
571  m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
572 
573  m_q.setIdentity(cols);
574  for(Index k = 0; k < size; ++k)
575  m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
576 
577  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
578 
579  m_isInitialized = true;
580 }
581 
582 template<typename MatrixType>
583 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
584 {
585  eigen_assert(m_isInitialized && "LU is not initialized.");
586  eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
587  return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
588 }
589 
593 template<typename MatrixType>
595 {
596  eigen_assert(m_isInitialized && "LU is not initialized.");
597  const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
598  // LU
599  MatrixType res(m_lu.rows(),m_lu.cols());
600  // FIXME the .toDenseMatrix() should not be needed...
601  res = m_lu.leftCols(smalldim)
602  .template triangularView<UnitLower>().toDenseMatrix()
603  * m_lu.topRows(smalldim)
604  .template triangularView<Upper>().toDenseMatrix();
605 
606  // P^{-1}(LU)
607  res = m_p.inverse() * res;
608 
609  // (P^{-1}LU)Q^{-1}
610  res = res * m_q.inverse();
611 
612  return res;
613 }
614 
615 /********* Implementation of kernel() **************************************************/
616 
617 namespace internal {
618 template<typename _MatrixType>
619 struct kernel_retval<FullPivLU<_MatrixType> >
620  : kernel_retval_base<FullPivLU<_MatrixType> >
621 {
622  EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
623 
624  enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
625  MatrixType::MaxColsAtCompileTime,
626  MatrixType::MaxRowsAtCompileTime)
627  };
628 
629  template<typename Dest> void evalTo(Dest& dst) const
630  {
631  using std::abs;
632  const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
633  if(dimker == 0)
634  {
635  // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
636  // avoid crashing/asserting as that depends on floating point calculations. Let's
637  // just return a single column vector filled with zeros.
638  dst.setZero();
639  return;
640  }
641 
642  /* Let us use the following lemma:
643  *
644  * Lemma: If the matrix A has the LU decomposition PAQ = LU,
645  * then Ker A = Q(Ker U).
646  *
647  * Proof: trivial: just keep in mind that P, Q, L are invertible.
648  */
649 
650  /* Thus, all we need to do is to compute Ker U, and then apply Q.
651  *
652  * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
653  * Thus, the diagonal of U ends with exactly
654  * dimKer zero's. Let us use that to construct dimKer linearly
655  * independent vectors in Ker U.
656  */
657 
659  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
660  Index p = 0;
661  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
662  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
663  pivots.coeffRef(p++) = i;
664  eigen_internal_assert(p == rank());
665 
666  // we construct a temporaty trapezoid matrix m, by taking the U matrix and
667  // permuting the rows and cols to bring the nonnegligible pivots to the top of
668  // the main diagonal. We need that to be able to apply our triangular solvers.
669  // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
670  Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
671  MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
672  m(dec().matrixLU().block(0, 0, rank(), cols));
673  for(Index i = 0; i < rank(); ++i)
674  {
675  if(i) m.row(i).head(i).setZero();
676  m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
677  }
678  m.block(0, 0, rank(), rank());
679  m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
680  for(Index i = 0; i < rank(); ++i)
681  m.col(i).swap(m.col(pivots.coeff(i)));
682 
683  // ok, we have our trapezoid matrix, we can apply the triangular solver.
684  // notice that the math behind this suggests that we should apply this to the
685  // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
686  m.topLeftCorner(rank(), rank())
687  .template triangularView<Upper>().solveInPlace(
688  m.topRightCorner(rank(), dimker)
689  );
690 
691  // now we must undo the column permutation that we had applied!
692  for(Index i = rank()-1; i >= 0; --i)
693  m.col(i).swap(m.col(pivots.coeff(i)));
694 
695  // see the negative sign in the next line, that's what we were talking about above.
696  for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
697  for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
698  for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
699  }
700 };
701 
702 /***** Implementation of image() *****************************************************/
703 
704 template<typename _MatrixType>
705 struct image_retval<FullPivLU<_MatrixType> >
706  : image_retval_base<FullPivLU<_MatrixType> >
707 {
708  EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
709 
710  enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
711  MatrixType::MaxColsAtCompileTime,
712  MatrixType::MaxRowsAtCompileTime)
713  };
714 
715  template<typename Dest> void evalTo(Dest& dst) const
716  {
717  using std::abs;
718  if(rank() == 0)
719  {
720  // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
721  // avoid crashing/asserting as that depends on floating point calculations. Let's
722  // just return a single column vector filled with zeros.
723  dst.setZero();
724  return;
725  }
726 
728  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
729  Index p = 0;
730  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
731  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
732  pivots.coeffRef(p++) = i;
733  eigen_internal_assert(p == rank());
734 
735  for(Index i = 0; i < rank(); ++i)
736  dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
737  }
738 };
739 
740 /***** Implementation of solve() *****************************************************/
741 
742 } // end namespace internal
743 
744 #ifndef EIGEN_PARSED_BY_DOXYGEN
745 template<typename _MatrixType>
746 template<typename RhsType, typename DstType>
747 void FullPivLU<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
748 {
749  /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
750  * So we proceed as follows:
751  * Step 1: compute c = P * rhs.
752  * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
753  * Step 3: replace c by the solution x to Ux = c. May or may not exist.
754  * Step 4: result = Q * c;
755  */
756 
757  const Index rows = this->rows(),
758  cols = this->cols(),
759  nonzero_pivots = this->rank();
760  eigen_assert(rhs.rows() == rows);
761  const Index smalldim = (std::min)(rows, cols);
762 
763  if(nonzero_pivots == 0)
764  {
765  dst.setZero();
766  return;
767  }
768 
769  typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
770 
771  // Step 1
772  c = permutationP() * rhs;
773 
774  // Step 2
775  m_lu.topLeftCorner(smalldim,smalldim)
776  .template triangularView<UnitLower>()
777  .solveInPlace(c.topRows(smalldim));
778  if(rows>cols)
779  c.bottomRows(rows-cols) -= m_lu.bottomRows(rows-cols) * c.topRows(cols);
780 
781  // Step 3
782  m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
783  .template triangularView<Upper>()
784  .solveInPlace(c.topRows(nonzero_pivots));
785 
786  // Step 4
787  for(Index i = 0; i < nonzero_pivots; ++i)
788  dst.row(permutationQ().indices().coeff(i)) = c.row(i);
789  for(Index i = nonzero_pivots; i < m_lu.cols(); ++i)
790  dst.row(permutationQ().indices().coeff(i)).setZero();
791 }
792 
793 template<typename _MatrixType>
794 template<bool Conjugate, typename RhsType, typename DstType>
795 void FullPivLU<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
796 {
797  /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1},
798  * and since permutations are real and unitary, we can write this
799  * as A^T = Q U^T L^T P,
800  * So we proceed as follows:
801  * Step 1: compute c = Q^T rhs.
802  * Step 2: replace c by the solution x to U^T x = c. May or may not exist.
803  * Step 3: replace c by the solution x to L^T x = c.
804  * Step 4: result = P^T c.
805  * If Conjugate is true, replace "^T" by "^*" above.
806  */
807 
808  const Index rows = this->rows(), cols = this->cols(),
809  nonzero_pivots = this->rank();
810  eigen_assert(rhs.rows() == cols);
811  const Index smalldim = (std::min)(rows, cols);
812 
813  if(nonzero_pivots == 0)
814  {
815  dst.setZero();
816  return;
817  }
818 
819  typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
820 
821  // Step 1
822  c = permutationQ().inverse() * rhs;
823 
824  if (Conjugate) {
825  // Step 2
826  m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
827  .template triangularView<Upper>()
828  .adjoint()
829  .solveInPlace(c.topRows(nonzero_pivots));
830  // Step 3
831  m_lu.topLeftCorner(smalldim, smalldim)
832  .template triangularView<UnitLower>()
833  .adjoint()
834  .solveInPlace(c.topRows(smalldim));
835  } else {
836  // Step 2
837  m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
838  .template triangularView<Upper>()
839  .transpose()
840  .solveInPlace(c.topRows(nonzero_pivots));
841  // Step 3
842  m_lu.topLeftCorner(smalldim, smalldim)
843  .template triangularView<UnitLower>()
844  .transpose()
845  .solveInPlace(c.topRows(smalldim));
846  }
847 
848  // Step 4
849  PermutationPType invp = permutationP().inverse().eval();
850  for(Index i = 0; i < smalldim; ++i)
851  dst.row(invp.indices().coeff(i)) = c.row(i);
852  for(Index i = smalldim; i < rows; ++i)
853  dst.row(invp.indices().coeff(i)).setZero();
854 }
855 
856 #endif
857 
858 namespace internal {
859 
860 
861 /***** Implementation of inverse() *****************************************************/
862 template<typename DstXprType, typename MatrixType>
863 struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivLU<MatrixType>::Scalar>, Dense2Dense>
864 {
865  typedef FullPivLU<MatrixType> LuType;
866  typedef Inverse<LuType> SrcXprType;
867  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename MatrixType::Scalar> &)
868  {
869  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
870  }
871 };
872 } // end namespace internal
873 
874 /******* MatrixBase methods *****************************************************************/
875 
882 template<typename Derived>
885 {
886  return FullPivLU<PlainObject>(eval());
887 }
888 
889 } // end namespace Eigen
890 
891 #endif // EIGEN_LU_H
const internal::image_retval< FullPivLU > image(const MatrixType &originalMatrix) const
Definition: FullPivLU.h:215
FullPivLU & setThreshold(const RealScalar &threshold)
Definition: FullPivLU.h:292
InverseReturnType inverse() const
Definition: PermutationMatrix.h:196
FullPivLU & compute(const EigenBase< InputType > &matrix)
Definition: FullPivLU.h:119
const PermutationQType & permutationQ() const
Definition: FullPivLU.h:169
Derived & applyTranspositionOnTheRight(Index i, Index j)
Definition: PermutationMatrix.h:185
RealScalar rcond() const
Definition: FullPivLU.h:252
Namespace containing all symbols from the Eigen library.
Definition: Core:303
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Derived & derived()
Definition: EigenBase.h:45
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:38
const FullPivLU< PlainObject > fullPivLu() const
Definition: FullPivLU.h:884
Index nonzeroPivots() const
Definition: FullPivLU.h:144
Definition: EigenBase.h:29
Index rank() const
Definition: FullPivLU.h:332
Expression of the inverse of another expression.
Definition: Inverse.h:43
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:183
RealScalar threshold() const
Definition: FullPivLU.h:317
internal::traits< MatrixType >::Scalar determinant() const
Definition: FullPivLU.h:583
Index size() const
Definition: EigenBase.h:66
const MatrixType & matrixLU() const
Definition: FullPivLU.h:131
const IndicesType & indices() const
Definition: PermutationMatrix.h:388
ConstTransposeReturnType transpose() const
Definition: SolverBase.h:90
FullPivLU & setThreshold(Default_t)
Definition: FullPivLU.h:307
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
const Inverse< FullPivLU > inverse() const
Definition: FullPivLU.h:400
bool isInjective() const
Definition: FullPivLU.h:362
bool isSurjective() const
Definition: FullPivLU.h:375
void setIdentity()
Definition: PermutationMatrix.h:142
MatrixType reconstructedMatrix() const
Definition: FullPivLU.h:594
bool isInvertible() const
Definition: FullPivLU.h:387
const PermutationPType & permutationP() const
Definition: FullPivLU.h:159
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
const internal::kernel_retval< FullPivLU > kernel() const
Definition: FullPivLU.h:189
const Solve< FullPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: FullPivLU.h:243
const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:160
LU decomposition of a matrix with complete pivoting, and related features.
Definition: ForwardDeclarations.h:249
FullPivLU()
Default Constructor.
Definition: FullPivLU.h:444
const int Dynamic
Definition: Constants.h:21
Pseudo expression representing a solving operation.
Definition: Solve.h:62
Index dimensionOfKernel() const
Definition: FullPivLU.h:349
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
AdjointReturnType adjoint() const
Definition: SolverBase.h:109
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
RealScalar maxPivot() const
Definition: FullPivLU.h:153