preCICE v3.1.1
Loading...
Searching...
No Matches
SVDFactorization.cpp
Go to the documentation of this file.
1#ifndef PRECICE_NO_MPI
2
4#include <Eigen/Core>
5#include <limits>
6#include <utility>
7
8#include "utils/IntraComm.hpp"
9
11
13 double eps,
14 PtrPreconditioner preconditioner)
15 : _preconditioner(std::move(preconditioner)),
16 _truncationEps(eps)
17{
18}
19
21 PtrParMatrixOps parOps,
22 int globalRows)
23{
24 _parMatrixOps = std::move(parOps);
25 _globalRows = globalRows;
26 _initialized = true;
27}
28
29/*
30void SVDFactorization::applyPreconditioner()
31{
32 PRECICE_TRACE();
33
34 if(_psi.size() > 0 && _phi.size() > 0){
35 // apply preconditioner: \psi_i * P_i, corresponds to Wtil_i * P_i, local!
36 _preconditioner->apply(_psi);
37 // apply preconditioner: \phi_i^T * P_i^{-1}, corresponds to Z_i * P_i^{-1}, local!
38 // here, \phi^T should be preconditioned from right with inv_weights, i.e., the columns
39 // of \phi^T are scaled. This is identical to scaling the rows of \phi, i.e., applying
40 // P^{-1} * \phi
41 _preconditioner->revert(_phi);
42 }
43 _preconditionerApplied = true;
44}
45
46void SVDFactorization::revertPreconditioner()
47{
48 PRECICE_TRACE();
49
50 if(_psi.size() > 0 && _phi.size() > 0){
51 // revert preconditioner: \psi_i * P_i^{-1}, corresponds to Wtil_i * P_i^{-1}, local!
52 _preconditioner->revert(_psi);
53 // revert preconditioner: \phi_i^T * P_i, corresponds to Z_i * P_i, local!
54 // here, \phi^T should be preconditioned from right with _weights, i.e., the columns
55 // of \phi^T are scaled. This is identical to scaling the rows of \phi, i.e., applying
56 // P * \phi
57 _preconditioner->apply(_phi);
58 }
59 _preconditionerApplied = false;
60}
61*/
62
64{
65 _psi.resize(0, 0);
66 _phi.resize(0, 0);
67 _sigma.resize(0);
69 _initialSVD = false;
70 _applyFilterQR = false;
71 _epsQR2 = 1e-3;
72}
73
75 Matrix const &A,
76 Matrix & Q,
77 Matrix & R)
78{
80
81 // if nothing is linear dependent, the dimensions stay like this
82 Q = Matrix::Zero(A.rows(), A.cols());
83 R = Matrix::Zero(A.cols(), A.cols());
84 int colsR = 0;
85 int rowsR = 0;
86 int colsQ = 0;
87
88 // magic params:
89 double omega = 0.;
90 double theta = std::sqrt(2);
91
92 // columns need to be inserted at the back, otherwise we would have to perform
93 // givens rotations, to re-establish the upper diagonal form of R
94 for (int colIndex = 0; colIndex < A.cols(); colIndex++) {
95
96 // invariants:
97 PRECICE_ASSERT(colsQ == rowsR, colsQ, rowsR);
98 PRECICE_ASSERT(colsQ <= colIndex, colsQ, colIndex);
99
100 Vector col = A.col(colIndex);
101
102 // if system is quadratic; discard
103 if (_globalRows == colIndex) {
104 PRECICE_WARN("The matrix that is about to be factorized is quadratic, i.e., the new column cannot be orthogonalized; discard.");
105 return;
106 }
107
111 Vector r = Vector::Zero(R.rows()); // gram-schmidt coefficients orthogonalization
112 Vector s = Vector::Zero(R.rows()); // gram-schmidt coefficients re-orthogonalization
113 Vector u = Vector::Zero(A.rows()); // sum of projections
114 double rho_orth = 0.;
115 double rho0 = utils::IntraComm::l2norm(col); // distributed l2norm;
116 double rho00 = rho0; // save norm of col for QR2 filter crit.
117
118 int its = 0;
119 bool termination = false;
120 bool orthogonalized = true;
121 // while col is not sufficiently orthogonal to Q
122 while (!termination) {
123
124 // take a gram-schmidt iteration
125 u = Vector::Zero(A.rows());
126 for (int j = 0; j < colsQ; j++) {
127
128 // dot-product <_Q(:,j), v >
129 Vector Qc = Q.col(j);
130 // dot product <_Q(:,j), v> =: r_ij
131 double r_ij = utils::IntraComm::dot(Qc, col);
132 // save r_ij in s(j) = column of R
133 s(j) = r_ij;
134 // u is the sum of projections r_ij * _Q(:,j) = _Q(:,j) * <_Q(:,j), v>
135 u += Q.col(j) * r_ij;
136 }
137 // add the gram-schmidt coefficients over all iterations of possible re-orthogonalizations
138 r += s;
139 // subtract projections from v, v is now orthogonal to columns of Q
140 col -= u;
141
142 // rho1 = norm of orthogonalized new column v_tilde (though not normalized)
143 rho_orth = utils::IntraComm::l2norm(col);
144 // t = norm of _r(:,j) with j = colNum-1
145 double norm_coefficients = utils::IntraComm::l2norm(s); // distributed l2norm
146
147 its++;
148
149 // if ||v_orth|| is nearly zero, col is not well orthogonalized; discard
150 if (rho_orth <= std::numeric_limits<double>::min()) {
151 PRECICE_DEBUG("The norm of v_orthogonal is almost zero, i.e., failed to orthogonalize column v; discard.");
152 orthogonalized = false;
153 termination = true;
154 }
155
168 if (rho_orth * theta <= rho0 + omega * norm_coefficients) {
169 // exit to fail if too many iterations
170 if (its >= 4) {
171 PRECICE_WARN("Matrix Q is not sufficiently orthogonal. Failed to orthogonalize new column after 4 iterations. New column will be discarded.");
172 orthogonalized = false;
173 termination = true;
174 }
175
176 // for re-orthogonalization
177 rho0 = rho_orth;
178
179 } else {
180 termination = true;
181 }
182 }
183
184 // if the QR2-filter crit. kicks in with threshold eps.
185 if (_applyFilterQR && orthogonalized && rho_orth <= _epsQR2 * rho00) {
186 orthogonalized = false;
187 }
188
189 // normalize col
190 double rho = orthogonalized ? rho_orth : 1.0;
191 col /= rho;
192 r(rowsR) = rho;
193
194 // as we always insert at the rightmost position, no need to shift
195 // entries of R or apply givens rotations on the QR-dec to maintain
196 // the upper triangular structure of R
197 Q.col(colsQ) = col; // insert orthogonalized column to the right in Q
198 R.col(colsR) = r; // insert gram-schmidt coefficients to the right in R
199
200 colsR++;
201 PRECICE_ASSERT(colsR <= R.cols(), colsR, R.cols());
202 rowsR++;
203 PRECICE_ASSERT(rowsR <= R.rows(), rowsR, R.rows());
204 colsQ++;
205
206 // failed to orthogonalize the column, i.e., it is linear dependent;
207 // modify the QR-dec such that it stays valid (column deleted) while
208 // also staying aligned with the dimension of A, MUST have the same
209 // number of cols (cannot delete from A)
210 if (not orthogonalized) {
211
212 colsQ--;
213 PRECICE_ASSERT(colsQ >= 0, colsQ);
214 rowsR--;
215 PRECICE_ASSERT(rowsR >= 0, rowsR);
216 // delete column that was just inserted (as it is not orthogonal to Q)
217 Q.col(colsQ) = Vector::Zero(A.rows());
218 // delete line in R that corresponds to the just inserted but not orthogonal column
219 // as we always insert to the right, no shifting/ application of givens roatations is
220 // necessary.
221 // Note: The corresponding column from R with index colIndex is not deleted: dimensions must align with A.
222 PRECICE_ASSERT(R(rowsR, colsR - 1) == 1.0, R(rowsR, colsR - 1));
223 R.row(rowsR) = Vector::Zero(A.cols());
224 }
225 }
226 // shrink matrices Q, R to actual size
227 Q.conservativeResize(A.rows(), colsQ);
228 R.conservativeResize(rowsR, colsR);
229}
230
235
240
245
250
252{
253 _applyFilterQR = b;
254 _epsQR2 = eps;
255}
256
257/*
258bool SVDFactorization::isPrecondApplied()
259{
260 return _preconditionerApplied;
261}
262*/
263
268
270{
271 _truncationEps = eps;
272}
273
275{
276 return _truncationEps;
277}
278
280{
281 int r = _waste;
282 _waste = 0;
283 return r;
284}
285
287{
288 return _cols;
289}
290
292{
293 return _rows;
294}
295
297{
298 return _cols;
299}
300
301} // namespace precice::acceleration::impl
302
303#endif // PRECICE_NO_MPI
#define PRECICE_WARN(...)
Definition LogMacros.hpp:11
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
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.
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
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)
int getWaste()
: returns the total number of truncated modes since last call to this method
static double l2norm(const Eigen::VectorXd &vec)
The l2 norm of a vector is calculated on distributed data.
Definition IntraComm.cpp:67
static double dot(const Eigen::VectorXd &vec1, const Eigen::VectorXd &vec2)
int Rank
Definition Types.hpp:37
STL namespace.
T sqrt(T... args)