QRMatrix.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "QRMatrix.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class MatrixType>
32 {
33  decompose(M);
34 }
35 
36 
37 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
38 
39 template<class MatrixType>
41 {
42  const label m = M.m();
43  const label n = M.n();
44 
45  // Initialise the R-matrix to M
46  R_ = M;
47 
48  // Initialise the Q-matrix to I
49  Q_.setSize(m);
50  Q_ = I;
51 
52  // Pre-allocate temporary storage for the Householder steps
53  QMatrixType Qk(m);
54  QMatrixType Rk(m);
55  Field<cmptType> uk(m);
56 
57  for (label k=0; k<n; k++)
58  {
59  // alpha = -|column k of Rk|
61  for (label bi=k; bi<m; bi++)
62  {
63  alpha += sqr(R_(bi, k));
64  }
65  alpha = sqrt(alpha);
66 
67  if (R_(k, k) > 0)
68  {
69  alpha = -alpha;
70  }
71 
72  // uk = column k of Rk - alpha*ek
73  // rSumSqrUk = 2/sum(sqr(uk))
74  uk[k] = R_(k, k) - alpha;
75  cmptType rSumSqrUk = sqr(uk[k]);
76  for (label bi=k+1; bi<m; bi++)
77  {
78  uk[bi] = R_(bi, k);
79  rSumSqrUk += sqr(uk[bi]);
80  }
81  rSumSqrUk = 2/rSumSqrUk;
82 
83  // Qk = I - 2*u*uT/sum(sqr(uk))
84  for (label bi=k; bi<m; bi++)
85  {
86  for (label bj=k; bj<m; bj++)
87  {
88  Qk(bi, bj) = -rSumSqrUk*uk[bi]*uk[bj];
89  }
90  Qk(bi, bi) += 1;
91  }
92 
93  // Rk = Qk*R
94  for (label bi=k; bi<m; bi++)
95  {
96  for (label bk=k; bk<n; bk++)
97  {
98  Rk(bi, bk) = Zero;
99  for (label bj=k; bj<n; bj++)
100  {
101  Rk(bi, bk) += Qk(bi, bj)*R_(bj, bk);
102  }
103  }
104  }
105 
106  // Ensure diagonal is positive
107  if (R_(k, k) < 0)
108  {
109  // R = -Rk
110  // Q = -Q
111  for (label bi=k; bi<m; bi++)
112  {
113  for (label bj=k; bj<n; bj++)
114  {
115  R_(bi, bj) = -Rk(bi, bj);
116  }
117  for (label bj=k; bj<m; bj++)
118  {
119  Q_(bi, bj) = -Q_(bi, bj);
120  }
121  }
122  }
123  else
124  {
125  // R = Rk
126  for (label bi=k; bi<m; bi++)
127  {
128  for (label bj=k; bj<n; bj++)
129  {
130  R_(bi, bj) = Rk(bi, bj);
131  }
132  }
133  }
134 
135  // Q = Q*Qk (using Rk as temporary storage)
136  for (label bi=0; bi<m; bi++)
137  {
138  for (label bk=k; bk<m; bk++)
139  {
140  Rk(bi, bk) = Zero;
141  for (label bj=k; bj<m; bj++)
142  {
143  Rk(bi, bk) += Q_(bi, bj)*Qk(bj, bk);
144  }
145  }
146  }
147  for (label bi=0; bi<m; bi++)
148  {
149  for (label bj=k; bj<n; bj++)
150  {
151  Q_(bi, bj) = Rk(bi, bj);
152  }
153  }
154  }
155 }
156 
157 
158 template<class MatrixType>
160 (
162 ) const
163 {
164  const label n = R_.n();
165 
166  for (int i=n-1; i>=0; i--)
167  {
168  cmptType sum = x[i];
169 
170  for (label j=i+1; j<n; j++)
171  {
172  sum -= x[j]*R_(i, j);
173  }
174 
175  if (mag(R_(i, i)) < small)
176  {
178  << "Back-substitution failed due to small diagonal"
179  << abort(FatalError);
180  }
181 
182  x[i] = sum/R_(i, i);
183  }
184 }
185 
186 
187 template<class MatrixType>
189 (
191  const Field<cmptType>& source
192 ) const
193 {
194  const label m = Q_.m();
195 
196  // x = Q_.T()*source;
197  for (label i=0; i<m; i++)
198  {
199  x[i] = 0;
200  for (label j=0; j<m; j++)
201  {
202  x[i] += Q_(j, i)*source[j];
203  }
204  }
205 
206  solvex(x);
207 }
208 
209 
210 template<class MatrixType>
213 (
214  const Field<cmptType>& source
215 ) const
216 {
217  tmp<Field<cmptType>> tx(new Field<cmptType>(Q_.m()));
218  Field<cmptType>& x = tx.ref();
219 
220  solve(x, source);
221 
222  return tx;
223 }
224 
225 
226 template<class MatrixType>
229 {
230  const label m = Q_.m();
231 
232  Field<cmptType> x(m);
233  QMatrixType inv(m);
234 
235  for (label i=0; i<m; i++)
236  {
237  for (label j=0; j<m; j++)
238  {
239  x[j] = Q_(i, j);
240  }
241  solvex(x);
242  inv.block(m, 1, 0, i) = x;
243  }
244 
245  return inv;
246 }
247 
248 
249 // ************************************************************************* //
label k
#define M(I)
label n
Pre-declare SubField and related Field type.
Definition: Field.H:83
label m() const
Return the number of rows.
Definition: MatrixI.H:57
Class templated on matrix type to perform the QR decomposition using Householder reflections on a squ...
Definition: QRMatrix.H:52
void decompose(const MatrixType &M)
Decompose given matrix.
Definition: QRMatrix.C:40
QRMatrix()
Construct null.
Definition: QRMatrixI.H:29
void solve(Field< cmptType > &x, const Field< cmptType > &source) const
Solve the linear system with the given source.
Definition: QRMatrix.C:189
MatrixType::cmptType cmptType
Definition: QRMatrix.H:56
QMatrixType inv() const
Return the inverse of a square matrix.
Definition: QRMatrix.C:228
A class for managing temporary objects.
Definition: tmp.H:55
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
static const zero Zero
Definition: zero.H:97
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static const Identity< scalar > I
Definition: Identity.H:93
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
error FatalError
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
rhoEqn solve()