SymGEigsSolver.h

Go to the documentation of this file.

// Copyright (C) 2016-2019 Yixuan Qiu <yixuan.qiu@cos.name>
//
// 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 https://mozilla.org/MPL/2.0/.
 
#ifndef SYM_GEIGS_SOLVER_H
#define SYM_GEIGS_SOLVER_H
 
#include "SymEigsBase.h"
#include "Util/GEigsMode.h"
#include "MatOp/internal/SymGEigsCholeskyOp.h"
#include "MatOp/internal/SymGEigsRegInvOp.h"
 
namespace Spectra {
 
 
 
// Empty class template
template <typename Scalar,
          int SelectionRule,
          typename OpType,
          typename BOpType,
          int GEigsMode>
class SymGEigsSolver
{};
 
 
// Partial specialization for GEigsMode = GEIGS_CHOLESKY
template <typename Scalar,
          int SelectionRule,
          typename OpType,
          typename BOpType>
class SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, GEIGS_CHOLESKY> :
    public SymEigsBase<Scalar, SelectionRule, SymGEigsCholeskyOp<Scalar, OpType, BOpType>, IdentityBOp>
{
private:
    typedef Eigen::Index Index;
    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
 
    BOpType* m_Bop;
 
public:
    SymGEigsSolver(OpType* op, BOpType* Bop, Index nev, Index ncv) :
        SymEigsBase<Scalar, SelectionRule, SymGEigsCholeskyOp<Scalar, OpType, BOpType>, IdentityBOp>(
            new SymGEigsCholeskyOp<Scalar, OpType, BOpType>(*op, *Bop), NULL, nev, ncv),
        m_Bop(Bop)
    {}
 
 
    ~SymGEigsSolver()
    {
        // m_op contains the constructed SymGEigsCholeskyOp object
        delete this->m_op;
    }
 
    Matrix eigenvectors(Index nvec) const
    {
        Matrix res = SymEigsBase<Scalar, SelectionRule, SymGEigsCholeskyOp<Scalar, OpType, BOpType>, IdentityBOp>::eigenvectors(nvec);
        Vector tmp(res.rows());
        const Index nconv = res.cols();
        for (Index i = 0; i < nconv; i++)
        {
            m_Bop->upper_triangular_solve(&res(0, i), tmp.data());
            res.col(i).noalias() = tmp;
        }
 
        return res;
    }
 
    Matrix eigenvectors() const
    {
        return SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, GEIGS_CHOLESKY>::eigenvectors(this->m_nev);
    }
 
};
 
 
// Partial specialization for GEigsMode = GEIGS_REGULAR_INVERSE
template <typename Scalar,
          int SelectionRule,
          typename OpType,
          typename BOpType>
class SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, GEIGS_REGULAR_INVERSE> :
    public SymEigsBase<Scalar, SelectionRule, SymGEigsRegInvOp<Scalar, OpType, BOpType>, BOpType>
{
private:
    typedef Eigen::Index Index;
 
public:
    SymGEigsSolver(OpType* op, BOpType* Bop, Index nev, Index ncv) :
        SymEigsBase<Scalar, SelectionRule, SymGEigsRegInvOp<Scalar, OpType, BOpType>, BOpType>(
            new SymGEigsRegInvOp<Scalar, OpType, BOpType>(*op, *Bop), Bop, nev, ncv)
    {}
 
    ~SymGEigsSolver()
    {
        // m_op contains the constructed SymGEigsRegInvOp object
        delete this->m_op;
    }
};
 
}  // namespace Spectra
 
#endif  // SYM_GEIGS_SOLVER_H