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 {
namespace acceleration {
namespace impl {
 
class SVDFactorization {
public:
  // Eigen
  typedef Eigen::MatrixXd Matrix;
  typedef Eigen::VectorXd Vector;
 
  SVDFactorization(
      double            eps,
      PtrPreconditioner preconditioner);
 
  virtual ~SVDFactorization() {}
 
  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() == _rows, A.rows(), _rows);
      PRECICE_ASSERT(B.rows() == _rows, B.rows(), _rows);
    } else {
      PRECICE_ASSERT(A.rows() == B.rows(), A.rows(), B.rows());
      PRECICE_ASSERT(A.cols() == B.cols(), A.cols(), B.cols());
      _rows  = A.rows();
      _cols  = 0;
      _psi   = Matrix::Zero(_rows, 0);
      _phi   = Matrix::Zero(_rows, 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(), _globalRows, (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, 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(), _globalRows, (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, 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(_rows, _psi.cols() + P.cols());
    Matrix rotRight(_rows, _phi.cols() + Q.cols());
 
    rotLeft.block(0, 0, _rows, _psi.cols())         = _psi;
    rotLeft.block(0, _psi.cols(), _rows, P.cols())  = P;
    rotRight.block(0, 0, _rows, _phi.cols())        = _phi;
    rotRight.block(0, _phi.cols(), _rows, 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(_rows, _cols);
    _phi.conservativeResize(_rows, _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 globalRows);
 
  void reset();
 
  Matrix &matrixPsi();
 
  Vector &singularValues();
 
  Matrix &matrixPhi();
 
  int cols();
 
  int rows();
 
  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, Matrix &Q, Matrix &R);
 
  logging::Logger _log{"acceleration::SVDFactorization"};
 
  PtrPreconditioner _preconditioner;
 
  PtrParMatrixOps _parMatrixOps = nullptr;
 
  Matrix _psi;
  Matrix _phi;
  Vector _sigma;
 
  int _rows = 0;
 
  int _cols = 0;
 
  int _globalRows = 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 impl
} // namespace acceleration
} // namespace precice
 
#endif /* PRECICE_NO_MPI */