Eigen libary for mbed
Diff: src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h
- Revision:
- 0:13a5d365ba16
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h Thu Oct 13 04:07:23 2016 +0000
@@ -0,0 +1,227 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
+#define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
+
+#include "./Tridiagonalization.h"
+
+namespace Eigen {
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+ *
+ *
+ * \class GeneralizedSelfAdjointEigenSolver
+ *
+ * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem
+ *
+ * \tparam _MatrixType the type of the matrix of which we are computing the
+ * eigendecomposition; this is expected to be an instantiation of the Matrix
+ * class template.
+ *
+ * This class solves the generalized eigenvalue problem
+ * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be
+ * selfadjoint and the matrix \f$ B \f$ should be positive definite.
+ *
+ * Only the \b lower \b triangular \b part of the input matrix is referenced.
+ *
+ * Call the function compute() to compute the eigenvalues and eigenvectors of
+ * a given matrix. Alternatively, you can use the
+ * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
+ * constructor which computes the eigenvalues and eigenvectors at construction time.
+ * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues()
+ * and eigenvectors() functions.
+ *
+ * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
+ * contains an example of the typical use of this class.
+ *
+ * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver
+ */
+template<typename _MatrixType>
+class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
+{
+ typedef SelfAdjointEigenSolver<_MatrixType> Base;
+ public:
+
+ typedef typename Base::Index Index;
+ typedef _MatrixType MatrixType;
+
+ /** \brief Default constructor for fixed-size matrices.
+ *
+ * The default constructor is useful in cases in which the user intends to
+ * perform decompositions via compute(). This constructor
+ * can only be used if \p _MatrixType is a fixed-size matrix; use
+ * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
+ */
+ GeneralizedSelfAdjointEigenSolver() : Base() {}
+
+ /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
+ *
+ * \param [in] size Positive integer, size of the matrix whose
+ * eigenvalues and eigenvectors will be computed.
+ *
+ * This constructor is useful for dynamic-size matrices, when the user
+ * intends to perform decompositions via compute(). The \p size
+ * parameter is only used as a hint. It is not an error to give a wrong
+ * \p size, but it may impair performance.
+ *
+ * \sa compute() for an example
+ */
+ GeneralizedSelfAdjointEigenSolver(Index size)
+ : Base(size)
+ {}
+
+ /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
+ *
+ * \param[in] matA Selfadjoint matrix in matrix pencil.
+ * Only the lower triangular part of the matrix is referenced.
+ * \param[in] matB Positive-definite matrix in matrix pencil.
+ * Only the lower triangular part of the matrix is referenced.
+ * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
+ * Default is #ComputeEigenvectors|#Ax_lBx.
+ *
+ * This constructor calls compute(const MatrixType&, const MatrixType&, int)
+ * to compute the eigenvalues and (if requested) the eigenvectors of the
+ * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the
+ * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix
+ * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property
+ * \f$ x^* B x = 1 \f$. The eigenvectors are computed if
+ * \a options contains ComputeEigenvectors.
+ *
+ * In addition, the two following variants can be solved via \p options:
+ * - \c ABx_lx: \f$ ABx = \lambda x \f$
+ * - \c BAx_lx: \f$ BAx = \lambda x \f$
+ *
+ * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
+ * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out
+ *
+ * \sa compute(const MatrixType&, const MatrixType&, int)
+ */
+ GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
+ int options = ComputeEigenvectors|Ax_lBx)
+ : Base(matA.cols())
+ {
+ compute(matA, matB, options);
+ }
+
+ /** \brief Computes generalized eigendecomposition of given matrix pencil.
+ *
+ * \param[in] matA Selfadjoint matrix in matrix pencil.
+ * Only the lower triangular part of the matrix is referenced.
+ * \param[in] matB Positive-definite matrix in matrix pencil.
+ * Only the lower triangular part of the matrix is referenced.
+ * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
+ * Default is #ComputeEigenvectors|#Ax_lBx.
+ *
+ * \returns Reference to \c *this
+ *
+ * Accoring to \p options, this function computes eigenvalues and (if requested)
+ * the eigenvectors of one of the following three generalized eigenproblems:
+ * - \c Ax_lBx: \f$ Ax = \lambda B x \f$
+ * - \c ABx_lx: \f$ ABx = \lambda x \f$
+ * - \c BAx_lx: \f$ BAx = \lambda x \f$
+ * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite
+ * matrix \f$ B \f$.
+ * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$.
+ *
+ * The eigenvalues() function can be used to retrieve
+ * the eigenvalues. If \p options contains ComputeEigenvectors, then the
+ * eigenvectors are also computed and can be retrieved by calling
+ * eigenvectors().
+ *
+ * The implementation uses LLT to compute the Cholesky decomposition
+ * \f$ B = LL^* \f$ and computes the classical eigendecomposition
+ * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
+ * and of \f$ L^{*} A L \f$ otherwise. This solves the
+ * generalized eigenproblem, because any solution of the generalized
+ * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution
+ * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the
+ * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements
+ * can be made for the two other variants.
+ *
+ * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
+ * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
+ *
+ * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
+ */
+ GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
+ int options = ComputeEigenvectors|Ax_lBx);
+
+ protected:
+
+};
+
+
+template<typename MatrixType>
+GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
+compute(const MatrixType& matA, const MatrixType& matB, int options)
+{
+ eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
+ eigen_assert((options&~(EigVecMask|GenEigMask))==0
+ && (options&EigVecMask)!=EigVecMask
+ && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
+ || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
+ && "invalid option parameter");
+
+ bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
+
+ // Compute the cholesky decomposition of matB = L L' = U'U
+ LLT<MatrixType> cholB(matB);
+
+ int type = (options&GenEigMask);
+ if(type==0)
+ type = Ax_lBx;
+
+ if(type==Ax_lBx)
+ {
+ // compute C = inv(L) A inv(L')
+ MatrixType matC = matA.template selfadjointView<Lower>();
+ cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
+ cholB.matrixU().template solveInPlace<OnTheRight>(matC);
+
+ Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
+
+ // transform back the eigen vectors: evecs = inv(U) * evecs
+ if(computeEigVecs)
+ cholB.matrixU().solveInPlace(Base::m_eivec);
+ }
+ else if(type==ABx_lx)
+ {
+ // compute C = L' A L
+ MatrixType matC = matA.template selfadjointView<Lower>();
+ matC = matC * cholB.matrixL();
+ matC = cholB.matrixU() * matC;
+
+ Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
+
+ // transform back the eigen vectors: evecs = inv(U) * evecs
+ if(computeEigVecs)
+ cholB.matrixU().solveInPlace(Base::m_eivec);
+ }
+ else if(type==BAx_lx)
+ {
+ // compute C = L' A L
+ MatrixType matC = matA.template selfadjointView<Lower>();
+ matC = matC * cholB.matrixL();
+ matC = cholB.matrixU() * matC;
+
+ Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
+
+ // transform back the eigen vectors: evecs = L * evecs
+ if(computeEigVecs)
+ Base::m_eivec = cholB.matrixL() * Base::m_eivec;
+ }
+
+ return *this;
+}
+
+} // end namespace Eigen
+
+#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
\ No newline at end of file