52 const label size = matrix.m();
54 pivotIndices.resize_nocopy(size);
56 List<scalar> vv(size);
59 for (label i = 0; i < size; ++i)
61 scalar largestCoeff = 0.0;
63 const scalar* __restrict__ matrixi = matrix[i];
65 for (label j = 0; j < size; ++j)
67 if ((temp =
mag(matrixi[j])) > largestCoeff)
73 if (largestCoeff == 0.0)
79 vv[i] = 1.0/largestCoeff;
82 for (label j = 0; j < size; ++j)
84 scalar* __restrict__ matrixj = matrix[j];
86 for (label i = 0; i < j; ++i)
88 scalar* __restrict__ matrixi = matrix[i];
90 scalar
sum = matrixi[j];
91 for (label
k = 0;
k < i; ++
k)
93 sum -= matrixi[
k]*matrix(
k, j);
100 scalar largestCoeff = 0.0;
101 for (label i = j; i < size; ++i)
103 scalar* __restrict__ matrixi = matrix[i];
104 scalar
sum = matrixi[j];
106 for (label
k = 0;
k < j; ++
k)
108 sum -= matrixi[
k]*matrix(
k, j);
114 if ((temp = vv[i]*
mag(
sum)) >= largestCoeff)
121 pivotIndices[j] = iMax;
125 scalar* __restrict__ matrixiMax = matrix[iMax];
127 for (label
k = 0;
k < size; ++
k)
129 std::swap(matrixj[
k], matrixiMax[
k]);
136 if (matrixj[j] == 0.0)
143 scalar rDiag = 1.0/matrixj[j];
145 for (label i = j + 1; i < size; ++i)
147 matrix(i, j) *= rDiag;
157 const label size = matrix.m();
160 for (label j = 0; j < size; ++j)
162 for (label
k = j + 1;
k < size; ++
k)
168 for (label j = 0; j < size; ++j)
172 for (label
k = 0;
k < j; ++
k)
176 for (label i = 0; i <
k; ++i)
178 s += matrix(i,
k)*matrix(i, j);
181 s = (matrix(j,
k) -
s)/matrix(
k,
k);
189 d = matrix(j, j) - d;
194 <<
"Matrix is not symmetric positive-definite. Unable to " 199 matrix(j, j) =
sqrt(d);
217 <<
"A and B must have identical inner dimensions but A.n = " 218 <<
A.n() <<
" and B.m = " <<
B.m()
225 <<
"B and C must have identical inner dimensions but B.n = " 226 <<
B.n() <<
" and C.m = " <<
C.m()
232 for (label i = 0; i <
A.m(); ++i)
234 for (label
g = 0;
g <
C.n(); ++
g)
236 for (label l = 0; l <
C.m(); ++l)
239 for (label j = 0; j <
A.n(); ++j)
241 ab +=
A(i, j)*
B(j, l);
243 ans(i,
g) +=
C(l,
g) * ab;
254 const DiagonalMatrix<scalar>&
B,
258 if (
A.n() !=
B.size())
261 <<
"A and B must have identical inner dimensions but A.n = " 262 <<
A.n() <<
" and B.m = " <<
B.size()
266 if (
B.size() !=
C.m())
269 <<
"B and C must have identical inner dimensions but B.n = " 270 <<
B.size() <<
" and C.m = " <<
C.m()
276 for (label i = 0; i <
A.m(); ++i)
278 for (label
g = 0;
g <
C.n(); ++
g)
280 for (label l = 0; l <
C.m(); ++l)
282 ans(i,
g) +=
C(l,
g) *
A(i, l)*
B[l];
293 const DiagonalMatrix<scalar>&
B,
297 if (
A.m() !=
B.size())
300 <<
"A and B must have identical dimensions but A.m = " 301 <<
A.m() <<
" and B.m = " <<
B.size()
305 if (
B.size() !=
C.m())
308 <<
"B and C must have identical dimensions but B.m = " 309 <<
B.size() <<
" and C.m = " <<
C.m()
313 const label size =
A.m();
317 for (label i = 0; i < size; ++i)
319 for (label
g = 0;
g < size; ++
g)
321 for (label l = 0; l < size; ++l)
323 ans(i,
g) +=
C(l,
g)*
A(i, l)*
B[l];
337 SVD svd(
A, minCondition);
338 return svd.VSinvUt();
dimensionedScalar sign(const dimensionedScalar &ds)
scalarRectangularMatrix SVDinv(const scalarRectangularMatrix &A, scalar minCondition=0)
Return the inverse of matrix A using SVD.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
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.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
label k
Boltzmann constant.
errorManip< error > abort(error &err)
const uniformDimensionedVectorField & g
SymmetricSquareMatrix< scalar > scalarSymmetricSquareMatrix
A class representing the concept of 0 (zero) that can be used to avoid manipulating objects known to ...
List< label > labelList
A List of labels.
void multiply(DimensionedField< Type, GeoMesh > &result, const DimensionedField< Type, GeoMesh > &f1, const DimensionedField< scalar, GeoMesh > &f2)
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
SquareMatrix< scalar > scalarSquareMatrix
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
RectangularMatrix< scalar > scalarRectangularMatrix
static constexpr const zero Zero
Global zero (0)