53 for (
label i=0; i<m; i++)
55 scalar largestCoeff = 0.0;
57 const scalar* __restrict__ matrixi = matrix[i];
59 for (
label j=0; j<m; j++)
61 if ((temp =
mag(matrixi[j])) > largestCoeff)
67 if (largestCoeff == 0.0)
73 vv[i] = 1.0/largestCoeff;
76 for (
label j=0; j<m; j++)
78 scalar* __restrict__ matrixj = matrix[j];
80 for (
label i=0; i<j; i++)
82 scalar* __restrict__ matrixi = matrix[i];
84 scalar
sum = matrixi[j];
87 sum -= matrixi[
k]*matrix(
k, j);
94 scalar largestCoeff = 0.0;
95 for (
label i=j; i<m; i++)
97 scalar* __restrict__ matrixi = matrix[i];
98 scalar
sum = matrixi[j];
102 sum -= matrixi[
k]*matrix(
k, j);
108 if ((temp = vv[i]*
mag(sum)) >= largestCoeff)
115 pivotIndices[j] = iMax;
119 scalar* __restrict__ matrixiMax = matrix[iMax];
123 Swap(matrixj[
k], matrixiMax[k]);
130 if (matrixj[j] == 0.0)
137 scalar rDiag = 1.0/matrixj[j];
139 for (
label i=j+1; i<m; i++)
141 matrix(i, j) *= rDiag;
154 for (
label j=0; j<size; j++)
162 for (
label j=0; j<size; j++)
172 s += matrix(i, k)*matrix(i, j);
175 s = (matrix(j, k) -
s)/matrix(k, k);
183 d = matrix(j, j) - d;
188 <<
"Matrix is not symmetric positive-definite. Unable to " 193 matrix(j, j) =
sqrt(d);
211 <<
"A and B must have identical inner dimensions but A.n = " 212 << A.
n() <<
" and B.m = " << B.
m()
219 <<
"B and C must have identical inner dimensions but B.n = " 220 << B.
n() <<
" and C.m = " << C.
m()
226 for (
label i=0; i<A.
m(); i++)
230 for (
label l=0; l<C.
m(); l++)
233 for (
label j=0; j<A.
n(); j++)
235 ab += A(i, j)*B(j, l);
237 ans(i,
g) += C(l,
g) * ab;
252 if (A.
n() != B.
size())
255 <<
"A and B must have identical inner dimensions but A.n = " 256 << A.
n() <<
" and B.m = " << B.
size()
260 if (B.
size() != C.
m())
263 <<
"B and C must have identical inner dimensions but B.n = " 264 << B.
size() <<
" and C.m = " << C.
m()
270 for (
label i=0; i<A.
m(); i++)
274 for (
label l=0; l<C.
m(); l++)
276 ans(i,
g) += C(l,
g) * A(i, l)*B[l];
291 if (A.
m() != B.
size())
294 <<
"A and B must have identical dimensions but A.m = " 295 << A.
m() <<
" and B.m = " << B.
size()
299 if (B.
size() != C.
m())
302 <<
"B and C must have identical dimensions but B.m = " 303 << B.
size() <<
" and C.m = " << C.
m()
311 for (
label i=0; i<size; i++)
315 for (
label l=0; l<size; l++)
317 ans(i,
g) += C(l,
g)*A(i, l)*B[l];
330 SVD svd(A, minCondition);
dimensionedScalar sign(const dimensionedScalar &ds)
scalarRectangularMatrix SVDinv(const scalarRectangularMatrix &A, scalar minCondition=0)
Return the inverse of matrix A using SVD.
label n() const
Return the number of columns.
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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
dimensionedScalar sqrt(const dimensionedScalar &ds)
label k
Boltzmann constant.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Singular value decomposition of a rectangular matrix.
errorManip< error > abort(error &err)
scalarRectangularMatrix VSinvUt() const
Return the matrix product V S^(-1) U^T (the pseudo inverse)
A templated 2D square symmetric matrix of objects of <T>, where the n x n matrix dimension is known a...
label m() const
Return the number of rows.
dimensioned< scalar > mag(const dimensioned< Type > &)
const dimensionedVector & g
SquareMatrix< scalar > scalarSquareMatrix
RectangularMatrix< scalar > scalarRectangularMatrix