aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/LU/FullPivLU.h45
-rw-r--r--Eigen/src/LU/PartialPivLU.h4
-rw-r--r--doc/snippets/class_FullPivLU.cpp6
-rw-r--r--test/lu.cpp17
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);
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-04 12:52:09 +0000