Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h
index a00bd5d..7409fca 100644
--- a/Eigen/src/SparseQR/SparseQR.h
+++ b/Eigen/src/SparseQR/SparseQR.h
@@ -21,8 +21,12 @@
template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> >
{
typedef typename SparseQRType::MatrixType ReturnType;
- typedef typename ReturnType::Index Index;
+ typedef typename ReturnType::StorageIndex StorageIndex;
typedef typename ReturnType::StorageKind StorageKind;
+ enum {
+ RowsAtCompileTime = Dynamic,
+ ColsAtCompileTime = Dynamic
+ };
};
template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
{
@@ -48,7 +52,7 @@
* rank-revealing permutations. Use colsPermutation() to get it.
*
* Q is the orthogonal matrix represented as products of Householder reflectors.
- * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
+ * Use matrixQ() to get an expression and matrixQ().adjoint() to get the adjoint.
* You can then apply it to a vector.
*
* R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient.
@@ -58,24 +62,37 @@
* \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module
* OrderingMethods \endlink module for the list of built-in and external ordering methods.
*
+ * \implsparsesolverconcept
+ *
* \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()).
+ * \warning For complex matrices matrixQ().transpose() will actually return the adjoint matrix.
*
*/
template<typename _MatrixType, typename _OrderingType>
-class SparseQR
+class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
{
+ protected:
+ typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base;
+ using Base::m_isInitialized;
public:
+ using Base::_solve_impl;
typedef _MatrixType MatrixType;
typedef _OrderingType OrderingType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
- typedef SparseMatrix<Scalar,ColMajor,Index> QRMatrixType;
- typedef Matrix<Index, Dynamic, 1> IndexVector;
+ typedef typename MatrixType::StorageIndex StorageIndex;
+ typedef SparseMatrix<Scalar,ColMajor,StorageIndex> QRMatrixType;
+ typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
- typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
+ typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
+
+ enum {
+ ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+ MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+ };
+
public:
- SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
+ SparseQR () : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
{ }
/** Construct a QR factorization of the matrix \a mat.
@@ -84,7 +101,7 @@
*
* \sa compute()
*/
- SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
+ explicit SparseQR(const MatrixType& mat) : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
{
compute(mat);
}
@@ -112,6 +129,17 @@
inline Index cols() const { return m_pmat.cols();}
/** \returns a const reference to the \b sparse upper triangular matrix R of the QR factorization.
+ * \warning The entries of the returned matrix are not sorted. This means that using it in algorithms
+ * expecting sorted entries will fail. This include random coefficient accesses (SpaseMatrix::coeff()),
+ * and coefficient-wise operations. Matrix products and triangular solves are fine though.
+ *
+ * To sort the entries, you can assign it to a row-major matrix, and if a column-major matrix
+ * is required, you can copy it again:
+ * \code
+ * SparseMatrix<double> R = qr.matrixR(); // column-major, not sorted!
+ * SparseMatrix<double,RowMajor> Rr = qr.matrixR(); // row-major, sorted
+ * SparseMatrix<double> Rc = Rr; // column-major, sorted
+ * \endcode
*/
const QRMatrixType& matrixR() const { return m_R; }
@@ -119,7 +147,7 @@
*
* \sa setPivotThreshold()
*/
- Index rank() const
+ Index rank() const
{
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
return m_nonzeropivots;
@@ -162,23 +190,23 @@
/** \internal */
template<typename Rhs, typename Dest>
- bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
+ bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
{
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
Index rank = this->rank();
- // Compute Q^T * b;
+ // Compute Q^* * b;
typename Dest::PlainObject y, b;
- y = this->matrixQ().transpose() * B;
+ y = this->matrixQ().adjoint() * B;
b = y;
// Solve with the triangular matrix R
- y.resize((std::max)(cols(),Index(y.rows())),y.cols());
+ y.resize((std::max<Index>)(cols(),y.rows()),y.cols());
y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
y.bottomRows(y.rows()-rank).setZero();
-
+
// Apply the column permutation
if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols());
else dest = y.topRows(cols());
@@ -186,7 +214,6 @@
m_info = Success;
return true;
}
-
/** Sets the threshold that is used to determine linearly dependent columns during the factorization.
*
@@ -204,18 +231,18 @@
* \sa compute()
*/
template<typename Rhs>
- inline const internal::solve_retval<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const
+ inline const Solve<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const
{
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
- return internal::solve_retval<SparseQR, Rhs>(*this, B.derived());
+ return Solve<SparseQR, Rhs>(*this, B.derived());
}
template<typename Rhs>
- inline const internal::sparse_solve_retval<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
+ inline const Solve<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
{
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
- return internal::sparse_solve_retval<SparseQR, Rhs>(*this, B.derived());
+ return Solve<SparseQR, Rhs>(*this, B.derived());
}
/** \brief Reports whether previous computation was successful.
@@ -232,8 +259,9 @@
return m_info;
}
- protected:
- inline void sort_matrix_Q()
+
+ /** \internal */
+ inline void _sort_matrix_Q()
{
if(this->m_isQSorted) return;
// The matrix Q is sorted during the transposition
@@ -244,7 +272,6 @@
protected:
- bool m_isInitialized;
bool m_analysisIsok;
bool m_factorizationIsok;
mutable ComputationInfo m_info;
@@ -258,14 +285,13 @@
PermutationType m_outputPerm_c; // The final column permutation
RealScalar m_threshold; // Threshold to determine null Householder reflections
bool m_useDefaultThreshold; // Use default threshold
- Index m_nonzeropivots; // Number of non zero pivots found
+ Index m_nonzeropivots; // Number of non zero pivots found
IndexVector m_etree; // Column elimination tree
IndexVector m_firstRowElt; // First element in each row
bool m_isQSorted; // whether Q is sorted or not
bool m_isEtreeOk; // whether the elimination tree match the initial input matrix
template <typename, typename > friend struct SparseQR_QProduct;
- template <typename > friend struct SparseQRMatrixQReturnType;
};
@@ -294,7 +320,7 @@
if (!m_perm_c.size())
{
m_perm_c.resize(n);
- m_perm_c.indices().setLinSpaced(n, 0,n-1);
+ m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
}
// Compute the column elimination tree of the permuted matrix
@@ -323,12 +349,11 @@
void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
{
using std::abs;
- using std::max;
eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
- Index m = mat.rows();
- Index n = mat.cols();
- Index diagSize = (std::min)(m,n);
+ StorageIndex m = StorageIndex(mat.rows());
+ StorageIndex n = StorageIndex(mat.cols());
+ StorageIndex diagSize = (std::min)(m,n);
IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes
IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q
Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q
@@ -353,7 +378,7 @@
// otherwise directly use the input matrix
//
IndexVector originalOuterIndicesCpy;
- const Index *originalOuterIndices = mat.outerIndexPtr();
+ const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
if(MatrixType::IsRowMajor)
{
originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
@@ -375,7 +400,7 @@
if(m_useDefaultThreshold)
{
RealScalar max2Norm = 0.0;
- for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm());
+ for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
if(max2Norm==RealScalar(0))
max2Norm = RealScalar(1);
pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
@@ -384,11 +409,11 @@
// Initialize the numerical permutation
m_pivotperm.setIdentity(n);
- Index nonzeroCol = 0; // Record the number of valid pivots
+ StorageIndex nonzeroCol = 0; // Record the number of valid pivots
m_Q.startVec(0);
// Left looking rank-revealing QR factorization: compute a column of R and Q at a time
- for (Index col = 0; col < n; ++col)
+ for (StorageIndex col = 0; col < n; ++col)
{
mark.setConstant(-1);
m_R.startVec(col);
@@ -404,12 +429,12 @@
// thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found.
for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
{
- Index curIdx = nonzeroCol;
- if(itp) curIdx = itp.row();
+ StorageIndex curIdx = nonzeroCol;
+ if(itp) curIdx = StorageIndex(itp.row());
if(curIdx == nonzeroCol) found_diag = true;
// Get the nonzeros indexes of the current column of R
- Index st = m_firstRowElt(curIdx); // The traversal of the etree starts here
+ StorageIndex st = m_firstRowElt(curIdx); // The traversal of the etree starts here
if (st < 0 )
{
m_lastError = "Empty row found during numerical factorization";
@@ -466,7 +491,7 @@
{
for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
{
- Index iQ = itq.row();
+ StorageIndex iQ = StorageIndex(itq.row());
if (mark(iQ) != col)
{
Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q,
@@ -476,7 +501,7 @@
}
} // End update current column
- Scalar tau = 0;
+ Scalar tau = RealScalar(0);
RealScalar beta = 0;
if(nonzeroCol < diagSize)
@@ -572,44 +597,15 @@
m_info = Success;
}
-namespace internal {
-
-template<typename _MatrixType, typename OrderingType, typename Rhs>
-struct solve_retval<SparseQR<_MatrixType,OrderingType>, Rhs>
- : solve_retval_base<SparseQR<_MatrixType,OrderingType>, Rhs>
-{
- typedef SparseQR<_MatrixType,OrderingType> Dec;
- EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- dec()._solve(rhs(),dst);
- }
-};
-template<typename _MatrixType, typename OrderingType, typename Rhs>
-struct sparse_solve_retval<SparseQR<_MatrixType, OrderingType>, Rhs>
- : sparse_solve_retval_base<SparseQR<_MatrixType, OrderingType>, Rhs>
-{
- typedef SparseQR<_MatrixType, OrderingType> Dec;
- EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec, Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- this->defaultEvalTo(dst);
- }
-};
-} // end namespace internal
-
template <typename SparseQRType, typename Derived>
struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
{
typedef typename SparseQRType::QRMatrixType MatrixType;
typedef typename SparseQRType::Scalar Scalar;
- typedef typename SparseQRType::Index Index;
// Get the references
SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) :
m_qr(qr),m_other(other),m_transpose(transpose) {}
- inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); }
+ inline Index rows() const { return m_qr.matrixQ().rows(); }
inline Index cols() const { return m_other.cols(); }
// Assign to a vector
@@ -637,7 +633,10 @@
}
else
{
- eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
+ eigen_assert(m_qr.matrixQ().cols() == m_other.rows() && "Non conforming object sizes");
+
+ res.conservativeResize(rows(), cols());
+
// Compute res = Q * other column by column
for(Index j = 0; j < res.cols(); j++)
{
@@ -646,7 +645,7 @@
Scalar tau = Scalar(0);
tau = m_qr.m_Q.col(k).dot(res.col(j));
if(tau==Scalar(0)) continue;
- tau = tau * m_qr.m_hcoeffs(k);
+ tau = tau * numext::conj(m_qr.m_hcoeffs(k));
res.col(j) -= tau * m_qr.m_Q.col(k);
}
}
@@ -655,52 +654,44 @@
const SparseQRType& m_qr;
const Derived& m_other;
- bool m_transpose;
+ bool m_transpose; // TODO this actually means adjoint
};
template<typename SparseQRType>
struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
{
- typedef typename SparseQRType::Index Index;
typedef typename SparseQRType::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
- SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
+ enum {
+ RowsAtCompileTime = Dynamic,
+ ColsAtCompileTime = Dynamic
+ };
+ explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
template<typename Derived>
SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
{
return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false);
}
+ // To use for operations with the adjoint of Q
SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
{
return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
}
inline Index rows() const { return m_qr.rows(); }
- inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); }
- // To use for operations with the transpose of Q
+ inline Index cols() const { return m_qr.rows(); }
+ // To use for operations with the transpose of Q FIXME this is the same as adjoint at the moment
SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
{
return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
}
- template<typename Dest> void evalTo(MatrixBase<Dest>& dest) const
- {
- dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows());
- }
- template<typename Dest> void evalTo(SparseMatrixBase<Dest>& dest) const
- {
- Dest idMat(m_qr.rows(), m_qr.rows());
- idMat.setIdentity();
- // Sort the sparse householder reflectors if needed
- const_cast<SparseQRType *>(&m_qr)->sort_matrix_Q();
- dest.derived() = SparseQR_QProduct<SparseQRType, Dest>(m_qr, idMat, false);
- }
-
const SparseQRType& m_qr;
};
+// TODO this actually represents the adjoint of Q
template<typename SparseQRType>
struct SparseQRMatrixQTransposeReturnType
{
- SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
+ explicit SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
template<typename Derived>
SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other)
{
@@ -709,6 +700,46 @@
const SparseQRType& m_qr;
};
+namespace internal {
+
+template<typename SparseQRType>
+struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
+{
+ typedef typename SparseQRType::MatrixType MatrixType;
+ typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
+ typedef SparseShape Shape;
+};
+
+template< typename DstXprType, typename SparseQRType>
+struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Sparse>
+{
+ typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
+ typedef typename DstXprType::Scalar Scalar;
+ typedef typename DstXprType::StorageIndex StorageIndex;
+ static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/)
+ {
+ typename DstXprType::PlainObject idMat(src.rows(), src.cols());
+ idMat.setIdentity();
+ // Sort the sparse householder reflectors if needed
+ const_cast<SparseQRType *>(&src.m_qr)->_sort_matrix_Q();
+ dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat, false);
+ }
+};
+
+template< typename DstXprType, typename SparseQRType>
+struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Dense>
+{
+ typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
+ typedef typename DstXprType::Scalar Scalar;
+ typedef typename DstXprType::StorageIndex StorageIndex;
+ static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/)
+ {
+ dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
+ }
+};
+
+} // end namespace internal
+
} // end namespace Eigen
#endif