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);
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];
143 else if (
sum != Type(Zero))
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];
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];
189 else if (
sum != Type(Zero))
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];
237template<
class Form,
class Type>
248 <<
"A and B must have identical inner dimensions but A.n = "
249 <<
A.n() <<
" and B.m = " <<
B.m()
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);
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
Various functions to operate on Lists.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
A templated (m x n) matrix of objects of <T>. The layout is (mRows x nCols) - row-major order:
label m() const noexcept
The number of rows.
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
const volScalarField & psi
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
void Swap(DynamicList< T, SizeMinA > &a, DynamicList< T, SizeMinB > &b)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
void multiply(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
errorManipArg< error, int > exit(error &err, const int errNo=1)
SolverPerformance< Type > solve(faMatrix< Type > &, const dictionary &)
Solve returning the solution statistics given convergence tolerance.