preCICE v3.1.1
Loading...
Searching...
No Matches
QRFactorization.cpp
Go to the documentation of this file.
1
2#include <Eigen/Core>
3#include <algorithm> // std::sort
4#include <cmath>
5#include <cstddef>
6#include <iostream>
7#include <memory>
8#include <utility>
9#include <vector>
10
13#include "com/Communication.hpp"
14#include "com/SharedPointer.hpp"
15#include "logging/LogMacros.hpp"
17#include "utils/IntraComm.hpp"
18#include "utils/assertion.hpp"
19
21
23 Eigen::MatrixXd Q,
24 Eigen::MatrixXd R,
25 int rows,
26 int cols,
27 int filter,
28 double omega,
29 double theta,
30 double sigma)
31
32 : _Q(std::move(Q)),
33 _R(std::move(R)),
34 _rows(rows),
35 _cols(cols),
36 _filter(filter),
37 _omega(omega),
38 _theta(theta),
39 _sigma(sigma),
40 _infostream(),
41 _fstream_set(false),
42 _globalRows(rows)
43{
44 PRECICE_ASSERT(_R.rows() == _cols, _R.rows(), _cols);
45 PRECICE_ASSERT(_R.cols() == _cols, _R.cols(), _cols);
46 PRECICE_ASSERT(_Q.cols() == _cols, _Q.cols(), _cols);
47 PRECICE_ASSERT(_Q.rows() == _rows, _Q.rows(), _rows);
48}
49
54 Eigen::MatrixXd A,
55 int filter,
56 double omega,
57 double theta,
58 double sigma)
59 : _Q(),
60 _R(),
61 _rows(A.rows()),
62 _cols(0),
63 _filter(filter),
64 _omega(omega),
65 _theta(theta),
66 _sigma(sigma),
67 _infostream(),
68 _fstream_set(false),
69 _globalRows(A.rows())
70{
71 int m = A.cols();
72 for (int k = 0; k < m; k++) {
73 Eigen::VectorXd v = A.col(k);
74 insertColumn(k, v);
75 }
76 //PRECICE_ASSERT(_R.rows() == _cols, _R.rows(), _cols);
77 PRECICE_ASSERT(_R.cols() == _cols, _R.cols(), _cols);
78 PRECICE_ASSERT(_Q.cols() == _cols, _Q.cols(), _cols);
79 PRECICE_ASSERT(_Q.rows() == _rows, _Q.rows(), _rows);
80 PRECICE_ASSERT(_cols == m, _cols, m);
81}
82
87 int filter,
88 double omega,
89 double theta,
90 double sigma)
91 : _Q(),
92 _R(),
93 _rows(0),
94 _cols(0),
95 _filter(filter),
96 _omega(omega),
97 _theta(theta),
98 _sigma(sigma),
99 _infostream(),
100 _fstream_set(false),
101 _globalRows(0)
102{
103}
104
105void QRFactorization::applyFilter(double singularityLimit, std::vector<int> &delIndices, Eigen::MatrixXd &V)
106{
108 delIndices.resize(0);
110 bool linearDependence = true;
111 std::vector<int> delFlag(_cols, 0);
112 int delCols = 0;
113 while (linearDependence) {
114 linearDependence = false;
115 int index = 0; // actual index of checked column, \in [0, _cols] and _cols is decreasing
116 if (_cols > 1) {
117 for (size_t i = 0; i < delFlag.size(); i++) {
118 // index is not incremented, if columns has been deleted in previous rounds
119 if (delFlag[i] > 0)
120 continue;
121
122 // QR1-filter
123 if (index >= cols())
124 break;
126 double factor = (_filter == Acceleration::QR1FILTER_ABS) ? 1.0 : _R.norm();
127 if (std::fabs(_R(index, index)) < singularityLimit * factor) {
128
129 linearDependence = true;
131 delFlag[i]++;
132 delIndices.push_back(i);
133 delCols++;
134 //break;
135 index--; // check same column index, as cols are shifted left
136 }
137 PRECICE_ASSERT(delCols + _cols == (int) delFlag.size(), (delCols + _cols), delFlag.size());
138 index++;
139 }
140 }
141 }
142 } else if (_filter == Acceleration::QR2FILTER) {
143 _Q.resize(0, 0);
144 _R.resize(0, 0);
145 _cols = 0;
146 _rows = V.rows();
147 // starting with the most recent input/output information, i.e., the latest column
148 // which is at position 0 in _matrixV (latest information is never filtered out!)
149 for (int k = 0; k < V.cols(); k++) {
150 Eigen::VectorXd v = V.col(k);
151 // this is the same as pushBack(v) as _cols grows within the insertion process
152 bool inserted = insertColumn(_cols, v, singularityLimit);
153 if (!inserted) {
154 delIndices.push_back(k);
155 }
156 }
157 }
158 std::sort(delIndices.begin(), delIndices.end());
159}
160
166{
167
169
170 PRECICE_ASSERT(k >= 0, k);
171 PRECICE_ASSERT(k < _cols, k, _cols);
172
173 // maintain decomposition and orthogonalization by application of givens rotations
174
175 for (int l = k; l < _cols - 1; l++) {
177 computeReflector(grot, _R(l, l + 1), _R(l + 1, l + 1));
178 Eigen::VectorXd Rr1 = _R.row(l);
179 Eigen::VectorXd Rr2 = _R.row(l + 1);
180 applyReflector(grot, l + 2, _cols, Rr1, Rr2);
181 _R.row(l) = Rr1;
182 _R.row(l + 1) = Rr2;
183 Eigen::VectorXd Qc1 = _Q.col(l);
184 Eigen::VectorXd Qc2 = _Q.col(l + 1);
185 applyReflector(grot, 0, _rows, Qc1, Qc2);
186 _Q.col(l) = Qc1;
187 _Q.col(l + 1) = Qc2;
188 }
189 // copy values and resize R and Q
190 for (int j = k; j < _cols - 1; j++) {
191 for (int i = 0; i <= j; i++) {
192 _R(i, j) = _R(i, j + 1);
193 }
194 }
195 _R.conservativeResize(_cols - 1, _cols - 1);
196 //_Q.conservativeResize(Eigen::NoChange_t, _cols-1);
197 _Q.conservativeResize(_rows, _cols - 1);
198 _cols--;
199
200 PRECICE_ASSERT(_Q.cols() == _cols, _Q.cols(), _cols);
201 PRECICE_ASSERT(_Q.rows() == _rows, _Q.rows(), _rows);
202 PRECICE_ASSERT(_R.cols() == _cols, _R.cols(), _cols);
203 //PRECICE_ASSERT(_R.rows() == _cols, _Q.rows(), _cols);
204}
205
206// ATTENTION: This method works on the memory of vector v, thus changes the vector v.
207bool QRFactorization::insertColumn(int k, const Eigen::VectorXd &vec, double singularityLimit)
208{
209 PRECICE_TRACE(k);
210
211 Eigen::VectorXd v(vec);
212
213 if (_cols == 0)
214 _rows = v.size();
215
216 bool applyFilter = (singularityLimit > 0.0);
217
218 PRECICE_ASSERT(k >= 0, k);
219 PRECICE_ASSERT(k <= _cols, k, _cols);
220 PRECICE_ASSERT(v.size() == _rows, v.size(), _rows);
221
222 _cols++;
223
224 // orthogonalize v to columns of Q
225 Eigen::VectorXd u(_cols);
226 double rho_orth = 0., rho0 = 0.;
227 if (applyFilter)
228 rho0 = utils::IntraComm::l2norm(v);
229
230 int err = orthogonalize(v, u, rho_orth, _cols - 1);
231
232 // on of the following is true
233 // - either ||v_orth|| / ||v|| <= 0.7 was true and the re-orthogonalization process failed 4 times
234 // - or the system is quadratic, thus no further column can be orthogonalized (||v_orth|| = rho_orth = 0)
235 // - or ||v_orth|| = rho_orth extremely small (almost zero), thus the column v is not very orthogonal to Q, discard
236 if (rho_orth <= std::numeric_limits<double>::min() || err < 0) {
237 PRECICE_DEBUG("The ratio ||v_orth|| / ||v|| is extremely small and either the orthogonalization process of column v failed or the system is quadratic.");
238
239 // necessary for applyFilter with the QR-2 filer. In this case, the new column is not inserted, but discarded.
240 _cols--;
241 return false;
242 }
243
244 // QR2-filter based on the new information added to the orthogonal system.
245 // if the new column incorporates less new information to the system than a
246 // prescribed threshold, the column is discarded
247 // rho_orth: the norm of the orthogonalized (but not normalized) column
248 // rho0: the norm of the initial column that is to be inserted
249 if (applyFilter && (rho0 * singularityLimit > rho_orth)) {
250 PRECICE_DEBUG("discarding column as it is filtered out by the QR2-filter: rho0*eps > rho_orth: {} > {}", rho0 * singularityLimit, rho_orth);
251 _cols--;
252 return false;
253 }
254
255 // resize R(1:m, 1:m) -> R(1:m+1, 1:m+1)
256 _R.conservativeResize(_cols, _cols);
257 _R.col(_cols - 1) = Eigen::VectorXd::Zero(_cols);
258 _R.row(_cols - 1) = Eigen::VectorXd::Zero(_cols);
259
260 for (int j = _cols - 2; j >= k; j--) {
261 for (int i = 0; i <= j; i++) {
262 _R(i, j + 1) = _R(i, j);
263 }
264 }
265
266 for (int j = k + 1; j < _cols; j++) {
267 _R(j, j) = 0.;
268 }
269
270 PRECICE_ASSERT(_R.cols() == _cols, _R.cols(), _cols);
271 //PRECICE_ASSERT(_R.rows() == _cols, _R.rows(), _cols);
272
273 // resize Q(1:n, 1:m) -> Q(1:n, 1:m+1)
274 _Q.conservativeResize(_rows, _cols);
275 _Q.col(_cols - 1) = v;
276
277 PRECICE_ASSERT(_Q.cols() == _cols, _Q.cols(), _cols);
278 PRECICE_ASSERT(_Q.rows() == _rows, _Q.rows(), _rows);
279
280 // maintain decomposition and orthogonalization by application of givens rotations
281 for (int l = _cols - 2; l >= k; l--) {
283 computeReflector(grot, u(l), u(l + 1));
284 Eigen::VectorXd Rr1 = _R.row(l);
285 Eigen::VectorXd Rr2 = _R.row(l + 1);
286 applyReflector(grot, l + 1, _cols, Rr1, Rr2);
287 _R.row(l) = Rr1;
288 _R.row(l + 1) = Rr2;
289 Eigen::VectorXd Qc1 = _Q.col(l);
290 Eigen::VectorXd Qc2 = _Q.col(l + 1);
291 applyReflector(grot, 0, _rows, Qc1, Qc2);
292 _Q.col(l) = Qc1;
293 _Q.col(l + 1) = Qc2;
294 }
295 for (int i = 0; i <= k; i++) {
296 _R(i, k) = u(i);
297 }
298
299 return true;
300}
301
314 Eigen::VectorXd &v,
315 Eigen::VectorXd &r,
316 double & rho,
317 int colNum)
318{
320
323 }
324
325 bool null = false;
326 bool termination = false;
327 double rho0 = 0., rho1 = 0.;
328 Eigen::VectorXd u = Eigen::VectorXd::Zero(_rows);
329 Eigen::VectorXd s = Eigen::VectorXd::Zero(colNum);
330 r = Eigen::VectorXd::Zero(_cols);
331
332 rho = utils::IntraComm::l2norm(v); // distributed l2norm
333 rho0 = rho;
334 int k = 0;
335 while (!termination) {
336
337 // take a gram-schmidt iteration
338 u = Eigen::VectorXd::Zero(_rows);
339 for (int j = 0; j < colNum; j++) {
340
341 // dot-product <_Q(:,j), v >
342 Eigen::VectorXd Qc = _Q.col(j);
343
344 // dot product <_Q(:,j), v> =: r_ij
345 double r_ij = utils::IntraComm::dot(Qc, v);
346 // save r_ij in s(j) = column of R
347 s(j) = r_ij;
348 // u is the sum of projections r_ij * _Q(:,j) = _Q(:,j) * <_Q(:,j), v>
349 u += _Q.col(j) * r_ij;
350 }
351 // add the furier coefficients over all orthogonalize iterations
352 for (int j = 0; j < colNum; j++) {
353 r(j) = r(j) + s(j);
354 }
355 // subtract projections from v, v is now orthogonal to columns of _Q
356 for (int i = 0; i < _rows; i++) {
357 v(i) = v(i) - u(i);
358 }
359 // rho1 = norm of orthogonalized new column v_tilde (though not normalized)
360 rho1 = utils::IntraComm::l2norm(v); // distributed l2norm
361
362 // t = norm of r_(:,j) with j = colNum-1
363 double norm_coefficients = utils::IntraComm::l2norm(s); // distributed l2norm
364 k++;
365
366 // treat the special case m=n
367 // Attention (intra-participant communication): Here, we need to compare the global _rows with colNum and NOT the local
368 // rows on the processor.
369 if (_globalRows == colNum) {
370 PRECICE_WARN("The least-squares system matrix is quadratic, i.e., the new column cannot be orthogonalized (and thus inserted) to the LS-system.\nOld columns need to be removed.");
371 v = Eigen::VectorXd::Zero(_rows);
372 rho = 0.;
373 return k;
374 }
375
376 // take correct action if v_orth is null
377 if (rho1 <= std::numeric_limits<double>::min()) {
378 PRECICE_DEBUG("The norm of v_orthogonal is almost zero, i.e., failed to orthogonalize column v; discard.");
379 null = true;
380 rho1 = 1;
381 termination = true;
382 }
383
396 if (rho1 * _theta <= rho0 + _omega * norm_coefficients) {
397 // exit to fail if too many iterations
398 if (k >= 4) {
399 PRECICE_WARN("Matrix Q is not sufficiently orthogonal. Failed to orthogonalize new column after 4 iterations. New column will be discarded. The least-squares system is very bad conditioned and the quasi-Newton will most probably fail to converge.");
400 return -1;
401 }
402 rho0 = rho1;
403
404 // termination, i.e., (rho0 + _omega * t < _theta *rho1)
405 } else {
406 termination = true;
407 }
408 }
409
410 // normalize v
411 v /= rho1;
412 rho = null ? 0 : rho1;
413 r(colNum) = rho;
414 return k;
415}
416
429 Eigen::VectorXd &v,
430 Eigen::VectorXd &r,
431 double & rho,
432 int colNum)
433{
435
436 // serial case
439 }
440
441 bool restart = false;
442 bool null = false;
443 bool termination = false;
444 double rho0 = 0., rho1 = 0.;
445 double t = 0;
446 Eigen::VectorXd u = Eigen::VectorXd::Zero(_rows);
447 Eigen::VectorXd s = Eigen::VectorXd::Zero(colNum);
448 r = Eigen::VectorXd::Zero(_cols);
449
450 rho = utils::IntraComm::l2norm(v); // distributed l2norm
451 rho0 = rho;
452 int k = 0;
453 while (!termination) {
454 // take a gram-schmidt iteration, ignoring r on later steps if previous v was null
455 u = Eigen::VectorXd::Zero(_rows);
456 for (int j = 0; j < colNum; j++) {
457
458 /*
459 * dot-product <_Q(:,j), v >
460 */
461 Eigen::VectorXd Qc = _Q.col(j);
462
463 // dot product <_Q(:,j), v> =: r_ij
464 double ss = utils::IntraComm::dot(Qc, v);
465 t = ss;
466 // save r_ij in s(j) = column of R
467 s(j) = t;
468 // u is the sum of projections r_ij * _Q(i,:) = _Q(i,:) * <_Q(:,j), v>
469 for (int i = 0; i < _rows; i++) {
470 u(i) = u(i) + _Q(i, j) * t;
471 }
472 }
473 if (!null) {
474 // add over all runs: r_ij = r_ij_prev + r_ij
475 for (int j = 0; j < colNum; j++) {
476 r(j) = r(j) + s(j);
477 }
478 }
479 // subtract projections from v, v is now orthogonal to columns of _Q
480 for (int i = 0; i < _rows; i++) {
481 v(i) = v(i) - u(i);
482 }
483 // rho1 = norm of orthogonalized new column v_tilde (though not normalized)
484 rho1 = utils::IntraComm::l2norm(v); // distributed l2norm
485
486 // t = norm of r_(:,j) with j = colNum-1
487 t = utils::IntraComm::l2norm(s); // distributed l2norm
488 k++;
489
490 // treat the special case m=n
491 // Attention (intra-participant communication): Here, we need to compare the global _rows with colNum and NOT the local
492 // rows on the processor.
493 if (_globalRows == colNum) {
494 PRECICE_WARN("The least-squares system matrix is quadratic, i.e., the new column cannot be orthogonalized (and thus inserted) to the LS-system.\nOld columns need to be removed.");
495 v = Eigen::VectorXd::Zero(_rows);
496 rho = 0.;
497 return k;
498 }
499
512 if (rho1 * _theta <= rho0 + _omega * t) {
513 // exit to fail if too many iterations
514 if (k >= 4) {
516 << "\ntoo many iterations in orthogonalize, termination failed\n";
517 PRECICE_WARN("Matrix Q is not sufficiently orthogonal. Failed to orthogonalize new column after 4 iterations. New column will be discarded. The least-squares system is very bad conditioned and the quasi-Newton will most probably fail to converge.");
518 return -1;
519 }
520
521 // if ||v_orth|| / ||v|| is extremely small (numeric limit)
522 // discard information from column, use any unit vector orthogonal to Q
523 if (!restart && rho1 <= rho * _sigma) {
524 PRECICE_WARN("The new column is in the range of Q, thus not possible to orthogonalize. Try to insert a unit vector that is orthogonal to the columns space of Q.");
525 //PRECICE_DEBUG("[QR-dec] - reorthogonalization");
526 if (_fstream_set)
527 (*_infostream) << "[QR-dec] - reorthogonalization\n";
528
529 restart = true;
530
538 u = Eigen::VectorXd::Zero(_rows);
539 for (int j = 0; j < colNum; j++) {
540 for (int i = 0; i < _rows; i++) {
541 u(i) = u(i) + _Q(i, j) * _Q(i, j);
542 }
543 }
544 t = 2;
545
546 // ATTENTION: maybe in the following is something wrong @todo: fix this comment
547
548 for (int i = 0; i < _rows; i++) {
549 if (u(i) < t) {
550 k = i;
551 t = u(k);
552 }
553 }
554
555 int global_k = k;
556 int local_k = 0;
557 double local_uk = 0.;
558 double global_uk = 0.;
559 int rank = 0;
560
563 utils::IntraComm::getCommunication()->send(u(k), 0);
564 }
565
567 global_uk = u(k);
568 for (Rank secondaryRank : utils::IntraComm::allSecondaryRanks()) {
569 utils::IntraComm::getCommunication()->receive(local_k, secondaryRank);
570 utils::IntraComm::getCommunication()->receive(local_uk, secondaryRank);
571 if (local_uk < global_uk) {
572 rank = secondaryRank;
573 global_uk = local_uk;
574 global_k = local_k;
575 }
576 }
577 if (_fstream_set)
578 (*_infostream) << " global u(k):" << global_uk << ", global k: " << global_k << ", rank: " << rank << '\n';
579 }
580
581 // take correct action if v is null
582 if (rho1 == 0) {
583 null = true;
584 rho1 = 1;
585 }
586 // reinitialize v and k
587 v = Eigen::VectorXd::Zero(_rows);
588
589 // insert rho1 at position k with smallest u(i) = Q(i,:) * Q(i,:)
591 v(k) = rho1;
592 } else {
593 if (utils::IntraComm::getRank() == rank)
594 v(global_k) = rho1;
595 }
596 k = 0;
597 }
598 rho0 = rho1;
599
600 // termination, i.e., (rho0 + _omega * t < _theta *rho1)
601 } else {
602 termination = true;
603 }
604 }
605
606 // normalize v
607 v /= rho1;
608 rho = null ? 0 : rho1;
609 r(colNum) = rho;
610 return k;
611}
612
618 double & x,
619 double & y)
620{
621 double u = x;
622 double v = y;
623 if (v == 0) {
624 grot.sigma = 0;
625 grot.gamma = 1;
626 } else {
627 double mu = std::max(std::fabs(u), std::fabs(v));
628 double t = mu * std::sqrt(std::pow(u / mu, 2) + std::pow(v / mu, 2));
629 t *= (u < 0) ? -1 : 1;
630 grot.gamma = u / t;
631 grot.sigma = v / t;
632 x = t;
633 y = 0;
634 }
635}
636
642 const QRFactorization::givensRot &grot,
643 int k,
644 int l,
645 Eigen::VectorXd & p,
646 Eigen::VectorXd & q)
647{
648 double nu = grot.sigma / (1. + grot.gamma);
649 for (int j = k; j < l; j++) {
650 double u = p(j);
651 double v = q(j);
652 double t = u * grot.gamma + v * grot.sigma;
653 p(j) = t;
654 q(j) = (t + u) * nu - v;
655 }
656}
657
659{
660 _globalRows = gr;
661}
662
663Eigen::MatrixXd &QRFactorization::matrixQ()
664{
665 return _Q;
666}
667
668Eigen::MatrixXd &QRFactorization::matrixR()
669{
670 return _R;
671}
672
674{
675 return _cols;
676}
677
679{
680 return _rows;
681}
682
684{
685 _Q.resize(0, 0);
686 _R.resize(0, 0);
687 _cols = 0;
688 _rows = 0;
689 _globalRows = 0;
690}
691
693 Eigen::MatrixXd const &Q,
694 Eigen::MatrixXd const &R,
695 int rows,
696 int cols,
697 double omega,
698 double theta,
699 double sigma)
700{
701 _Q = Q;
702 _R = R;
703 _rows = rows;
704 _cols = cols;
705 _omega = omega;
706 _theta = theta;
707 _sigma = sigma;
709 PRECICE_ASSERT(_R.rows() == _cols, _R.rows(), _cols);
710 PRECICE_ASSERT(_R.cols() == _cols, _R.cols(), _cols);
711 PRECICE_ASSERT(_Q.cols() == _cols, _Q.cols(), _cols);
712 PRECICE_ASSERT(_Q.rows() == _rows, _Q.rows(), _rows);
713}
714
716 Eigen::MatrixXd const &A,
717 int globalRows,
718 double omega,
719 double theta,
720 double sigma)
721{
723 _Q.resize(0, 0);
724 _R.resize(0, 0);
725 _cols = 0;
726 _rows = A.rows();
727 _omega = omega;
728 _theta = theta;
729 _sigma = sigma;
730 _globalRows = globalRows;
731
732 int m = A.cols();
733 int col = 0, k = 0;
734 for (; col < m; k++, col++) {
735 Eigen::VectorXd v = A.col(col);
736 bool inserted = insertColumn(k, v);
737 if (not inserted) {
738 k--;
739 PRECICE_DEBUG("column {} has not been inserted in the QR-factorization, failed to orthogonalize.", col);
740 }
741 }
742 PRECICE_ASSERT(_R.rows() == _cols, _R.rows(), _cols);
743 PRECICE_ASSERT(_R.cols() == _cols, _R.cols(), _cols);
744 PRECICE_ASSERT(_Q.cols() == _cols, _Q.cols(), _cols);
745 PRECICE_ASSERT(_Q.rows() == _rows, _Q.rows(), _rows);
746 PRECICE_ASSERT(_cols == m, _cols, m);
747}
748
749void QRFactorization::pushFront(const Eigen::VectorXd &v)
750{
751 insertColumn(0, v);
752}
753
754void QRFactorization::pushBack(const Eigen::VectorXd &v)
755{
757}
758
763
765{
766 deleteColumn(_cols - 1);
767}
768
770{
771 _infostream = stream;
772 _fstream_set = true;
773}
774
776{
777 _filter = filter;
778}
779
780} // namespace precice::acceleration::impl
#define PRECICE_WARN(...)
Definition LogMacros.hpp:11
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
unsigned int index
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T begin(T... args)
int orthogonalize_stable(Eigen::VectorXd &v, Eigen::VectorXd &r, double &rho, int colNum)
assuming Q(1:n,1:m) has nearly orthonormal columns, this procedure orthogonlizes v(1:n) to the column...
bool insertColumn(int k, const Eigen::VectorXd &v, double singularityLimit=0)
inserts a new column at arbitrary position and updates the QR factorization This function works on th...
int orthogonalize(Eigen::VectorXd &v, Eigen::VectorXd &r, double &rho, int colNum)
assuming Q(1:n,1:m) has nearly orthonormal columns, this procedure orthogonlizes v(1:n) to the column...
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 computeReflector(givensRot &grot, double &x, double &y)
computes parameters for givens matrix G for which (x,y)G = (z,0). replaces (x,y) by (z,...
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)
void popBack()
deletes the column at position _cols-1, i.e., deletes the last column and updates the QR factorizatio...
void pushFront(const Eigen::VectorXd &v)
inserts a new column at position 0, i.e., shifts right and inserts at first position and updates the ...
QRFactorization(int filter=0, double omega=0, double theta=1./0.7, double sigma=std::numeric_limits< double >::min())
Constructor.
void applyReflector(const givensRot &grot, int k, int l, Eigen::VectorXd &p, Eigen::VectorXd &q)
this procedure replaces the two column matrix [p(k:l-1), q(k:l-1)] by [p(k:l), q(k:l)]*G,...
Eigen::MatrixXd & matrixQ()
returns a matrix representation of the orthogonal matrix Q
void popFront()
deletes the column at position 0, i.e., deletes and shifts columns to the left and updates the QR fac...
void deleteColumn(int k)
updates the factorization A=Q[1:n,1:m]R[1:m,1:n] when the kth column of A is deleted....
static double l2norm(const Eigen::VectorXd &vec)
The l2 norm of a vector is calculated on distributed data.
Definition IntraComm.cpp:67
static Rank getRank()
Current rank.
Definition IntraComm.cpp:42
static double dot(const Eigen::VectorXd &vec1, const Eigen::VectorXd &vec2)
static bool isPrimary()
True if this process is running the primary rank.
Definition IntraComm.cpp:52
static auto allSecondaryRanks()
Returns an iterable range over salve ranks [1, _size)
Definition IntraComm.hpp:37
static bool isParallel()
True if this process is running in parallel.
Definition IntraComm.cpp:62
static bool isSecondary()
True if this process is running a secondary rank.
Definition IntraComm.cpp:57
static com::PtrCommunication & getCommunication()
Intra-participant communication.
Definition IntraComm.hpp:31
T end(T... args)
T fabs(T... args)
T max(T... args)
int Rank
Definition Types.hpp:37
STL namespace.
T pow(T... args)
T push_back(T... args)
T resize(T... args)
T size(T... args)
T sort(T... args)
T sqrt(T... args)