// 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>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.

#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