GeneralizedSelfAdjointEigenSolver.h

// 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/.
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
#define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
 
#include "./Tridiagonalization.h"
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {

template <typename MatrixType_>
class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<MatrixType_> {
  using Base = SelfAdjointEigenSolver<MatrixType_>;
 
 public:
  using MatrixType = MatrixType_;

  GeneralizedSelfAdjointEigenSolver() : Base(), m_cholB(), m_matC() {}

  explicit GeneralizedSelfAdjointEigenSolver(Index size) : Base(size), m_cholB(size), m_matC(size, size) {}

  template <typename InputTypeA, typename InputTypeB>
  GeneralizedSelfAdjointEigenSolver(const EigenBase<InputTypeA>& matA, const EigenBase<InputTypeB>& matB,
                                    int options = ComputeEigenvectors | Ax_lBx)
      : Base(matA.cols()) {
    compute(matA.derived(), matB.derived(), options);
  }

  template <typename InputTypeA, typename InputTypeB, bool IsRef = internal::is_ref<MatrixType>::value,
            std::enable_if_t<IsRef, int> = 0>
  GeneralizedSelfAdjointEigenSolver(EigenBase<InputTypeA>& matA, EigenBase<InputTypeB>& matB,
                                    int options = ComputeEigenvectors | Ax_lBx)
      : Base(matA, typename Base::BindStorageTag()), m_cholB(matB.derived()), m_matC(matA.derived()) {
    // Complete the upper triangle from the referenced lower one.
    m_matC = m_matC.template selfadjointView<Lower>();
    computeInPlace(options);
  }

  template <typename InputTypeA, typename InputTypeB>
  GeneralizedSelfAdjointEigenSolver& compute(const EigenBase<InputTypeA>& matA, const EigenBase<InputTypeB>& matB,
                                             int options = ComputeEigenvectors | Ax_lBx);
 
 protected:
  // Reused across compute() calls so that a solver constructed with its size,
  // or computed with once, does not allocate again.
  LLT<MatrixType> m_cholB;
  MatrixType m_matC;
 
 private:
  GeneralizedSelfAdjointEigenSolver& computeInPlace(int options);
};
 
template <typename MatrixType>
template <typename InputTypeA, typename InputTypeB>
GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::compute(
    const EigenBase<InputTypeA>& matA, const EigenBase<InputTypeB>& matB, int options) {
  eigen_assert(matA.cols() == matA.rows() && matB.rows() == matA.rows() && matB.cols() == matB.rows());
 
  // Compute the cholesky decomposition of matB = L L' = U'U
  m_cholB.compute(matB.derived());
  m_matC = matA.derived().template selfadjointView<Lower>();
  return computeInPlace(options);
}

template <typename MatrixType>
GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::computeInPlace(
    int options) {
  eigen_assert(m_matC.cols() == m_matC.rows() && m_cholB.rows() == m_matC.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);
 
  int type = (options & GenEigMask);
  if (type == 0) type = Ax_lBx;
 
  // In the inplace decomposition the eigenvector storage is the matrix A itself, i.e. m_matC: the products of the
  // ABx_lx and BAx_lx forms, which cannot run in place, then go through a temporary instead.
  const bool eivecIsMatC = internal::is_same_dense(Base::m_eivec, m_matC);
 
  if (type == Ax_lBx) {
    // compute C = inv(L) A inv(L')
    m_cholB.matrixL().template solveInPlace<OnTheLeft>(m_matC);
    m_cholB.matrixU().template solveInPlace<OnTheRight>(m_matC);
 
    Base::compute(m_matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
 
    // transform back the eigen vectors: evecs = inv(U) * evecs
    if (computeEigVecs) m_cholB.matrixU().solveInPlace(Base::m_eivec);
  } else if (type == ABx_lx) {
    // compute C = L' A L, using Base::m_eivec for the intermediate product: Base::compute()
    // overwrites it before reading it.
    if (eivecIsMatC) {
      typename Base::PlainMatrixType tmp = m_matC * m_cholB.matrixL();
      m_matC.noalias() = m_cholB.matrixU() * tmp;
    } else {
      Base::m_eivec.noalias() = m_matC * m_cholB.matrixL();
      m_matC.noalias() = m_cholB.matrixU() * Base::m_eivec;
    }
 
    Base::compute(m_matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
 
    // transform back the eigen vectors: evecs = inv(U) * evecs
    if (computeEigVecs) m_cholB.matrixU().solveInPlace(Base::m_eivec);
  } else if (type == BAx_lx) {
    // compute C = L' A L
    if (eivecIsMatC) {
      typename Base::PlainMatrixType tmp = m_matC * m_cholB.matrixL();
      m_matC.noalias() = m_cholB.matrixU() * tmp;
    } else {
      Base::m_eivec.noalias() = m_matC * m_cholB.matrixL();
      m_matC.noalias() = m_cholB.matrixU() * Base::m_eivec;
    }
 
    Base::compute(m_matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
 
    // transform back the eigen vectors: evecs = L * evecs, using m_matC as the
    // intermediate: Base::compute() has consumed it.
    if (computeEigVecs) {
      if (eivecIsMatC) {
        Base::m_eivec = m_cholB.matrixL() * Base::m_eivec;
      } else {
        m_matC.noalias() = m_cholB.matrixL() * Base::m_eivec;
        Base::m_eivec = m_matC;
      }
    }
  }
 
  // The reduction above is only valid for a positive definite B.
  if (m_cholB.info() != Success) Base::m_info = NumericalIssue;
 
  return *this;
}
 
}  // end namespace Eigen
 
#endif  // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H