aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Core/MatrixBase.h2
-rw-r--r--Eigen/src/Core/Part.h1
-rw-r--r--Eigen/src/Core/Swap.h25
-rw-r--r--Eigen/src/Core/util/Constants.h5
-rw-r--r--Eigen/src/LU/LU.h108
5 files changed, 73 insertions, 68 deletions
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index f5ffbed56..0243faaed 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -535,7 +535,7 @@ template<typename Derived> class MatrixBase
/////////// LU module ///////////
- const LU<EvalType> lu(int pivoting) const;
+ const LU<EvalType> lu() const;
const EvalType inverse() const;
void computeInverse(EvalType *result) const;
Scalar determinant() const;
diff --git a/Eigen/src/Core/Part.h b/Eigen/src/Core/Part.h
index 871e27427..4d39c4c08 100644
--- a/Eigen/src/Core/Part.h
+++ b/Eigen/src/Core/Part.h
@@ -179,6 +179,7 @@ struct ei_part_assignment_impl
}
else
{
+ ei_assert(Mode == Upper || Mode == Lower || Mode == StrictlyUpper || Mode == StrictlyLower);
if((Mode == Upper && row <= col)
|| (Mode == Lower && row >= col)
|| (Mode == StrictlyUpper && row < col)
diff --git a/Eigen/src/Core/Swap.h b/Eigen/src/Core/Swap.h
index 0ee57017e..99ca5f72c 100644
--- a/Eigen/src/Core/Swap.h
+++ b/Eigen/src/Core/Swap.h
@@ -27,14 +27,9 @@
/** \class SwapWrapper
*
- * \brief Expression which must be nested by value
+ * \internal
*
- * \param ExpressionType the type of the object of which we are requiring nesting-by-value
- *
- * This class is the return type of MatrixBase::nestByValue()
- * and most of the time this is the only way it is used.
- *
- * \sa MatrixBase::nestByValue()
+ * \brief Internal helper class for swapping two expressions
*/
template<typename ExpressionType>
struct ei_traits<SwapWrapper<ExpressionType> >
@@ -116,12 +111,10 @@ template<typename ExpressionType> class SwapWrapper
/** swaps *this with the expression \a other.
*
- * \note \a other is only marked const because I couln't find another way
- * to get g++ (4.2 and 4.3) to accept that template parameter resolution.
- * The problem seems to be that when swapping expressions as in
- * m.row(i).swap(m.row(j)); the Row object returned by row(j) is a temporary
- * and g++ doesn't dare to pass it by non-constant reference.
- * It gets const_cast'd of course. TODO: get rid of const here.
+ * \note \a other is only marked for internal reasons, but of course
+ * it gets const-casted. One reason is that one will often call swap
+ * on temporary objects (hence non-const references are forbidden).
+ * Another reason is that lazyAssign takes a const argument anyway.
*/
template<typename Derived>
template<typename OtherDerived>
@@ -131,3 +124,9 @@ void MatrixBase<Derived>::swap(const MatrixBase<OtherDerived>& other)
}
#endif // EIGEN_SWAP_H
+
+
+
+
+
+
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index dd854e9ac..24c653e2e 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -209,11 +209,6 @@ enum {
HasDirectAccess = DirectAccessBit
};
-enum {
- PartialPivoting,
- CompletePivoting
-};
-
const int FullyCoherentAccessPattern = 0x1;
const int InnerCoherentAccessPattern = 0x2 | FullyCoherentAccessPattern;
const int OuterCoherentAccessPattern = 0x4 | InnerCoherentAccessPattern;
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index a707ce8d7..b891b2fbf 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -54,7 +54,7 @@ template<typename MatrixType> class LU
typedef Matrix<int, MatrixType::ColsAtCompileTime, 1> IntRowVectorType;
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
- LU(const MatrixType& matrix, int pivoting = CompletePivoting);
+ LU(const MatrixType& matrix);
inline const MatrixType& matrixLU() const
{
@@ -81,6 +81,10 @@ template<typename MatrixType> class LU
return m_q;
}
+ inline const Matrix<Scalar, MatrixType::RowsAtCompileTime, Dynamic,
+ MatrixType::MaxRowsAtCompileTime,
+ MatrixType::MaxColsAtCompileTime> kernel() const;
+
template<typename OtherDerived>
typename ProductReturnType<Transpose<MatrixType>, OtherDerived>::Type::Eval
solve(const MatrixBase<MatrixType> &b) const;
@@ -107,11 +111,21 @@ template<typename MatrixType> class LU
return m_lu.cols() - m_rank;
}
- inline bool isInvertible() const
+ inline bool isInjective() const
{
return m_rank == m_lu.cols();
}
+ inline bool isSurjective() const
+ {
+ return m_rank == m_lu.rows();
+ }
+
+ inline bool isInvertible() const
+ {
+ return isInjective() && isSurjective();
+ }
+
protected:
MatrixType m_lu;
IntColVectorType m_p;
@@ -119,61 +133,48 @@ template<typename MatrixType> class LU
int m_det_pq;
Scalar m_biggest_eigenvalue_of_u;
int m_rank;
- int m_pivoting;
};
template<typename MatrixType>
-LU<MatrixType>::LU(const MatrixType& matrix, int pivoting)
+LU<MatrixType>::LU(const MatrixType& matrix)
: m_lu(matrix),
m_p(matrix.rows()),
- m_q(matrix.cols()),
- m_pivoting(pivoting)
+ m_q(matrix.cols())
{
const int size = matrix.diagonal().size();
const int rows = matrix.rows();
const int cols = matrix.cols();
- ei_assert(pivoting == PartialPivoting || pivoting == CompletePivoting);
IntColVectorType rows_transpositions(matrix.rows());
IntRowVectorType cols_transpositions(matrix.cols());
int number_of_transpositions = 0;
+ RealScalar biggest;
for(int k = 0; k < size; k++)
{
- int row_of_biggest, col_of_biggest;
- Scalar biggest;
- if(m_pivoting == CompletePivoting)
- {
- biggest = m_lu.corner(Eigen::BottomRight, rows-k, cols-k)
- .cwise().abs()
- .maxCoeff(&row_of_biggest, &col_of_biggest);
- row_of_biggest += k;
- col_of_biggest += k;
- rows_transpositions.coeffRef(k) = row_of_biggest;
- cols_transpositions.coeffRef(k) = col_of_biggest;
- if(k != row_of_biggest) {
- m_lu.row(k).swap(m_lu.row(row_of_biggest));
- number_of_transpositions++;
- }
- if(k != col_of_biggest) {
- m_lu.col(k).swap(m_lu.col(col_of_biggest));
- number_of_transpositions++;
- }
+ int row_of_biggest_in_corner, col_of_biggest_in_corner;
+ RealScalar biggest_in_corner;
+
+ biggest_in_corner = m_lu.corner(Eigen::BottomRight, rows-k, cols-k)
+ .cwise().abs()
+ .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
+ row_of_biggest_in_corner += k;
+ col_of_biggest_in_corner += k;
+ rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
+ cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
+ if(k != row_of_biggest_in_corner) {
+ m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
+ number_of_transpositions++;
}
- else // partial pivoting
- {
- biggest = m_lu.col(k).end(rows-k)
- .cwise().abs()
- .maxCoeff(&row_of_biggest);
- row_of_biggest += k;
- rows_transpositions.coeffRef(k) = row_of_biggest;
- if(k != row_of_biggest) {
- m_lu.row(k).swap(m_lu.row(row_of_biggest));
- number_of_transpositions++;
- }
+ if(k != col_of_biggest_in_corner) {
+ m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
+ number_of_transpositions++;
}
+
+ if(k==0) biggest = biggest_in_corner;
const Scalar lu_k_k = m_lu.coeff(k,k);
- if(ei_isMuchSmallerThan(lu_k_k, biggest)) continue;
+ std::cout << lu_k_k << " " << biggest << std::endl;
+ if(ei_isMuchSmallerThan(lu_k_k, biggest)) { std::cout << "hello" << std::endl; continue; }
if(k<rows-1)
m_lu.col(k).end(rows-k-1) /= lu_k_k;
if(k<size-1)
@@ -185,18 +186,15 @@ LU<MatrixType>::LU(const MatrixType& matrix, int pivoting)
for(int k = size-1; k >= 0; k--)
std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
- if(pivoting == CompletePivoting)
- {
- for(int k = 0; k < matrix.cols(); k++) m_q.coeffRef(k) = k;
- for(int k = 0; k < size; k++)
- std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
- }
+ for(int k = 0; k < matrix.cols(); k++) m_q.coeffRef(k) = k;
+ for(int k = 0; k < size; k++)
+ std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
- int index_of_biggest;
- m_lu.diagonal().cwise().abs().maxCoeff(&index_of_biggest);
- m_biggest_eigenvalue_of_u = m_lu.diagonal().coeff(index_of_biggest);
+ int index_of_biggest_in_diagonal;
+ m_lu.diagonal().cwise().abs().maxCoeff(&index_of_biggest_in_diagonal);
+ m_biggest_eigenvalue_of_u = m_lu.diagonal().coeff(index_of_biggest_in_diagonal);
m_rank = 0;
for(int k = 0; k < size; k++)
@@ -212,15 +210,27 @@ typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
return res;
}
+#if 0
+template<typename MatrixType>
+inline const Matrix<Scalar, RowsAtCompileTime, Dynamic,
+ MaxRowsAtCompileTime, MaxColsAtCompileTime> kernel() const
+{
+ Matrix<Scalar, RowsAtCompileTime, Dynamic,
+ MaxRowsAtCompileTime, MaxColsAtCompileTime>
+ result(m_lu.rows(), dimensionOfKernel());
+
+}
+#endif
+
/** \return the LU decomposition of \c *this.
*
* \sa class LU
*/
template<typename Derived>
const LU<typename MatrixBase<Derived>::EvalType>
-MatrixBase<Derived>::lu(int pivoting = CompletePivoting) const
+MatrixBase<Derived>::lu() const
{
- return LU<typename MatrixBase<Derived>::EvalType>(eval(), pivoting);
+ return eval();
}
#endif // EIGEN_LU_H
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-09 01:09:24 +0000