40 label m = tmpMatrix.
m();
43 for (label i = 0; i < m; ++i)
46 scalar largestCoeff =
mag(tmpMatrix(iMax, i));
49 for (label j = i + 1; j < m; ++j)
51 if (
mag(tmpMatrix(j, i)) > largestCoeff)
54 largestCoeff =
mag(tmpMatrix(iMax, i));
60 for (label
k = i;
k < m; ++
k)
62 Swap(tmpMatrix(i,
k), tmpMatrix(iMax,
k));
64 Swap(sourceSol[i], sourceSol[iMax]);
68 if (
mag(tmpMatrix(i, i)) < 1
e-20)
76 for (label j = i + 1; j < m; ++j)
78 sourceSol[j] -= sourceSol[i]*(tmpMatrix(j, i)/tmpMatrix(i, i));
80 for (label
k = m - 1;
k >= i; --
k)
83 tmpMatrix(i,
k)*tmpMatrix(j, i)/tmpMatrix(i, i);
89 for (label j = m - 1; j >= 0; --j)
93 for (label
k = j + 1;
k < m; ++
k)
95 ntempvec += tmpMatrix(j,
k)*sourceSol[
k];
98 sourceSol[j] = (sourceSol[j] - ntempvec)/tmpMatrix(j, j);
108 const List<Type>& source
122 List<Type>& sourceSol
125 label m = luMatrix.m();
129 for (label i = 0; i < m; ++i)
131 label ip = pivotIndices[i];
132 Type
sum = sourceSol[ip];
133 sourceSol[ip] = sourceSol[i];
134 const scalar* __restrict__ luMatrixi = luMatrix[i];
138 for (label j = ii - 1; j < i; ++j)
140 sum -= luMatrixi[j]*sourceSol[j];
151 for (label i = m - 1; i >= 0; --i)
153 Type
sum = sourceSol[i];
154 const scalar* __restrict__ luMatrixi = luMatrix[i];
156 for (label j = i + 1; j < m; ++j)
158 sum -= luMatrixi[j]*sourceSol[j];
161 sourceSol[i] =
sum/luMatrixi[i];
170 List<Type>& sourceSol
173 label m = luMatrix.m();
177 for (label i = 0; i < m; ++i)
179 Type
sum = sourceSol[i];
180 const scalar* __restrict__ luMatrixi = luMatrix[i];
184 for (label j = ii - 1; j < i; ++j)
186 sum -= luMatrixi[j]*sourceSol[j];
194 sourceSol[i] =
sum/luMatrixi[i];
197 for (label i = m - 1; i >= 0; --i)
199 Type
sum = sourceSol[i];
200 const scalar* __restrict__ luMatrixi = luMatrix[i];
202 for (label j = i + 1; j < m; ++j)
204 sum -= luMatrixi[j]*sourceSol[j];
207 sourceSol[i] =
sum/luMatrixi[i];
216 List<Type>& sourceSol
229 List<Type>& sourceSol
237 template<
class Form,
class Type>
240 Matrix<Form, Type>& ans,
241 const Matrix<Form, Type>&
A,
242 const Matrix<Form, Type>&
B 248 <<
"A and B must have identical inner dimensions but A.n = " 249 <<
A.n() <<
" and B.m = " <<
B.m()
253 ans = Matrix<Form, Type>(
A.m(),
B.n(),
Zero);
255 for (label i = 0; i <
A.m(); ++i)
257 for (label j = 0; j <
B.n(); ++j)
259 for (label l = 0; l <
B.m(); ++l)
261 ans(i, j) +=
A(i, l)*
B(l, j);
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
SolverPerformance< Type > solve(faMatrix< Type > &, const dictionary &)
Solve returning the solution statistics given convergence tolerance.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
label k
Boltzmann constant.
Various functions to operate on Lists.
const dimensionedScalar e
Elementary charge.
errorManip< error > abort(error &err)
label m() const noexcept
The number of rows.
SymmetricSquareMatrix< scalar > scalarSymmetricSquareMatrix
void Swap(DynamicList< T, SizeMinA > &a, DynamicList< T, SizeMinB > &b)
const volScalarField & psi
List< label > labelList
A List of labels.
void multiply(DimensionedField< Type, GeoMesh > &result, const DimensionedField< Type, GeoMesh > &f1, const DimensionedField< scalar, GeoMesh > &f2)
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
SquareMatrix< scalar > scalarSquareMatrix
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Solve the matrix using LU decomposition with pivoting returning the LU form of the matrix and the sol...
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution in the source.
static constexpr const zero Zero
Global zero (0)