Eigen  3.3.3
FullPivLU.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_LU_H
00011 #define EIGEN_LU_H
00012 
00013 namespace Eigen {
00014 
00015 namespace internal {
00016 template<typename _MatrixType> struct traits<FullPivLU<_MatrixType> >
00017  : traits<_MatrixType>
00018 {
00019   typedef MatrixXpr XprKind;
00020   typedef SolverStorage StorageKind;
00021   enum { Flags = 0 };
00022 };
00023 
00024 } // end namespace internal
00025 
00059 template<typename _MatrixType> class FullPivLU
00060   : public SolverBase<FullPivLU<_MatrixType> >
00061 {
00062   public:
00063     typedef _MatrixType MatrixType;
00064     typedef SolverBase<FullPivLU> Base;
00065 
00066     EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivLU)
00067     // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int
00068     enum {
00069       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00070       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00071     };
00072     typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType;
00073     typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType;
00074     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
00075     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
00076     typedef typename MatrixType::PlainObject PlainObject;
00077 
00084     FullPivLU();
00085 
00092     FullPivLU(Index rows, Index cols);
00093 
00099     template<typename InputType>
00100     explicit FullPivLU(const EigenBase<InputType>& matrix);
00101 
00108     template<typename InputType>
00109     explicit FullPivLU(EigenBase<InputType>& matrix);
00110 
00118     template<typename InputType>
00119     FullPivLU& compute(const EigenBase<InputType>& matrix) {
00120       m_lu = matrix.derived();
00121       computeInPlace();
00122       return *this;
00123     }
00124 
00131     inline const MatrixType& matrixLU() const
00132     {
00133       eigen_assert(m_isInitialized && "LU is not initialized.");
00134       return m_lu;
00135     }
00136 
00144     inline Index nonzeroPivots() const
00145     {
00146       eigen_assert(m_isInitialized && "LU is not initialized.");
00147       return m_nonzero_pivots;
00148     }
00149 
00153     RealScalar maxPivot() const { return m_maxpivot; }
00154 
00159     EIGEN_DEVICE_FUNC inline const PermutationPType& permutationP() const
00160     {
00161       eigen_assert(m_isInitialized && "LU is not initialized.");
00162       return m_p;
00163     }
00164 
00169     inline const PermutationQType& permutationQ() const
00170     {
00171       eigen_assert(m_isInitialized && "LU is not initialized.");
00172       return m_q;
00173     }
00174 
00189     inline const internal::kernel_retval<FullPivLU> kernel() const
00190     {
00191       eigen_assert(m_isInitialized && "LU is not initialized.");
00192       return internal::kernel_retval<FullPivLU>(*this);
00193     }
00194 
00214     inline const internal::image_retval<FullPivLU>
00215       image(const MatrixType& originalMatrix) const
00216     {
00217       eigen_assert(m_isInitialized && "LU is not initialized.");
00218       return internal::image_retval<FullPivLU>(*this, originalMatrix);
00219     }
00220 
00240     // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion.
00241     template<typename Rhs>
00242     inline const Solve<FullPivLU, Rhs>
00243     solve(const MatrixBase<Rhs>& b) const
00244     {
00245       eigen_assert(m_isInitialized && "LU is not initialized.");
00246       return Solve<FullPivLU, Rhs>(*this, b.derived());
00247     }
00248 
00252     inline RealScalar rcond() const
00253     {
00254       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00255       return internal::rcond_estimate_helper(m_l1_norm, *this);
00256     }
00257 
00273     typename internal::traits<MatrixType>::Scalar determinant() const;
00274 
00292     FullPivLU& setThreshold(const RealScalar& threshold)
00293     {
00294       m_usePrescribedThreshold = true;
00295       m_prescribedThreshold = threshold;
00296       return *this;
00297     }
00298 
00307     FullPivLU& setThreshold(Default_t)
00308     {
00309       m_usePrescribedThreshold = false;
00310       return *this;
00311     }
00312 
00317     RealScalar threshold() const
00318     {
00319       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00320       return m_usePrescribedThreshold ? m_prescribedThreshold
00321       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00322       // and turns out to be identical to Higham's formula used already in LDLt.
00323                                       : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
00324     }
00325 
00332     inline Index rank() const
00333     {
00334       using std::abs;
00335       eigen_assert(m_isInitialized && "LU is not initialized.");
00336       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
00337       Index result = 0;
00338       for(Index i = 0; i < m_nonzero_pivots; ++i)
00339         result += (abs(m_lu.coeff(i,i)) > premultiplied_threshold);
00340       return result;
00341     }
00342 
00349     inline Index dimensionOfKernel() const
00350     {
00351       eigen_assert(m_isInitialized && "LU is not initialized.");
00352       return cols() - rank();
00353     }
00354 
00362     inline bool isInjective() const
00363     {
00364       eigen_assert(m_isInitialized && "LU is not initialized.");
00365       return rank() == cols();
00366     }
00367 
00375     inline bool isSurjective() const
00376     {
00377       eigen_assert(m_isInitialized && "LU is not initialized.");
00378       return rank() == rows();
00379     }
00380 
00387     inline bool isInvertible() const
00388     {
00389       eigen_assert(m_isInitialized && "LU is not initialized.");
00390       return isInjective() && (m_lu.rows() == m_lu.cols());
00391     }
00392 
00400     inline const Inverse<FullPivLU> inverse() const
00401     {
00402       eigen_assert(m_isInitialized && "LU is not initialized.");
00403       eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
00404       return Inverse<FullPivLU>(*this);
00405     }
00406 
00407     MatrixType reconstructedMatrix() const;
00408 
00409     EIGEN_DEVICE_FUNC inline Index rows() const { return m_lu.rows(); }
00410     EIGEN_DEVICE_FUNC inline Index cols() const { return m_lu.cols(); }
00411 
00412     #ifndef EIGEN_PARSED_BY_DOXYGEN
00413     template<typename RhsType, typename DstType>
00414     EIGEN_DEVICE_FUNC
00415     void _solve_impl(const RhsType &rhs, DstType &dst) const;
00416 
00417     template<bool Conjugate, typename RhsType, typename DstType>
00418     EIGEN_DEVICE_FUNC
00419     void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
00420     #endif
00421 
00422   protected:
00423 
00424     static void check_template_parameters()
00425     {
00426       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00427     }
00428 
00429     void computeInPlace();
00430 
00431     MatrixType m_lu;
00432     PermutationPType m_p;
00433     PermutationQType m_q;
00434     IntColVectorType m_rowsTranspositions;
00435     IntRowVectorType m_colsTranspositions;
00436     Index m_nonzero_pivots;
00437     RealScalar m_l1_norm;
00438     RealScalar m_maxpivot, m_prescribedThreshold;
00439     signed char m_det_pq;
00440     bool m_isInitialized, m_usePrescribedThreshold;
00441 };
00442 
00443 template<typename MatrixType>
00444 FullPivLU<MatrixType>::FullPivLU()
00445   : m_isInitialized(false), m_usePrescribedThreshold(false)
00446 {
00447 }
00448 
00449 template<typename MatrixType>
00450 FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
00451   : m_lu(rows, cols),
00452     m_p(rows),
00453     m_q(cols),
00454     m_rowsTranspositions(rows),
00455     m_colsTranspositions(cols),
00456     m_isInitialized(false),
00457     m_usePrescribedThreshold(false)
00458 {
00459 }
00460 
00461 template<typename MatrixType>
00462 template<typename InputType>
00463 FullPivLU<MatrixType>::FullPivLU(const EigenBase<InputType>& matrix)
00464   : m_lu(matrix.rows(), matrix.cols()),
00465     m_p(matrix.rows()),
00466     m_q(matrix.cols()),
00467     m_rowsTranspositions(matrix.rows()),
00468     m_colsTranspositions(matrix.cols()),
00469     m_isInitialized(false),
00470     m_usePrescribedThreshold(false)
00471 {
00472   compute(matrix.derived());
00473 }
00474 
00475 template<typename MatrixType>
00476 template<typename InputType>
00477 FullPivLU<MatrixType>::FullPivLU(EigenBase<InputType>& matrix)
00478   : m_lu(matrix.derived()),
00479     m_p(matrix.rows()),
00480     m_q(matrix.cols()),
00481     m_rowsTranspositions(matrix.rows()),
00482     m_colsTranspositions(matrix.cols()),
00483     m_isInitialized(false),
00484     m_usePrescribedThreshold(false)
00485 {
00486   computeInPlace();
00487 }
00488 
00489 template<typename MatrixType>
00490 void FullPivLU<MatrixType>::computeInPlace()
00491 {
00492   check_template_parameters();
00493 
00494   // the permutations are stored as int indices, so just to be sure:
00495   eigen_assert(m_lu.rows()<=NumTraits<int>::highest() && m_lu.cols()<=NumTraits<int>::highest());
00496 
00497   m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
00498 
00499   const Index size = m_lu.diagonalSize();
00500   const Index rows = m_lu.rows();
00501   const Index cols = m_lu.cols();
00502 
00503   // will store the transpositions, before we accumulate them at the end.
00504   // can't accumulate on-the-fly because that will be done in reverse order for the rows.
00505   m_rowsTranspositions.resize(m_lu.rows());
00506   m_colsTranspositions.resize(m_lu.cols());
00507   Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
00508 
00509   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00510   m_maxpivot = RealScalar(0);
00511 
00512   for(Index k = 0; k < size; ++k)
00513   {
00514     // First, we need to find the pivot.
00515 
00516     // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
00517     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00518     typedef internal::scalar_score_coeff_op<Scalar> Scoring;
00519     typedef typename Scoring::result_type Score;
00520     Score biggest_in_corner;
00521     biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
00522                         .unaryExpr(Scoring())
00523                         .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00524     row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
00525     col_of_biggest_in_corner += k; // need to add k to them.
00526 
00527     if(biggest_in_corner==Score(0))
00528     {
00529       // before exiting, make sure to initialize the still uninitialized transpositions
00530       // in a sane state without destroying what we already have.
00531       m_nonzero_pivots = k;
00532       for(Index i = k; i < size; ++i)
00533       {
00534         m_rowsTranspositions.coeffRef(i) = i;
00535         m_colsTranspositions.coeffRef(i) = i;
00536       }
00537       break;
00538     }
00539 
00540     RealScalar abs_pivot = internal::abs_knowing_score<Scalar>()(m_lu(row_of_biggest_in_corner, col_of_biggest_in_corner), biggest_in_corner);
00541     if(abs_pivot > m_maxpivot) m_maxpivot = abs_pivot;
00542 
00543     // Now that we've found the pivot, we need to apply the row/col swaps to
00544     // bring it to the location (k,k).
00545 
00546     m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
00547     m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
00548     if(k != row_of_biggest_in_corner) {
00549       m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
00550       ++number_of_transpositions;
00551     }
00552     if(k != col_of_biggest_in_corner) {
00553       m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
00554       ++number_of_transpositions;
00555     }
00556 
00557     // Now that the pivot is at the right location, we update the remaining
00558     // bottom-right corner by Gaussian elimination.
00559 
00560     if(k<rows-1)
00561       m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
00562     if(k<size-1)
00563       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);
00564   }
00565 
00566   // the main loop is over, we still have to accumulate the transpositions to find the
00567   // permutations P and Q
00568 
00569   m_p.setIdentity(rows);
00570   for(Index k = size-1; k >= 0; --k)
00571     m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
00572 
00573   m_q.setIdentity(cols);
00574   for(Index k = 0; k < size; ++k)
00575     m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00576 
00577   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00578 
00579   m_isInitialized = true;
00580 }
00581 
00582 template<typename MatrixType>
00583 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
00584 {
00585   eigen_assert(m_isInitialized && "LU is not initialized.");
00586   eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
00587   return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
00588 }
00589 
00593 template<typename MatrixType>
00594 MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
00595 {
00596   eigen_assert(m_isInitialized && "LU is not initialized.");
00597   const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
00598   // LU
00599   MatrixType res(m_lu.rows(),m_lu.cols());
00600   // FIXME the .toDenseMatrix() should not be needed...
00601   res = m_lu.leftCols(smalldim)
00602             .template triangularView<UnitLower>().toDenseMatrix()
00603       * m_lu.topRows(smalldim)
00604             .template triangularView<Upper>().toDenseMatrix();
00605 
00606   // P^{-1}(LU)
00607   res = m_p.inverse() * res;
00608 
00609   // (P^{-1}LU)Q^{-1}
00610   res = res * m_q.inverse();
00611 
00612   return res;
00613 }
00614 
00615 /********* Implementation of kernel() **************************************************/
00616 
00617 namespace internal {
00618 template<typename _MatrixType>
00619 struct kernel_retval<FullPivLU<_MatrixType> >
00620   : kernel_retval_base<FullPivLU<_MatrixType> >
00621 {
00622   EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
00623 
00624   enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00625             MatrixType::MaxColsAtCompileTime,
00626             MatrixType::MaxRowsAtCompileTime)
00627   };
00628 
00629   template<typename Dest> void evalTo(Dest& dst) const
00630   {
00631     using std::abs;
00632     const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
00633     if(dimker == 0)
00634     {
00635       // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
00636       // avoid crashing/asserting as that depends on floating point calculations. Let's
00637       // just return a single column vector filled with zeros.
00638       dst.setZero();
00639       return;
00640     }
00641 
00642     /* Let us use the following lemma:
00643       *
00644       * Lemma: If the matrix A has the LU decomposition PAQ = LU,
00645       * then Ker A = Q(Ker U).
00646       *
00647       * Proof: trivial: just keep in mind that P, Q, L are invertible.
00648       */
00649 
00650     /* Thus, all we need to do is to compute Ker U, and then apply Q.
00651       *
00652       * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
00653       * Thus, the diagonal of U ends with exactly
00654       * dimKer zero's. Let us use that to construct dimKer linearly
00655       * independent vectors in Ker U.
00656       */
00657 
00658     Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00659     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00660     Index p = 0;
00661     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00662       if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00663         pivots.coeffRef(p++) = i;
00664     eigen_internal_assert(p == rank());
00665 
00666     // we construct a temporaty trapezoid matrix m, by taking the U matrix and
00667     // permuting the rows and cols to bring the nonnegligible pivots to the top of
00668     // the main diagonal. We need that to be able to apply our triangular solvers.
00669     // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
00670     Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
00671            MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
00672       m(dec().matrixLU().block(0, 0, rank(), cols));
00673     for(Index i = 0; i < rank(); ++i)
00674     {
00675       if(i) m.row(i).head(i).setZero();
00676       m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
00677     }
00678     m.block(0, 0, rank(), rank());
00679     m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
00680     for(Index i = 0; i < rank(); ++i)
00681       m.col(i).swap(m.col(pivots.coeff(i)));
00682 
00683     // ok, we have our trapezoid matrix, we can apply the triangular solver.
00684     // notice that the math behind this suggests that we should apply this to the
00685     // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
00686     m.topLeftCorner(rank(), rank())
00687      .template triangularView<Upper>().solveInPlace(
00688         m.topRightCorner(rank(), dimker)
00689       );
00690 
00691     // now we must undo the column permutation that we had applied!
00692     for(Index i = rank()-1; i >= 0; --i)
00693       m.col(i).swap(m.col(pivots.coeff(i)));
00694 
00695     // see the negative sign in the next line, that's what we were talking about above.
00696     for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
00697     for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00698     for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
00699   }
00700 };
00701 
00702 /***** Implementation of image() *****************************************************/
00703 
00704 template<typename _MatrixType>
00705 struct image_retval<FullPivLU<_MatrixType> >
00706   : image_retval_base<FullPivLU<_MatrixType> >
00707 {
00708   EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
00709 
00710   enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00711             MatrixType::MaxColsAtCompileTime,
00712             MatrixType::MaxRowsAtCompileTime)
00713   };
00714 
00715   template<typename Dest> void evalTo(Dest& dst) const
00716   {
00717     using std::abs;
00718     if(rank() == 0)
00719     {
00720       // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
00721       // avoid crashing/asserting as that depends on floating point calculations. Let's
00722       // just return a single column vector filled with zeros.
00723       dst.setZero();
00724       return;
00725     }
00726 
00727     Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00728     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00729     Index p = 0;
00730     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00731       if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00732         pivots.coeffRef(p++) = i;
00733     eigen_internal_assert(p == rank());
00734 
00735     for(Index i = 0; i < rank(); ++i)
00736       dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
00737   }
00738 };
00739 
00740 /***** Implementation of solve() *****************************************************/
00741 
00742 } // end namespace internal
00743 
00744 #ifndef EIGEN_PARSED_BY_DOXYGEN
00745 template<typename _MatrixType>
00746 template<typename RhsType, typename DstType>
00747 void FullPivLU<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
00748 {
00749   /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
00750   * So we proceed as follows:
00751   * Step 1: compute c = P * rhs.
00752   * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
00753   * Step 3: replace c by the solution x to Ux = c. May or may not exist.
00754   * Step 4: result = Q * c;
00755   */
00756 
00757   const Index rows = this->rows(),
00758               cols = this->cols(),
00759               nonzero_pivots = this->rank();
00760   eigen_assert(rhs.rows() == rows);
00761   const Index smalldim = (std::min)(rows, cols);
00762 
00763   if(nonzero_pivots == 0)
00764   {
00765     dst.setZero();
00766     return;
00767   }
00768 
00769   typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
00770 
00771   // Step 1
00772   c = permutationP() * rhs;
00773 
00774   // Step 2
00775   m_lu.topLeftCorner(smalldim,smalldim)
00776       .template triangularView<UnitLower>()
00777       .solveInPlace(c.topRows(smalldim));
00778   if(rows>cols)
00779     c.bottomRows(rows-cols) -= m_lu.bottomRows(rows-cols) * c.topRows(cols);
00780 
00781   // Step 3
00782   m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
00783       .template triangularView<Upper>()
00784       .solveInPlace(c.topRows(nonzero_pivots));
00785 
00786   // Step 4
00787   for(Index i = 0; i < nonzero_pivots; ++i)
00788     dst.row(permutationQ().indices().coeff(i)) = c.row(i);
00789   for(Index i = nonzero_pivots; i < m_lu.cols(); ++i)
00790     dst.row(permutationQ().indices().coeff(i)).setZero();
00791 }
00792 
00793 template<typename _MatrixType>
00794 template<bool Conjugate, typename RhsType, typename DstType>
00795 void FullPivLU<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
00796 {
00797   /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1},
00798    * and since permutations are real and unitary, we can write this
00799    * as   A^T = Q U^T L^T P,
00800    * So we proceed as follows:
00801    * Step 1: compute c = Q^T rhs.
00802    * Step 2: replace c by the solution x to U^T x = c. May or may not exist.
00803    * Step 3: replace c by the solution x to L^T x = c.
00804    * Step 4: result = P^T c.
00805    * If Conjugate is true, replace "^T" by "^*" above.
00806    */
00807 
00808   const Index rows = this->rows(), cols = this->cols(),
00809     nonzero_pivots = this->rank();
00810    eigen_assert(rhs.rows() == cols);
00811   const Index smalldim = (std::min)(rows, cols);
00812 
00813   if(nonzero_pivots == 0)
00814   {
00815     dst.setZero();
00816     return;
00817   }
00818 
00819   typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
00820 
00821   // Step 1
00822   c = permutationQ().inverse() * rhs;
00823 
00824   if (Conjugate) {
00825     // Step 2
00826     m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
00827         .template triangularView<Upper>()
00828         .adjoint()
00829         .solveInPlace(c.topRows(nonzero_pivots));
00830     // Step 3
00831     m_lu.topLeftCorner(smalldim, smalldim)
00832         .template triangularView<UnitLower>()
00833         .adjoint()
00834         .solveInPlace(c.topRows(smalldim));
00835   } else {
00836     // Step 2
00837     m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
00838         .template triangularView<Upper>()
00839         .transpose()
00840         .solveInPlace(c.topRows(nonzero_pivots));
00841     // Step 3
00842     m_lu.topLeftCorner(smalldim, smalldim)
00843         .template triangularView<UnitLower>()
00844         .transpose()
00845         .solveInPlace(c.topRows(smalldim));
00846   }
00847 
00848   // Step 4
00849   PermutationPType invp = permutationP().inverse().eval();
00850   for(Index i = 0; i < smalldim; ++i)
00851     dst.row(invp.indices().coeff(i)) = c.row(i);
00852   for(Index i = smalldim; i < rows; ++i)
00853     dst.row(invp.indices().coeff(i)).setZero();
00854 }
00855 
00856 #endif
00857 
00858 namespace internal {
00859 
00860 
00861 /***** Implementation of inverse() *****************************************************/
00862 template<typename DstXprType, typename MatrixType>
00863 struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivLU<MatrixType>::Scalar>, Dense2Dense>
00864 {
00865   typedef FullPivLU<MatrixType> LuType;
00866   typedef Inverse<LuType> SrcXprType;
00867   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename MatrixType::Scalar> &)
00868   {
00869     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
00870   }
00871 };
00872 } // end namespace internal
00873 
00874 /******* MatrixBase methods *****************************************************************/
00875 
00882 template<typename Derived>
00883 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
00884 MatrixBase<Derived>::fullPivLu() const
00885 {
00886   return FullPivLU<PlainObject>(eval());
00887 }
00888 
00889 } // end namespace Eigen
00890 
00891 #endif // EIGEN_LU_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends