QRMatrix.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2016 OpenFOAM Foundation
9  Copyright (C) 2019-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "QRMatrix.H"
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 template<class MatrixType>
35 (
36  MatrixType& A
37 )
38 {
39  const label nIter = min(A.m() - 1, A.n());
40 
41  // Loop through all subcolumns of which the diagonal elem is the first elem
42  for (label k = 0; k < nIter; ++k)
43  {
44  const RMatrix u(householderReflector(A.subColumn(k, k)));
45 
46  applyHouseholder(A, u, k);
47  }
48 
49  if (outputType_ == outputTypes::REDUCED_R)
50  {
51  A.resize(A.n(), A.n());
52  }
53 }
54 
55 
56 template<class MatrixType>
58 (
59  MatrixType& A
60 )
61 {
62  const label nCols = A.n();
63  const label nIter = min(A.m() - 1, A.n());
64 
65  // Initialise permutation vector, and column norm vector
66  P_ = identity(nCols);
67 
68  // Initialise vector norms of each column of A
69  List<scalar> colNorms(nCols);
70  for (label k = 0; k < nCols; ++k)
71  {
72  colNorms[k] = A.columnNorm(k, true);
73  }
74 
75  // Loop through all subcolumns of which the diagonal elem is the first elem
76  for (label k = 0; k < nIter; ++k)
77  {
78  const labelRange colRange(k, nCols);
79  const SubList<scalar> subColNorms(colNorms, colRange);
80 
81  // Column pivoting
82  const label maxColNormi =
84  (
85  subColNorms.cbegin(),
86  std::max_element(subColNorms.cbegin(), subColNorms.cend())
87  );
88 
89  // Swap R_, P_ and colNorms_ according to pivot column if the current
90  // column is not the max norm column by avoiding any column swap where
91  // the leading elem is very small
92  if (maxColNormi != 0 && SMALL < mag(A(k, k + maxColNormi)))
93  {
94  const RMatrix R1(A.subColumn(k));
95  const RMatrix R2(A.subColumn(maxColNormi + k));
96  A.subColumn(k) = R2;
97  A.subColumn(maxColNormi + k) = R1;
98 
99  Swap(P_[k], P_[maxColNormi + k]);
100  Swap(colNorms[k], colNorms[maxColNormi + k]);
101  }
102 
103  {
104  const RMatrix u(householderReflector(A.subColumn(k, k)));
105 
106  applyHouseholder(A, u, k);
107  }
108 
109  // Update norms
110  if (k < nIter - 1)
111  {
112  label q = k + 1;
113  for (const auto& val : RMatrix(A.subRow(k, q)))
114  {
115  colNorms[q] -= magSqr(val);
116  ++q;
117  }
118  }
119  }
120 
121  if (outputType_ == outputTypes::REDUCED_R)
122  {
123  A.resize(A.n(), A.n());
124  }
125 }
126 
127 
128 template<class MatrixType>
129 template<template<typename> class ListContainer>
131 (
132  ListContainer<cmptType>& x
133 ) const
134 {
135  const label n = R_.n();
136 
137  for (label i = n - 1; 0 <= i; --i)
138  {
139  cmptType sum = x[i];
140 
141  for (label j = i + 1; j < n; ++j)
142  {
143  sum -= x[j]*R_(i, j);
144  }
145 
146  if (mag(R_(i, i)) < SMALL)
147  {
149  << "Back-substitution failed due to small diagonal"
150  << abort(FatalError);
151  }
152 
153  x[i] = sum/R_(i, i);
154  }
155 }
156 
157 
158 template<class MatrixType>
159 template<template<typename> class ListContainer>
161 (
162  List<cmptType>& x,
163  const ListContainer<cmptType>& source
164 ) const
165 {
166  // Assert (&x != &source) ?
167  const label m = Q_.m();
168 
169  // x = Q_.T()*source; i.e., Q_.Tmul(source)
170  for (label i = 0; i < m; ++i)
171  {
172  x[i] = 0;
173  for (label j = 0; j < m; ++j)
174  {
175  x[i] += Q_(j, i)*source[j];
176  }
177  }
178 
179  solvex(x);
180 }
181 
182 
183 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
184 
185 template<class MatrixType>
187 :
188  outputType_(),
189  storeMethod_(),
190  colPivot_()
191 {}
192 
193 
194 template<class MatrixType>
196 (
197  const outputTypes outputType,
198  const storeMethods storeMethod,
199  const colPivoting colPivot
200 )
201 :
202  outputType_(outputType),
203  storeMethod_(storeMethod),
204  colPivot_(colPivot),
205  Q_(),
206  R_(),
207  P_()
208 {}
209 
210 
211 template<class MatrixType>
213 (
214  MatrixType& A,
215  const outputTypes outputType,
216  const storeMethods storeMethod,
217  const colPivoting colPivot
218 )
219 :
220  outputType_(outputType),
221  storeMethod_(storeMethod),
222  colPivot_(colPivot),
223  Q_(),
224  R_(),
225  P_()
226 {
227  decompose(A);
228 }
229 
230 
231 template<class MatrixType>
233 (
234  const MatrixType& A,
235  const outputTypes outputType,
236  const colPivoting colPivot
237 )
238 :
239  outputType_(outputType),
240  storeMethod_(storeMethods::OUT_OF_PLACE),
241  colPivot_(colPivot),
242  Q_(),
243  R_(),
244  P_()
245 {
246  decompose(A);
247 }
248 
249 
250 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
251 
252 template<class MatrixType>
254 (
255  MatrixType& A
256 )
257 {
258  // Check whether settings and input are consistent for reduced QRMatrix
259  if (A.m() <= A.n() && outputType_ == outputTypes::REDUCED_R)
260  {
261  outputType_ = outputTypes::FULL_R;
262 
263  #if FULLDEBUG
265  << "Reduced QR decomposition is by definition limited to matrices"
266  << " wherein rows > columns, thus computing FULL decomposition."
267  << nl;
268  #endif
269  }
270 
271  // Allocate resources for Q_, if need be
272  if (outputType_ == outputTypes::FULL_QR)
273  {
274  Q_ = MatrixType({A.m(), A.m()}, I);
275  }
276 
277  // Allocate resources for R_ and execute decomposition
278  switch (storeMethod_)
279  {
280  case storeMethods::IN_PLACE:
281  {
282  if (colPivot_)
283  {
284  qrPivot(A);
285  }
286  else
287  {
288  qr(A);
289  }
290  break;
291  }
292 
293  case storeMethods::OUT_OF_PLACE:
294  {
295  R_ = A;
296 
297  if (colPivot_)
298  {
299  qrPivot(R_);
300  }
301  else
302  {
303  qr(R_);
304  }
305  break;
306  }
307  }
308 }
309 
310 
311 template<class MatrixType>
313 (
314  const MatrixType& A
315 )
316 {
317  if (storeMethod_ == storeMethods::IN_PLACE)
318  {
320  << "const type qualifier invalidates storeMethods::IN_PLACE." << nl;
321  }
322 
323  // Check whether settings and input are consistent for reduced QRMatrix
324  if (A.m() <= A.n() && outputType_ == outputTypes::REDUCED_R)
325  {
326  outputType_ = outputTypes::FULL_R;
327 
328  #if FULLDEBUG
330  << "Reduced QR decomposition is by definition limited to matrices"
331  << " wherein rows > columns, thus computing FULL decomposition."
332  << nl;
333  #endif
334  }
335 
336  // Allocate resources for Q_, if need be
337  if (outputType_ == outputTypes::FULL_QR)
338  {
339  Q_ = MatrixType({A.m(), A.m()}, I);
340  }
341 
342  // Allocate resources for R_ and execute decomposition
343  R_ = A;
344 
345  if (colPivot_)
346  {
347  qrPivot(R_);
348  }
349  else
350  {
351  qr(R_);
352  }
353 }
354 
355 
356 template<class MatrixType>
358 (
359  List<cmptType>& x,
360  const UList<cmptType>& source
361 ) const
362 {
363  solveImpl(x, source);
364 }
365 
366 
367 template<class MatrixType>
368 template<class Addr>
370 (
371  List<cmptType>& x,
373 ) const
374 {
375  solveImpl(x, source);
376 }
377 
378 
379 template<class MatrixType>
382 (
383  const UList<cmptType>& source
384 ) const
385 {
386  auto tresult(Q_.Tmul(source));
387 
388  solvex(tresult.ref());
389 
390  return tresult;
391 }
392 
393 
394 template<class MatrixType>
395 template<class Addr>
398 (
400 ) const
401 {
402  auto tresult(Q_.Tmul(source));
403 
404  solvex(tresult.ref());
405 
406  return tresult;
407 }
408 
409 
410 template<class MatrixType>
413 {
414  const label m = Q_.m();
415 
416  Field<cmptType> x(m);
417  SMatrix inv(m);
418 
419  for (label i = 0; i < m; ++i)
420  {
421  for (label j = 0; j < m; ++j)
422  {
423  x[j] = Q_(i, j);
424  }
425  solvex(x);
426  inv.subColumn(i) = x;
427  }
428 
429  return inv;
430 }
431 
432 
433 template<class MatrixType>
436 (
437  const SMatrix& A,
438  const RMatrix& rhs
439 )
440 {
441  const label mRows = A.m();
442  const label nCols = A.n();
443  const label pCols = rhs.n();
444 
445  #ifdef FULLDEBUG
446  {
447  const label qRows = rhs.m();
448  if (mRows != qRows)
449  {
451  << "Linear system is not solvable since the number of rows of"
452  << "A and rhs are not equal:" << tab << mRows << "vs." << qRows
453  << abort(FatalError);
454  }
455 
456  const List<cmptType> diag(A.diag());
457  for (const auto& val : diag)
458  {
459  if (mag(val) < SMALL)
460  {
462  << "SMALL-valued diagonal elem in back-substitution." << nl;
463  }
464  }
465  }
466  #endif
467 
468  RMatrix X({nCols, pCols}, Zero);
469 
470  for (label i = 0; i < pCols; ++i)
471  {
472  for (label j = mRows - 1; -1 < j; --j)
473  {
474  cmptType alpha(rhs(j, i));
475 
476  for (label k = j + 1; k < mRows; ++k)
477  {
478  alpha -= X(k, i)*A(j, k);
479  }
480 
481  X(j, i) = alpha/A(j, j);
482  }
483  }
484 
485  return X;
486 }
487 
488 
489 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
490 
491 template<class MatrixType>
492 MatrixType Foam::pinv
493 (
494  const MatrixType& A,
495  const scalar tolerance
496 )
497 {
498  scalar tol = tolerance;
499  typedef typename MatrixType::cmptType cmptType;
500 
501  if (A.empty())
502  {
504  << "Empty matrix found."
505  << abort(FatalError);
506  }
507 
508  if (A.size() == 1)
509  {
510  if (A(0,0) == cmptType(0))
511  {
512  return MatrixType({1, 1}, cmptType(0));
513  }
514 
515  return MatrixType({1, 1}, cmptType(1)/(A(0,0) + cmptType(VSMALL)));
516  }
517 
519  (
520  A,
523  );
524  const MatrixType& R(QRM.R());
525  const MatrixType& Q(QRM.Q());
526 
527  // R1 (KP:p. 648)
528  // Find the first diagonal element index with (almost) zero value in R
529  label firstZeroElemi = 0;
530  label i = 0;
531  while (i < 2)
532  {
533  const List<cmptType> diag(R.diag());
534 
535  auto lessThan = [=](const cmptType& x) { return tol > mag(x); };
536 
537  firstZeroElemi =
539  (
540  diag.cbegin(),
541  std::find_if(diag.cbegin(), diag.cend(), lessThan)
542  );
543 
544  if (firstZeroElemi == 0)
545  {
546  if (i == 0)
547  {
549  << "The largest diagonal element magnitude is nearly zero. "
550  << "Tightening the tolerance."
551  << endl;
552 
553  tol = 1e-13;
554  ++i;
555  }
556  else
557  {
559  << "The largest diagonal element magnitude is nearly zero. "
560  << "Returning a zero matrix."
561  << endl;
562  ++i;
563 
564  return MatrixType({A.n(), A.m()}, Zero);
565  }
566  }
567  else
568  {
569  i += 2;
570  }
571  }
572 
573  // Exclude the first (almost) zero diagonal row and the rows below
574  // since R diagonal is already descending due to the QR column pivoting
575  const RectangularMatrix<cmptType> R1(R.subMatrix(0, 0, firstZeroElemi));
576 
577  // R2 (KP:p. 648)
578  if (R1.n() < R1.m())
579  {
580  const SquareMatrix<cmptType> C(R1&R1);
581 
583  (
584  C,
585  QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
586  );
587 
589  (
590  QRSolve.backSubstitution
591  (
592  QRSolve.R(),
593  QRSolve.Q() & (R1.T())
594  )
595  );
596 
597  // R3 (KP:p. 648)
598  R2.resize(R.m(), R.n());
599 
600  return MatrixType((QRM.P()^R2)^Q);
601  }
602  else
603  {
604  const SquareMatrix<cmptType> C(R1^R1);
605 
607  (
608  C,
609  QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
610  );
611 
613  (
614  QRSolve.backSubstitution
615  (
616  QRSolve.R(),
617  QRSolve.Q() & R1
618  )
619  );
620 
621  // R3
622  R2.resize(R.m(), R.n());
623 
624  return MatrixType((QRM.P()^R2)^Q);
625  }
626 }
627 
628 
629 // ************************************************************************* //
Foam::QRMatrix::solve
void solve(List< cmptType > &x, const UList< cmptType > &source) const
Definition: QRMatrix.C:358
Foam::Swap
void Swap(DynamicList< T, SizeMinA > &a, DynamicList< T, SizeMinB > &b)
Definition: DynamicList.H:429
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Matrix::n
label n() const noexcept
The number of columns.
Definition: MatrixI.H:103
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::diag
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
Definition: pointPatchFieldFunctions.H:287
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Definition: readThermalProperties.H:212
Foam::pinv
MatrixType pinv(const MatrixType &A, const scalar tol=1e-5)
Definition: QRMatrix.C:493
QRMatrix.H
Foam::QRMatrix::cmptType
MatrixType::cmptType cmptType
Definition: QRMatrix.H:157
A
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::QRMatrix::inv
SMatrix inv() const
Return the inverse of (Q*R), so that solving x = (Q*R).inv()*source.
Definition: QRMatrix.C:412
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::QRMatrix::Q
const MatrixType & Q() const
Return the unitary similarity matrix.
Definition: QRMatrixI.H:155
C
volScalarField & C
Definition: readThermalProperties.H:102
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::QRMatrix
QRMatrix (i.e. QR decomposition, QR factorisation or orthogonal-triangular decomposition) decomposes ...
Definition: QRMatrix.H:152
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::Matrix::resize
void resize(const label m, const label n)
Change Matrix dimensions, preserving the elements.
Definition: Matrix.C:324
R
#define R(A, B, C, D, E, F, K, M)
Foam::QRMatrix::backSubstitution
RMatrix backSubstitution(const SMatrix &A, const RMatrix &b)
Solve a row-echelon-form linear system starting from the bottom.
Definition: QRMatrix.C:436
Foam::Field
Generic templated field type.
Definition: Field.H:63
Foam::inv
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:73
Foam::RectangularMatrix
A templated (M x N) rectangular matrix of objects of <Type>, containing M*N elements,...
Definition: RectangularMatrix.H:57
Foam::QRMatrix::R
const MatrixType & R() const
Return the upper triangular matrix.
Definition: QRMatrixI.H:164
Foam::FatalError
error FatalError
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::distance
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
Foam::QRMatrix::colPivoting
colPivoting
Options for the computation of column pivoting.
Definition: QRMatrix.H:177
Foam::SquareMatrix< cmptType >
Foam::QRMatrix::storeMethods
storeMethods
Options where to store R.
Definition: QRMatrix.H:170
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::tab
constexpr char tab
Definition: Ostream.H:403
Foam::Matrix::m
label m() const noexcept
The number of rows.
Definition: MatrixI.H:96
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::QRMatrix::QRMatrix
QRMatrix()
Construct null.
Definition: QRMatrix.C:186
Foam::List< cmptType >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::UList
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:103
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::identity
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
Definition: labelList.C:38
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::IndirectListBase
Base for lists with indirect addressing, templated on the list contents type and the addressing type....
Definition: IndirectListBase.H:56
Foam::QRMatrix::outputTypes
outputTypes
Options for the output matrix forms of QRMatrix.
Definition: QRMatrix.H:162
Foam::C
Graphite solid properties.
Definition: C.H:50
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::I
static const Identity< scalar > I
Definition: Identity.H:95
Foam::QRMatrix::decompose
void decompose(MatrixType &A)
Compute QR decomposition according to constructor settings.
Definition: QRMatrix.C:254
Foam::lessThan
static bool lessThan(const scalar &val, const scalar &upper)
Definition: foamVtkSeriesWriter.C:81
Foam::QRMatrix::P
SMatrix P() const
Create and return the permutation matrix.
Definition: QRMatrixI.H:179