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::acceleration::impl {

class Preconditioner {
public:
  Preconditioner(int maxNonConstTimeWindows)
      : _maxNonConstTimeWindows(maxNonConstTimeWindows)
  {
  }

  virtual ~Preconditioner() = default;

  virtual void initialize(std::vector<size_t> &svs)
  {
    PRECICE_TRACE();
 
    _subVectorSizes = svs;
 
    // Compute offsets of each subvector
    _subVectorOffsets.resize(_subVectorSizes.size(), 0);
    std::partial_sum(_subVectorSizes.begin(), --_subVectorSizes.end(), ++_subVectorOffsets.begin());
 
    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_DEBUG_IF((int) _weights.size() != M.cols(), "The number of columns of the matrix {} and weights size {} mismatched.", M.rows(), _weights.size());
 
      int validCols = std::min(static_cast<int>(M.cols()), (int) _weights.size());
      for (int i = 0; i < validCols; i++) {
        for (int j = 0; j < M.rows(); j++) {
          M(j, i) *= _weights[i];
        }
      }
    } else {
      PRECICE_DEBUG_IF((int) _weights.size() != M.rows(), "The number of rows of the matrix {} and weights size {} mismatched.", M.rows(), _weights.size());
 
      int validRows = std::min(static_cast<int>(M.rows()), (int) _weights.size());
      for (int i = 0; i < M.cols(); i++) {
        for (int j = 0; j < validRows; j++) {
          M(j, i) *= _weights[j];
        }
      }
    }
  }

  void revert(Eigen::MatrixXd &M, bool transpose)
  {
    PRECICE_TRACE();
    // PRECICE_ASSERT(_needsGlobalWeights);
    if (transpose) {
      PRECICE_DEBUG_IF((int) _weights.size() != M.cols(), "The number of columns of the matrix {} and weights size {} mismatched.", M.cols(), _weights.size());
 
      int validCols = std::min(static_cast<int>(M.cols()), (int) _weights.size());
      for (int i = 0; i < validCols; i++) {
        for (int j = 0; j < M.rows(); j++) {
          M(j, i) *= _invWeights[i];
        }
      }
    } else {
      PRECICE_DEBUG_IF((int) _weights.size() != M.rows(), "The number of rows of the matrix {} and weights size {} mismatched.", M.rows(), _weights.size());
 
      int validRows = std::min(static_cast<int>(M.rows()), (int) _weights.size());
      for (int i = 0; i < M.cols(); i++) {
        for (int j = 0; j < validRows; j++) {
          M(j, i) *= _invWeights[j];
        }
      }
    }
  }

  void apply(Eigen::MatrixXd &M)
  {
    PRECICE_TRACE();
    PRECICE_DEBUG_IF((int) _weights.size() != M.rows(), "The number of rows of the matrix {} and weights size {} mismatched.", M.rows(), _weights.size());
 
    // scale matrix M
    int validRows = std::min(static_cast<int>(M.rows()), (int) _weights.size());
    for (int i = 0; i < M.cols(); i++) {
      for (int j = 0; j < validRows; j++) {
        M(j, i) *= _weights[j];
      }
    }
  }

  void apply(Eigen::VectorXd &v)
  {
    PRECICE_TRACE();
    PRECICE_DEBUG_IF((int) _weights.size() != v.size(), "The vector size {} and weights size {} mismatched.", v.size(), _weights.size());
 
    // scale vector
    int validSize = std::min(static_cast<int>(v.size()), (int) _weights.size());
    for (int j = 0; j < validSize; j++) {
      v[j] *= _weights[j];
    }
  }

  void revert(Eigen::MatrixXd &M)
  {
    PRECICE_TRACE();
    PRECICE_DEBUG_IF((int) _weights.size() != M.rows(), "The number of rows of the matrix {} and weights size {} mismatched.", M.rows(), _weights.size());
 
    // scale matrix M
    int validRows = std::min(static_cast<int>(M.rows()), (int) _weights.size());
    for (int i = 0; i < M.cols(); i++) {
      for (int j = 0; j < validRows; j++) {
        M(j, i) *= _invWeights[j];
      }
    }
  }

  void revert(Eigen::VectorXd &v)
  {
    PRECICE_TRACE();
    PRECICE_DEBUG_IF((int) _weights.size() != v.size(), "The vector size {} and weights size {} mismatched.", v.size(), _weights.size());
 
    // revert vector scaling
    int validSize = std::min(static_cast<int>(v.size()), (int) _weights.size());
    for (int j = 0; j < validSize; 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;

  std::vector<size_t> _subVectorOffsets;

  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 precice::acceleration::impl