GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

port the QR module to PermutationMatrix

Browse Source
This commit is contained in:
Benoit Jacob 2009-11-17 08:14:54 -05:00
parent 30b610a10f
commit 9f21e2aab7
4 changed files with 17 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Eigen/src/LU/FullPivLU.h
View File

@ -594,7 +594,8 @@ struct ei_image_retval<FullPivLU<_MatrixType> >
pivots.coeffRef(p++) = i; pivots.coeffRef(p++) = i;
ei_internal_assert(p == rank()); ei_internal_assert(p == rank());
dst.block(0,0,dst.rows(),rank()) = (originalMatrix()*dec().permutationQ()).block(0,0,dst.rows(),rank()); for(int i = 0; i < rank(); ++i)
dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
} }
}; };
@ -630,8 +631,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs>
typename Rhs::PlainMatrixType c(rhs().rows(), rhs().cols()); typename Rhs::PlainMatrixType c(rhs().rows(), rhs().cols());
// Step 1 // Step 1
for(int i = 0; i < rows; ++i) c = dec().permutationP() * rhs();
c.row(dec().permutationP().indices().coeff(i)) = rhs().row(i);
// Step 2 // Step 2
dec().matrixLU() dec().matrixLU()

15
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -57,10 +57,9 @@ template<typename _MatrixType> class ColPivHouseholderQR
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType; typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
typedef PermutationMatrix<ColsAtCompileTime> PermutationType;
typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType; typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType; typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType;
typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
@ -117,7 +116,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
ColPivHouseholderQR& compute(const MatrixType& matrix); ColPivHouseholderQR& compute(const MatrixType& matrix);
const IntRowVectorType& colsPermutation() const const PermutationType& colsPermutation() const
{ {
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_cols_permutation; return m_cols_permutation;
@ -230,7 +229,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
protected: protected:
MatrixType m_qr; MatrixType m_qr;
HCoeffsType m_hCoeffs; HCoeffsType m_hCoeffs;
IntRowVectorType m_cols_permutation; PermutationType m_cols_permutation;
bool m_isInitialized; bool m_isInitialized;
RealScalar m_precision; RealScalar m_precision;
int m_rank; int m_rank;
@ -271,7 +270,6 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
m_precision = epsilon<Scalar>() * size; m_precision = epsilon<Scalar>() * size;
IntRowVectorType cols_transpositions(matrix.cols()); IntRowVectorType cols_transpositions(matrix.cols());
m_cols_permutation.resize(matrix.cols());
int number_of_transpositions = 0; int number_of_transpositions = 0;
RealRowVectorType colSqNorms(cols); RealRowVectorType colSqNorms(cols);
@ -314,9 +312,9 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2(); colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2();
} }
for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k; m_cols_permutation.setIdentity(cols);
for(int k = 0; k < size; ++k) for(int k = 0; k < size; ++k)
std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k))); m_cols_permutation.applyTranspositionOnTheLeft(k, cols_transpositions.coeff(k));
m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true; m_isInitialized = true;
@ -368,8 +366,7 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
.template triangularView<UpperTriangular>() .template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, dec().rank(), c.cols())); .solveInPlace(c.corner(TopLeft, dec().rank(), c.cols()));
for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().coeff(i)) = c.row(i); dst = dec().colsPermutation() * c;
for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero();
} }
}; };

14
Eigen/src/QR/FullPivHouseholderQR.h
View File

@ -58,6 +58,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType; typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType; typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType;
typedef PermutationMatrix<ColsAtCompileTime> PermutationType;
typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType; typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
@ -112,7 +113,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
FullPivHouseholderQR& compute(const MatrixType& matrix); FullPivHouseholderQR& compute(const MatrixType& matrix);
const IntRowVectorType& colsPermutation() const const PermutationType& colsPermutation() const
{ {
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return m_cols_permutation; return m_cols_permutation;
@ -231,7 +232,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
MatrixType m_qr; MatrixType m_qr;
HCoeffsType m_hCoeffs; HCoeffsType m_hCoeffs;
IntColVectorType m_rows_transpositions; IntColVectorType m_rows_transpositions;
IntRowVectorType m_cols_permutation; PermutationType m_cols_permutation;
bool m_isInitialized; bool m_isInitialized;
RealScalar m_precision; RealScalar m_precision;
int m_rank; int m_rank;
@ -273,7 +274,6 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
m_rows_transpositions.resize(matrix.rows()); m_rows_transpositions.resize(matrix.rows());
IntRowVectorType cols_transpositions(matrix.cols()); IntRowVectorType cols_transpositions(matrix.cols());
m_cols_permutation.resize(matrix.cols());
int number_of_transpositions = 0; int number_of_transpositions = 0;
RealScalar biggest(0); RealScalar biggest(0);
@ -322,9 +322,9 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
.applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
} }
for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k; m_cols_permutation.setIdentity(cols);
for(int k = 0; k < size; ++k) for(int k = 0; k < size; ++k)
std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k))); m_cols_permutation.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k));
m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true; m_isInitialized = true;
@ -379,8 +379,8 @@ struct ei_solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
.template triangularView<UpperTriangular>() .template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, dec().rank(), c.cols())); .solveInPlace(c.corner(TopLeft, dec().rank(), c.cols()));
for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().coeff(i)) = c.row(i); for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero(); for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
} }
}; };

5
test/qr_fullpivoting.cpp
View File

@ -51,11 +51,8 @@ template<typename MatrixType> void qr()
// FIXME need better way to construct trapezoid // FIXME need better way to construct trapezoid
for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
MatrixType b = qr.matrixQ() * r; MatrixType c = qr.matrixQ() * r * qr.colsPermutation().inverse();
MatrixType c = MatrixType::Zero(rows,cols);
for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i);
VERIFY_IS_APPROX(m1, c); VERIFY_IS_APPROX(m1, c);
MatrixType m2 = MatrixType::Random(cols,cols2); MatrixType m2 = MatrixType::Random(cols,cols2);