scalarMatrices.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) 2011-2020 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 "scalarMatrices.H"
27 #include "SVD.H"
28 
29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 
34 defineCompoundTypeName(RectangularMatrix<scalar>, scalarRectangularMatrix);
36 (
37  RectangularMatrix<scalar>,
39 );
40 
41 defineCompoundTypeName(SquareMatrix<scalar>, scalarSquareMatrix);
43 
45 (
46  SymmetricSquareMatrix<scalar>,
48 );
50 (
51  SymmetricSquareMatrix<scalar>,
53 );
54 
55 defineCompoundTypeName(DiagonalMatrix<scalar>, scalarDiagonalMatrix);
57 (
58  DiagonalMatrix<scalar>,
60 );
61 
62 }
63 
64 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
65 
67 (
68  scalarSquareMatrix& matrix,
69  labelList& pivotIndices
70 )
71 {
72  label sign;
73  LUDecompose(matrix, pivotIndices, sign);
74 }
75 
76 
78 (
79  scalarSquareMatrix& matrix,
80  labelList& pivotIndices,
81  label& sign
82 )
83 {
84  label m = matrix.m();
85  scalar vv[m];
86  sign = 1;
87 
88  for (label i=0; i<m; i++)
89  {
90  scalar largestCoeff = 0.0;
91  scalar temp;
92  const scalar* __restrict__ matrixi = matrix[i];
93 
94  for (label j=0; j<m; j++)
95  {
96  if ((temp = mag(matrixi[j])) > largestCoeff)
97  {
98  largestCoeff = temp;
99  }
100  }
101 
102  if (largestCoeff == 0.0)
103  {
105  << "Singular matrix" << exit(FatalError);
106  }
107 
108  vv[i] = 1.0/largestCoeff;
109  }
110 
111  for (label j=0; j<m; j++)
112  {
113  scalar* __restrict__ matrixj = matrix[j];
114 
115  for (label i=0; i<j; i++)
116  {
117  scalar* __restrict__ matrixi = matrix[i];
118 
119  scalar sum = matrixi[j];
120  for (label k=0; k<i; k++)
121  {
122  sum -= matrixi[k]*matrix(k, j);
123  }
124  matrixi[j] = sum;
125  }
126 
127  label iMax = 0;
128 
129  scalar largestCoeff = 0.0;
130  for (label i=j; i<m; i++)
131  {
132  scalar* __restrict__ matrixi = matrix[i];
133  scalar sum = matrixi[j];
134 
135  for (label k=0; k<j; k++)
136  {
137  sum -= matrixi[k]*matrix(k, j);
138  }
139 
140  matrixi[j] = sum;
141 
142  scalar temp;
143  if ((temp = vv[i]*mag(sum)) >= largestCoeff)
144  {
145  largestCoeff = temp;
146  iMax = i;
147  }
148  }
149 
150  pivotIndices[j] = iMax;
151 
152  if (j != iMax)
153  {
154  scalar* __restrict__ matrixiMax = matrix[iMax];
155 
156  for (label k=0; k<m; k++)
157  {
158  Swap(matrixj[k], matrixiMax[k]);
159  }
160 
161  sign *= -1;
162  vv[iMax] = vv[j];
163  }
164 
165  if (matrixj[j] == 0.0)
166  {
167  matrixj[j] = small;
168  }
169 
170  if (j != m-1)
171  {
172  scalar rDiag = 1.0/matrixj[j];
173 
174  for (label i=j+1; i<m; i++)
175  {
176  matrix(i, j) *= rDiag;
177  }
178  }
179  }
180 }
181 
182 
184 {
185  // Store result in upper triangular part of matrix
186  label size = matrix.m();
187 
188  // Set upper triangular parts to zero.
189  for (label j=0; j<size; j++)
190  {
191  for (label k=j + 1; k<size; k++)
192  {
193  matrix(j, k) = 0.0;
194  }
195  }
196 
197  for (label j=0; j<size; j++)
198  {
199  scalar d = 0.0;
200 
201  for (label k=0; k<j; k++)
202  {
203  scalar s = 0.0;
204 
205  for (label i=0; i<k; i++)
206  {
207  s += matrix(i, k)*matrix(i, j);
208  }
209 
210  s = (matrix(j, k) - s)/matrix(k, k);
211 
212  matrix(k, j) = s;
213  matrix(j, k) = s;
214 
215  d += sqr(s);
216  }
217 
218  d = matrix(j, j) - d;
219 
220  if (d < 0.0)
221  {
223  << "Matrix is not symmetric positive-definite. Unable to "
224  << "decompose."
225  << abort(FatalError);
226  }
227 
228  matrix(j, j) = sqrt(d);
229  }
230 }
231 
232 
233 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
234 
235 void Foam::multiply
236 (
237  scalarRectangularMatrix& ans, // value changed in return
238  const scalarRectangularMatrix& A,
239  const scalarRectangularMatrix& B,
240  const scalarRectangularMatrix& C
241 )
242 {
243  if (A.n() != B.m())
244  {
246  << "A and B must have identical inner dimensions but A.n = "
247  << A.n() << " and B.m = " << B.m()
248  << abort(FatalError);
249  }
250 
251  if (B.n() != C.m())
252  {
254  << "B and C must have identical inner dimensions but B.n = "
255  << B.n() << " and C.m = " << C.m()
256  << abort(FatalError);
257  }
258 
259  ans = scalarRectangularMatrix(A.m(), C.n(), scalar(0));
260 
261  for (label i=0; i<A.m(); i++)
262  {
263  for (label g = 0; g < C.n(); g++)
264  {
265  for (label l=0; l<C.m(); l++)
266  {
267  scalar ab = 0;
268  for (label j=0; j<A.n(); j++)
269  {
270  ab += A(i, j)*B(j, l);
271  }
272  ans(i, g) += C(l, g) * ab;
273  }
274  }
275  }
276 }
277 
278 
279 void Foam::multiply
280 (
281  scalarRectangularMatrix& ans, // value changed in return
282  const scalarRectangularMatrix& A,
283  const DiagonalMatrix<scalar>& B,
284  const scalarRectangularMatrix& C
285 )
286 {
287  if (A.n() != B.size())
288  {
290  << "A and B must have identical inner dimensions but A.n = "
291  << A.n() << " and B.m = " << B.size()
292  << abort(FatalError);
293  }
294 
295  if (B.size() != C.m())
296  {
298  << "B and C must have identical inner dimensions but B.n = "
299  << B.size() << " and C.m = " << C.m()
300  << abort(FatalError);
301  }
302 
303  ans = scalarRectangularMatrix(A.m(), C.n(), scalar(0));
304 
305  for (label i=0; i<A.m(); i++)
306  {
307  for (label g=0; g<C.n(); g++)
308  {
309  for (label l=0; l<C.m(); l++)
310  {
311  ans(i, g) += C(l, g) * A(i, l)*B[l];
312  }
313  }
314  }
315 }
316 
317 
318 void Foam::multiply
319 (
320  scalarSquareMatrix& ans, // value changed in return
321  const scalarSquareMatrix& A,
322  const DiagonalMatrix<scalar>& B,
323  const scalarSquareMatrix& C
324 )
325 {
326  if (A.m() != B.size())
327  {
329  << "A and B must have identical dimensions but A.m = "
330  << A.m() << " and B.m = " << B.size()
331  << abort(FatalError);
332  }
333 
334  if (B.size() != C.m())
335  {
337  << "B and C must have identical dimensions but B.m = "
338  << B.size() << " and C.m = " << C.m()
339  << abort(FatalError);
340  }
341 
342  const label size = A.m();
343 
344  ans = scalarSquareMatrix(size, Zero);
345 
346  for (label i=0; i<size; i++)
347  {
348  for (label g=0; g<size; g++)
349  {
350  for (label l=0; l<size; l++)
351  {
352  ans(i, g) += C(l, g)*A(i, l)*B[l];
353  }
354  }
355  }
356 }
357 
358 
360 (
361  const scalarRectangularMatrix& A,
362  scalar minCondition
363 )
364 {
365  SVD svd(A, minCondition);
366  return svd.VSinvUt();
367 }
368 
369 
370 // ************************************************************************* //
dimensionedScalar sign(const dimensionedScalar &ds)
scalarRectangularMatrix SVDinv(const scalarRectangularMatrix &A, scalar minCondition=0)
Return the inverse of matrix A using SVD.
label n() const
Return the number of columns.
Definition: MatrixI.H:64
DiagonalMatrix< scalar > scalarDiagonalMatrix
void multiply(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
addCompoundToRunTimeSelectionTable(List< complex >, complexList)
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
defineCompoundTypeName(List< complex >, complexList)
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
dimensionedScalar sqrt(const dimensionedScalar &ds)
label k
Boltzmann constant.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
void Swap(T &a, T &b)
Definition: Swap.H:43
Singular value decomposition of a rectangular matrix.
Definition: SVD.H:51
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
SymmetricSquareMatrix< scalar > scalarSymmetricSquareMatrix
scalarRectangularMatrix VSinvUt() const
Return the matrix product V S^(-1) U^T (the pseudo inverse)
Definition: SVD.C:385
A templated 2D square symmetric matrix of objects of <T>, where the n x n matrix dimension is known a...
label m() const
Return the number of rows.
Definition: MatrixI.H:57
dimensioned< scalar > mag(const dimensioned< Type > &)
const dimensionedVector & g
SquareMatrix< scalar > scalarSquareMatrix
Namespace for OpenFOAM.
RectangularMatrix< scalar > scalarRectangularMatrix