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