diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-11-16 17:05:12 -0500 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-11-16 17:05:12 -0500 |
commit | b90744dc053d176f3ae00f67fbaac7f3b083b55e (patch) | |
tree | d2b7f049337f60f8385eaa6de0e3f2fe242db533 | |
parent | 76c614f9bd8e80a530a4aa685f1ae7506b64b8dd (diff) |
Port FullPivLU to PermutationMatrix
-rw-r--r-- | Eigen/src/LU/FullPivLU.h | 45 | ||||
-rw-r--r-- | Eigen/src/LU/PartialPivLU.h | 4 | ||||
-rw-r--r-- | doc/snippets/class_FullPivLU.cpp | 6 | ||||
-rw-r--r-- | test/lu.cpp | 17 |
4 files changed, 38 insertions, 34 deletions
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index d8975b0b6..62bdfba56 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -63,8 +63,8 @@ template<typename _MatrixType> class FullPivLU typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; - typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType; + typedef PermutationMatrix<MatrixType::ColsAtCompileTime> PermutationQType; + typedef PermutationMatrix<MatrixType::RowsAtCompileTime> PermutationPType; /** * \brief Default Constructor. @@ -120,25 +120,21 @@ template<typename _MatrixType> class FullPivLU */ RealScalar maxPivot() const { return m_maxpivot; } - /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, - * representing the P permutation i.e. the permutation of the rows. For its precise meaning, - * see the examples given in the documentation of class FullPivLU. + /** \returns the permutation matrix P * * \sa permutationQ() */ - inline const IntColVectorType& permutationP() const + inline const PermutationPType& permutationP() const { ei_assert(m_isInitialized && "LU is not initialized."); return m_p; } - /** \returns a vector of integers, whose size is the number of columns of the matrix being - * decomposed, representing the Q permutation i.e. the permutation of the columns. - * For its precise meaning, see the examples given in the documentation of class FullPivLU. + /** \returns the permutation matrix Q * * \sa permutationP() */ - inline const IntRowVectorType& permutationQ() const + inline const PermutationQType& permutationQ() const { ei_assert(m_isInitialized && "LU is not initialized."); return m_q; @@ -368,8 +364,8 @@ template<typename _MatrixType> class FullPivLU protected: MatrixType m_lu; - IntColVectorType m_p; - IntRowVectorType m_q; + PermutationPType m_p; + PermutationQType m_q; int m_det_pq, m_nonzero_pivots; RealScalar m_maxpivot, m_prescribedThreshold; bool m_isInitialized, m_usePrescribedThreshold; @@ -463,13 +459,13 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix) // the main loop is over, we still have to accumulate the transpositions to find the // permutations P and Q - for(int k = 0; k < matrix.rows(); ++k) m_p.coeffRef(k) = k; + for(int k = 0; k < matrix.rows(); ++k) m_p.indices().coeffRef(k) = k; for(int k = size-1; k >= 0; --k) - std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k))); + std::swap(m_p.indices().coeffRef(k), m_p.indices().coeffRef(rows_transpositions.coeff(k))); - for(int k = 0; k < matrix.cols(); ++k) m_q.coeffRef(k) = k; + for(int k = 0; k < matrix.cols(); ++k) m_q.indices().coeffRef(k) = k; for(int k = 0; k < size; ++k) - std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k))); + std::swap(m_q.indices().coeffRef(k), m_q.indices().coeffRef(cols_transpositions.coeff(k))); m_det_pq = (number_of_transpositions%2) ? -1 : 1; return *this; @@ -562,9 +558,9 @@ struct ei_kernel_retval<FullPivLU<_MatrixType> > m.col(i).swap(m.col(pivots.coeff(i))); // see the negative sign in the next line, that's what we were talking about above. - for(int i = 0; i < rank(); ++i) dst.row(dec().permutationQ().coeff(i)) = -m.row(i).end(dimker); - for(int i = rank(); i < cols; ++i) dst.row(dec().permutationQ().coeff(i)).setZero(); - for(int k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().coeff(rank()+k), k) = Scalar(1); + for(int i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).end(dimker); + for(int i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero(); + for(int k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1); } }; @@ -599,9 +595,8 @@ struct ei_image_retval<FullPivLU<_MatrixType> > if(ei_abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold) pivots.coeffRef(p++) = i; ei_internal_assert(p == rank()); - - for(int i = 0; i < rank(); ++i) - dst.col(i) = originalMatrix().col(dec().permutationQ().coeff(pivots.coeff(i))); + + dst.block(0,0,dst.rows(),rank()) = (originalMatrix()*dec().permutationQ()).block(0,0,dst.rows(),rank()); } }; @@ -638,7 +633,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs> // Step 1 for(int i = 0; i < rows; ++i) - c.row(dec().permutationP().coeff(i)) = rhs().row(i); + c.row(dec().permutationP().indices().coeff(i)) = rhs().row(i); // Step 2 dec().matrixLU() @@ -660,9 +655,9 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs> // Step 4 for(int i = 0; i < nonzero_pivots; ++i) - dst.row(dec().permutationQ().coeff(i)) = c.row(i); + dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i); for(int i = nonzero_pivots; i < dec().matrixLU().cols(); ++i) - dst.row(dec().permutationQ().coeff(i)).setZero(); + dst.row(dec().permutationQ().indices().coeff(i)).setZero(); } }; diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 847c895d7..4eb1162c1 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -103,9 +103,7 @@ template<typename _MatrixType> class PartialPivLU return m_lu; } - /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, - * representing the P permutation i.e. the permutation of the rows. For its precise meaning, - * see the examples given in the documentation of class FullPivLU. + /** \returns the permutation matrix P. */ inline const PermutationType& permutationP() const { diff --git a/doc/snippets/class_FullPivLU.cpp b/doc/snippets/class_FullPivLU.cpp index 40d76e8e6..855782d15 100644 --- a/doc/snippets/class_FullPivLU.cpp +++ b/doc/snippets/class_FullPivLU.cpp @@ -13,8 +13,4 @@ cout << "Here is the U part:" << endl; Matrix5x3 u = lu.matrixLU().part<UpperTriangular>(); cout << u << endl; cout << "Let us now reconstruct the original matrix m:" << endl; -Matrix5x3 x = l * u; -Matrix5x3 y; -for(int i = 0; i < 5; i++) for(int j = 0; j < 3; j++) - y(i, lu.permutationQ()[j]) = x(lu.permutationP()[i], j); -cout << y << endl; // should be equal to the original matrix m +cout << lu.permutationP().inverse() * l * u * lu.permutationQ().inverse() << endl; diff --git a/test/lu.cpp b/test/lu.cpp index 9dcebbeaa..c2237febf 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -24,9 +24,11 @@ #include "main.h" #include <Eigen/LU> +using namespace std; template<typename MatrixType> void lu_non_invertible() { + typedef typename MatrixType::Scalar Scalar; /* this test covers the following files: LU.h */ @@ -65,7 +67,20 @@ template<typename MatrixType> void lu_non_invertible() createRandomMatrixOfRank(rank, rows, cols, m1); FullPivLU<MatrixType> lu(m1); - std::cout << lu.kernel().rows() << " " << lu.kernel().cols() << std::endl; + // FIXME need better way to construct trapezoid matrices. extend triangularView to support rectangular. + DynamicMatrixType u(rows,cols); + for(int i = 0; i < rows; i++) + for(int j = 0; j < cols; j++) + u(i,j) = i>j ? Scalar(0) : lu.matrixLU()(i,j); + DynamicMatrixType l(rows,rows); + for(int i = 0; i < rows; i++) + for(int j = 0; j < rows; j++) + l(i,j) = (i<j || j>=cols) ? Scalar(0) + : i==j ? Scalar(1) + : lu.matrixLU()(i,j); + + VERIFY_IS_APPROX(lu.permutationP() * m1 * lu.permutationQ(), l*u); + KernelMatrixType m1kernel = lu.kernel(); ImageMatrixType m1image = lu.image(m1); |