LOBPCG.h

Go to the documentation of this file.

#ifndef LOBPCG_H
#define LOBPCG_H
 
 
#include "MUQ/Modeling/LinearAlgebra/LinearOperator.h"
#include "MUQ/Modeling/LinearAlgebra/GeneralizedEigenSolver.h"
 
#include <boost/property_tree/ptree_fwd.hpp>
#include <Eigen/Core>
 
namespace muq{
namespace Modeling{
 
  class LOBPCG : public GeneralizedEigenSolver
  {
  public:
 
    LOBPCG(int    numEigsIn,
           double eigRelTolIn=0.0,
           double eigAbsTolIn=0.0,
           int    blockSizeIn=-1,
           double solverTolIn=-1,
           int    maxItsIn=-1,
           int    verbosityIn=0);
 
    LOBPCG(boost::property_tree::ptree const& options);
 
    virtual ~LOBPCG() = default;
 
    LOBPCG& compute(std::shared_ptr<LinearOperator> const& A,
                    std::shared_ptr<LinearOperator>        B=nullptr,
                    std::shared_ptr<LinearOperator>        M=nullptr){return compute(A,Eigen::MatrixXd(),B,M);};
 
    LOBPCG& compute(std::shared_ptr<LinearOperator> const& A,
                    Eigen::MatrixXd                 const& constMat,
                    std::shared_ptr<LinearOperator>        B = nullptr,
                    std::shared_ptr<LinearOperator>        M = nullptr);
 
 
    void InitializeVectors(Eigen::MatrixXd const& vecs);
 
    LOBPCG& reset(int dim);
 
  private:
 
    std::pair<Eigen::VectorXd, Eigen::MatrixXd> ComputeBlock(std::shared_ptr<LinearOperator>   const& A,
                                                             Eigen::Ref<const Eigen::MatrixXd> const& X0,
                                                             Eigen::Ref<const Eigen::MatrixXd> const& constMat,
                                                             std::shared_ptr<LinearOperator>          B,
                                                             std::shared_ptr<LinearOperator>          M);
 
    class Orthonormalizer{
    public:
      Orthonormalizer(std::shared_ptr<LinearOperator> const& Bin) : B(Bin){};
 
      bool ComputeInPlace(Eigen::Ref<Eigen::MatrixXd> V);
      bool ComputeInPlace(Eigen::Ref<Eigen::MatrixXd> V, Eigen::Ref<const Eigen::MatrixXd> const& BVin);
 
      Eigen::MatrixXd Compute(Eigen::Ref<const Eigen::MatrixXd> const& V);
      Eigen::MatrixXd Compute(Eigen::Ref<const Eigen::MatrixXd> const& V, Eigen::Ref<const Eigen::MatrixXd> const& BVin);
 
      Eigen::MatrixXd InverseVBV() const;
      Eigen::MatrixXd const& GetBV() const{return BV;};
      Eigen::MatrixXd& GetBV(){return BV;};
 
      int vDim;
      std::shared_ptr<LinearOperator> B;
      Eigen::MatrixXd BV;
      Eigen::MatrixXd VBV_Chol;
    };
 
    class Constraints{
    public:
      Constraints(std::shared_ptr<LinearOperator>   const& B,
                  Eigen::Ref<const Eigen::MatrixXd> const& constVec);
 
      void ApplyInPlace(Eigen::Ref<Eigen::MatrixXd> x);
 
      int size() const{return BY.cols();};
 
    private:
      Eigen::MatrixXd BY;
      Eigen::Ref<const Eigen::MatrixXd> const& Y;
      Eigen::LLT<Eigen::MatrixXd> YBY_llt;
 
    }; // class LOBPCG::Constraints
 
    int numEigs;
 
    int blockSize;
 
    double solverTol;
 
    double eigRelTol, eigAbsTol;
 
    int maxIts;
 
    int verbosity;
 
    //Eigen::VectorXd eigVals;
    //Eigen::MatrixXd eigVecs;
 
    //std::shared_ptr<LinearOperator> A, B, M;
    double Anorm, Bnorm;
 
  }; // class LOBPCG
}
}
 
 
 
#endif // LOBPCG_H