scalarMatricesTemplates.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 | Copyright (C) 2011-2016 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 "Swap.H"
28 #include "ListOps.H"
29 
30 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
31 
32 template<class Type>
33 void Foam::solve
34 (
35  scalarSquareMatrix& tmpMatrix,
36  List<Type>& sourceSol
37 )
38 {
39  label m = tmpMatrix.m();
40 
41  // Elimination
42  for (label i=0; i<m; i++)
43  {
44  label iMax = i;
45  scalar largestCoeff = mag(tmpMatrix(iMax, i));
46 
47  // Swap elements around to find a good pivot
48  for (label j=i+1; j<m; j++)
49  {
50  if (mag(tmpMatrix(j, i)) > largestCoeff)
51  {
52  iMax = j;
53  largestCoeff = mag(tmpMatrix(iMax, i));
54  }
55  }
56 
57  if (i != iMax)
58  {
59  for (label k=i; k<m; k++)
60  {
61  Swap(tmpMatrix(i, k), tmpMatrix(iMax, k));
62  }
63  Swap(sourceSol[i], sourceSol[iMax]);
64  }
65 
66  // Check that the system of equations isn't singular
67  if (mag(tmpMatrix(i, i)) < 1e-20)
68  {
70  << "Singular Matrix"
71  << exit(FatalError);
72  }
73 
74  // Reduce to upper triangular form
75  for (label j=i+1; j<m; j++)
76  {
77  sourceSol[j] -= sourceSol[i]*(tmpMatrix(j, i)/tmpMatrix(i, i));
78 
79  for (label k=m-1; k>=i; k--)
80  {
81  tmpMatrix(j, k) -=
82  tmpMatrix(i, k)*tmpMatrix(j, i)/tmpMatrix(i, i);
83  }
84  }
85  }
86 
87  // Back-substitution
88  for (label j=m-1; j>=0; j--)
89  {
90  Type ntempvec = Zero;
91 
92  for (label k=j+1; k<m; k++)
93  {
94  ntempvec += tmpMatrix(j, k)*sourceSol[k];
95  }
96 
97  sourceSol[j] = (sourceSol[j] - ntempvec)/tmpMatrix(j, j);
98  }
99 }
100 
101 
102 template<class Type>
103 void Foam::solve
104 (
105  List<Type>& psi,
106  const scalarSquareMatrix& matrix,
107  const List<Type>& source
108 )
109 {
110  scalarSquareMatrix tmpMatrix = matrix;
111  psi = source;
112  solve(tmpMatrix, psi);
113 }
114 
115 
116 template<class Type>
118 (
119  const scalarSquareMatrix& luMatrix,
120  const labelList& pivotIndices,
121  List<Type>& sourceSol
122 )
123 {
124  label m = luMatrix.m();
125 
126  label ii = 0;
127 
128  for (label i=0; i<m; i++)
129  {
130  label ip = pivotIndices[i];
131  Type sum = sourceSol[ip];
132  sourceSol[ip] = sourceSol[i];
133  const scalar* __restrict__ luMatrixi = luMatrix[i];
134 
135  if (ii != 0)
136  {
137  for (label j=ii-1; j<i; j++)
138  {
139  sum -= luMatrixi[j]*sourceSol[j];
140  }
141  }
142  else if (sum != pTraits<Type>::zero)
143  {
144  ii = i+1;
145  }
146 
147  sourceSol[i] = sum;
148  }
149 
150  for (label i=m-1; i>=0; i--)
151  {
152  Type sum = sourceSol[i];
153  const scalar* __restrict__ luMatrixi = luMatrix[i];
154 
155  for (label j=i+1; j<m; j++)
156  {
157  sum -= luMatrixi[j]*sourceSol[j];
158  }
159 
160  sourceSol[i] = sum/luMatrixi[i];
161  }
162 }
163 
164 
165 template<class Type>
167 (
168  const scalarSymmetricSquareMatrix& luMatrix,
169  List<Type>& sourceSol
170 )
171 {
172  label m = luMatrix.m();
173 
174  label ii = 0;
175 
176  for (label i=0; i<m; i++)
177  {
178  Type sum = sourceSol[i];
179  const scalar* __restrict__ luMatrixi = luMatrix[i];
180 
181  if (ii != 0)
182  {
183  for (label j=ii-1; j<i; j++)
184  {
185  sum -= luMatrixi[j]*sourceSol[j];
186  }
187  }
188  else if (sum != pTraits<Type>::zero)
189  {
190  ii = i+1;
191  }
192 
193  sourceSol[i] = sum/luMatrixi[i];
194  }
195 
196  for (label i=m-1; i>=0; i--)
197  {
198  Type sum = sourceSol[i];
199  const scalar* __restrict__ luMatrixi = luMatrix[i];
200 
201  for (label j=i+1; j<m; j++)
202  {
203  sum -= luMatrixi[j]*sourceSol[j];
204  }
205 
206  sourceSol[i] = sum/luMatrixi[i];
207  }
208 }
209 
210 
211 template<class Type>
212 void Foam::LUsolve
213 (
214  scalarSquareMatrix& matrix,
215  List<Type>& sourceSol
216 )
217 {
218  labelList pivotIndices(matrix.m());
219  LUDecompose(matrix, pivotIndices);
220  LUBacksubstitute(matrix, pivotIndices, sourceSol);
221 }
222 
223 
224 template<class Type>
225 void Foam::LUsolve
226 (
228  List<Type>& sourceSol
229 )
230 {
231  LUDecompose(matrix);
232  LUBacksubstitute(matrix, sourceSol);
233 }
234 
235 
236 template<class Form, class Type>
237 void Foam::multiply
238 (
239  Matrix<Form, Type>& ans, // value changed in return
240  const Matrix<Form, Type>& A,
241  const Matrix<Form, Type>& B
242 )
243 {
244  if (A.n() != B.m())
245  {
247  << "A and B must have identical inner dimensions but A.n = "
248  << A.n() << " and B.m = " << B.m()
249  << abort(FatalError);
250  }
251 
252  ans = Matrix<Form, Type>(A.m(), B.n(), scalar(0));
253 
254  for (label i=0; i<A.m(); i++)
255  {
256  for (label j=0; j<B.n(); j++)
257  {
258  for (label l=0; l<B.m(); l++)
259  {
260  ans(i, j) += A(i, l)*B(l, j);
261  }
262  }
263  }
264 }
265 
266 
267 // ************************************************************************* //
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.
Traits class for primitives.
Definition: pTraits.H:50
label k
Boltzmann constant.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Various functions to operate on Lists.
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
void Swap(T &a, T &b)
Definition: Swap.H:43
SolverPerformance< Type > solve(fvMatrix< Type > &, const dictionary &)
Solve returning the solution statistics given convergence tolerance.
rhoEqn solve()
static const zero Zero
Definition: zero.H:91
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A templated (m x n) matrix of objects of <T>.
A templated 2D square symmetric matrix of objects of <T>, where the n x n matrix dimension is known a...
Swap its arguments.
label m() const
Return the number of rows.
Definition: MatrixI.H:57
dimensioned< scalar > mag(const dimensioned< Type > &)
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Solve the matrix using LU decomposition with pivoting.
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution.