aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/LU/FullPivLU.h6
-rw-r--r--Eigen/src/QR/ColPivHouseholderQR.h15
-rw-r--r--Eigen/src/QR/FullPivHouseholderQR.h14
-rw-r--r--test/qr_fullpivoting.cpp5
4 files changed, 17 insertions, 23 deletions
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 607fe261f..28dc0cd47 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -594,7 +594,8 @@ struct ei_image_retval<FullPivLU<_MatrixType> >
pivots.coeffRef(p++) = i;
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());
// Step 1
- for(int i = 0; i < rows; ++i)
- c.row(dec().permutationP().indices().coeff(i)) = rhs().row(i);
+ c = dec().permutationP() * rhs();
// Step 2
dec().matrixLU()
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index d3c6b4c7c..e764af529 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -57,10 +57,9 @@ template<typename _MatrixType> class ColPivHouseholderQR
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
+ typedef PermutationMatrix<ColsAtCompileTime> PermutationType;
typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType;
- typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
- typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType;
typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
@@ -117,7 +116,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
ColPivHouseholderQR& compute(const MatrixType& matrix);
- const IntRowVectorType& colsPermutation() const
+ const PermutationType& colsPermutation() const
{
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_cols_permutation;
@@ -230,7 +229,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
- IntRowVectorType m_cols_permutation;
+ PermutationType m_cols_permutation;
bool m_isInitialized;
RealScalar m_precision;
int m_rank;
@@ -271,7 +270,6 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
m_precision = epsilon<Scalar>() * size;
IntRowVectorType cols_transpositions(matrix.cols());
- m_cols_permutation.resize(matrix.cols());
int number_of_transpositions = 0;
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();
}
- 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)
- 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_isInitialized = true;
@@ -368,8 +366,7 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
.template triangularView<UpperTriangular>()
.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 = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero();
+ dst = dec().colsPermutation() * c;
}
};
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index 07df5ed6b..ac0108046 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -58,6 +58,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType;
+ typedef PermutationMatrix<ColsAtCompileTime> PermutationType;
typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
@@ -112,7 +113,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
FullPivHouseholderQR& compute(const MatrixType& matrix);
- const IntRowVectorType& colsPermutation() const
+ const PermutationType& colsPermutation() const
{
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return m_cols_permutation;
@@ -231,7 +232,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
MatrixType m_qr;
HCoeffsType m_hCoeffs;
IntColVectorType m_rows_transpositions;
- IntRowVectorType m_cols_permutation;
+ PermutationType m_cols_permutation;
bool m_isInitialized;
RealScalar m_precision;
int m_rank;
@@ -273,7 +274,6 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
m_rows_transpositions.resize(matrix.rows());
IntRowVectorType cols_transpositions(matrix.cols());
- m_cols_permutation.resize(matrix.cols());
int number_of_transpositions = 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));
}
- 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)
- 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_isInitialized = true;
@@ -379,8 +379,8 @@ struct ei_solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
.template triangularView<UpperTriangular>()
.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 = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero();
+ 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().indices().coeff(i)).setZero();
}
};
diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp
index 65d9a071f..c82ba4c7e 100644
--- a/test/qr_fullpivoting.cpp
+++ b/test/qr_fullpivoting.cpp
@@ -51,11 +51,8 @@ template<typename MatrixType> void qr()
// 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);
- 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);
MatrixType m2 = MatrixType::Random(cols,cols2);
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-08 23:37:47 +0000