45 scalar largestCoeff =
mag(tmpMatrix[iMax][i]);
48 for (
label j=i+1; j<
n; j++)
50 if (
mag(tmpMatrix[j][i]) > largestCoeff)
53 largestCoeff =
mag(tmpMatrix[iMax][i]);
63 Swap(tmpMatrix[i][
k], tmpMatrix[iMax][k]);
65 Swap(sourceSol[i], sourceSol[iMax]);
69 if (
mag(tmpMatrix[i][i]) < 1
e-20)
71 FatalErrorIn(
"solve(scalarSquareMatrix&, Field<Type>& sourceSol)")
77 for (
label j=i+1; j<
n; j++)
79 sourceSol[j] -= sourceSol[i]*(tmpMatrix[j][i]/tmpMatrix[i][i]);
84 tmpMatrix[i][
k]*tmpMatrix[j][i]/tmpMatrix[i][i];
90 for (
label j=n-1; j>=0; j--)
96 ntempvec += tmpMatrix[j][
k]*sourceSol[
k];
99 sourceSol[j] = (sourceSol[j] - ntempvec)/tmpMatrix[j][j];
114 solve(tmpMatrix, psi);
132 label ip = pivotIndices[i];
133 Type
sum = sourceSol[ip];
134 sourceSol[ip] = sourceSol[i];
135 const scalar* __restrict__ luMatrixi = luMatrix[i];
139 for (
label j=ii-1; j<i; j++)
141 sum -= luMatrixi[j]*sourceSol[j];
152 for (
label i=n-1; i>=0; i--)
154 Type
sum = sourceSol[i];
155 const scalar* __restrict__ luMatrixi = luMatrix[i];
157 for (
label j=i+1; j<
n; j++)
159 sum -= luMatrixi[j]*sourceSol[j];
162 sourceSol[i] = sum/luMatrixi[i];
180 Type
sum = sourceSol[i];
181 const scalar* __restrict__ luMatrixi = luMatrix[i];
185 for (
label j=ii-1; j<i; j++)
187 sum -= luMatrixi[j]*sourceSol[j];
195 sourceSol[i] = sum/luMatrixi[i];
198 for (
label i=n-1; i>=0; i--)
200 Type
sum = sourceSol[i];
201 const scalar* __restrict__ luMatrixi = luMatrix[i];
203 for (
label j=i+1; j<
n; j++)
205 sum -= luMatrixi[j]*sourceSol[j];
208 sourceSol[i] = sum/luMatrixi[i];
238 template<
class Form,
class Type>
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()
261 for (
label i = 0; i < A.
n(); i++)
263 for (
label j = 0; j < B.
m(); j++)
265 for (
label l = 0; l < B.
n(); l++)
267 ans[i][j] += A[i][l]*B[l][j];
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.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n() const
Return the number of rows.
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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Various functions to operate on Lists.
const dimensionedScalar e
Elementary charge.
label k
Boltzmann constant.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
errorManip< error > abort(error &err)
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.
solverPerformance solve(fvMatrix< Type > &, const dictionary &)
Solve returning the solution statistics given convergence tolerance.
Traits class for primitives.
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.