improve ComplexShur api and doc

This commit is contained in:
Gael Guennebaud 2009-09-16 14:34:08 +02:00
parent a4fd0aa25b
commit 77f858f6ab
2 changed files with 43 additions and 17 deletions

View File

@ -31,8 +31,15 @@
* *
* \class ComplexShur * \class ComplexShur
* *
* \brief Performs a complex Shur decomposition of a real or complex square matrix * \brief Performs a complex Schur decomposition of a real or complex square matrix
* *
* Given a real or complex square matrix A, this class computes the Schur decomposition:
* \f$ A = U T U^*\f$ where U is a unitary complex matrix, and T is a complex upper
* triangular matrix.
*
* The diagonal of the matrix T corresponds to the eigenvalues of the matrix A.
*
* \sa class RealSchur, class EigenSolver
*/ */
template<typename _MatrixType> class ComplexSchur template<typename _MatrixType> class ComplexSchur
{ {
@ -42,41 +49,56 @@ template<typename _MatrixType> class ComplexSchur
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> Complex; typedef std::complex<RealScalar> Complex;
typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> ComplexMatrixType; typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> ComplexMatrixType;
enum {
Size = MatrixType::RowsAtCompileTime
};
/** /** \brief Default Constructor.
* \brief Default Constructor.
* *
* The default constructor is useful in cases in which the user intends to * The default constructor is useful in cases in which the user intends to
* perform decompositions via ComplexSchur::compute(const MatrixType&). * perform decompositions via ComplexSchur::compute().
*/ */
ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false) ComplexSchur(int size = Size==Dynamic ? 0 : Size)
: m_matT(size,size), m_matU(size,size), m_isInitialized(false), m_matUisUptodate(false)
{} {}
ComplexSchur(const MatrixType& matrix) /** Constructor computing the Schur decomposition of the matrix \a matrix.
* If \a skipU is true, then the matrix U is not computed. */
ComplexSchur(const MatrixType& matrix, bool skipU = false)
: 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_isInitialized(false) m_isInitialized(false),
m_matUisUptodate(false)
{ {
compute(matrix); compute(matrix, skipU);
} }
ComplexMatrixType matrixU() const /** \returns a const reference to the matrix U of the respective Schur decomposition. */
const ComplexMatrixType& matrixU() const
{ {
ei_assert(m_isInitialized && "ComplexSchur is not initialized."); ei_assert(m_isInitialized && "ComplexSchur is not initialized.");
ei_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
return m_matU; return m_matU;
} }
ComplexMatrixType matrixT() const /** \returns a const reference to the matrix T of the respective Schur decomposition.
* Note that this function returns a plain square matrix. If you want to reference
* only the upper triangular part, use:
* \code schur.matrixT().triangularView<Upper>() \endcode. */
const ComplexMatrixType& matrixT() const
{ {
ei_assert(m_isInitialized && "ComplexShur is not initialized."); ei_assert(m_isInitialized && "ComplexShur is not initialized.");
return m_matT; return m_matT;
} }
void compute(const MatrixType& matrix); /** Computes the Schur decomposition of the matrix \a matrix.
* If \a skipU is true, then the matrix U is not computed. */
void compute(const MatrixType& matrix, bool skipU = false);
protected: protected:
ComplexMatrixType m_matT, m_matU; ComplexMatrixType m_matT, m_matU;
bool m_isInitialized; bool m_isInitialized;
bool m_matUisUptodate;
}; };
/** Computes the principal value of the square root of the complex \a z. */ /** Computes the principal value of the square root of the complex \a z. */
@ -117,17 +139,20 @@ std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
} }
template<typename MatrixType> template<typename MatrixType>
void ComplexSchur<MatrixType>::compute(const MatrixType& matrix) void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
{ {
// this code is inspired from Jampack // this code is inspired from Jampack
m_matUisUptodate = false;
assert(matrix.cols() == matrix.rows()); assert(matrix.cols() == matrix.rows());
int n = matrix.cols(); int n = matrix.cols();
// Reduce to Hessenberg form // Reduce to Hessenberg form
// TODO skip Q if skipU = true
HessenbergDecomposition<MatrixType> hess(matrix); HessenbergDecomposition<MatrixType> hess(matrix);
m_matT = hess.matrixH(); m_matT = hess.matrixH();
m_matU = hess.matrixQ(); if(!skipU) m_matU = hess.matrixQ();
int iu = m_matT.cols() - 1; int iu = m_matT.cols() - 1;
int il; int il;
@ -206,7 +231,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
{ {
m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint()); m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint());
m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot); m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot);
m_matU.applyOnTheRight(i, i+1, rot); if(!skipU) m_matU.applyOnTheRight(i, i+1, rot);
if(i != iu-1) if(i != iu-1)
{ {
@ -232,6 +257,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
*/ */
m_isInitialized = true; m_isInitialized = true;
m_matUisUptodate = !skipU;
} }
#endif // EIGEN_COMPLEX_SCHUR_H #endif // EIGEN_COMPLEX_SCHUR_H

View File

@ -88,14 +88,14 @@ template<typename _MatrixType> class HessenbergDecomposition
_compute(m_matrix, m_hCoeffs); _compute(m_matrix, m_hCoeffs);
} }
/** \returns the householder coefficients allowing to /** \returns a const reference to the householder coefficients allowing to
* reconstruct the matrix Q from the packed data. * reconstruct the matrix Q from the packed data.
* *
* \sa packedMatrix() * \sa packedMatrix()
*/ */
CoeffVectorType householderCoefficients() const { return m_hCoeffs; } const CoeffVectorType& householderCoefficients() const { return m_hCoeffs; }
/** \returns the internal result of the decomposition. /** \returns a const reference to the internal representation of the decomposition.
* *
* The returned matrix contains the following information: * The returned matrix contains the following information:
* - the upper part and lower sub-diagonal represent the Hessenberg matrix H * - the upper part and lower sub-diagonal represent the Hessenberg matrix H