preCICE v3.1.1
Loading...
Searching...
No Matches
IQNIMVJAcceleration.cpp
Go to the documentation of this file.
1#ifndef PRECICE_NO_MPI
2
3#include <Eigen/Core>
4#include <algorithm>
5#include <fstream>
6#include <map>
7#include <memory>
8#include <string>
9#include <utility>
10
15#include "com/Communication.hpp"
19#include "logging/LogMacros.hpp"
21#include "profiling/Event.hpp"
23#include "utils/IntraComm.hpp"
24#include "utils/assertion.hpp"
25
27
28namespace precice::acceleration {
29
30// ==================================================================================
32 double initialRelaxation,
33 bool forceInitialRelaxation,
34 int maxIterationsUsed,
35 int pastTimeWindowsReused,
36 int filter,
37 double singularityLimit,
38 std::vector<int> dataIDs,
39 const impl::PtrPreconditioner &preconditioner,
40 bool alwaysBuildJacobian,
41 int imvjRestartType,
42 int chunkSize,
43 int RSLSreusedTimeWindows,
44 double RSSVDtruncationEps)
45 : BaseQNAcceleration(initialRelaxation, forceInitialRelaxation, maxIterationsUsed, pastTimeWindowsReused,
46 filter, singularityLimit, std::move(dataIDs), preconditioner),
47 // _secondaryOldXTildes(),
48 _invJacobian(),
49 _oldInvJacobian(),
50 _Wtil(),
51 _WtilChunk(),
52 _pseudoInverseChunk(),
53 _matrixV_RSLS(),
54 _matrixW_RSLS(),
55 _matrixCols_RSLS(),
56 _parMatrixOps(nullptr),
57 _svdJ(RSSVDtruncationEps, preconditioner),
58 _alwaysBuildJacobian(alwaysBuildJacobian),
59 _imvjRestartType(imvjRestartType),
60 _imvjRestart(false),
61 _chunkSize(chunkSize),
62 _RSLSreusedTimeWindows(RSLSreusedTimeWindows),
63 _nbRestarts(0),
64 _avgRank(0)
65{
66}
67
68// ==================================================================================
70
71// ==================================================================================
73 const DataMap &cplData)
74{
76
77 // do common QN acceleration initialization
79
80 if (_imvjRestartType > 0)
81 _imvjRestart = true;
82
83 // initialize parallel matrix-matrix operation module
84 _parMatrixOps = std::make_shared<impl::ParallelMatrixOperations>();
85 _parMatrixOps->initialize(not _imvjRestart);
87
88 int entries = _residuals.size();
89 int global_n = 0;
90
92 global_n = entries;
93 } else {
94 global_n = _dimOffsets.back();
95 }
96
97 if (not _imvjRestart) {
98 // only need memory for Jacobain of not in restart mode
99 _invJacobian = Eigen::MatrixXd::Zero(global_n, entries);
100 _oldInvJacobian = Eigen::MatrixXd::Zero(global_n, entries);
101 }
102 // initialize V, W matrices for the LS restart
103 if (_imvjRestartType == RS_LS) {
105 _matrixV_RSLS = Eigen::MatrixXd::Zero(entries, 0);
106 _matrixW_RSLS = Eigen::MatrixXd::Zero(entries, 0);
107 }
108 _Wtil = Eigen::MatrixXd::Zero(entries, 0);
109
111 _infostringstream << " IMVJ restart mode: " << _imvjRestart << "\n chunk size: " << _chunkSize << "\n trunc eps: " << _svdJ.getThreshold() << "\n R_RS: " << _RSLSreusedTimeWindows << "\n--------\n"
112 << '\n';
113 }
114}
115
116// ==================================================================================
118 const DataMap &cplData)
119{
120 // Perform underrelaxation with initial relaxation factor for secondary data
121 for (int id : _secondaryDataIDs) {
122 PtrCouplingData data = cplData.at(id);
123 Eigen::VectorXd &values = data->values();
124 values *= _initialRelaxation; // new * omg
125 Eigen::VectorXd &secResiduals = _secondaryResiduals[id];
126 secResiduals = data->previousIteration();
127 secResiduals *= 1.0 - _initialRelaxation; // (1-omg) * old
128 values += secResiduals; // (1-omg) * old + new * omg
129 }
130}
131
132// ==================================================================================
134 const DataMap &cplData)
135{
142
143 // call the base method for common update of V, W matrices
144 // important that base method is called before updating _Wtil
146
147 // update _Wtil if the efficient computation of the quasi-Newton update is used
148 // or update current Wtil if the restart mode of imvj is used
151 // do nothing: constant relaxation
152 } else {
153 if (not _firstIteration) {
154 // Update matrix _Wtil = (W - J_prev*V) with newest information
155
156 Eigen::VectorXd v = _matrixV.col(0);
157 Eigen::VectorXd w = _matrixW.col(0);
158
159 // here, we check for _Wtil.cols() as the matrices V, W need to be updated before hand
160 // and thus getLSSystemCols() does not yield the correct result.
161 bool columnLimitReached = _Wtil.cols() == _maxIterationsUsed;
162 bool overdetermined = _Wtil.cols() <= getLSSystemRows();
163
164 Eigen::VectorXd wtil = Eigen::VectorXd::Zero(_matrixV.rows());
165
166 // add column: Wtil(:,0) = W(:,0) - sum_q [ Wtil^q * ( Z^q * V(:,0)) ]
167 // |--- J_prev ---|
168 // iterate over all stored Wtil and Z matrices in current chunk
169 if (_imvjRestart) {
170 for (int i = 0; i < static_cast<int>(_WtilChunk.size()); i++) {
171 int colsLSSystemBackThen = _pseudoInverseChunk[i].rows();
172 PRECICE_ASSERT(colsLSSystemBackThen == _WtilChunk[i].cols(), colsLSSystemBackThen, _WtilChunk[i].cols());
173 Eigen::VectorXd Zv = Eigen::VectorXd::Zero(colsLSSystemBackThen);
174 // multiply: Zv := Z^q * V(:,0) of size (m x 1)
175 _parMatrixOps->multiply(_pseudoInverseChunk[i], v, Zv, colsLSSystemBackThen, getLSSystemRows(), 1);
176 // multiply: Wtil^q * Zv dimensions: (n x m) * (m x 1), fully local
177 wtil += _WtilChunk[i] * Zv;
178 }
179
180 // store columns if restart mode = RS-LS
181 if (_imvjRestartType == RS_LS) {
185 }
186
187 // imvj without restart is used, but efficient update, i.e. no Jacobian assembly in each iteration
188 // add column: Wtil(:,0) = W(:,0) - J_prev * V(:,0)
189 } else {
190 // compute J_prev * V(0) := wtil the new column in _Wtil of dimension: (n x n) * (n x 1) = (n x 1),
191 // parallel: (n_global x n_local) * (n_local x 1) = (n_local x 1)
192 _parMatrixOps->multiply(_oldInvJacobian, v, wtil, _dimOffsets, getLSSystemRows(), getLSSystemRows(), 1, false);
193 }
194 wtil *= -1;
195 wtil += w;
196
197 if (not columnLimitReached && overdetermined) {
199 } else {
201 }
202 }
203 }
204 }
205}
206
207// ==================================================================================
209 const DataMap & cplData,
210 Eigen::VectorXd &xUpdate)
211{
230 computeNewtonUpdate(cplData, xUpdate);
231 } else {
232 computeNewtonUpdateEfficient(cplData, xUpdate);
233 }
234}
235
236// ==================================================================================
238 Eigen::MatrixXd &pseudoInverse)
239{
245 auto Q = _qrV.matrixQ();
246 auto R = _qrV.matrixR();
247
250
251 Eigen::VectorXd yVec(pseudoInverse.rows());
252
253 // assertions for the case of processors with no vertices
256 PRECICE_ASSERT(_qrV.rows() == 0, _qrV.rows());
257 PRECICE_ASSERT(Q.size() == 0, Q.size());
258 }
259
260 // backsubstitution
261 for (int i = 0; i < Q.rows(); i++) {
262 Eigen::VectorXd Qrow = Q.row(i);
263 yVec = R.triangularView<Eigen::Upper>().solve<Eigen::OnTheLeft>(Qrow);
264 pseudoInverse.col(i) = yVec;
265 } // ----------------
266
267 // scale pseudo inverse back Z := Z' * P,
268 // Z' is scaled pseudo inverse i.e, Z' = R^-1 * Q^T * P^-1
269 _preconditioner->apply(pseudoInverse, true);
270 // e.stop(true);
271}
272
273// ==================================================================================
275{
280
281 PRECICE_ASSERT(_matrixV.rows() == _qrV.rows(), _matrixV.rows(), _qrV.rows());
283
284 _Wtil = Eigen::MatrixXd::Zero(_qrV.rows(), _qrV.cols());
285
286 // imvj restart mode: re-compute Wtil: Wtil = W - sum_q [ Wtil^q * (Z^q*V) ]
287 // |--- J_prev ---|
288 // iterate over all stored Wtil and Z matrices in current chunk
289 if (_imvjRestart) {
290 for (int i = 0; i < static_cast<int>(_WtilChunk.size()); i++) {
291 int colsLSSystemBackThen = _pseudoInverseChunk[i].rows();
292 PRECICE_ASSERT(colsLSSystemBackThen == _WtilChunk[i].cols(), colsLSSystemBackThen, _WtilChunk[i].cols());
293 Eigen::MatrixXd ZV = Eigen::MatrixXd::Zero(colsLSSystemBackThen, _qrV.cols());
294 // multiply: ZV := Z^q * V of size (m x m) with m=#cols, stored on each proc.
295 _parMatrixOps->multiply(_pseudoInverseChunk[i], _matrixV, ZV, colsLSSystemBackThen, getLSSystemRows(), _qrV.cols());
296 // multiply: Wtil^q * ZV dimensions: (n x m) * (m x m), fully local and embarrassingly parallel
297 _Wtil += _WtilChunk[i] * ZV;
298 }
299
300 // imvj without restart is used, i.e., recompute Wtil: Wtil = W - J_prev * V
301 } else {
302 // multiply J_prev * V = W_til of dimension: (n x n) * (n x m) = (n x m),
303 // parallel: (n_global x n_local) * (n_local x m) = (n_local x m)
305 }
306
307 // W_til = (W-J_inv_n*V) = (W-V_tilde)
308 _Wtil *= -1.;
309 _Wtil = _Wtil + _matrixW;
310
311 _resetLS = false;
312 // e.stop(true);
313}
314
315// ==================================================================================
317{
329 Eigen::MatrixXd Z(_qrV.cols(), _qrV.rows());
330 pseudoInverse(Z);
331
335 PRECICE_ASSERT(_matrixV.rows() == _qrV.rows(), _matrixV.rows(), _qrV.rows());
337 if (_resetLS) {
338 buildWtil();
339 PRECICE_WARN(" ATTENTION, in buildJacobian call for buildWtill() - this should not be the case except the coupling did only one iteration");
340 }
341
348 // --------
349
350 // update Jacobian
352}
353
354// ==================================================================================
356 const DataMap & cplData,
357 Eigen::VectorXd &xUpdate)
358{
360
381 Eigen::MatrixXd Z(_qrV.cols(), _qrV.rows());
382 pseudoInverse(Z);
383
387 PRECICE_ASSERT(_matrixV.rows() == _qrV.rows(), _matrixV.rows(), _qrV.rows());
389
390 // rebuild matrix Wtil if V changes substantially.
391 if (_resetLS) {
392 buildWtil();
393 }
394
406 Eigen::VectorXd negativeResiduals = -_residuals;
407 Eigen::VectorXd r_til = Eigen::VectorXd::Zero(getLSSystemCols());
408 _parMatrixOps->multiply(Z, negativeResiduals, r_til, getLSSystemCols(), getLSSystemRows(), 1); // --------
409
418 Eigen::VectorXd xUptmp(_residuals.size());
419 xUpdate = Eigen::VectorXd::Zero(_residuals.size());
420 xUptmp = _Wtil * r_til; // local product, result is naturally distributed.
421
428 if (_imvjRestart) {
429 for (int i = 0; i < static_cast<int>(_WtilChunk.size()); i++) {
430 int colsLSSystemBackThen = _pseudoInverseChunk[i].rows();
431 PRECICE_ASSERT(colsLSSystemBackThen == _WtilChunk[i].cols(), colsLSSystemBackThen, _WtilChunk[i].cols());
432 r_til = Eigen::VectorXd::Zero(colsLSSystemBackThen);
433 // multiply: r_til := Z^q * (-res) of size (m x 1) with m=#cols of LS at that time, result stored on each proc.
434 _parMatrixOps->multiply(_pseudoInverseChunk[i], negativeResiduals, r_til, colsLSSystemBackThen, getLSSystemRows(), 1);
435 // multiply: Wtil^q * r_til dimensions: (n x m) * (m x 1), fully local and embarrassingly parallel
436 xUpdate += _WtilChunk[i] * r_til;
437 }
438
439 // imvj without restart is used, i.e., compute directly J_prev * (-res)
440 } else {
441 _parMatrixOps->multiply(_oldInvJacobian, negativeResiduals, xUpdate, _dimOffsets, getLSSystemRows(), getLSSystemRows(), 1, false);
442 PRECICE_DEBUG("Mult J*V DONE");
443 }
444
445 xUpdate += xUptmp;
446
447 // pending deletion: delete Wtil
449 _Wtil.conservativeResize(0, 0);
450 _resetLS = true;
451 }
452}
453
454// ==================================================================================
455void IQNIMVJAcceleration::computeNewtonUpdate(const DataMap &cplData, Eigen::VectorXd &xUpdate)
456{
458
466 Eigen::MatrixXd Z(_qrV.cols(), _qrV.rows());
467 pseudoInverse(Z);
468
471 buildWtil();
472
479
480 // update Jacobian
482
485 Eigen::VectorXd negativeResiduals = -_residuals;
486
487 // multiply J_inv * (-res) = x_Update of dimension: (n x n) * (n x 1) = (n x 1),
488 // parallel: (n_global x n_local) * (n_local x 1) = (n_local x 1)
489 _parMatrixOps->multiply(_invJacobian, negativeResiduals, xUpdate, _dimOffsets, getLSSystemRows(), getLSSystemRows(), 1, false); // --------
490}
491
492// ==================================================================================
494{
496
497 // int used_storage = 0;
498 // int theoreticalJ_storage = 2*getLSSystemRows()*_residuals.size() + 3*_residuals.size()*getLSSystemCols() + _residuals.size()*_residuals.size();
499 // ------------ RESTART SVD ------------
501
502 // we need to compute the updated SVD of the scaled Jacobian matrix
503 // |= APPLY PRECONDITIONING J_prev = Wtil^q, Z^q ===|
504 for (int i = 0; i < static_cast<int>(_WtilChunk.size()); i++) {
505 _preconditioner->apply(_WtilChunk[i]);
506 _preconditioner->revert(_pseudoInverseChunk[i], true);
507 }
508 // |=================== ===|
509
510 Rank rankBefore = _svdJ.isSVDinitialized() ? _svdJ.rank() : 0;
511
512 // if it is the first time window, there is no initial SVD, so take all Wtil, Z matrices
513 // otherwise, the first element of each container holds the decomposition of the current
514 // truncated SVD, i.e., Wtil^0 = \phi, Z^0 = S\psi^T, this should not be added to the SVD.
515 int q = _svdJ.isSVDinitialized() ? 1 : 0;
516
517 // perform M-1 rank-1 updates of the truncated SVD-dec of the Jacobian
518 for (; q < static_cast<int>(_WtilChunk.size()); q++) {
519 // update SVD, i.e., PSI * SIGMA * PHI^T <-- PSI * SIGMA * PHI^T + Wtil^q * Z^q
520 _svdJ.update(_WtilChunk[q], _pseudoInverseChunk[q].transpose());
521 // used_storage += 2*_WtilChunk.size();
522 }
523 // int m = _WtilChunk[q].cols(), n = _WtilChunk[q].rows();
524 // used_storage += 2*rankBefore*m + 4*m*n + 2*m*m + (rankBefore+m)*(rankBefore+m) + 2*n*(rankBefore+m);
525
526 // drop all stored Wtil^q, Z^q matrices
529
530 auto &psi = _svdJ.matrixPsi();
531 auto &sigma = _svdJ.singularValues();
532 auto &phi = _svdJ.matrixPhi();
533
534 // multiply sigma * phi^T, phi is distributed block-row wise, phi^T is distributed block-column wise
535 // sigma is stored local on each proc, thus, the multiplication is fully local, no communication.
536 // Z = sigma * phi^T
537 Eigen::MatrixXd Z(phi.cols(), phi.rows());
538 for (int i = 0; i < static_cast<int>(Z.rows()); i++)
539 for (int j = 0; j < static_cast<int>(Z.cols()); j++)
540 Z(i, j) = phi(j, i) * sigma[i];
541
542 Rank rankAfter = _svdJ.rank();
543 int waste = _svdJ.getWaste();
544 _avgRank += rankAfter;
545
546 // store factorized truncated SVD of J
549
550 // |= REVERT PRECONDITIONING J_prev = Wtil^0, Z^0 ==|
553 // |=================== ==|
554
555 PRECICE_DEBUG("MVJ-RESTART, mode=SVD. Rank of truncated SVD of Jacobian {}, new modes: {}, truncated modes: {} avg rank: {}", rankAfter, rankAfter - rankBefore, waste, _avgRank / _nbRestarts);
556
557 // double percentage = 100.0*used_storage/(double)theoreticalJ_storage;
559 _infostringstream << " - MVJ-RESTART " << _nbRestarts << ", mode= SVD -\n new modes: " << rankAfter - rankBefore << "\n rank svd: " << rankAfter << "\n avg rank: " << _avgRank / _nbRestarts << "\n truncated modes: " << waste << "\n"
560 << '\n';
561 }
562
563 // ------------ RESTART LEAST SQUARES ------------
565 // drop all stored Wtil^q, Z^q matrices
568
569 if (_matrixV_RSLS.cols() > 0) {
570 // avoid that the system is getting too squared
571 while (_matrixV_RSLS.cols() * 2 >= getLSSystemRows()) {
573 }
574
575 // preconditioning
576 // V needs to be sclaed to compute the pseudo inverse
577 // W only needs to be scaled, as the design requires to store scaled
578 // matrices Wtil^0 and Z^0 as initial guess after the restart
581
584 // for QR2-filter, the QR-dec is computed in qr-applyFilter()
586 for (int i = 0; i < static_cast<int>(_matrixV_RSLS.cols()); i++) {
587 Eigen::VectorXd v = _matrixV_RSLS.col(i);
588 qr.pushBack(v); // same order as matrix V_RSLS
589 }
590 }
591
592 // apply filter
594 std::vector<int> delIndices(0);
596 // start with largest index (as V,W matrices are shrunk and shifted
597 for (int i = delIndices.size() - 1; i >= 0; i--) {
598 removeMatrixColumnRSLS(delIndices[i]);
599 }
600 PRECICE_ASSERT(_matrixV_RSLS.cols() == qr.cols(), _matrixV_RSLS.cols(), qr.cols());
601 }
602
607 auto Q = qr.matrixQ();
608 auto R = qr.matrixR();
609 Eigen::MatrixXd pseudoInverse(qr.cols(), qr.rows());
610 Eigen::VectorXd yVec(pseudoInverse.rows());
611
612 // backsubstitution
613 for (int i = 0; i < Q.rows(); i++) {
614 Eigen::VectorXd Qrow = Q.row(i);
615 yVec = R.triangularView<Eigen::Upper>().solve<Eigen::OnTheLeft>(Qrow);
616 pseudoInverse.col(i) = yVec;
617 }
618
619 // scale pseudo inverse back Z := Z' * P,
620 // Z' is scaled pseudo inverse i.e, Z' = R^-1 * Q^T * P^-1
621 //_preconditioner->apply(pseudoInverse, true, false);
622
623 // store factorization of least-squares initial guess for Jacobian
626
627 // |= REVERT PRECONDITIONING J_prev = Wtil^0, Z^0 ==|
632 // |=================== ==|
633 }
634
635 PRECICE_DEBUG("MVJ-RESTART, mode=LS. Restart with {} columns from {} time windows.", _matrixV_RSLS.cols(), _RSLSreusedTimeWindows);
637 _infostringstream << " - MVJ-RESTART" << _nbRestarts << ", mode= LS -\n used cols: " << _matrixV_RSLS.cols() << "\n R_RS: " << _RSLSreusedTimeWindows << "\n"
638 << '\n';
639 }
640
641 // ------------ RESTART ZERO ------------
643 // drop all stored Wtil^q, Z^q matrices
646
647 PRECICE_DEBUG("MVJ-RESTART, mode=Zero");
648
650
651 // re-compute Wtil -- compensate for dropping of Wtil_0 and Z_0:
652 // Wtil_q <-- Wtil_q + Wtil^0 * (Z^0*V_q)
653 for (int i = static_cast<int>(_WtilChunk.size()) - 1; i >= 1; i--) {
654
655 int colsLSSystemBackThen = _pseudoInverseChunk.front().rows();
656 PRECICE_ASSERT(colsLSSystemBackThen == _WtilChunk.front().cols(), colsLSSystemBackThen, _WtilChunk.front().cols());
657 Eigen::MatrixXd ZV = Eigen::MatrixXd::Zero(colsLSSystemBackThen, _qrV.cols());
658 // multiply: ZV := Z^q * V of size (m x m) with m=#cols, stored on each proc.
659 _parMatrixOps->multiply(_pseudoInverseChunk.front(), _matrixV, ZV, colsLSSystemBackThen, getLSSystemRows(), _qrV.cols());
660 // multiply: Wtil^0 * (Z_0*V) dimensions: (n x m) * (m x m), fully local and embarrassingly parallel
661 Eigen::MatrixXd tmp = Eigen::MatrixXd::Zero(_qrV.rows(), _qrV.cols());
662 tmp = _WtilChunk.front() * ZV;
663 _WtilChunk[i] += tmp;
664
665 // drop oldest pair Wtil_0 and Z_0
670 }
671
673 PRECICE_ASSERT(false); // should not happen, in this case _imvjRestart=false
674 } else {
675 PRECICE_ASSERT(false);
676 }
677}
678
679// ==================================================================================
681 const DataMap &cplData)
682{
684
685 // truncate V_RSLS and W_RSLS matrices according to _RSLSreusedTimeWindows
686 if (_imvjRestartType == RS_LS) {
687 if (_matrixCols_RSLS.front() == 0) { // Did only one iteration
689 }
690 if (_RSLSreusedTimeWindows == 0) {
691 _matrixV_RSLS.resize(0, 0);
692 _matrixW_RSLS.resize(0, 0);
694 } else if (static_cast<int>(_matrixCols_RSLS.size()) > _RSLSreusedTimeWindows) {
695 int toRemove = _matrixCols_RSLS.back();
696 PRECICE_ASSERT(toRemove > 0, toRemove);
697 if (_matrixV_RSLS.size() > 0) {
698 PRECICE_ASSERT(_matrixV_RSLS.cols() > toRemove, _matrixV_RSLS.cols(), toRemove);
699 }
700
701 // remove columns
702 for (int i = 0; i < toRemove; i++) {
705 }
707 }
709 }
710
711 //_info2<<'\n';
712
713 // if efficient update of imvj is enabled
715 // need to apply the preconditioner, as all data structures are reverted after
716 // call to computeQNUpdate. Need to call this before the preconditioner is updated.
717
718 // |= REBUILD QR-dec if needed ============|
719 // apply scaling to V, V' := P * V (only needed to reset the QR-dec of V)
721
722 if (_preconditioner->requireNewQR()) {
723 if (not(_filter == Acceleration::QR2FILTER)) { // for QR2 filter, there is no need to do this twice
725 }
726 _preconditioner->newQRfulfilled();
727 }
728 // apply the configured filter to the LS system
729 // as it changed in BaseQNAcceleration::iterationsConverged()
731 _preconditioner->revert(_matrixV);
732 // |=================== ============|
733
734 // ------- RESTART/ JACOBIAN ASSEMBLY -------
735 if (_imvjRestart) {
736
737 // add the matrices Wtil and Z of the converged configuration to the storage containers
738 Eigen::MatrixXd Z(_qrV.cols(), _qrV.rows());
739 // compute pseudo inverse using QR factorization and back-substitution
740 // also compensates for the scaling of V, i.e.,
741 // reverts Z' = R^-1 * Q^T * P^-1 as Z := Z' * P
742 pseudoInverse(Z);
743
744 // push back unscaled pseudo Inverse, Wtil is also unscaled.
745 // all objects in Wtil chunk and Z chunk are NOT PRECONDITIONED
748
752 if (static_cast<int>(_WtilChunk.size()) >= _chunkSize + 1) {
753
754 // < RESTART >
755 _nbRestarts++;
756 profiling::Event restartUpdate("IMVJRestart");
757 restartIMVJ();
758 restartUpdate.stop();
759 }
760
761 // only in imvj normal mode with efficient update:
762 } else {
763
764 // compute explicit representation of Jacobian
766 }
767
773 //_Wtil.conservativeResize(0, 0);
774 _resetLS = true;
775 }
776 }
777
778 // only store Jacobian if imvj is in normal mode, i.e., the Jacobian is build
779 if (not _imvjRestart) {
780 // store inverse Jacobian from converged time window. NOT SCALED with preconditioner
782 }
783}
784
785// ==================================================================================
787 int columnIndex)
788{
789 PRECICE_TRACE(columnIndex, _matrixV.cols());
790 PRECICE_ASSERT(_matrixV.cols() > 1, _matrixV.cols());
791 PRECICE_ASSERT(_Wtil.cols() > 1);
792
793 // remove column from matrix _Wtil
794 if (not _resetLS && not _alwaysBuildJacobian)
796
798}
799
800// ==================================================================================
802 int columnIndex)
803{
804 PRECICE_TRACE(columnIndex, _matrixV_RSLS.cols());
805 PRECICE_ASSERT(_matrixV_RSLS.cols() > 1);
806
809
810 // Reduce column count
811 std::deque<int>::iterator iter = _matrixCols_RSLS.begin();
812 int cols = 0;
813 while (iter != _matrixCols_RSLS.end()) {
814 cols += *iter;
815 if (cols > columnIndex) {
816 PRECICE_ASSERT(*iter > 0);
817 *iter -= 1;
818 if (*iter == 0) {
820 }
821 break;
822 }
823 iter++;
824 }
825}
826
827} // namespace precice::acceleration
828
829#endif
#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
T at(T... args)
T back(T... args)
T begin(T... args)
Base Class for quasi-Newton acceleration schemes.
virtual void updateDifferenceMatrices(const DataMap &cplData)
Updates the V, W matrices (as well as the matrices for the secondary data)
int _filter
filter method that is used to maintain good conditioning of the least-squares system Either of two ty...
double _initialRelaxation
Constant relaxation factor used for first iteration.
int _timeWindowsReused
Maximum number of old time windows (with data values) kept.
virtual int getLSSystemCols() const
: computes number of cols in least squares system, i.e, number of cols in _matrixV,...
virtual void applyFilter()
Applies the filter method for the least-squares system, defined in the configuration.
std::ostringstream _infostringstream
write some debug/acceleration info to file
bool _resetLS
If true, the LS system has been modified (reset or recomputed) in such a way, that mere updating of m...
Eigen::MatrixXd _matrixV
Stores residual deltas.
std::vector< int > _dimOffsets
Stores the local dimensions, i.e., the offsets in _invJacobian for all processors.
Eigen::VectorXd _residuals
Current iteration residuals of IQN data. Temporary.
virtual void removeMatrixColumn(int columnIndex)
Removes one iteration from V,W matrices and adapts _matrixCols.
Eigen::MatrixXd _matrixW
Stores x tilde deltas, where x tilde are values computed by solvers.
double _singularityLimit
Determines sensitivity when two matrix columns are considered equal.
bool _firstIteration
Indicates the first iteration, where constant relaxation is used.
std::map< int, Eigen::VectorXd > _secondaryResiduals
Current iteration residuals of secondary data.
impl::PtrPreconditioner _preconditioner
Preconditioner for least-squares system if vectorial system is used.
impl::QRFactorization _qrV
Stores the current QR decomposition ov _matrixV, can be updated via deletion/insertion of columns.
int _maxIterationsUsed
Maximum number of old data iterations kept.
std::vector< int > _secondaryDataIDs
Data IDs of data not involved in IQN coefficient computation.
virtual void initialize(const DataMap &cplData)
Initializes the acceleration.
int _RSLSreusedTimeWindows
: Number of reused time windows at restart if restart-mode = RS-LS
void computeNewtonUpdateEfficient(const DataMap &cplData, Eigen::VectorXd &update)
: computes the quasi-Newton update vector based on the same numerics as above. However,...
bool _imvjRestart
: If true, the imvj method is used with the restart chunk based approach that avoids to explicitly bu...
virtual void updateDifferenceMatrices(const DataMap &cplData)
: updates the V, W matrices (as well as the matrices for the secondary data)
void pseudoInverse(Eigen::MatrixXd &pseudoInverse)
: computes the pseudo inverse of V multiplied with V^T, i.e., Z = (V^TV)^-1V^T via QR-dec
virtual void initialize(const DataMap &cplData)
Initializes the acceleration.
virtual void computeQNUpdate(const DataMap &cplData, Eigen::VectorXd &xUpdate)
: computes the IQNIMVJ update using QR decomposition of V, furthermore it updates the inverse of the ...
Eigen::MatrixXd _matrixW_RSLS
stores columns from previous _RSLSreusedTimeWindows time windows if RS-LS restart-mode is active
IQNIMVJAcceleration(double initialRelaxation, bool forceInitialRelaxation, int maxIterationsUsed, int pastTimeWindowsReused, int filter, double singularityLimit, std::vector< int > dataIDs, const impl::PtrPreconditioner &preconditioner, bool alwaysBuildJacobian, int imvjRestartType, int chunkSize, int RSLSreusedTimeWindows, double RSSVDtruncationEps)
Constructor.
int _nbRestarts
tracks the number of restarts of IMVJ
int _imvjRestartType
: Indicates the type of the imvj restart-mode:
Eigen::MatrixXd _invJacobian
stores the approximation of the inverse Jacobian of the system at current time window.
virtual void computeUnderrelaxationSecondaryData(const DataMap &cplData)
: computes underrelaxation for the secondary data
virtual void removeMatrixColumn(int columnIndex)
: Removes one iteration from V,W matrices and adapts _matrixCols.
impl::PtrParMatrixOps _parMatrixOps
encapsulates matrix-matrix and matrix-vector multiplications for serial and parallel execution
int _chunkSize
: Number of time windows between restarts for the imvj method in restart mode
Eigen::MatrixXd _Wtil
stores the sub result (W-J_prev*V) for the current iteration
void removeMatrixColumnRSLS(int columnINdex)
: Removes one column form the V_RSLS and W_RSLS matrices and adapts _matrixCols_RSLS
void buildWtil()
: re-computes the matrix _Wtil = ( W - J_prev * V) instead of updating it according to V
bool _alwaysBuildJacobian
If true, the less efficient method to compute the quasi-Newton update is used, that explicitly builds...
virtual void specializedIterationsConverged(const DataMap &cplData)
Marks a iteration sequence as converged.
void restartIMVJ()
: restarts the imvj method, i.e., drops all stored matrices Wtil and Z and computes a initial guess o...
void buildJacobian()
: computes a explicit representation of the Jacobian, i.e., n x n matrix
std::vector< Eigen::MatrixXd > _pseudoInverseChunk
stores all pseudo inverses within the current chunk of the imvj restart mode, disabled if _imvjRestar...
Eigen::MatrixXd _oldInvJacobian
stores the approximation of the inverse Jacobian from the previous time window.
Eigen::MatrixXd _matrixV_RSLS
stores columns from previous _RSLSreusedTimeWindows time windows if RS-LS restart-mode is active
virtual ~IQNIMVJAcceleration()
Destructor, empty.
std::deque< int > _matrixCols_RSLS
number of cols per time window
std::vector< Eigen::MatrixXd > _WtilChunk
stores all Wtil matrices within the current chunk of the imvj restart mode, disabled if _imvjRestart ...
impl::SVDFactorization _svdJ
holds and maintains a truncated SVD decomposition of the Jacobian matrix
void computeNewtonUpdate(const DataMap &cplData, Eigen::VectorXd &update)
: computes the quasi-Newton update vector based on the matrices V and W using a QR decomposition of V...
Class that provides functionality for a dynamic QR-decomposition, that can be updated in O(mn) flops ...
Eigen::MatrixXd & matrixR()
returns a matrix representation of the upper triangular matrix R
void applyFilter(double singularityLimit, std::vector< int > &delIndices, Eigen::MatrixXd &V)
filters the least squares system, i.e., the decomposition Q*R = V according to the defined filter tec...
void pushBack(const Eigen::VectorXd &v)
inserts a new column at position _cols-1, i.e., appends a column at the end and updates the QR factor...
void reset()
resets the QR factorization zo zero Q(0:0, 0:0)R(0:0, 0:0)
Eigen::MatrixXd & matrixQ()
returns a matrix representation of the orthogonal matrix Q
double getThreshold()
: returns the truncation threshold for the SVD
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
Rank rank()
: returns the rank of the truncated SVD factorization
Matrix & matrixPsi()
: returns a matrix representation of the orthogonal matrix Psi, A = Psi * Sigma * Phi^bs
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 getWaste()
: returns the total number of truncated modes since last call to this method
void stop()
Stops a running event.
Definition Event.cpp:37
static bool isPrimary()
True if this process is running the primary rank.
Definition IntraComm.cpp:52
static bool isParallel()
True if this process is running in parallel.
Definition IntraComm.cpp:62
T clear(T... args)
T empty(T... args)
T end(T... args)
T erase(T... args)
T front(T... args)
contains implementations of acceleration schemes.
void removeColumnFromMatrix(Eigen::MatrixXd &A, int col)
void appendFront(Eigen::MatrixXd &A, Eigen::VectorXd &v)
void shiftSetFirst(Eigen::MatrixXd &A, const Eigen::VectorXd &v)
int Rank
Definition Types.hpp:37
STL namespace.
T pop_back(T... args)
T pop_front(T... args)
T push_back(T... args)
T push_front(T... args)
T size(T... args)