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-2018 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 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
30 
32 (
33  scalarSquareMatrix& matrix,
34  labelList& pivotIndices
35 )
36 {
37  label sign;
38  LUDecompose(matrix, pivotIndices, sign);
39 }
40 
41 
43 (
44  scalarSquareMatrix& matrix,
45  labelList& pivotIndices,
46  label& sign
47 )
48 {
49  label m = matrix.m();
50  scalar vv[m];
51  sign = 1;
52 
53  for (label i=0; i<m; i++)
54  {
55  scalar largestCoeff = 0.0;
56  scalar temp;
57  const scalar* __restrict__ matrixi = matrix[i];
58 
59  for (label j=0; j<m; j++)
60  {
61  if ((temp = mag(matrixi[j])) > largestCoeff)
62  {
63  largestCoeff = temp;
64  }
65  }
66 
67  if (largestCoeff == 0.0)
68  {
70  << "Singular matrix" << exit(FatalError);
71  }
72 
73  vv[i] = 1.0/largestCoeff;
74  }
75 
76  for (label j=0; j<m; j++)
77  {
78  scalar* __restrict__ matrixj = matrix[j];
79 
80  for (label i=0; i<j; i++)
81  {
82  scalar* __restrict__ matrixi = matrix[i];
83 
84  scalar sum = matrixi[j];
85  for (label k=0; k<i; k++)
86  {
87  sum -= matrixi[k]*matrix(k, j);
88  }
89  matrixi[j] = sum;
90  }
91 
92  label iMax = 0;
93 
94  scalar largestCoeff = 0.0;
95  for (label i=j; i<m; i++)
96  {
97  scalar* __restrict__ matrixi = matrix[i];
98  scalar sum = matrixi[j];
99 
100  for (label k=0; k<j; k++)
101  {
102  sum -= matrixi[k]*matrix(k, j);
103  }
104 
105  matrixi[j] = sum;
106 
107  scalar temp;
108  if ((temp = vv[i]*mag(sum)) >= largestCoeff)
109  {
110  largestCoeff = temp;
111  iMax = i;
112  }
113  }
114 
115  pivotIndices[j] = iMax;
116 
117  if (j != iMax)
118  {
119  scalar* __restrict__ matrixiMax = matrix[iMax];
120 
121  for (label k=0; k<m; k++)
122  {
123  Swap(matrixj[k], matrixiMax[k]);
124  }
125 
126  sign *= -1;
127  vv[iMax] = vv[j];
128  }
129 
130  if (matrixj[j] == 0.0)
131  {
132  matrixj[j] = small;
133  }
134 
135  if (j != m-1)
136  {
137  scalar rDiag = 1.0/matrixj[j];
138 
139  for (label i=j+1; i<m; i++)
140  {
141  matrix(i, j) *= rDiag;
142  }
143  }
144  }
145 }
146 
147 
149 {
150  // Store result in upper triangular part of matrix
151  label size = matrix.m();
152 
153  // Set upper triangular parts to zero.
154  for (label j=0; j<size; j++)
155  {
156  for (label k=j + 1; k<size; k++)
157  {
158  matrix(j, k) = 0.0;
159  }
160  }
161 
162  for (label j=0; j<size; j++)
163  {
164  scalar d = 0.0;
165 
166  for (label k=0; k<j; k++)
167  {
168  scalar s = 0.0;
169 
170  for (label i=0; i<k; i++)
171  {
172  s += matrix(i, k)*matrix(i, j);
173  }
174 
175  s = (matrix(j, k) - s)/matrix(k, k);
176 
177  matrix(k, j) = s;
178  matrix(j, k) = s;
179 
180  d += sqr(s);
181  }
182 
183  d = matrix(j, j) - d;
184 
185  if (d < 0.0)
186  {
188  << "Matrix is not symmetric positive-definite. Unable to "
189  << "decompose."
190  << abort(FatalError);
191  }
192 
193  matrix(j, j) = sqrt(d);
194  }
195 }
196 
197 
198 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
199 
200 void Foam::multiply
201 (
202  scalarRectangularMatrix& ans, // value changed in return
203  const scalarRectangularMatrix& A,
204  const scalarRectangularMatrix& B,
205  const scalarRectangularMatrix& C
206 )
207 {
208  if (A.n() != B.m())
209  {
211  << "A and B must have identical inner dimensions but A.n = "
212  << A.n() << " and B.m = " << B.m()
213  << abort(FatalError);
214  }
215 
216  if (B.n() != C.m())
217  {
219  << "B and C must have identical inner dimensions but B.n = "
220  << B.n() << " and C.m = " << C.m()
221  << abort(FatalError);
222  }
223 
224  ans = scalarRectangularMatrix(A.m(), C.n(), scalar(0));
225 
226  for (label i=0; i<A.m(); i++)
227  {
228  for (label g = 0; g < C.n(); g++)
229  {
230  for (label l=0; l<C.m(); l++)
231  {
232  scalar ab = 0;
233  for (label j=0; j<A.n(); j++)
234  {
235  ab += A(i, j)*B(j, l);
236  }
237  ans(i, g) += C(l, g) * ab;
238  }
239  }
240  }
241 }
242 
243 
244 void Foam::multiply
245 (
246  scalarRectangularMatrix& ans, // value changed in return
247  const scalarRectangularMatrix& A,
248  const DiagonalMatrix<scalar>& B,
249  const scalarRectangularMatrix& C
250 )
251 {
252  if (A.n() != B.size())
253  {
255  << "A and B must have identical inner dimensions but A.n = "
256  << A.n() << " and B.m = " << B.size()
257  << abort(FatalError);
258  }
259 
260  if (B.size() != C.m())
261  {
263  << "B and C must have identical inner dimensions but B.n = "
264  << B.size() << " and C.m = " << C.m()
265  << abort(FatalError);
266  }
267 
268  ans = scalarRectangularMatrix(A.m(), C.n(), scalar(0));
269 
270  for (label i=0; i<A.m(); i++)
271  {
272  for (label g=0; g<C.n(); g++)
273  {
274  for (label l=0; l<C.m(); l++)
275  {
276  ans(i, g) += C(l, g) * A(i, l)*B[l];
277  }
278  }
279  }
280 }
281 
282 
283 void Foam::multiply
284 (
285  scalarSquareMatrix& ans, // value changed in return
286  const scalarSquareMatrix& A,
287  const DiagonalMatrix<scalar>& B,
288  const scalarSquareMatrix& C
289 )
290 {
291  if (A.m() != B.size())
292  {
294  << "A and B must have identical dimensions but A.m = "
295  << A.m() << " and B.m = " << B.size()
296  << abort(FatalError);
297  }
298 
299  if (B.size() != C.m())
300  {
302  << "B and C must have identical dimensions but B.m = "
303  << B.size() << " and C.m = " << C.m()
304  << abort(FatalError);
305  }
306 
307  const label size = A.m();
308 
309  ans = scalarSquareMatrix(size, Zero);
310 
311  for (label i=0; i<size; i++)
312  {
313  for (label g=0; g<size; g++)
314  {
315  for (label l=0; l<size; l++)
316  {
317  ans(i, g) += C(l, g)*A(i, l)*B[l];
318  }
319  }
320  }
321 }
322 
323 
325 (
326  const scalarRectangularMatrix& A,
327  scalar minCondition
328 )
329 {
330  SVD svd(A, minCondition);
331  return svd.VSinvUt();
332 }
333 
334 
335 // ************************************************************************* //
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
void multiply(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
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:163
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
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
RectangularMatrix< scalar > scalarRectangularMatrix