QRFactorization.hpp

Go to the documentation of this file.

#pragma once
 
#include <Eigen/Core>
#include <fstream>
#include <limits>
#include <string>
#include <vector>
#include "logging/Logger.hpp"
#include "mesh/SharedPointer.hpp"
 
namespace precice {
namespace acceleration {
namespace impl {
 
class QRFactorization {
public:
  QRFactorization(
      int    filter = 0,
      double omega  = 0,
      double theta  = 1. / 0.7,
      double sigma  = std::numeric_limits<double>::min());
 
  QRFactorization(
      Eigen::MatrixXd A,
      int             filter,
      double          omega = 0,
      double          theta = 1. / 0.7,
      double          sigma = std::numeric_limits<double>::min());
 
  QRFactorization(
      Eigen::MatrixXd Q,
      Eigen::MatrixXd R,
      int             rows,
      int             cols,
      int             filter,
      double          omega = 0,
      double          theta = 1. / 0.7,
      double          sigma = std::numeric_limits<double>::min());
 
  virtual ~QRFactorization() {}
 
  void reset();
 
  void reset(
      Eigen::MatrixXd const &Q,
      Eigen::MatrixXd const &R,
      int                    rows,
      int                    cols,
      double                 omega = 0,
      double                 theta = 1. / 0.7,
      double                 sigma = std::numeric_limits<double>::min());
 
  void reset(
      Eigen::MatrixXd const &A,
      int                    globalRows,
      double                 omega = 0,
      double                 theta = 1. / 0.7,
      double                 sigma = std::numeric_limits<double>::min());
 
  bool insertColumn(int k, const Eigen::VectorXd &v, double singularityLimit = 0);
 
  void deleteColumn(int k);
 
  void pushFront(const Eigen::VectorXd &v);
 
  void pushBack(const Eigen::VectorXd &v);
 
  void popFront();
 
  void popBack();
 
  void applyFilter(double singularityLimit, std::vector<int> &delIndices, Eigen::MatrixXd &V);
 
  Eigen::MatrixXd &matrixQ();
 
  Eigen::MatrixXd &matrixR();
 
  // @brief returns the number of columns in the QR-decomposition
  int cols() const;
  // @brief returns the number of rows in the QR-decomposition
  int rows() const;
 
  // @brief optional file-stream for logging output
  void setfstream(std::fstream *stream);
 
  // @brief set number of global rows
  void setGlobalRows(int gr);
 
  // @brief sets the filtering technique to maintain good conditioning of the least squares system
  void setFilter(int filter);
 
private:
  struct givensRot {
    int    i, j;
    double sigma, gamma;
  };
 
  int orthogonalize_stable(Eigen::VectorXd &v, Eigen::VectorXd &r, double &rho, int colNum);
 
  int orthogonalize(Eigen::VectorXd &v, Eigen::VectorXd &r, double &rho, int colNum);
 
  void computeReflector(givensRot &grot, double &x, double &y);
 
  void applyReflector(const givensRot &grot, int k, int l, Eigen::VectorXd &p, Eigen::VectorXd &q);
 
  logging::Logger _log{"acceleration::QRFactorization"};
 
  Eigen::MatrixXd _Q;
  Eigen::MatrixXd _R;
 
  int _rows;
  int _cols;
 
  int    _filter;
  double _omega;
  double _theta;
  double _sigma;
 
  // @brief optional infostream that writes information to file
  std::fstream *_infostream;
  bool          _fstream_set;
 
  int _globalRows;
};
 
} // namespace impl
} // namespace acceleration
} // namespace precice