- add MatrixBase::eigenvalues() convenience method

- add MatrixBase::matrixNorm(); in the non-selfadjoint case, we reduce to the
  selfadjoint case by using the "C*-identity" a.k.a.
	norm of x = sqrt(norm of x * x.adjoint())
This commit is contained in:
Benoit Jacob 2008-06-02 04:42:45 +00:00
parent 92b7e2d6a1
commit 0444e3601a
3 changed files with 53 additions and 3 deletions

View File

@ -202,6 +202,7 @@ template<typename Derived> class MatrixBase : public ArrayBase<Derived>
/** the return type of MatrixBase::adjoint() */ /** the return type of MatrixBase::adjoint() */
typedef Transpose<NestByValue<typename ei_unref<ConjugateReturnType>::type> > typedef Transpose<NestByValue<typename ei_unref<ConjugateReturnType>::type> >
AdjointReturnType; AdjointReturnType;
typedef Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1> EigenvaluesReturnType;
//@} //@}
/// \name Copying and initialization /// \name Copying and initialization
@ -605,6 +606,9 @@ template<typename Derived> class MatrixBase : public ArrayBase<Derived>
*/ */
//@{ //@{
const QR<typename ei_eval<Derived>::type> qr() const; const QR<typename ei_eval<Derived>::type> qr() const;
EigenvaluesReturnType eigenvalues() const;
RealScalar matrixNorm() const;
//@} //@}
}; };

View File

@ -269,7 +269,7 @@ template<typename Derived>
const Inverse<typename ei_eval<Derived>::type, true> const Inverse<typename ei_eval<Derived>::type, true>
MatrixBase<Derived>::inverse() const MatrixBase<Derived>::inverse() const
{ {
return Inverse<typename ei_eval<Derived>::type, true>(eval()); return Inverse<typename Derived::Eval, true>(eval());
} }
/** \return the matrix inverse of \c *this, which is assumed to exist. /** \return the matrix inverse of \c *this, which is assumed to exist.
@ -283,7 +283,7 @@ template<typename Derived>
const Inverse<typename ei_eval<Derived>::type, false> const Inverse<typename ei_eval<Derived>::type, false>
MatrixBase<Derived>::quickInverse() const MatrixBase<Derived>::quickInverse() const
{ {
return Inverse<typename ei_eval<Derived>::type, false>(eval()); return Inverse<typename Derived::Eval, false>(eval());
} }
#endif // EIGEN_INVERSE_H #endif // EIGEN_INVERSE_H

View File

@ -31,6 +31,8 @@
* *
* \param MatrixType the type of the matrix of which we are computing the eigen decomposition * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
* *
* \note MatrixType must be an actual Matrix type, it can't be an expression type.
*
* \sa MatrixBase::eigenvalues(), class EigenSolver * \sa MatrixBase::eigenvalues(), class EigenSolver
*/ */
template<typename _MatrixType> class SelfAdjointEigenSolver template<typename _MatrixType> class SelfAdjointEigenSolver
@ -58,7 +60,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
RealVectorType eigenvalues(void) const { return m_eivalues; } RealVectorType eigenvalues(void) const { return m_eivalues; }
protected: protected:
MatrixType m_eivec; MatrixType m_eivec;
RealVectorType m_eivalues; RealVectorType m_eivalues;
@ -169,4 +170,49 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix)
std::cout << "ei values = " << m_eivalues.transpose() << "\n\n"; std::cout << "ei values = " << m_eivalues.transpose() << "\n\n";
} }
template<typename Derived>
inline Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1>
MatrixBase<Derived>::eigenvalues() const
{
ei_assert(Flags&SelfAdjointBit);
return SelfAdjointEigenSolver<typename Derived::Eval>(eval()).eigenvalues();
}
template<typename Derived, bool IsSelfAdjoint>
struct ei_matrixNorm_selector
{
static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
matrixNorm(const MatrixBase<Derived>& m)
{
// FIXME if it is really guaranteed that the eigenvalues are already sorted,
// then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
return m.eigenvalues().cwiseAbs().maxCoeff();
}
};
template<typename Derived> struct ei_matrixNorm_selector<Derived, false>
{
static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
matrixNorm(const MatrixBase<Derived>& m)
{
// FIXME if it is really guaranteed that the eigenvalues are already sorted,
// then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
return ei_sqrt(
(m*m.adjoint())
.template marked<SelfAdjoint>()
.eigenvalues()
.cwiseAbs()
.maxCoeff()
);
}
};
template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
MatrixBase<Derived>::matrixNorm() const
{
return ei_matrixNorm_selector<Derived, Flags&SelfAdjointBit>
::matrixNorm(derived());
}
#endif // EIGEN_SELFADJOINTEIGENSOLVER_H #endif // EIGEN_SELFADJOINTEIGENSOLVER_H