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

Change skipU argument to computeU - this reverses the meaning.

Browse Source 
See "skipXxx / computeXxx parameters in Eigenvalues module" on mailing list.
This commit is contained in:
Jitse Niesen 2010-05-31 16:48:41 +01:00
parent c21390a611
commit 609941380a
3 changed files with 35 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/src/Array/VectorwiseOp.h
View File

@ -380,8 +380,8 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
/** \returns a matrix expression /** \returns a matrix expression
* where each column (or row) are reversed. * where each column (or row) are reversed.
* *
* Example: \include VectorWise_reverse.cpp * Example: \include Vectorwise_reverse.cpp
* Output: \verbinclude VectorWise_reverse.out * Output: \verbinclude Vectorwise_reverse.out
* *
* \sa DenseBase::reverse() */ * \sa DenseBase::reverse() */
const Reverse<ExpressionType, Direction> reverse() const const Reverse<ExpressionType, Direction> reverse() const

57
Eigen/src/Eigenvalues/ComplexSchur.h
View File

@ -110,21 +110,21 @@ template<typename _MatrixType> class ComplexSchur
/** \brief Constructor; computes Schur decomposition of given matrix. /** \brief Constructor; computes Schur decomposition of given matrix.
* *
* \param[in] matrix Square matrix whose Schur decomposition is to be computed. * \param[in] matrix Square matrix whose Schur decomposition is to be computed.
* \param[in] skipU If true, then the unitary matrix U in the decomposition is not computed. * \param[in] computeU If true, both T and U are computed; if false, only T is computed.
* *
* This constructor calls compute() to compute the Schur decomposition. * This constructor calls compute() to compute the Schur decomposition.
* *
* \sa matrixT() and matrixU() for examples. * \sa matrixT() and matrixU() for examples.
*/ */
ComplexSchur(const MatrixType& matrix, bool skipU = false) ComplexSchur(const MatrixType& matrix, bool computeU = true)
: m_matT(matrix.rows(),matrix.cols()), : m_matT(matrix.rows(),matrix.cols()),
m_matU(matrix.rows(),matrix.cols()), m_matU(matrix.rows(),matrix.cols()),
m_hess(matrix.rows()), m_hess(matrix.rows()),
m_isInitialized(false), m_isInitialized(false),
m_matUisUptodate(false) m_matUisUptodate(false)
{ {
compute(matrix, skipU); compute(matrix, computeU);
} }
/** \brief Returns the unitary matrix in the Schur decomposition. /** \brief Returns the unitary matrix in the Schur decomposition.
@ -132,10 +132,10 @@ template<typename _MatrixType> class ComplexSchur
* \returns A const reference to the matrix U. * \returns A const reference to the matrix U.
* *
* It is assumed that either the constructor * It is assumed that either the constructor
* ComplexSchur(const MatrixType& matrix, bool skipU) or the * ComplexSchur(const MatrixType& matrix, bool computeU) or the
* member function compute(const MatrixType& matrix, bool skipU) * member function compute(const MatrixType& matrix, bool computeU)
* has been called before to compute the Schur decomposition of a * has been called before to compute the Schur decomposition of a
* matrix, and that \p skipU was set to false (the default * matrix, and that \p computeU was set to true (the default
* value). * value).
* *
* Example: \include ComplexSchur_matrixU.cpp * Example: \include ComplexSchur_matrixU.cpp
@ -153,8 +153,8 @@ template<typename _MatrixType> class ComplexSchur
* \returns A const reference to the matrix T. * \returns A const reference to the matrix T.
* *
* It is assumed that either the constructor * It is assumed that either the constructor
* ComplexSchur(const MatrixType& matrix, bool skipU) or the * ComplexSchur(const MatrixType& matrix, bool computeU) or the
* member function compute(const MatrixType& matrix, bool skipU) * member function compute(const MatrixType& matrix, bool computeU)
* has been called before to compute the Schur decomposition of a * has been called before to compute the Schur decomposition of a
* matrix. * matrix.
* *
@ -174,7 +174,7 @@ template<typename _MatrixType> class ComplexSchur
/** \brief Computes Schur decomposition of given matrix. /** \brief Computes Schur decomposition of given matrix.
* *
* \param[in] matrix Square matrix whose Schur decomposition is to be computed. * \param[in] matrix Square matrix whose Schur decomposition is to be computed.
* \param[in] skipU If true, then the unitary matrix U in the decomposition is not computed. * \param[in] computeU If true, both T and U are computed; if false, only T is computed.
* *
* The Schur decomposition is computed by first reducing the * The Schur decomposition is computed by first reducing the
* matrix to Hessenberg form using the class * matrix to Hessenberg form using the class
@ -182,13 +182,14 @@ template<typename _MatrixType> class ComplexSchur
* to triangular form by performing QR iterations with a single * to triangular form by performing QR iterations with a single
* shift. The cost of computing the Schur decomposition depends * shift. The cost of computing the Schur decomposition depends
* on the number of iterations; as a rough guide, it may be taken * on the number of iterations; as a rough guide, it may be taken
* on the number of iterations; as a rough guide, it may be taken
* to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops
* if \a skipU is true. * if \a computeU is false.
* *
* Example: \include ComplexSchur_compute.cpp * Example: \include ComplexSchur_compute.cpp
* Output: \verbinclude ComplexSchur_compute.out * Output: \verbinclude ComplexSchur_compute.out
*/ */
void compute(const MatrixType& matrix, bool skipU = false); void compute(const MatrixType& matrix, bool computeU = true);
protected: protected:
ComplexMatrixType m_matT, m_matU; ComplexMatrixType m_matT, m_matU;
@ -199,7 +200,7 @@ template<typename _MatrixType> class ComplexSchur
private: private:
bool subdiagonalEntryIsNeglegible(Index i); bool subdiagonalEntryIsNeglegible(Index i);
ComplexScalar computeShift(Index iu, Index iter); ComplexScalar computeShift(Index iu, Index iter);
void reduceToTriangularForm(bool skipU); void reduceToTriangularForm(bool computeU);
friend struct ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>; friend struct ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
}; };
@ -295,22 +296,22 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
template<typename MatrixType> template<typename MatrixType>
void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
{ {
m_matUisUptodate = false; m_matUisUptodate = false;
ei_assert(matrix.cols() == matrix.rows()); ei_assert(matrix.cols() == matrix.rows());
if(matrix.cols() == 1) if(matrix.cols() == 1)
{ {
m_matU = ComplexMatrixType::Identity(1,1); m_matT = matrix.template cast<ComplexScalar>();
if(!skipU) m_matT = matrix.template cast<ComplexScalar>(); if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
m_isInitialized = true; m_isInitialized = true;
m_matUisUptodate = !skipU; m_matUisUptodate = computeU;
return; return;
} }
ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, skipU); ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
reduceToTriangularForm(skipU); reduceToTriangularForm(computeU);
} }
/* Reduce given matrix to Hessenberg form */ /* Reduce given matrix to Hessenberg form */
@ -318,28 +319,26 @@ template<typename MatrixType, bool IsComplex>
struct ei_complex_schur_reduce_to_hessenberg struct ei_complex_schur_reduce_to_hessenberg
{ {
// this is the implementation for the case IsComplex = true // this is the implementation for the case IsComplex = true
static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool skipU) static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
{ {
// TODO skip Q if skipU = true
_this.m_hess.compute(matrix); _this.m_hess.compute(matrix);
_this.m_matT = _this.m_hess.matrixH(); _this.m_matT = _this.m_hess.matrixH();
if(!skipU) _this.m_matU = _this.m_hess.matrixQ(); if(computeU) _this.m_matU = _this.m_hess.matrixQ();
} }
}; };
template<typename MatrixType> template<typename MatrixType>
struct ei_complex_schur_reduce_to_hessenberg<MatrixType, false> struct ei_complex_schur_reduce_to_hessenberg<MatrixType, false>
{ {
static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool skipU) static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
{ {
typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar; typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType; typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
// Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
// TODO skip Q if skipU = true
_this.m_hess.compute(matrix); _this.m_hess.compute(matrix);
_this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>(); _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
if(!skipU) if(computeU)
{ {
// This may cause an allocation which seems to be avoidable // This may cause an allocation which seems to be avoidable
MatrixType Q = _this.m_hess.matrixQ(); MatrixType Q = _this.m_hess.matrixQ();
@ -350,7 +349,7 @@ struct ei_complex_schur_reduce_to_hessenberg<MatrixType, false>
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration. // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
template<typename MatrixType> template<typename MatrixType>
void ComplexSchur<MatrixType>::reduceToTriangularForm(bool skipU) void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
{ {
// The matrix m_matT is divided in three parts. // The matrix m_matT is divided in three parts.
// Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
@ -393,7 +392,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool skipU)
rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il)); rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint()); m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
m_matT.topRows(std::min(il+2,iu)+1).applyOnTheRight(il, il+1, rot); m_matT.topRows(std::min(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
if(!skipU) m_matU.applyOnTheRight(il, il+1, rot); if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
for(Index i=il+1 ; i<iu ; i++) for(Index i=il+1 ; i<iu ; i++)
{ {
@ -401,7 +400,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool skipU)
m_matT.coeffRef(i+1,i-1) = ComplexScalar(0); m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint()); m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
m_matT.topRows(std::min(i+2,iu)+1).applyOnTheRight(i, i+1, rot); m_matT.topRows(std::min(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
if(!skipU) m_matU.applyOnTheRight(i, i+1, rot); if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
} }
} }
@ -413,7 +412,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool skipU)
} }
m_isInitialized = true; m_isInitialized = true;
m_matUisUptodate = !skipU; m_matUisUptodate = computeU;
} }
#endif // EIGEN_COMPLEX_SCHUR_H #endif // EIGEN_COMPLEX_SCHUR_H

5
test/schur_complex.cpp
View File

@ -56,6 +56,11 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
ComplexSchur<MatrixType> cs2(A); ComplexSchur<MatrixType> cs2(A);
VERIFY_IS_EQUAL(cs1.matrixT(), cs2.matrixT()); VERIFY_IS_EQUAL(cs1.matrixT(), cs2.matrixT());
VERIFY_IS_EQUAL(cs1.matrixU(), cs2.matrixU()); VERIFY_IS_EQUAL(cs1.matrixU(), cs2.matrixU());
// Test computation of only T, not U
ComplexSchur<MatrixType> csOnlyT(A, false);
VERIFY_IS_EQUAL(cs1.matrixT(), csOnlyT.matrixT());
VERIFY_RAISES_ASSERT(csOnlyT.matrixU());
} }
void test_schur_complex() void test_schur_complex()