53 for (
label i=0; i<Un; i++)
63 scale +=
mag(U_(
k, i));
71 s += U_(
k, i)*U_(
k, i);
75 g = -sign(
sqrt(s), f);
79 for (
label j=l-1; j<Un; j++)
84 s += U_(
k, i)*U_(
k, j);
90 U_(
k, j) += f*U_(
k, i);
105 if (i+1 <= Um && i != Un)
109 scale +=
mag(U_(i,
k));
117 s += U_(i,
k)*U_(i,
k);
120 scalar
f = U_(i, l-1);
121 g = -sign(
sqrt(s), f);
130 for (
label j=l-1; j<Um; j++)
135 s += U_(j,
k)*U_(i,
k);
140 U_(j,
k) += s*rv1[
k];
150 anorm =
max(anorm,
mag(S_[i]) +
mag(rv1[i]));
155 for (
label i=Un-1; i >= 0; i--)
161 for (
label j=l; j<Un; j++)
163 V_(j, i) = (U_(i, j)/U_(i, l))/g;
166 for (
label j=l; j<Un; j++)
171 s += U_(i,
k)*V_(
k, j);
176 V_(
k, j) += s*V_(
k, i);
181 for (
label j=l; j<Un;j++)
183 V_(i, j) = V_(j, i) = 0;
192 for (
label i=
min(Un, Um) - 1; i>=0; i--)
197 for (
label j=l; j<Un; j++)
206 for (
label j=l; j<Un; j++)
211 s += U_(
k, i)*U_(
k, j);
214 scalar
f = (s/U_(i, i))*g;
218 U_(
k, j) += f*U_(
k, i);
222 for (
label j=i; j<Um; j++)
229 for (
label j=i; j<Um; j++)
240 for (
label its = 0; its < 30; its++)
245 for (l =
k; l >= 0; l--)
249 if (l == 0 ||
mag(rv1[l]) <= anorm)
255 if (
mag(S_[mn]) <= anorm)
265 for (
label i=l; i<
k+1; i++)
282 for (
label j=0; j<Um; j++)
284 scalar
y = U_(j, mn);
286 U_(j, mn) = y*c + z*
s;
287 U_(j, i) = z*c - y*
s;
299 for (
label j=0; j<Un; j++)
301 V_(j,
k) = -V_(j,
k);
317 scalar
f = ((y - z)*(y + z) + (g -
h)*(g + h))/(2*h*y);
319 f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) -
h))/
x;
323 for (
label j=l; j <= mn; j++)
339 for (
label jj = 0; jj < Un; jj++)
343 V_(jj, j) = x*c + z*
s;
344 V_(jj, i) = z*c - x*
s;
358 for (
label jj=0; jj < Um; jj++)
362 U_(jj, j) = y*c + z*
s;
363 U_(jj, i) = z*c - y*
s;
373 const scalar minS = minCondition*S_[
findMax(S_)];
label n() const
Return the number of columns.
#define forAll(list, i)
Loop across all elements in list.
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.
SVD(const scalarRectangularMatrix &A, const scalar minCondition=0)
Construct from a rectangular Matrix.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Form T() const
Return the transpose of the matrix.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
label k
Boltzmann constant.
const dimensionedScalar & c
Speed of light in a vacuum.
Various functions to operate on Lists.
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))
Scalar sqrtSumSqr(const Scalar a, const Scalar b)
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
scalarRectangularMatrix VSinvUt() const
Return the matrix product V S^(-1) U^T (the pseudo inverse)
label m() const
Return the number of rows.
dimensioned< scalar > mag(const dimensioned< Type > &)
const dimensionedScalar & h
Planck constant.
const dimensionedVector & g