36 template<
class cmptType>
39 for (label j = 0; j < n_; ++j)
41 EValsRe_[j] = EVecs_[n_ - 1][j];
45 for (label i = n_ - 1; i > 0; --i)
48 cmptType scale =
Zero;
51 for (label
k = 0;
k < i; ++
k)
53 scale = scale +
mag(EValsRe_[
k]);
58 EValsIm_[i] = EValsRe_[i - 1];
60 for (label j = 0; j < i; ++j)
62 EValsRe_[j] = EVecs_[i - 1][j];
70 for (label
k = 0;
k < i; ++
k)
73 h += EValsRe_[
k]*EValsRe_[
k];
76 cmptType
f = EValsRe_[i - 1];
84 EValsIm_[i] = scale*
g;
86 EValsRe_[i - 1] =
f -
g;
88 for (label j = 0; j < i; ++j)
94 for (label j = 0; j < i; ++j)
98 g = EValsIm_[j] + EVecs_[j][j]*
f;
100 for (label
k = j + 1;
k <= i - 1; ++
k)
102 g += EVecs_[
k][j]*EValsRe_[
k];
103 EValsIm_[
k] += EVecs_[
k][j]*
f;
111 for (label j = 0; j < i; ++j)
114 f += EValsIm_[j]*EValsRe_[j];
117 cmptType hh =
f/(2.0*
h);
119 for (label j = 0; j < i; ++j)
121 EValsIm_[j] -= hh*EValsRe_[j];
124 for (label j = 0; j < i; ++j)
129 for (label
k = j;
k <= i - 1; ++
k)
132 (
f*EValsIm_[
k] +
g*EValsRe_[
k]);
135 EValsRe_[j] = EVecs_[i - 1][j];
144 for (label i = 0; i < n_ - 1; ++i)
146 EVecs_[n_ - 1][i] = EVecs_[i][i];
147 EVecs_[i][i] = cmptType(1);
148 cmptType
h = EValsRe_[i + 1];
152 for (label
k = 0;
k <= i; ++
k)
154 EValsRe_[
k] = EVecs_[
k][i + 1]/
h;
157 for (label j = 0; j <= i; ++j)
161 for (label
k = 0;
k <= i; ++
k)
163 g += EVecs_[
k][i + 1]*EVecs_[
k][j];
166 for (label
k = 0;
k <= i; ++
k)
168 EVecs_[
k][j] -=
g*EValsRe_[
k];
173 for (label
k = 0;
k <= i; ++
k)
175 EVecs_[
k][i + 1] =
Zero;
179 for (label j = 0; j < n_; ++j)
181 EValsRe_[j] = EVecs_[n_ - 1][j];
182 EVecs_[n_ - 1][j] =
Zero;
185 EVecs_[n_ - 1][n_ - 1] = cmptType(1);
190 template<
class cmptType>
193 for (label i = 1; i < n_; ++i)
195 EValsIm_[i - 1] = EValsIm_[i];
198 EValsIm_[n_-1] =
Zero;
201 cmptType tst1 =
Zero;
202 cmptType eps =
pow(2.0, -52.0);
204 for (label l = 0; l < n_; l++)
207 tst1 =
max(tst1,
mag(EValsRe_[l]) +
mag(EValsIm_[l]));
213 if (
mag(EValsIm_[m]) <= eps*tst1)
227 cmptType
g = EValsRe_[l];
228 cmptType
p = (EValsRe_[l + 1] -
g)/(2.0*EValsIm_[l]);
236 EValsRe_[l] = EValsIm_[l]/(
p + r);
237 EValsRe_[l + 1] = EValsIm_[l]*(
p + r);
238 cmptType dl1 = EValsRe_[l + 1];
239 cmptType
h =
g - EValsRe_[l];
241 for (label i = l + 2; i < n_; ++i)
250 cmptType
c = cmptType(1);
253 cmptType el1 = EValsIm_[l + 1];
257 for (label i = m - 1; i >= l; --i)
265 EValsIm_[i + 1] =
s*r;
268 p =
c*EValsRe_[i] -
s*
g;
269 EValsRe_[i + 1] =
h +
s*(
c*
g +
s*EValsRe_[i]);
272 for (label
k = 0;
k < n_; ++
k)
274 h = EVecs_[
k][i + 1];
275 EVecs_[
k][i + 1] =
s*EVecs_[
k][i] +
c*
h;
276 EVecs_[
k][i] =
c*EVecs_[
k][i] -
s*
h;
280 p = -
s*s2*c3*el1*EValsIm_[l]/dl1;
284 while (
mag(EValsIm_[l]) > eps*tst1);
287 EValsRe_[l] = EValsRe_[l] +
f;
292 for (label i = 0; i < n_ - 1; ++i)
295 cmptType
p = EValsRe_[i];
297 for (label j = i + 1; j < n_; ++j)
308 EValsRe_[
k] = EValsRe_[i];
311 for (label j = 0; j < n_; ++j)
314 EVecs_[j][i] = EVecs_[j][
k];
322 template<
class cmptType>
325 DiagonalMatrix<cmptType> orth_(n_);
330 for (label m = low + 1; m <= high - 1; ++m)
333 cmptType scale =
Zero;
335 for (label i = m; i <= high; ++i)
337 scale = scale +
mag(H_[i][m - 1]);
345 for (label i = high; i >= m; --i)
347 orth_[i] = H_[i][m - 1]/scale;
348 h += orth_[i]*orth_[i];
363 for (label j = m; j < n_; ++j)
367 for (label i = high; i >= m; --i)
369 f += orth_[i]*H_[i][j];
374 for (label i = m; i <= high; ++i)
376 H_[i][j] -=
f*orth_[i];
380 for (label i = 0; i <= high; ++i)
384 for (label j = high; j >= m; --j)
386 f += orth_[j]*H_[i][j];
391 for (label j = m; j <= high; ++j)
393 H_[i][j] -=
f*orth_[j];
397 orth_[m] = scale*orth_[m];
398 H_[m][m-1] = scale*
g;
403 for (label i = 0; i < n_; ++i)
405 for (label j = 0; j < n_; ++j)
407 EVecs_[i][j] = (i == j ? 1.0 : 0.0);
411 for (label m = high - 1; m >= low + 1; --m)
413 if (H_[m][m - 1] != 0.0)
415 for (label i = m + 1; i <= high; ++i)
417 orth_[i] = H_[i][m - 1];
420 for (label j = m; j <= high; ++j)
424 for (label i = m; i <= high; ++i)
426 g += orth_[i]*EVecs_[i][j];
430 g = (
g/orth_[m])/H_[m][m - 1];
432 for (label i = m; i <= high; ++i)
434 EVecs_[i][j] +=
g*orth_[i];
442 template<
class cmptType>
450 cmptType eps =
pow(2.0, -52.0);
451 cmptType exshift =
Zero;
452 cmptType
p = 0, q = 0, r = 0,
s = 0, z = 0, t, w,
x,
y;
455 cmptType norm =
Zero;
457 for (label i = 0; i < nn; ++i)
459 if ((i < low) || (i > high))
461 EValsRe_[i] = H_[i][i];
465 for (label j =
max(i - 1, 0); j < nn; ++j)
467 norm +=
mag(H_[i][j]);
480 s =
mag(H_[l - 1][l - 1]) +
mag(H_[l][l]);
487 if (
mag(H_[l][l - 1]) < eps*
s)
498 H_[
n][
n] = H_[
n][
n] + exshift;
499 EValsRe_[
n] = H_[
n][
n];
506 w = H_[
n][
n - 1]*H_[
n - 1][
n];
507 p = (H_[
n - 1][
n - 1] - H_[
n][
n])/2.0;
511 H_[
n - 1][
n - 1] += exshift;
526 EValsRe_[
n - 1] =
x + z;
527 EValsRe_[
n] = EValsRe_[
n - 1];
531 EValsRe_[
n] =
x - w/z;
534 EValsIm_[
n - 1] =
Zero;
545 for (label j =
n - 1; j < nn; ++j)
548 H_[
n - 1][j] = q*z +
p*H_[
n][j];
549 H_[
n][j] = q*H_[
n][j] -
p*z;
553 for (label i = 0; i <=
n; ++i)
556 H_[i][
n - 1] = q*z +
p*H_[i][
n];
557 H_[i][
n] = q*H_[i][
n] -
p*z;
561 for (label i = low; i <= high; ++i)
563 z = EVecs_[i][
n - 1];
564 EVecs_[i][
n - 1] = q*z +
p*EVecs_[i][
n];
565 EVecs_[i][
n] = q*EVecs_[i][
n] -
p*z;
570 EValsRe_[
n - 1] =
x +
p;
588 y = H_[
n - 1][
n - 1];
589 w = H_[
n][
n - 1]*H_[
n - 1][
n];
596 for (label i = low; i <=
n; ++i)
622 s =
x - w/((
y -
x)/2.0 +
s);
624 for (label i = low; i <=
n; ++i)
646 p = (r*
s - w)/H_[m + 1][m] + H_[m][m + 1];
647 q = H_[m + 1][m + 1] - z - r -
s;
648 r = H_[m + 2][m + 1];
662 < eps*(
mag(
p)*(
mag(H_[m - 1][m - 1])
663 +
mag(z) +
mag(H_[m + 1][m + 1])))
672 for (label i = m + 2; i <=
n; ++i)
683 for (label
k = m;
k <=
n - 1; ++
k)
685 label notlast = (
k !=
n - 1);
690 q = H_[
k + 1][
k - 1];
691 r = (notlast ? H_[
k + 2][
k - 1] : 0.0);
722 H_[
k][
k - 1] = -H_[
k][
k - 1];
733 for (label j =
k; j < nn; ++j)
735 p = H_[
k][j] + q*H_[
k + 1][j];
748 for (label i = 0; i <=
min(
n,
k + 3); ++i)
750 p =
x*H_[i][
k] +
y*H_[i][
k + 1];
763 for (label i = low; i <= high; ++i)
765 p =
x*EVecs_[i][
k] +
y*EVecs_[i][
k + 1];
769 p += z*EVecs_[i][
k + 2];
770 EVecs_[i][
k + 2] -=
p*r;
774 EVecs_[i][
k + 1] -=
p*q;
787 for (
n = nn-1;
n >= 0; --
n)
796 H_[
n][
n] = cmptType(1);
798 for (label i =
n-1; i >= 0; --i)
803 for (label j = l; j <=
n; ++j)
805 r += H_[i][j]*H_[j][
n];
808 if (EValsIm_[i] < 0.0)
817 if (EValsIm_[i] == 0.0)
825 H_[i][
n] = -r/(eps*norm);
832 q = (EValsRe_[i] -
p)*(EValsRe_[i] -
p)
833 + EValsIm_[i]*EValsIm_[i];
840 H_[i + 1][
n] = (-r - w*t)/
x;
844 H_[i + 1][
n] = (-
s -
y*t)/z;
853 for (label j = i; j <=
n; ++j)
868 H_[
n - 1][
n - 1] = q/H_[
n][
n - 1];
869 H_[
n - 1][
n] = -(H_[
n][
n] -
p)/H_[
n][
n - 1];
876 H_[
n - 1][
n - 1] = cDiv.real();
877 H_[
n - 1][
n] = cDiv.imag();
881 H_[
n][
n] = cmptType(1);
883 for (label i =
n - 2; i >= 0; --i)
885 cmptType ra, sa, vr, vi;
889 for (label j = l; j <=
n; ++j)
891 ra += H_[i][j]*H_[j][
n-1];
892 sa += H_[i][j]*H_[j][
n];
897 if (EValsIm_[i] < 0.0)
907 if (EValsIm_[i] == 0)
910 H_[i][
n - 1] = cDiv.real();
911 H_[i][
n] = cDiv.imag();
918 vr = (EValsRe_[i] -
p)*(EValsRe_[i] -
p)
919 + EValsIm_[i]*EValsIm_[i] - q*q;
921 vi = 2.0*(EValsRe_[i] -
p)*q;
923 if ((vr == 0.0) && (vi == 0.0))
933 H_[i][
n - 1] = cDiv.real();
934 H_[i][
n] = cDiv.imag();
938 H_[i + 1][
n - 1] = (-ra - w*H_[i][
n - 1]
941 H_[i + 1][
n] = (-sa - w*H_[i][
n]
949 H_[i + 1][
n - 1] = cDiv.real();
950 H_[i + 1][
n] = cDiv.imag();
959 for (label j = i; j <=
n; ++j)
971 for (label i = 0; i < nn; ++i)
973 if (i < low || i > high)
975 for (label j = i; j < nn; ++j)
977 EVecs_[i][j] = H_[i][j];
983 for (label j = nn - 1; j >= low; --j)
985 for (label i = low; i <= high; ++i)
989 for (label
k = low;
k <=
min(j, high); ++
k)
991 z = z + EVecs_[i][
k]*H_[
k][j];
1002 template<
class cmptType>
1014 <<
"Input matrix has zero size." 1021 tridiagonaliseSymmMatrix();
1022 symmTridiagonalQL();
1033 template<
class cmptType>
1036 const SquareMatrix<cmptType>&
A,
1049 <<
"Input matrix has zero size." 1056 tridiagonaliseSymmMatrix();
1057 symmTridiagonalQL();
1070 template<
class cmptType>
1081 [&](
const cmptType&
x){
return complex(
x); }
1084 for (label i = 0; i < EValsIm_.size(); ++i)
1086 if (
mag(EValsIm_[i]) > VSMALL)
1088 for (label j = 0; j < EVecs.m(); ++j)
1090 EVecs(j, i) =
complex(EVecs_(j, i), EVecs_(j, i+1));
1091 EVecs(j, i+1) = EVecs(j, i).conjugate();
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
EigenMatrix()=delete
No default construct.
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.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
dimensionedScalar sqrt(const dimensionedScalar &ds)
label k
Boltzmann constant.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
const uniformDimensionedVectorField & g
const SquareMatrix< complex > complexEVecs() const
Return right eigenvectors in unpacked form.
const dimensionedScalar h
Planck constant.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const dimensionedScalar c
Speed of light in a vacuum.
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))
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
EigenMatrix (i.e. eigendecomposition or spectral decomposition) decomposes a diagonalisable nonsymmet...
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix...
static constexpr const zero Zero
Global zero (0)