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::acceleration::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() = default;

  void resetFilter(double singularityLimit, std::vector<int> &delIndices, Eigen::MatrixXd &V);

  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());

  [[nodiscard]] std::vector<int> 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);
 
  // @brief returns the number of times the QR2 filter step was performed
  size_t getResetFilterCounter() const;
 
  // @brief marks that the prescaling coefficients aren't constant and the QR3 filter needs to fallback to QR2
  void requireQR3Fallback();
 
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;
 
  bool   _computeQR2Filter   = true; // flag to indicate if the QR2 filter step should be performed
  size_t _resetFilterCounter = 0;    // counter for the number of times the QR2 filter step was performed
 
  // @brief optional infostream that writes information to file
  std::fstream *_infostream;
  bool          _fstream_set;
 
  int _globalRows;
};
 
} // namespace precice::acceleration::impl