preCICE v3.1.1
Loading...
Searching...
No Matches
SVDFactorization.hpp
Go to the documentation of this file.
1#pragma once
2/*
3 * SVDFactorization.hpp
4 *
5 * Created on: Feb 5, 2016
6 * Author: Klaudius Scheufele
7 */
8
9#ifndef PRECICE_NO_MPI
10
11#include <Eigen/Core>
12#include <Eigen/Dense>
13#include <string>
14
18#include "logging/Logger.hpp"
20
21// ------- CLASS DEFINITION
22
23namespace precice {
24namespace acceleration {
25namespace impl {
26
32public:
33 // Eigen
34 typedef Eigen::MatrixXd Matrix;
35 typedef Eigen::VectorXd Vector;
36
41 double eps,
42 PtrPreconditioner preconditioner);
43
47 virtual ~SVDFactorization() {}
48
54 template <typename Derived1, typename Derived2>
55 void update(
56 const Eigen::MatrixBase<Derived1> &A,
57 const Eigen::MatrixBase<Derived2> &B)
58 {
66 if (_initialSVD) {
67 PRECICE_ASSERT(A.rows() == _rows, A.rows(), _rows);
68 PRECICE_ASSERT(B.rows() == _rows, B.rows(), _rows);
69 } else {
70 PRECICE_ASSERT(A.rows() == B.rows(), A.rows(), B.rows());
71 PRECICE_ASSERT(A.cols() == B.cols(), A.cols(), B.cols());
72 _rows = A.rows();
73 _cols = 0;
74 _psi = Matrix::Zero(_rows, 0);
75 _phi = Matrix::Zero(_rows, 0);
76 _sigma = Vector::Zero(0);
77 }
78
81 Matrix Atil(_psi.cols(), A.cols()); // Atil is of size (K_bar x m)
82
83 // Atil := \psi^T *A
84 // local computation of \psi^T * A and allreduce_sum to Atil (global), stored local on each proc
85 _parMatrixOps->multiply(_psi.transpose(), A, Atil, (int) _psi.cols(), _globalRows, (int) A.cols());
86
87 // Ptil := (I-\psi\psi^T)A
88 // Atil is local on each proc, thus fully local computation, embarrassingly parallel
89 Matrix Ptil = A - _psi * Atil;
90
91 // compute orthogonal basis P of Ptil, i.e., QR-dec (P, R_A) = QR(Ptil)
92 Matrix P, R_A;
93 computeQRdecomposition(Ptil, P, R_A);
94
97 Matrix Btil(_phi.cols(), B.cols()); // Btil is of size (K_bar x m)
98 // Btil := \phi^T *B
99 _parMatrixOps->multiply(_phi.transpose(), B, Btil, (int) _phi.cols(), _globalRows, (int) B.cols());
100 // Qtil := (I-\phi\phi^T)B
101 Matrix Qtil = B - _phi * Btil;
102
103 // compute orthogonal basis Q of Qtil, i.e., QR-dec (Q, R_B) = QR(Qtil)
104 Matrix Q, R_B;
105 computeQRdecomposition(Qtil, Q, R_B);
106
107 // e_orthModes.stop(true);
108
117 Matrix K = Matrix::Zero(_psi.cols() + R_A.rows(), _psi.cols() + R_B.rows());
118 Matrix K_A(_psi.cols() + R_A.rows(), Atil.cols());
119 Matrix K_B(_phi.cols() + R_B.rows(), Btil.cols());
120
121 for (int i = 0; i < _sigma.size(); i++)
122 K(i, i) = _sigma(i);
123
124 K_A.block(0, 0, Atil.rows(), Atil.cols()) = Atil;
125 K_A.block(Atil.rows(), 0, R_A.rows(), R_A.cols()) = R_A;
126 K_B.block(0, 0, Btil.rows(), Btil.cols()) = Btil;
127 K_B.block(Btil.rows(), 0, R_B.rows(), R_B.cols()) = R_B;
128 K += K_A * K_B.transpose();
129
130 // compute svd of K
131 Eigen::JacobiSVD<Matrix> svd(K, Eigen::ComputeThinU | Eigen::ComputeThinV);
132 _sigma = svd.singularValues();
133 auto &psiPrime = svd.matrixU();
134 auto &phiPrime = svd.matrixV();
135 // e_matK.stop(true);
136
139 Matrix rotLeft(_rows, _psi.cols() + P.cols());
140 Matrix rotRight(_rows, _phi.cols() + Q.cols());
141
142 rotLeft.block(0, 0, _rows, _psi.cols()) = _psi;
143 rotLeft.block(0, _psi.cols(), _rows, P.cols()) = P;
144 rotRight.block(0, 0, _rows, _phi.cols()) = _phi;
145 rotRight.block(0, _phi.cols(), _rows, Q.cols()) = Q;
146
147 // [\psi,P] is distributed block-row wise, but \psiPrime is local on each proc, hence local mult.
148 _psi = rotLeft * psiPrime;
149 _phi = rotRight * phiPrime;
150
151 // e_rot.stop(true);
152
155 _cols = _sigma.size();
156
157 int waste = 0;
158 for (int i = 0; i < (int) _sigma.size(); i++) {
159 if (_sigma(i) < _sigma(0) * _truncationEps) {
160 _cols = i;
161 waste = _sigma.size() - i;
162 break;
163 }
164 }
165 _waste += waste;
166
167 _psi.conservativeResize(_rows, _cols);
168 _phi.conservativeResize(_rows, _cols);
169 _sigma.conservativeResize(_cols);
170 PRECICE_ASSERT(_sigma(0) >= 0.0);
171 PRECICE_DEBUG("SVD factorization of Jacobian is truncated to {} DOFs. Cut off {} DOFs", _cols, waste);
172
173 _initialSVD = true;
174 }
175
180 void initialize(PtrParMatrixOps parMatOps, int globalRows);
181
185 void reset();
186
190 Matrix &matrixPsi();
191
196
200 Matrix &matrixPhi();
201
203 int cols();
204
206 int rows();
207
209 Rank rank();
210
212 int getWaste();
213
215 void setThreshold(double eps);
216
218 double getThreshold();
219
221 // void applyPreconditioner();
222
224 // void revertPreconditioner();
225
226 void setPrecondApplied(bool b);
227
229 void setApplyFilterQR(bool b, double eps = 1e-3);
230
231 // bool isPrecondApplied();
232
233 bool isSVDinitialized();
234
237
238private:
251 void computeQRdecomposition(Matrix const &A, Matrix &Q, Matrix &R);
252
253 logging::Logger _log{"acceleration::SVDFactorization"};
254
257
260
265
267 int _rows = 0;
268
270 int _cols = 0;
271
273 int _globalRows = 0;
274
275 // Total number of truncated modes after last call to method getWaste()
276 int _waste = 0;
277
280
282 double _epsQR2 = 1e-3;
283
286
288 bool _initialized = false;
289
291 bool _initialSVD = false;
292
293 bool _applyFilterQR = false;
294};
295
296} // namespace impl
297} // namespace acceleration
298} // namespace precice
299
300#endif /* PRECICE_NO_MPI */
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
Class that provides functionality to maintain a SVD decomposition of a matrix via successive rank-1 u...
int _cols
Number of columns, i.e., rank of the truncated svd.
PtrParMatrixOps _parMatrixOps
: object for parallel matrix operations, i.e., parallel mat-mat/ mat-vec multiplications
void setApplyFilterQR(bool b, double eps=1e-3)
: enables or disables an additional QR-2 filter for the QR-decomposition
int _rows
Number of rows (on each proc, i.e., local)
double _truncationEps
Truncation parameter for the updated SVD decomposition.
bool _initialSVD
true, if at least one update has been made, i.e., the number of rows is known and a initial rank is g...
double getThreshold()
: returns the truncation threshold for the SVD
bool _preconditionerApplied
true if the preconditioner has been applied appropriate to the updated SVD decomposition
void setThreshold(double eps)
: sets the threshold for the truncation of the SVD factorization
int _globalRows
Number of global rows, i.e., sum of _rows for all procs.
Matrix _psi
: SVD factorization of the matrix J = _psi * _sigma * _phi^T
void setPrecondApplied(bool b)
: applies the preconditioner to the factorized and truncated representation of the Jacobian matrix
Matrix & matrixPhi()
: returns a matrix representation of the orthogonal matrix Phi, A = Psi * Sigma * Phi^T
void initialize(PtrParMatrixOps parMatOps, int globalRows)
: initializes the updated SVD factorization, i.e., sets the object for parallel matrix-matrix operati...
Vector & singularValues()
: returns a matrix representation of the orthogonal matrix Sigma, A = Psi * Sigma * Phi^T
void reset()
: resets the SVD factorization
Rank rank()
: returns the rank of the truncated SVD factorization
double _epsQR2
Threshold for the QR2 filter for the QR decomposition.
void setfstream(std::fstream *stream)
Optional file-stream for logging output.
Matrix & matrixPsi()
: returns a matrix representation of the orthogonal matrix Psi, A = Psi * Sigma * Phi^bs
bool _initialized
true, if ParallelMatrixOperations object is set, i.e., initialized
void update(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
: updates the SVD decomposition with the rank-1 update A*B^T, i.e., _psi * _sigma * _phi^T + A*B^T an...
int cols()
: returns the number of columns in the QR-decomposition
SVDFactorization(double eps, PtrPreconditioner preconditioner)
Constructor.
int rows()
: returns the number of rows in the QR-decomposition
void computeQRdecomposition(Matrix const &A, Matrix &Q, Matrix &R)
: computes the QR decomposition of a matrix A of type A = PSI^T*A \in R^(rank x n)
PtrPreconditioner _preconditioner
: preconditioner for least-squares system if vectorial system is used.
int getWaste()
: returns the total number of truncated modes since last call to this method
This class provides a lightweight logger.
Definition Logger.hpp:16
Main namespace of the precice library.
int Rank
Definition Types.hpp:37
static std::unique_ptr< precice::Participant > impl
Definition preciceC.cpp:21