Preconditioner.hpp

Go to the documentation of this file.

#pragma once
 
#include <Eigen/Core>
#include <numeric>
#include <vector>
 
#include "cplscheme/SharedPointer.hpp"
#include "logging/LogMacros.hpp"
#include "logging/Logger.hpp"
#include "utils/assertion.hpp"
 
namespace precice {
namespace acceleration {
namespace impl {
 
class Preconditioner {
public:
  Preconditioner(int maxNonConstTimeWindows)
      : _maxNonConstTimeWindows(maxNonConstTimeWindows)
  {
  }
 
  virtual ~Preconditioner() {}
 
  virtual void initialize(std::vector<size_t> &svs)
  {
    PRECICE_TRACE();
 
    PRECICE_ASSERT(_weights.empty());
    _subVectorSizes = svs;
 
    size_t N = std::accumulate(_subVectorSizes.begin(), _subVectorSizes.end(), static_cast<std::size_t>(0));
 
    // cannot do this already in the constructor as the size is unknown at that point
    _weights.resize(N, 1.0);
    _invWeights.resize(N, 1.0);
  }
 
  void apply(Eigen::MatrixXd &M, bool transpose)
  {
    PRECICE_TRACE();
    if (transpose) {
      PRECICE_ASSERT(M.cols() == (int) _weights.size(), M.cols(), _weights.size());
      for (int i = 0; i < M.cols(); i++) {
        for (int j = 0; j < M.rows(); j++) {
          M(j, i) *= _weights[i];
        }
      }
    } else {
      PRECICE_ASSERT(M.rows() == (int) _weights.size(), M.rows(), (int) _weights.size());
      for (int i = 0; i < M.cols(); i++) {
        for (int j = 0; j < M.rows(); j++) {
          M(j, i) *= _weights[j];
        }
      }
    }
  }
 
  void revert(Eigen::MatrixXd &M, bool transpose)
  {
    PRECICE_TRACE();
    //PRECICE_ASSERT(_needsGlobalWeights);
    if (transpose) {
      PRECICE_ASSERT(M.cols() == (int) _invWeights.size());
      for (int i = 0; i < M.cols(); i++) {
        for (int j = 0; j < M.rows(); j++) {
          M(j, i) *= _invWeights[i];
        }
      }
    } else {
      PRECICE_ASSERT(M.rows() == (int) _invWeights.size(), M.rows(), (int) _invWeights.size());
      for (int i = 0; i < M.cols(); i++) {
        for (int j = 0; j < M.rows(); j++) {
          M(j, i) *= _invWeights[j];
        }
      }
    }
  }
 
  void apply(Eigen::MatrixXd &M)
  {
    PRECICE_TRACE();
    PRECICE_ASSERT(M.rows() == (int) _weights.size(), M.rows(), (int) _weights.size());
 
    // scale matrix M
    for (int i = 0; i < M.cols(); i++) {
      for (int j = 0; j < M.rows(); j++) {
        M(j, i) *= _weights[j];
      }
    }
  }
 
  void apply(Eigen::VectorXd &v)
  {
    PRECICE_TRACE();
 
    PRECICE_ASSERT(v.size() == (int) _weights.size());
 
    // scale residual
    for (int j = 0; j < v.size(); j++) {
      v[j] *= _weights[j];
    }
  }
 
  void revert(Eigen::MatrixXd &M)
  {
    PRECICE_TRACE();
 
    PRECICE_ASSERT(M.rows() == (int) _weights.size());
 
    // scale matrix M
    for (int i = 0; i < M.cols(); i++) {
      for (int j = 0; j < M.rows(); j++) {
        M(j, i) *= _invWeights[j];
      }
    }
  }
 
  void revert(Eigen::VectorXd &v)
  {
    PRECICE_TRACE();
 
    PRECICE_ASSERT(v.size() == (int) _weights.size());
 
    // scale residual
    for (int j = 0; j < v.size(); j++) {
      v[j] *= _invWeights[j];
    }
  }
 
  void update(bool timeWindowComplete, const Eigen::VectorXd &oldValues, const Eigen::VectorXd &res)
  {
    PRECICE_TRACE(_nbNonConstTimeWindows, _frozen);
 
    // if number of allowed non-const time windows is exceeded, do not update weights
    if (_frozen)
      return;
 
    // increment number of time windows that has been scaled with changing preconditioning weights
    if (timeWindowComplete) {
      _nbNonConstTimeWindows++;
      if (_nbNonConstTimeWindows >= _maxNonConstTimeWindows && _maxNonConstTimeWindows > 0)
        _frozen = true;
    }
 
    // type specific update functionality
    _update_(timeWindowComplete, oldValues, res);
  }
 
  bool requireNewQR()
  {
    PRECICE_TRACE(_requireNewQR);
    return _requireNewQR;
  }
 
  void newQRfulfilled()
  {
    _requireNewQR = false;
  }
 
  std::vector<double> &getWeights()
  {
    return _weights;
  }
 
  bool isConst()
  {
    return _frozen;
  }
 
protected:
  std::vector<double> _weights;
 
  std::vector<double> _invWeights;
 
  std::vector<size_t> _subVectorSizes;
 
  int _maxNonConstTimeWindows;
 
  int _nbNonConstTimeWindows = 0;
 
  bool _requireNewQR = false;
 
  bool _frozen = false;
 
  virtual void _update_(bool timeWindowComplete, const Eigen::VectorXd &oldValues, const Eigen::VectorXd &res) = 0;
 
private:
  logging::Logger _log{"acceleration::Preconditioner"};
};
 
} // namespace impl
} // namespace acceleration
} // namespace precice