SVDFactorization.hpp

Go to the documentation of this file.

#pragma once
/*
 * SVDFactorization.hpp
 *
 *  Created on: Feb 5, 2016
 *      Author: Klaudius Scheufele
 */
 
#ifndef PRECICE_NO_MPI
 
#include <Eigen/Core>
#include <Eigen/Dense>
#include <string>
 
#include "acceleration/impl/ParallelMatrixOperations.hpp"
#include "acceleration/impl/Preconditioner.hpp"
#include "acceleration/impl/SharedPointer.hpp"
#include "logging/Logger.hpp"
#include "precice/impl/Types.hpp"
 
// ------- CLASS DEFINITION
 
namespace precice::acceleration::impl {

class SVDFactorization {
public:
  // Eigen
  typedef Eigen::MatrixXd Matrix;
  typedef Eigen::VectorXd Vector;

  SVDFactorization(
      double            eps,
      PtrPreconditioner preconditioner);

  virtual ~SVDFactorization() = default;

  template <typename Derived1, typename Derived2>
  void update(
      const Eigen::MatrixBase<Derived1> &A,
      const Eigen::MatrixBase<Derived2> &B)
  {
    PRECICE_TRACE();
    PRECICE_ASSERT(_initialized);
    if (_initialSVD) {
      PRECICE_ASSERT(A.rows() == _rowsA, A.rows(), _rowsA);
      PRECICE_ASSERT(B.rows() == _rowsB, B.rows(), _rowsB);
    } else {
      PRECICE_ASSERT(A.cols() == B.cols(), A.cols(), B.cols());
      _rowsA = A.rows();
      _rowsB = B.rows();
      _cols  = 0;
      _psi   = Matrix::Zero(_rowsA, 0);
      _phi   = Matrix::Zero(_rowsB, 0);
      _sigma = Vector::Zero(0);
    }

    Matrix Atil(_psi.cols(), A.cols()); // Atil is of size (K_bar x m)
 
    // Atil := \psi^T *A
    // local computation of \psi^T * A and allreduce_sum to Atil (global), stored local on each proc
    _parMatrixOps->multiply(_psi.transpose(), A, Atil, (int) _psi.cols(), _globalRowsA, (int) A.cols());
 
    // Ptil := (I-\psi\psi^T)A
    // Atil is local on each proc, thus fully local computation, embarrassingly parallel
    Matrix Ptil = A - _psi * Atil;
 
    // compute orthogonal basis P of Ptil, i.e., QR-dec (P, R_A) = QR(Ptil)
    Matrix P, R_A;
    computeQRdecomposition(Ptil, _globalRowsA, P, R_A);

    Matrix Btil(_phi.cols(), B.cols()); // Btil is of size (K_bar x m)
    // Btil := \phi^T *B
    _parMatrixOps->multiply(_phi.transpose(), B, Btil, (int) _phi.cols(), _globalRowsB, (int) B.cols());
    // Qtil := (I-\phi\phi^T)B
    Matrix Qtil = B - _phi * Btil;
 
    // compute orthogonal basis Q of Qtil, i.e., QR-dec (Q, R_B) = QR(Qtil)
    Matrix Q, R_B;
    computeQRdecomposition(Qtil, _globalRowsB, Q, R_B);
 
    //     e_orthModes.stop(true);

    Matrix K = Matrix::Zero(_psi.cols() + R_A.rows(), _psi.cols() + R_B.rows());
    Matrix K_A(_psi.cols() + R_A.rows(), Atil.cols());
    Matrix K_B(_phi.cols() + R_B.rows(), Btil.cols());
 
    for (int i = 0; i < _sigma.size(); i++)
      K(i, i) = _sigma(i);
 
    K_A.block(0, 0, Atil.rows(), Atil.cols())         = Atil;
    K_A.block(Atil.rows(), 0, R_A.rows(), R_A.cols()) = R_A;
    K_B.block(0, 0, Btil.rows(), Btil.cols())         = Btil;
    K_B.block(Btil.rows(), 0, R_B.rows(), R_B.cols()) = R_B;
    K += K_A * K_B.transpose();
 
    // compute svd of K
    Eigen::JacobiSVD<Matrix> svd(K, Eigen::ComputeThinU | Eigen::ComputeThinV);
    _sigma         = svd.singularValues();
    auto &psiPrime = svd.matrixU();
    auto &phiPrime = svd.matrixV();
    //     e_matK.stop(true);

    Matrix rotLeft(_rowsA, _psi.cols() + P.cols());
    Matrix rotRight(_rowsB, _phi.cols() + Q.cols());
 
    rotLeft.block(0, 0, _rowsA, _psi.cols())         = _psi;
    rotLeft.block(0, _psi.cols(), _rowsA, P.cols())  = P;
    rotRight.block(0, 0, _rowsB, _phi.cols())        = _phi;
    rotRight.block(0, _phi.cols(), _rowsB, Q.cols()) = Q;
 
    // [\psi,P] is distributed block-row wise, but \psiPrime is local on each proc, hence local mult.
    _psi = rotLeft * psiPrime;
    _phi = rotRight * phiPrime;
 
    //     e_rot.stop(true);

    _cols = _sigma.size();
 
    int waste = 0;
    for (int i = 0; i < (int) _sigma.size(); i++) {
      if (_sigma(i) < _sigma(0) * _truncationEps) {
        _cols = i;
        waste = _sigma.size() - i;
        break;
      }
    }
    _waste += waste;
 
    _psi.conservativeResize(_rowsA, _cols);
    _phi.conservativeResize(_rowsB, _cols);
    _sigma.conservativeResize(_cols);
    PRECICE_ASSERT(_sigma(0) >= 0.0);
    PRECICE_DEBUG("SVD factorization of Jacobian is truncated to {} DOFs. Cut off {} DOFs", _cols, waste);
 
    _initialSVD = true;
  }

  void initialize(PtrParMatrixOps parMatOps, int globalRowsA, int globalRowsB);

  void reset();

  Matrix &matrixPsi();

  Vector &singularValues();

  Matrix &matrixPhi();

  int cols();

  Rank rank();

  int getWaste();

  void setThreshold(double eps);

  double getThreshold();

  // void applyPreconditioner();

  // void revertPreconditioner();
 
  void setPrecondApplied(bool b);

  void setApplyFilterQR(bool b, double eps = 1e-3);
 
  // bool isPrecondApplied();
 
  bool isSVDinitialized();

  void setfstream(std::fstream *stream);
 
private:
  void computeQRdecomposition(Matrix const &A, int globalRows, Matrix &Q, Matrix &R);
 
  logging::Logger _log{"acceleration::SVDFactorization"};

  PtrPreconditioner _preconditioner;

  PtrParMatrixOps _parMatrixOps = nullptr;

  Matrix _psi;
  Matrix _phi;
  Vector _sigma;

  int _rowsA = 0;

  int _rowsB = 0;

  int _cols = 0;

  int _globalRowsA = 0;

  int _globalRowsB = 0;
 
  // Total number of truncated modes after last call to method getWaste()
  int _waste = 0;

  double _truncationEps;

  double _epsQR2 = 1e-3;

  bool _preconditionerApplied = false;

  bool _initialized = false;

  bool _initialSVD = false;
 
  bool _applyFilterQR = false;
};
 
} // namespace precice::acceleration::impl
 
#endif /* PRECICE_NO_MPI */