45 const label Un = U_.
n();
46 const label Um = U_.
m();
55 for (label i=0; i<Un; i++)
63 for (label
k=i;
k<Um;
k++)
65 scale +=
mag(U_(
k, i));
70 for (label
k=i;
k<Um;
k++)
73 s += U_(
k, i)*U_(
k, i);
81 for (label j=l-1; j<Un; j++)
84 for (label
k=i;
k<Um;
k++)
86 s += U_(
k, i)*U_(
k, j);
90 for (label
k=i;
k<
A.m();
k++)
92 U_(
k, j) +=
f*U_(
k, i);
96 for (label
k=i;
k<Um;
k++)
107 if (i+1 <= Um && i != Un)
109 for (label
k=l-1;
k<Un;
k++)
111 scale +=
mag(U_(i,
k));
116 for (label
k=l-1;
k<Un;
k++)
119 s += U_(i,
k)*U_(i,
k);
122 scalar
f = U_(i, l-1);
127 for (label
k=l-1;
k<Un;
k++)
132 for (label j=l-1; j<Um; j++)
135 for (label
k=l-1;
k<Un;
k++)
137 s += U_(j,
k)*U_(i,
k);
140 for (label
k=l-1;
k<Un;
k++)
142 U_(j,
k) +=
s*rv1[
k];
145 for (label
k=l-1;
k<Un;
k++)
152 anorm =
max(anorm,
mag(S_[i]) +
mag(rv1[i]));
157 for (label i=Un-1; i >= 0; i--)
163 for (label j=l; j<Un; j++)
165 V_(j, i) = (U_(i, j)/U_(i, l))/
g;
168 for (label j=l; j<Un; j++)
171 for (label
k=l;
k<Un;
k++)
173 s += U_(i,
k)*V_(
k, j);
176 for (label
k=l;
k<Un;
k++)
178 V_(
k, j) +=
s*V_(
k, i);
183 for (label j=l; j<Un;j++)
185 V_(i, j) = V_(j, i) = 0;
194 for (label i=
min(Un, Um) - 1; i>=0; i--)
199 for (label j=l; j<Un; j++)
208 for (label j=l; j<Un; j++)
211 for (label
k=l;
k<Um;
k++)
213 s += U_(
k, i)*U_(
k, j);
216 scalar
f = (
s/U_(i, i))*
g;
218 for (label
k=i;
k<Um;
k++)
220 U_(
k, j) +=
f*U_(
k, i);
224 for (label j=i; j<Um; j++)
231 for (label j=i; j<Um; j++)
240 for (label
k=Un-1;
k >= 0;
k--)
242 for (label its = 0; its < 30; its++)
247 for (l =
k; l >= 0; l--)
251 if (l == 0 ||
mag(rv1[l]) <= anorm)
257 if (
mag(S_[mn]) <= anorm)
267 for (label i=l; i<
k+1; i++)
284 for (label j=0; j<Um; j++)
286 scalar
y = U_(j, mn);
288 U_(j, mn) =
y*
c + z*
s;
289 U_(j, i) = z*
c -
y*
s;
301 for (label j=0; j<Un; j++)
303 V_(j,
k) = -V_(j,
k);
319 scalar
f = ((
y - z)*(
y + z) + (
g -
h)*(
g +
h))/(2*
h*
y);
321 f = ((
x - z)*(
x + z) +
h*((
y/(
f + sign(
g,
f))) -
h))/
x;
325 for (label j=l; j <= mn; j++)
341 for (label jj = 0; jj < Un; jj++)
345 V_(jj, j) =
x*
c + z*
s;
346 V_(jj, i) = z*
c -
x*
s;
360 for (label jj=0; jj < Um; jj++)
364 U_(jj, j) =
y*
c + z*
s;
365 U_(jj, i) = z*
c -
y*
s;
375 const scalar minS = minCondition*S_[
findMax(S_)];
label findMax(const ListType &input, label start=0)
Linear search for the index of the max element, similar to std::max_element but for lists and returns...
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
label k
Boltzmann constant.
Various functions to operate on Lists.
#define forAll(list, i)
Loop across all elements in list.
label n() const noexcept
The number of columns.
Scalar sqrtSumSqr(const Scalar a, const Scalar b)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
const uniformDimensionedVectorField & g
dimensioned< Type > T() const
Return transpose.
scalarRectangularMatrix VSinvUt() const
Return the matrix product V S^(-1) U^T (the pseudo inverse)
label m() const noexcept
The number of rows.
const dimensionedScalar h
Planck constant.
const dimensionedScalar c
Speed of light in a vacuum.
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)
RectangularMatrix< scalar > scalarRectangularMatrix