cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

MatrixBaseEigenvalues.h (5575B)


      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
      5 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
      6 //
      7 // This Source Code Form is subject to the terms of the Mozilla
      8 // Public License v. 2.0. If a copy of the MPL was not distributed
      9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
     10 
     11 #ifndef EIGEN_MATRIXBASEEIGENVALUES_H
     12 #define EIGEN_MATRIXBASEEIGENVALUES_H
     13 
     14 namespace Eigen { 
     15 
     16 namespace internal {
     17 
     18 template<typename Derived, bool IsComplex>
     19 struct eigenvalues_selector
     20 {
     21   // this is the implementation for the case IsComplex = true
     22   static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
     23   run(const MatrixBase<Derived>& m)
     24   {
     25     typedef typename Derived::PlainObject PlainObject;
     26     PlainObject m_eval(m);
     27     return ComplexEigenSolver<PlainObject>(m_eval, false).eigenvalues();
     28   }
     29 };
     30 
     31 template<typename Derived>
     32 struct eigenvalues_selector<Derived, false>
     33 {
     34   static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
     35   run(const MatrixBase<Derived>& m)
     36   {
     37     typedef typename Derived::PlainObject PlainObject;
     38     PlainObject m_eval(m);
     39     return EigenSolver<PlainObject>(m_eval, false).eigenvalues();
     40   }
     41 };
     42 
     43 } // end namespace internal
     44 
     45 /** \brief Computes the eigenvalues of a matrix 
     46   * \returns Column vector containing the eigenvalues.
     47   *
     48   * \eigenvalues_module
     49   * This function computes the eigenvalues with the help of the EigenSolver
     50   * class (for real matrices) or the ComplexEigenSolver class (for complex
     51   * matrices). 
     52   *
     53   * The eigenvalues are repeated according to their algebraic multiplicity,
     54   * so there are as many eigenvalues as rows in the matrix.
     55   *
     56   * The SelfAdjointView class provides a better algorithm for selfadjoint
     57   * matrices.
     58   *
     59   * Example: \include MatrixBase_eigenvalues.cpp
     60   * Output: \verbinclude MatrixBase_eigenvalues.out
     61   *
     62   * \sa EigenSolver::eigenvalues(), ComplexEigenSolver::eigenvalues(),
     63   *     SelfAdjointView::eigenvalues()
     64   */
     65 template<typename Derived>
     66 inline typename MatrixBase<Derived>::EigenvaluesReturnType
     67 MatrixBase<Derived>::eigenvalues() const
     68 {
     69   return internal::eigenvalues_selector<Derived, NumTraits<Scalar>::IsComplex>::run(derived());
     70 }
     71 
     72 /** \brief Computes the eigenvalues of a matrix
     73   * \returns Column vector containing the eigenvalues.
     74   *
     75   * \eigenvalues_module
     76   * This function computes the eigenvalues with the help of the
     77   * SelfAdjointEigenSolver class.  The eigenvalues are repeated according to
     78   * their algebraic multiplicity, so there are as many eigenvalues as rows in
     79   * the matrix.
     80   *
     81   * Example: \include SelfAdjointView_eigenvalues.cpp
     82   * Output: \verbinclude SelfAdjointView_eigenvalues.out
     83   *
     84   * \sa SelfAdjointEigenSolver::eigenvalues(), MatrixBase::eigenvalues()
     85   */
     86 template<typename MatrixType, unsigned int UpLo> 
     87 EIGEN_DEVICE_FUNC inline typename SelfAdjointView<MatrixType, UpLo>::EigenvaluesReturnType
     88 SelfAdjointView<MatrixType, UpLo>::eigenvalues() const
     89 {
     90   PlainObject thisAsMatrix(*this);
     91   return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix, false).eigenvalues();
     92 }
     93 
     94 
     95 
     96 /** \brief Computes the L2 operator norm
     97   * \returns Operator norm of the matrix.
     98   *
     99   * \eigenvalues_module
    100   * This function computes the L2 operator norm of a matrix, which is also
    101   * known as the spectral norm. The norm of a matrix \f$ A \f$ is defined to be
    102   * \f[ \|A\|_2 = \max_x \frac{\|Ax\|_2}{\|x\|_2} \f]
    103   * where the maximum is over all vectors and the norm on the right is the
    104   * Euclidean vector norm. The norm equals the largest singular value, which is
    105   * the square root of the largest eigenvalue of the positive semi-definite
    106   * matrix \f$ A^*A \f$.
    107   *
    108   * The current implementation uses the eigenvalues of \f$ A^*A \f$, as computed
    109   * by SelfAdjointView::eigenvalues(), to compute the operator norm of a
    110   * matrix.  The SelfAdjointView class provides a better algorithm for
    111   * selfadjoint matrices.
    112   *
    113   * Example: \include MatrixBase_operatorNorm.cpp
    114   * Output: \verbinclude MatrixBase_operatorNorm.out
    115   *
    116   * \sa SelfAdjointView::eigenvalues(), SelfAdjointView::operatorNorm()
    117   */
    118 template<typename Derived>
    119 inline typename MatrixBase<Derived>::RealScalar
    120 MatrixBase<Derived>::operatorNorm() const
    121 {
    122   using std::sqrt;
    123   typename Derived::PlainObject m_eval(derived());
    124   // FIXME if it is really guaranteed that the eigenvalues are already sorted,
    125   // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
    126   return sqrt((m_eval*m_eval.adjoint())
    127                  .eval()
    128 		 .template selfadjointView<Lower>()
    129 		 .eigenvalues()
    130 		 .maxCoeff()
    131 		 );
    132 }
    133 
    134 /** \brief Computes the L2 operator norm
    135   * \returns Operator norm of the matrix.
    136   *
    137   * \eigenvalues_module
    138   * This function computes the L2 operator norm of a self-adjoint matrix. For a
    139   * self-adjoint matrix, the operator norm is the largest eigenvalue.
    140   *
    141   * The current implementation uses the eigenvalues of the matrix, as computed
    142   * by eigenvalues(), to compute the operator norm of the matrix.
    143   *
    144   * Example: \include SelfAdjointView_operatorNorm.cpp
    145   * Output: \verbinclude SelfAdjointView_operatorNorm.out
    146   *
    147   * \sa eigenvalues(), MatrixBase::operatorNorm()
    148   */
    149 template<typename MatrixType, unsigned int UpLo>
    150 EIGEN_DEVICE_FUNC inline typename SelfAdjointView<MatrixType, UpLo>::RealScalar
    151 SelfAdjointView<MatrixType, UpLo>::operatorNorm() const
    152 {
    153   return eigenvalues().cwiseAbs().maxCoeff();
    154 }
    155 
    156 } // end namespace Eigen
    157 
    158 #endif