| [eigen-commits] commit/eigen: 4 new changesets |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen-commits Archives
]
4 new commits in eigen:
https://bitbucket.org/eigen/eigen/changeset/7419a8d76ac7/
changeset: 7419a8d76ac7
user: adolfo.rodriguez
date: 2011-10-31 04:55:16
summary: Bug 206 - part 2: For HouseholderSequence objects, added non-allocating versions of evalTo() and applyThisOnTheRight/Left that take additional working vector parameters.
affected #: 1 file
diff -r 4b14b9286929eaf945c524235e9ba7154caa31a6 -r 7419a8d76ac708a4c964aa3e799225d6ba255733 Eigen/src/Householder/HouseholderSequence.h
--- a/Eigen/src/Householder/HouseholderSequence.h
+++ b/Eigen/src/Householder/HouseholderSequence.h
@@ -237,13 +237,20 @@
ConjugateReturnType inverse() const { return adjoint(); }
/** \internal */
- template<typename DestType> void evalTo(DestType& dst) const
+ template<typename DestType> inline void evalTo(DestType& dst) const
{
+ Matrix<Scalar, DestType::RowsAtCompileTime, 1,
+ AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
+ evalTo(dst, workspace);
+ }
+
+ /** \internal */
+ template<typename Dest, typename Workspace>
+ void evalTo(Dest& dst, Workspace& workspace) const
+ {
+ workspace.resize(rows());
Index vecs = m_length;
- // FIXME find a way to pass this temporary if the user wants to
- Matrix<Scalar, DestType::RowsAtCompileTime, 1,
- AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> temp(rows());
- if( internal::is_same<typename internal::remove_all<VectorsType>::type,DestType>::value
+ if( internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value
&& internal::extract_data(dst) == internal::extract_data(m_vectors))
{
// in-place
@@ -254,10 +261,10 @@
Index cornerSize = rows() - k - m_shift;
if(m_trans)
dst.bottomRightCorner(cornerSize, cornerSize)
- .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
+ .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
else
dst.bottomRightCorner(cornerSize, cornerSize)
- .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
+ .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
// clear the off diagonal vector
dst.col(k).tail(rows()-k-1).setZero();
@@ -274,10 +281,10 @@
Index cornerSize = rows() - k - m_shift;
if(m_trans)
dst.bottomRightCorner(cornerSize, cornerSize)
- .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
+ .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
else
dst.bottomRightCorner(cornerSize, cornerSize)
- .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
+ .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
}
}
}
@@ -285,24 +292,40 @@
/** \internal */
template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
{
- Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
+ Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
+ applyThisOnTheRight(dst, workspace);
+ }
+
+ /** \internal */
+ template<typename Dest, typename Workspace>
+ inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
+ {
+ workspace.resize(dst.rows());
for(Index k = 0; k < m_length; ++k)
{
Index actual_k = m_trans ? m_length-k-1 : k;
dst.rightCols(rows()-m_shift-actual_k)
- .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs..coeff(actual_k), &temp.coeffRef(0));
+ .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs..coeff(actual_k), workspace.data());
}
}
/** \internal */
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
{
- Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
+ Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
+ applyThisOnTheLeft(dst, workspace);
+ }
+
+ /** \internal */
+ template<typename Dest, typename Workspace>
+ inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
+ {
+ workspace.resize(dst.cols());
for(Index k = 0; k < m_length; ++k)
{
Index actual_k = m_trans ? k : m_length-k-1;
dst.bottomRows(rows()-m_shift-actual_k)
- .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
+ .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
}
}
https://bitbucket.org/eigen/eigen/changeset/053b0c623e42/
changeset: 053b0c623e42
user: adolfo.rodriguez
date: 2011-03-08 19:04:31
summary: Bug 206 - part 3: Reimplement FullPivHouseholderQR<T>::matrixQ() using ReturnByValue
affected #: 1 file
diff -r 7419a8d76ac708a4c964aa3e799225d6ba255733 -r 053b0c623e42578cb5c8b2bfc7a73822c38cc2ea Eigen/src/QR/FullPivHouseholderQR.h
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -26,6 +26,18 @@
#ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
#define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
+namespace internal {
+
+template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
+
+template<typename MatrixType>
+struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
+{
+ typedef typename MatrixType::PlainObject ReturnType;
+};
+
+}
+
/** \ingroup QR_Module
*
* \class FullPivHouseholderQR
@@ -62,7 +74,7 @@
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
- typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
+ typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
@@ -139,7 +151,9 @@
return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
}
- MatrixQType matrixQ(void) const;
+ /** \returns Expression object representing the matrix Q
+ */
+ MatrixQReturnType matrixQ(void) const;
/** \returns a reference to the matrix where the Householder QR decomposition is stored
*/
@@ -508,28 +522,73 @@
}
};
+/** \ingroup QR_Module
+ *
+ * \brief Expression type for return value of FullPivHouseholderQR::matrixQ()
+ *
+ * \tparam MatrixType type of underlying dense matrix
+ */
+template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
+ : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
+{
+public:
+ typedef typename MatrixType::Index Index;
+ typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
+ typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
+ typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
+ MatrixType::MaxRowsAtCompileTime> WorkVectorType;
+
+ FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
+ const HCoeffsType& hCoeffs,
+ const IntColVectorType& rowsTranspositions)
+ : m_qr(qr),
+ m_hCoeffs(hCoeffs),
+ m_rowsTranspositions(rowsTranspositions)
+ {}
+
+ template <typename ResultType>
+ void evalTo(ResultType& result) const
+ {
+ const Index rows = m_qr.rows();
+ WorkVectorType workspace(rows);
+ evalTo(result, workspace);
+ }
+
+ template <typename ResultType>
+ void evalTo(ResultType& result, WorkVectorType& workspace) const
+ {
+ // compute the product H'_0 H'_1 ... H'_n-1,
+ // where H_k is the k-th Householder transformation I - h_k v_k v_k'
+ // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
+ const Index rows = m_qr.rows();
+ const Index cols = m_qr.cols();
+ const Index size = std::min(rows, cols);
+ workspace.resize(rows);
+ result.setIdentity(rows, rows);
+ for (Index k = size-1; k >= 0; k--)
+ {
+ result.block(k, k, rows-k, rows-k)
+ .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
+ result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
+ }
+ }
+
+ Index rows() const { return m_qr.rows(); }
+ Index cols() const { return m_qr.rows(); }
+
+protected:
+ const typename MatrixType::Nested m_qr;
+ const typename HCoeffsType::Nested m_hCoeffs;
+ const typename IntColVectorType::Nested m_rowsTranspositions;
+};
+
} // end namespace internal
-/** \returns the matrix Q */
template<typename MatrixType>
-typename FullPivHouseholderQR<MatrixType>::MatrixQType FullPivHouseholderQR<MatrixType>::matrixQ() const
+inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
{
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized..");
- // compute the product H'_0 H'_1 ... H'_n-1,
- // where H_k is the k-th Householder transformation I - h_k v_k v_k'
- // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ....]
- Index rows = m_qr.rows();
- Index cols = m_qr.cols();
- Index size = (std::min)(rows,cols);
- MatrixQType res = MatrixQType::Identity(rows, rows);
- Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
- for (Index k = size-1; k >= 0; k--)
- {
- res.block(k, k, rows-k, rows-k)
- .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
- res.row(k).swap(res.row(m_rows_transpositions.coeff(k)));
- }
- return res;
+ return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
}
/** \return the full-pivoting Householder QR decomposition of \c *this.
https://bitbucket.org/eigen/eigen/changeset/0fa573497f3c/
changeset: 0fa573497f3c
user: adolfo.rodriguez
date: 2011-10-31 04:55:20
summary: Bug 206 - part 4: Removes heap allocations from JacobiSVD and its preconditioners
affected #: 2 files
diff -r 053b0c623e42578cb5c8b2bfc7a73822c38cc2ea -r 0fa573497f3c6aba3ccf579fb951e4f089596298 Eigen/src/QR/FullPivHouseholderQR.h
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -562,7 +562,7 @@
// and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
const Index rows = m_qr.rows();
const Index cols = m_qr.cols();
- const Index size = std::min(rows, cols);
+ const Index size = (std::min)(rows, cols);
workspace.resize(rows);
result.setIdentity(rows, rows);
for (Index k = size-1; k >= 0; k--)
diff -r 053b0c623e42578cb5c8b2bfc7a73822c38cc2ea -r 0fa573497f3c6aba3ccf579fb951e4f089596298 Eigen/src/SVD/JacobiSVD.h
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -61,9 +61,12 @@
> struct qr_preconditioner_impl {};
template<typename MatrixType, int QRPreconditioner, int Case>
-struct qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
+class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
{
- static bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
+public:
+ typedef typename MatrixType::Index Index;
+ void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {}
+ bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
{
return false;
}
@@ -72,134 +75,279 @@
/*** preconditioner using FullPivHouseholderQR ***/
template<typename MatrixType>
-struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
+class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
{
- static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
+public:
+ typedef typename MatrixType::Index Index;
+ typedef typename MatrixType::Scalar Scalar;
+ enum
+ {
+ RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+ MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
+ };
+ typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
+
+ void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
+ {
+ if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
+ {
+ m_qr = FullPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
+ }
+ if (svd.m_computeFullU) m_workspace.resize(svd.rows());
+ }
+
+ bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
{
if(matrix.rows() > matrix.cols())
{
- FullPivHouseholderQR<MatrixType> qr(matrix);
- svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
- if(svd.m_computeFullU) svd.m_matrixU = qr.matrixQ();
- if(svd.computeV()) svd.m_matrixV = qr.colsPermutation();
+ m_qr.compute(matrix);
+ svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
+ if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
+ if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
return true;
}
return false;
}
+private:
+ FullPivHouseholderQR<MatrixType> m_qr;
+ WorkspaceType m_workspace;
};
template<typename MatrixType>
-struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
+class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
{
- static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
+public:
+ typedef typename MatrixType::Index Index;
+ typedef typename MatrixType::Scalar Scalar;
+ enum
+ {
+ RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+ ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+ MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+ MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
+ Options = MatrixType::Options
+ };
+ typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
+ TransposeTypeWithSameStorageOrder;
+
+ void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
+ {
+ if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
+ {
+ m_qr = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd..cols(), svd.rows());
+ }
+ m_adjoint.resize(svd.cols(), svd.rows());
+ if (svd.m_computeFullV) m_workspace.resize(svd.cols());
+ }
+
+ bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
{
if(matrix.cols() > matrix.rows())
{
- typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
- MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
- TransposeTypeWithSameStorageOrder;
- FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
- svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
- if(svd.m_computeFullV) svd.m_matrixV = qr.matrixQ();
- if(svd.computeU()) svd.m_matrixU = qr.colsPermutation();
+ m_adjoint = matrix.adjoint();
+ m_qr.compute(m_adjoint);
+ svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
+ if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
+ if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
return true;
}
else return false;
}
+private:
+ FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
+ TransposeTypeWithSameStorageOrder m_adjoint;
+ typename internal::plain_row_type<MatrixType>::type m_workspace;
};
/*** preconditioner using ColPivHouseholderQR ***/
template<typename MatrixType>
-struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
+class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
{
- static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
+public:
+ typedef typename MatrixType::Index Index;
+
+ void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
+ {
+ if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
+ {
+ m_qr = ColPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
+ }
+ if (svd.m_computeFullU) m_workspace.resize(svd.rows());
+ else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
+ }
+
+ bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
{
if(matrix.rows() > matrix.cols())
{
- ColPivHouseholderQR<MatrixType> qr(matrix);
- svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
- if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ();
- else if(svd.m_computeThinU) {
+ m_qr.compute(matrix);
+ svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
+ if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
+ else if(svd.m_computeThinU)
+ {
svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
- qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
+ m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
}
- if(svd.computeV()) svd.m_matrixV = qr.colsPermutation();
+ if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
return true;
}
return false;
}
+
+private:
+ ColPivHouseholderQR<MatrixType> m_qr;
+ typename internal::plain_col_type<MatrixType>::type m_workspace;
};
template<typename MatrixType>
-struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
+class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
{
- static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
+public:
+ typedef typename MatrixType::Index Index;
+ typedef typename MatrixType::Scalar Scalar;
+ enum
+ {
+ RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+ ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+ MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+ MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
+ Options = MatrixType::Options
+ };
+
+ typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
+ TransposeTypeWithSameStorageOrder;
+
+ void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
+ {
+ if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
+ {
+ m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
+ }
+ if (svd.m_computeFullV) m_workspace.resize(svd.cols());
+ else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
+ m_adjoint.resize(svd.cols(), svd.rows());
+ }
+
+ bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
{
if(matrix.cols() > matrix.rows())
{
- typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
- MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
- TransposeTypeWithSameStorageOrder;
- ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
- svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
- if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ();
- else if(svd.m_computeThinV) {
+ m_adjoint = matrix.adjoint();
+ m_qr.compute(m_adjoint);
+
+ svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
+ if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
+ else if(svd.m_computeThinV)
+ {
svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
- qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
+ m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
}
- if(svd.computeU()) svd.m_matrixU = qr.colsPermutation();
+ if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
return true;
}
else return false;
}
+
+private:
+ ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
+ TransposeTypeWithSameStorageOrder m_adjoint;
+ typename internal::plain_row_type<MatrixType>::type m_workspace;
};
/*** preconditioner using HouseholderQR ***/
template<typename MatrixType>
-struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
+class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
{
- static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
+public:
+ typedef typename MatrixType::Index Index;
+
+ void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
+ {
+ if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
+ {
+ m_qr = HouseholderQR<MatrixType>(svd.rows(), svd.cols());
+ }
+ if (svd.m_computeFullU) m_workspace.resize(svd.rows());
+ else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
+ }
+
+ bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
{
if(matrix.rows() > matrix.cols())
{
- HouseholderQR<MatrixType> qr(matrix);
- svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
- if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ();
- else if(svd.m_computeThinU) {
+ m_qr.compute(matrix);
+ svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
+ if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
+ else if(svd.m_computeThinU)
+ {
svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
- qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
+ m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
}
if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
return true;
}
return false;
}
+private:
+ HouseholderQR<MatrixType> m_qr;
+ typename internal::plain_col_type<MatrixType>::type m_workspace;
};
template<typename MatrixType>
-struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
+class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
{
- static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
+public:
+ typedef typename MatrixType::Index Index;
+ typedef typename MatrixType::Scalar Scalar;
+ enum
+ {
+ RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+ ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+ MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+ MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
+ Options = MatrixType::Options
+ };
+
+ typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
+ TransposeTypeWithSameStorageOrder;
+
+ void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
+ {
+ if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
+ {
+ m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
+ }
+ if (svd.m_computeFullV) m_workspace.resize(svd.cols());
+ else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
+ m_adjoint.resize(svd.cols(), svd.rows());
+ }
+
+ bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
{
if(matrix.cols() > matrix.rows())
{
- typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
- MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
- TransposeTypeWithSameStorageOrder;
- HouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
- svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
- if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ();
- else if(svd.m_computeThinV) {
+ m_adjoint = matrix.adjoint();
+ m_qr.compute(m_adjoint);
+
+ svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
+ if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
+ else if(svd.m_computeThinV)
+ {
svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
- qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
+ m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
}
if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
return true;
}
else return false;
}
+
+private:
+ HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
+ TransposeTypeWithSameStorageOrder m_adjoint;
+ typename internal::plain_row_type<MatrixType>::type m_workspace;
};
/*** 2x2 SVD implementation
@@ -316,7 +464,7 @@
* Here's an example demonstrating basic usage:
* \include JacobiSVD_basic.cpp
* Output: \verbinclude JacobiSVD_basic.out
- *
+ *
* This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than
* bidiagonalizing SVD algorithms for large square matrices; however its complexity is still \f$ O(n^2p) \f$ where \a n is the smaller dimension and
* \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms.
@@ -324,7 +472,7 @@
*
* If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to
* terminate in finite (and reasonable) time.
- *
+ *
* The possible values for QRPreconditioner are:
* \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR.
* \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR.
@@ -494,7 +642,7 @@
* \param b the right-hand-side of the equation to solve.
*
* \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V.
- *
+ *
* \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving.
* In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$.
*/
@@ -535,6 +683,9 @@
friend struct internal::svd_precondition_2x2_block_to_be_real;
template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
friend struct internal::qr_preconditioner_impl;
+
+ internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
+ internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
};
template<typename MatrixType, int QRPreconditioner>
@@ -578,6 +729,9 @@
: m_computeThinV ? m_diagSize
: 0);
m_workMatrix.resize(m_diagSize, m_diagSize);
+
+ m_qr_precond_morecols.allocate(*this);
+ m_qr_precond_morerows.allocate(*this);
}
template<typename MatrixType, int QRPreconditioner>
@@ -595,8 +749,7 @@
/*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
- if(!internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>::run(*this, matrix)
- && !internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols>::run(*this, matrix))
+ if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix))
{
m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
https://bitbucket.org/eigen/eigen/changeset/13342df27e02/
changeset: 13342df27e02
user: bjacob
date: 2011-10-31 04:55:20
summary: Fix some unused-variable warnings with GCC 4.6
affected #: 5 files
diff -r 0fa573497f3c6aba3ccf579fb951e4f089596298 -r 13342df27e02bb7df69e151318355fbf513d7db8 test/basicstuff.cpp
--- a/test/basicstuff.cpp
+++ b/test/basicstuff.cpp
@@ -42,11 +42,8 @@
m2 = MatrixType::Random(rows, cols),
m3(rows, cols),
mzero = MatrixType::Zero(rows, cols),
- identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
- ::Identity(rows, rows),
square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>::Random(rows, rows);
VectorType v1 = VectorType::Random(rows),
- v2 = VectorType::Random(rows),
vzero = VectorType::Zero(rows);
SquareMatrixType sm1 = SquareMatrixType::Random(rows,rows), sm2(rows,rows);
diff -r 0fa573497f3c6aba3ccf579fb951e4f089596298 -r 13342df27e02bb7df69e151318355fbf513d7db8 test/eigensolver_complex.cpp
--- a/test/eigensolver_complex.cpp
+++ b/test/eigensolver_complex.cpp
@@ -125,4 +125,6 @@
// Test problem size constructors
CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf>(s));
+
+ EIGEN_UNUSED_VARIABLE(s)
}
diff -r 0fa573497f3c6aba3ccf579fb951e4f089596298 -r 13342df27e02bb7df69e151318355fbf513d7db8 test/eigensolver_generic.cpp
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -118,4 +118,6 @@
// Test problem size constructors
CALL_SUBTEST_5(EigenSolver<MatrixXf>(s));
+
+ EIGEN_UNUSED_VARIABLE(s)
}
diff -r 0fa573497f3c6aba3ccf579fb951e4f089596298 -r 13342df27e02bb7df69e151318355fbf513d7db8 test/eigensolver_selfadjoint.cpp
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -204,5 +204,7 @@
s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf>(s));
CALL_SUBTEST_8(Tridiagonalization<MatrixXf>(s));
+
+ EIGEN_UNUSED_VARIABLE(s)
}
diff -r 0fa573497f3c6aba3ccf579fb951e4f089596298 -r 13342df27e02bb7df69e151318355fbf513d7db8 test/linearstructure.cpp
--- a/test/linearstructure.cpp
+++ b/test/linearstructure.cpp
@@ -39,8 +39,7 @@
// to test it, hence I consider that we will have tested Random.h
MatrixType m1 = MatrixType::Random(rows, cols),
m2 = MatrixType::Random(rows, cols),
- m3(rows, cols),
- mzero = MatrixType::Zero(rows, cols);
+ m3(rows, cols);
Scalar s1 = internal::random<Scalar>();
while (internal::abs(s1)<1e-3) s1 = internal::random<Scalar>();
Repository URL: https://bitbucket.org/eigen/eigen/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.