GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-03 10:14:04 +08:00
Code Issues Packages Projects Releases Wiki Activity

Change return type of matrixH() method to HouseholderSequence.

Browse Source 
This method is a member of Tridiagonalization and HessenbergDecomposition.
This commit is contained in:
Jitse Niesen 2010-05-24 17:35:54 +01:00
parent 76dd0e5314
commit eb3ca68684
5 changed files with 92 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -3,6 +3,7 @@
//
// Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -26,6 +27,8 @@
#ifndef EIGEN_COMPLEX_SCHUR_H
#define EIGEN_COMPLEX_SCHUR_H
template<typename MatrixType, bool IsComplex> struct ei_complex_schur_reduce_to_hessenberg;
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
@ -50,6 +53,8 @@
* decomposition is computed, you can use the matrixU() and matrixT()
* functions to retrieve the matrices U and V in the decomposition.
*
* \note This code is inspired from Jampack
*
* \sa class RealSchur, class EigenSolver, class ComplexEigenSolver
*/
template<typename _MatrixType> class ComplexSchur
@ -194,6 +199,8 @@ template<typename _MatrixType> class ComplexSchur
private:
bool subdiagonalEntryIsNeglegible(int i);
ComplexScalar computeShift(int iu, int iter);
void reduceToTriangularForm(bool skipU);
friend struct ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
};
/** Computes the principal value of the square root of the complex \a z. */
@ -290,12 +297,10 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
template<typename MatrixType>
void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
{
// this code is inspired from Jampack
m_matUisUptodate = false;
ei_assert(matrix.cols() == matrix.rows());
int n = matrix.cols();
if(n==1)
if(matrix.cols() == 1)
{
m_matU = ComplexMatrixType::Identity(1,1);
if(!skipU) m_matT = matrix.template cast<ComplexScalar>();
@ -304,15 +309,49 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
return;
}
// Reduce to Hessenberg form
// TODO skip Q if skipU = true
m_hess.compute(matrix);
ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, skipU);
reduceToTriangularForm(skipU);
}
m_matT = m_hess.matrixH().template cast<ComplexScalar>();
if(!skipU) m_matU = m_hess.matrixQ().template cast<ComplexScalar>();
/* Reduce given matrix to Hessenberg form */
template<typename MatrixType, bool IsComplex>
struct ei_complex_schur_reduce_to_hessenberg
{
// this is the implementation for the case IsComplex = true
static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool skipU)
{
// TODO skip Q if skipU = true
_this.m_hess.compute(matrix);
_this.m_matT = _this.m_hess.matrixH();
if(!skipU) _this.m_matU = _this.m_hess.matrixQ();
}
};
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
template<typename MatrixType>
struct ei_complex_schur_reduce_to_hessenberg<MatrixType, false>
{
static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool skipU)
{
typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
// 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_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
if(!skipU)
{
// This may cause an allocation which seems to be avoidable
MatrixType Q = _this.m_hess.matrixQ();
_this.m_matU = Q.template cast<ComplexScalar>();
}
}
};
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
template<typename MatrixType>
void ComplexSchur<MatrixType>::reduceToTriangularForm(bool skipU)
{
// 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 il,...,iu is the part we are working on (the active submatrix).
@ -352,7 +391,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
ComplexScalar shift = computeShift(iu, iter);
PlanarRotation<ComplexScalar> rot;
rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
m_matT.rightCols(n-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);
if(!skipU) m_matU.applyOnTheRight(il, il+1, rot);
@ -360,7 +399,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
{
rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
m_matT.rightCols(n-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);
if(!skipU) m_matU.applyOnTheRight(i, i+1, rot);
}

42
Eigen/src/Eigenvalues/HessenbergDecomposition.h
View File

@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -59,7 +60,9 @@ template<typename _MatrixType> class HessenbergDecomposition
{
public:
/** \brief Synonym for the template parameter \p _MatrixType. */
typedef _MatrixType MatrixType;
enum {
Size = MatrixType::RowsAtCompileTime,
SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
@ -68,17 +71,20 @@ template<typename _MatrixType> class HessenbergDecomposition
MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
};
/** \brief Scalar type for matrices of type \p _MatrixType. */
/** \brief Scalar type for matrices of type #MatrixType. */
typedef typename MatrixType::Scalar Scalar;
/** \brief Type for vector of Householder coefficients.
*
* This is column vector with entries of type #Scalar. The length of the
* vector is one less than the size of \p _MatrixType, if it is a
* fixed-side type.
* vector is one less than the size of #MatrixType, if it is a fixed-side
* type.
*/
typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
/** \brief Return type of matrixQ() */
typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
/** \brief Default constructor; the decomposition will be computed later.
*
* \param [in] size The size of the matrix whose Hessenberg decomposition will be computed.
@ -190,19 +196,22 @@ template<typename _MatrixType> class HessenbergDecomposition
/** \brief Reconstructs the orthogonal matrix Q in the decomposition
*
* \returns the matrix Q
* \returns object representing the matrix Q
*
* \pre Either the constructor HessenbergDecomposition(const MatrixType&)
* or the member function compute(const MatrixType&) has been called
* before to compute the Hessenberg decomposition of a matrix.
*
* This function reconstructs the matrix Q from the Householder
* coefficients and the packed matrix stored internally. This
* reconstruction requires \f$ 4n^3 / 3 \f$ flops.
* This function returns a light-weight object of template class
* HouseholderSequence. You can either apply it directly to a matrix or
* you can convert it to a matrix of type #MatrixType.
*
* \sa matrixH() for an example
* \sa matrixH() for an example, class HouseholderSequence
*/
MatrixType matrixQ() const;
HouseholderSequenceType matrixQ() const
{
return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
}
/** \brief Constructs the Hessenberg matrix H in the decomposition
*
@ -281,21 +290,6 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
}
}
template<typename MatrixType>
typename HessenbergDecomposition<MatrixType>::MatrixType
HessenbergDecomposition<MatrixType>::matrixQ() const
{
int n = m_matrix.rows();
MatrixType matQ = MatrixType::Identity(n,n);
VectorType temp(n);
for (int i = n-2; i>=0; i--)
{
matQ.bottomRightCorner(n-i-1,n-i-1)
.applyHouseholderOnTheLeft(m_matrix.col(i).tail(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &temp.coeffRef(0,0));
}
return matQ;
}
#endif // EIGEN_HIDE_HEAVY_CODE
template<typename MatrixType>

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

@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -61,7 +62,9 @@ template<typename _MatrixType> class Tridiagonalization
{
public:
/** \brief Synonym for the template parameter \p _MatrixType. */
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
@ -89,6 +92,9 @@ template<typename _MatrixType> class Tridiagonalization
Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 >
>::ret SubDiagonalReturnType;
/** \brief Return type of matrixQ() */
typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
/** \brief Default constructor.
*
* \param [in] size Positive integer, size of the matrix whose tridiagonal
@ -195,29 +201,25 @@ template<typename _MatrixType> class Tridiagonalization
*/
inline const MatrixType& packedMatrix() const { return m_matrix; }
/** \brief Reconstructs the unitary matrix Q in the decomposition
/** \brief Returns the unitary matrix Q in the decomposition
*
* \returns the matrix Q
* \returns object representing the matrix Q
*
* \pre Either the constructor Tridiagonalization(const MatrixType&) or
* the member function compute(const MatrixType&) has been called before
* to compute the tridiagonal decomposition of a matrix.
*
* This function reconstructs the matrix Q from the Householder
* coefficients and the packed matrix stored internally. This
* reconstruction requires \f$ 4n^3 / 3 \f$ flops.
* This function returns a light-weight object of template class
* HouseholderSequence. You can either apply it directly to a matrix or
* you can convert it to a matrix of type #MatrixType.
*
* \sa Tridiagonalization(const MatrixType&) for an example,
* matrixT(), matrixQInPlace()
* matrixT(), class HouseholderSequence
*/
MatrixType matrixQ() const;
/** \brief Reconstructs the unitary matrix Q in the decomposition
*
* This is an in-place variant of matrixQ() which avoids the copy.
* This function will probably be deleted soon.
*/
template<typename QDerived> void matrixQInPlace(MatrixBase<QDerived>* q) const;
HouseholderSequenceType matrixQ() const
{
return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
}
/** \brief Constructs the tridiagonal matrix T in the decomposition
*
@ -386,31 +388,6 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
}
}
template<typename MatrixType>
typename Tridiagonalization<MatrixType>::MatrixType
Tridiagonalization<MatrixType>::matrixQ() const
{
MatrixType matQ;
matrixQInPlace(&matQ);
return matQ;
}
template<typename MatrixType>
template<typename QDerived>
void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) const
{
QDerived& matQ = q->derived();
int n = m_matrix.rows();
matQ = MatrixType::Identity(n,n);
typedef typename ei_plain_row_type<MatrixType>::type RowVectorType;
RowVectorType aux(n);
for (int i = n-2; i>=0; i--)
{
matQ.bottomRightCorner(n-i-1,n-i-1)
.applyHouseholderOnTheLeft(m_matrix.col(i).tail(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &aux.coeffRef(0,0));
}
}
template<typename MatrixType>
void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
{
@ -426,7 +403,7 @@ void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalT
diag = tridiag.diagonal();
subdiag = tridiag.subDiagonal();
if (extractQ)
tridiag.matrixQInPlace(&mat);
mat = tridiag.matrixQ();
}
}

6
doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp
View File

@ -1,11 +1,9 @@
MatrixXd X = MatrixXd::Random(5,5);
MatrixXd A = X + X.transpose();
cout << "Here is a random symmetric 5x5 matrix:" << endl << A << endl << endl;
Tridiagonalization<MatrixXd> triOfA(A);
cout << "The orthogonal matrix Q is:" << endl << triOfA.matrixQ() << endl;
cout << "The tridiagonal matrix T is:" << endl << triOfA.matrixT() << endl << endl;
MatrixXd Q = triOfA.matrixQ();
cout << "The orthogonal matrix Q is:" << endl << Q << endl;
MatrixXd T = triOfA.matrixT();
cout << "The tridiagonal matrix T is:" << endl << T << endl << endl;
cout << "Q * T * Q^T = " << endl << Q * T * Q.transpose() << endl;

7
test/hessenberg.cpp
View File

@ -34,8 +34,9 @@ template<typename Scalar,int Size> void hessenberg(int size = Size)
for(int counter = 0; counter < g_repeat; ++counter) {
MatrixType m = MatrixType::Random(size,size);
HessenbergDecomposition<MatrixType> hess(m);
VERIFY_IS_APPROX(m, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
MatrixType Q = hess.matrixQ();
MatrixType H = hess.matrixH();
VERIFY_IS_APPROX(m, Q * H * Q.adjoint());
for(int row = 2; row < size; ++row) {
for(int col = 0; col < row-1; ++col) {
VERIFY(H(row,col) == (typename MatrixType::Scalar)0);
@ -48,8 +49,10 @@ template<typename Scalar,int Size> void hessenberg(int size = Size)
HessenbergDecomposition<MatrixType> cs1;
cs1.compute(A);
HessenbergDecomposition<MatrixType> cs2(A);
VERIFY_IS_EQUAL(cs1.matrixQ(), cs2.matrixQ());
VERIFY_IS_EQUAL(cs1.matrixH(), cs2.matrixH());
MatrixType cs1Q = cs1.matrixQ();
MatrixType cs2Q = cs2.matrixQ();
VERIFY_IS_EQUAL(cs1Q, cs2Q);
// TODO: Add tests for packedMatrix() and householderCoefficients()
}