diff --git a/applications/test/matrices/EigenMatrix/Make/files b/applications/test/matrices/EigenMatrix/Make/files
new file mode 100644
index 0000000000..1bd5722a5f
--- /dev/null
+++ b/applications/test/matrices/EigenMatrix/Make/files
@@ -0,0 +1,3 @@
+Test-EigenMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-EigenMatrix
diff --git a/applications/test/matrices/EigenMatrix/Make/options b/applications/test/matrices/EigenMatrix/Make/options
new file mode 100644
index 0000000000..b23590707c
--- /dev/null
+++ b/applications/test/matrices/EigenMatrix/Make/options
@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/matrices/EigenMatrix/Test-EigenMatrix.C b/applications/test/matrices/EigenMatrix/Test-EigenMatrix.C
new file mode 100644
index 0000000000..bae8bc7630
--- /dev/null
+++ b/applications/test/matrices/EigenMatrix/Test-EigenMatrix.C
@@ -0,0 +1,556 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is derivative work of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Application
+ Test-EigenMatrix
+
+Description
+ Tests for \c EigenMatrix constructors, and member functions
+ using \c floatScalar, and \c doubleScalar base types.
+
+ Cross-checks were obtained from 'NumPy 1.15.1' if no theoretical
+ cross-check exists (like eigendecomposition relations), and
+ were hard-coded for elementwise comparisons.
+
+\*---------------------------------------------------------------------------*/
+
+#include "scalarMatrices.H"
+#include "RectangularMatrix.H"
+#include "SquareMatrix.H"
+#include "complex.H"
+#include "IOmanip.H"
+#include "EigenMatrix.H"
+#include "TestTools.H"
+#include
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Return the absolute tolerance value for bitwise comparisons of floatScalars
+floatScalar getTol(floatScalar)
+{
+ return 1e-2;
+}
+
+
+// Return the absolute tolerance value for bitwise comparisons of doubleScalars
+doubleScalar getTol(doubleScalar)
+{
+ return 1e-10;
+}
+
+
+// Create each constructor of EigenMatrix, and print output
+template
+void test_constructors(Type)
+{
+ {
+ Info<< "# Construct from a SquareMatrix" << nl;
+ const SquareMatrix A(5, Zero);
+ const EigenMatrix EM(A);
+ }
+ {
+ Info<< "# Construct from a SquareMatrix and symmetry flag" << nl;
+ const SquareMatrix A(5, Zero);
+ const EigenMatrix EM(A, true);
+ }
+}
+
+
+// Execute each member function of EigenMatrix, and print output
+template
+void test_member_funcs(Type)
+{
+ SquareMatrix A(3, Zero);
+ assignMatrix
+ (
+ A,
+ {
+ Type(1), Type(2), Type(3),
+ Type(4), Type(5), Type(6),
+ Type(7), Type(8), Type(9)
+ }
+ );
+
+ Info<< "# Operand: " << nl
+ << " SquareMatrix = " << A << endl;
+
+
+ // Since eigenvalues are unique and eigenvectors are not unique,
+ // the bitwise comparisons are limited to eigenvalue computations.
+ // Here, only the execution of the functions is tested, rather than
+ // the verification of the eigendecomposition through theoretical relations.
+ {
+ const EigenMatrix EM(A);
+ const DiagonalMatrix EValsRe = EM.EValsRe();
+ const DiagonalMatrix& EValsIm = EM.EValsIm();
+ const SquareMatrix& EVecs = EM.EVecs();
+
+ cmp
+ (
+ " Return real eigenvalues or real part of complex eigenvalues = ",
+ EValsRe,
+ List
+ ({
+ Type(16.116844),
+ Type(-1.116844),
+ Type(0)
+ }),
+ getTol(Type(0)),
+ 1e-6
+ );
+ cmp
+ (
+ " Return zero-matrix for real eigenvalues "
+ "or imaginary part of complex eigenvalues = ",
+ EValsIm,
+ List(3, Zero)
+ );
+ Info<< " Return eigenvectors matrix = " << EVecs << endl;
+ }
+}
+
+
+// Test the relation: "sum(eigenvalues) = trace(A)"
+// w.wiki/4zs (Retrieved: 16-06-19) # Item-1
+template
+void test_eigenvalues_sum
+(
+ const SquareMatrix& A,
+ const DiagonalMatrix& EValsRe
+)
+{
+ const Type trace = A.trace();
+ // Imaginary part of complex conjugates cancel each other
+ const Type EValsSum = sum(EValsRe);
+ Info<< " # A.mRows = " << A.m() << nl;
+ cmp
+ (
+ " # sum(eigenvalues) = trace(A) = ",
+ EValsSum,
+ trace,
+ getTol(Type(0))
+ );
+}
+
+
+// Test the relation: "prod(eigenvalues) = det(A)"
+// w.wiki/4zs (Retrieved: 16-06-19) # Item-2
+// Note that the determinant computation may fail
+// which is not a suggestion that eigendecomposition fails
+template
+void test_eigenvalues_prod
+(
+ const SquareMatrix& A,
+ const DiagonalMatrix& EValsRe,
+ const DiagonalMatrix& EValsIm
+)
+{
+ const Type determinant = mag(det(A));
+ Type EValsProd = Type(1);
+
+ if (EValsIm.empty())
+ {
+ for (label i = 0; i < EValsRe.size(); ++i)
+ {
+ EValsProd *= Foam::sqrt(sqr(EValsRe[i]));
+ }
+ }
+ else
+ {
+ for (label i = 0; i < EValsRe.size(); ++i)
+ {
+ EValsProd *= Foam::sqrt(sqr(EValsRe[i]) + sqr(EValsIm[i]));
+ }
+ }
+ cmp
+ (
+ " # prod(eigenvalues) = det(A) = ",
+ EValsProd,
+ determinant,
+ getTol(Type(0))
+ );
+}
+
+
+// Test eigenvalues in eigendecomposition relations
+// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307)
+template
+void test_eigenvalues(Type)
+{
+ Random rndGen(1234);
+ const label numberOfTests = 20;
+
+ // Non-symmetric
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(100, 200);
+ const labelPair m(mRows, mRows);
+
+ const SquareMatrix A
+ (
+ makeRandomMatrix>(m, rndGen)
+ );
+
+ const EigenMatrix EM(A);
+ const DiagonalMatrix& EValsRe = EM.EValsRe();
+
+ test_eigenvalues_sum(A, EValsRe);
+
+ // LUDecompose does not work with floatScalar at the time of writing,
+ // hence det function. Once LUDecompose is refactored, comment out below
+ // const DiagonalMatrix& EValsIm = EM.EValsIm();
+ // test_eigenvalues_prod(A, EValsRe, EValsIm);
+ }
+
+ // Symmetric
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(100, 200);
+ const labelPair m(mRows, mRows);
+
+ SquareMatrix A
+ (
+ makeRandomMatrix>(m, rndGen)
+ );
+ // Symmetrise with noise
+ for (label n = 0; n < A.n() - 1; ++n)
+ {
+ for (label m = A.m() - 1; m > n; --m)
+ {
+ A(n, m) = A(m, n) + SMALL;
+ }
+ }
+
+ const bool symmetric = true;
+ const EigenMatrix EM(A, symmetric);
+ const DiagonalMatrix& EValsRe = EM.EValsRe();
+
+ test_eigenvalues_sum(A, EValsRe);
+ }
+}
+
+
+// Test the relation: "(A & EVec - EVal*EVec) = 0"
+template
+void test_characteristic_eq
+(
+ const SquareMatrix& Aorig,
+ const DiagonalMatrix& EValsRe,
+ const DiagonalMatrix& EValsIm,
+ const SquareMatrix& EVecs
+)
+{
+ SquareMatrix A(Aorig.m());
+ auto convertToComplex = [&](const scalar& val) { return complex(val); };
+ std::transform
+ (
+ Aorig.cbegin(),
+ Aorig.cend(),
+ A.begin(),
+ convertToComplex
+ );
+
+ for (label i = 0; i < A.m(); ++i)
+ {
+ const RectangularMatrix& EVec(EVecs.subColumn(i));
+ const complex EVal(EValsRe[i], EValsIm[i]);
+ const RectangularMatrix leftSide(A*EVec);
+ const RectangularMatrix rightSide(EVal*EVec);
+
+ cmp
+ (
+ " # (A & EVec - EVal*EVec) = 0:",
+ flt(leftSide),
+ flt(rightSide),
+ getTol(Type(0))
+ );
+ }
+}
+
+
+// Test eigenvectors in eigendecomposition relations
+template
+void test_eigenvectors(Type)
+{
+ Random rndGen(1234);
+ const label numberOfTests = 20;
+
+ // Non-symmetric
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(100, 200);
+ const labelPair m(mRows, mRows);
+
+ const SquareMatrix A
+ (
+ makeRandomMatrix>(m, rndGen)
+ );
+
+ const EigenMatrix EM(A);
+ const DiagonalMatrix& EValsRe = EM.EValsRe();
+ const DiagonalMatrix& EValsIm = EM.EValsIm();
+ const SquareMatrix EVecs(EM.complexEVecs());
+
+ test_characteristic_eq(A, EValsRe, EValsIm, EVecs);
+ }
+
+ // Symmetric
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(100, 200);
+ const labelPair m(mRows, mRows);
+
+ SquareMatrix A
+ (
+ makeRandomMatrix>(m, rndGen)
+ );
+ // Symmetrise with noise
+ for (label n = 0; n < A.n() - 1; ++n)
+ {
+ for (label m = A.m() - 1; m > n; --m)
+ {
+ A(n, m) = A(m, n) + SMALL;
+ }
+ }
+
+ const bool symmetric = true;
+ const EigenMatrix EM(A, symmetric);
+ const DiagonalMatrix& EValsRe = EM.EValsRe();
+ const DiagonalMatrix& EValsIm = EM.EValsIm();
+ const SquareMatrix EVecs(EM.complexEVecs());
+
+ test_characteristic_eq(A, EValsRe, EValsIm, EVecs);
+ }
+}
+
+
+// Do compile-time recursion over the given types
+template
+inline typename std::enable_if::type
+run_tests(const std::tuple& types, const List& typeID){}
+
+
+template
+inline typename std::enable_if::type
+run_tests(const std::tuple& types, const List& typeID)
+{
+ Info<< nl << " ## Test constructors: "<< typeID[I] <<" ##" << nl;
+ test_constructors(std::get(types));
+
+ Info<< nl << " ## Test member functions: "<< typeID[I] <<" ##" << nl;
+ test_member_funcs(std::get(types));
+
+ Info<< nl << " ## Test eigenvalues: "<< typeID[I] <<" ##" << nl;
+ test_eigenvalues(std::get(types));
+
+ Info<< nl << " ## Test eigenvectors: "<< typeID[I] <<" ##" << nl;
+ test_eigenvectors(std::get(types));
+
+ run_tests(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
+
+int main()
+{
+ Info<< setprecision(15);
+
+ const std::tuple types
+ (
+ std::make_tuple(Zero, Zero)
+ );
+
+ const List typeID
+ ({
+ "SquareMatrix",
+ "SquareMatrix"
+ });
+
+ run_tests(types, typeID);
+
+
+ Info<< nl << " ## Test corner cases ##" << endl;
+ {
+ Info<< nl << " ## Rosser et al. (1951) matrix: ##" << nl;
+ // Rosser, J. B., Lanczos, C., Hestenes, M. R., & Karush, W. (1951).
+ // Separation of close eigenvalues of a real symmetric matrix.
+ // Jour. Research of the National Bureau of Standards, 47(4), 291-297.
+ // DOI:10.6028/jres.047.037
+ //
+ // 8x8 symmetric square matrix consisting of close real eigenvalues
+ // ibid, p. 294
+ // {
+ // 1020.04901843, 1020.000, 1019.90195136, 1000.000,
+ // 1000.000, 0.09804864072, 0.000, -1020.0490
+ // }
+ // Note that prod(eigenvalues) != determinant(A) for this matrix
+ // via the LAPACK routine z/dgetrf
+
+ SquareMatrix A(8, Zero);
+ assignMatrix
+ (
+ A,
+ {
+ 611, 196, -192, 407, -8, -52, -49, 29,
+ 196, 899, 113, -192, -71, -43, -8, -44,
+ -192, 113, 899, 196, 61, 49, 8, 52,
+ 407, -192, 196, 611, 8, 44, 59, -23,
+ -8, -71, 61, 8, 411, -599, 208, 208,
+ -52, -43, 49, 44, -599, 411, 208, 208,
+ -49, -8, 8, 59, 208, 208, 99, -911,
+ 29, -44, 52, -23, 208, 208, -911, 99
+ }
+ );
+ const EigenMatrix EM(A);
+ const DiagonalMatrix& EValsRe = EM.EValsRe();
+ const DiagonalMatrix& EValsIm = EM.EValsIm();
+
+ test_eigenvalues_sum(A, EValsRe);
+
+ cmp
+ (
+ " # Rosser et al. (1951) case, EValsRe = ",
+ EValsRe,
+ List // theoretical EValsRe
+ ({
+ -1020.0490, 0.000, 0.09804864072, 1000.000,
+ 1000.000, 1019.90195136, 1020.000, 1020.04901843
+ }),
+ 1e-3
+ );
+ cmp
+ (
+ " # Rosser et al. (1951) case, EValsIm = ",
+ EValsIm,
+ List(8, Zero)
+ );
+ }
+ {
+ Info<< nl << " ## Test eigenvector unpacking: ##" << nl;
+
+ SquareMatrix A(3, Zero);
+ assignMatrix
+ (
+ A,
+ {
+ 1, 2, 3,
+ -4, -5, 6,
+ 7, -8, 9
+ }
+ );
+ const EigenMatrix EM(A);
+ const SquareMatrix complexEVecs(EM.complexEVecs());
+
+ SquareMatrix B(3, Zero);
+ assignMatrix
+ (
+ B,
+ {
+ complex(-0.373220280),
+ complex(0.417439996, 0.642691344),
+ complex(0.417439996, -0.642691344),
+ complex(-0.263919251),
+ complex(-1.165275867, 0.685068715),
+ complex(-1.165275867, -0.685068715),
+ complex(-0.889411744),
+ complex(-0.89990601, -0.3672785281),
+ complex(-0.89990601, 0.3672785281),
+ }
+ );
+
+ cmp
+ (
+ " # ",
+ flt(complexEVecs),
+ flt(B)
+ );
+ }
+ {
+ Info<< nl << " ## Test matrices with small values: ##" << nl;
+
+ const List epsilons
+ ({
+ 0, SMALL, Foam::sqrt(SMALL), sqr(SMALL), Foam::cbrt(SMALL),
+ -SMALL, -Foam::sqrt(SMALL), -sqr(SMALL), -Foam::cbrt(SMALL)
+ });
+
+ Random rndGen(1234);
+ const label numberOfTests = 20;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(100, 200);
+
+ for (const auto& eps : epsilons)
+ {
+ const SquareMatrix A(mRows, eps);
+
+ const EigenMatrix EM(A);
+ const DiagonalMatrix& EValsRe = EM.EValsRe();
+ const DiagonalMatrix& EValsIm = EM.EValsIm();
+ const SquareMatrix EVecs(EM.complexEVecs());
+
+ test_eigenvalues_sum(A, EValsRe);
+ test_characteristic_eq(A, EValsRe, EValsIm, EVecs);
+ }
+ }
+ }
+ {
+ Info<< nl << " ## Test matrices with repeating eigenvalues: ##" << nl;
+ SquareMatrix A(3, Zero);
+ assignMatrix
+ (
+ A,
+ {
+ 0, 1, 1,
+ 1, 0, 1,
+ 1, 1, 0
+ }
+ );
+ const EigenMatrix EM(A);
+ const DiagonalMatrix& EValsRe = EM.EValsRe();
+ const DiagonalMatrix& EValsIm = EM.EValsIm();
+ const SquareMatrix EVecs(EM.complexEVecs());
+
+ test_eigenvalues_sum(A, EValsRe);
+ test_characteristic_eq(A, EValsRe, EValsIm, EVecs);
+ }
+
+
+ if (nFail_)
+ {
+ Info<< nl << " #### "
+ << "Failed in " << nFail_ << " tests "
+ << "out of total " << nTest_ << " tests "
+ << "####\n" << endl;
+ return 1;
+ }
+
+ Info<< nl << " #### Passed all " << nTest_ <<" tests ####\n" << endl;
+ return 0;
+}
diff --git a/src/OpenFOAM/matrices/EigenMatrix/EigenMatrix.C b/src/OpenFOAM/matrices/EigenMatrix/EigenMatrix.C
new file mode 100644
index 0000000000..170a9777bd
--- /dev/null
+++ b/src/OpenFOAM/matrices/EigenMatrix/EigenMatrix.C
@@ -0,0 +1,1105 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Code created 2014-2018 by Alberto Passalacqua
+ Contributed 2018-07-31 to the OpenFOAM Foundation
+ Copyright (C) 2018 OpenFOAM Foundation
+ Copyright (C) 2019-2020 Alberto Passalacqua
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+\*---------------------------------------------------------------------------*/
+
+#include "EigenMatrix.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template
+void Foam::EigenMatrix::tridiagonaliseSymmMatrix()
+{
+ for (label j = 0; j < n_; ++j)
+ {
+ EValsRe_[j] = EVecs_[n_ - 1][j];
+ }
+
+ // Householder reduction to tridiagonal form
+ for (label i = n_ - 1; i > 0; --i)
+ {
+ // Scale to avoid under/overflow.
+ cmptType scale = Zero;
+ cmptType h = Zero;
+
+ for (label k = 0; k < i; ++k)
+ {
+ scale = scale + mag(EValsRe_[k]);
+ }
+
+ if (scale == 0.0)
+ {
+ EValsIm_[i] = EValsRe_[i - 1];
+
+ for (label j = 0; j < i; ++j)
+ {
+ EValsRe_[j] = EVecs_[i - 1][j];
+ EVecs_[i][j] = Zero;
+ EVecs_[j][i] = Zero;
+ }
+ }
+ else
+ {
+ // Generate Householder vector
+ for (label k = 0; k < i; ++k)
+ {
+ EValsRe_[k] /= scale;
+ h += EValsRe_[k]*EValsRe_[k];
+ }
+
+ cmptType f = EValsRe_[i - 1];
+ cmptType g = Foam::sqrt(h);
+
+ if (f > 0)
+ {
+ g = -g;
+ }
+
+ EValsIm_[i] = scale*g;
+ h -= f*g;
+ EValsRe_[i - 1] = f - g;
+
+ for (label j = 0; j < i; ++j)
+ {
+ EValsIm_[j] = Zero;
+ }
+
+ // Apply similarity transformation to remaining columns
+ for (label j = 0; j < i; ++j)
+ {
+ f = EValsRe_[j];
+ EVecs_[j][i] = f;
+ g = EValsIm_[j] + EVecs_[j][j]*f;
+
+ for (label k = j + 1; k <= i - 1; ++k)
+ {
+ g += EVecs_[k][j]*EValsRe_[k];
+ EValsIm_[k] += EVecs_[k][j]*f;
+ }
+
+ EValsIm_[j] = g;
+ }
+
+ f = Zero;
+
+ for (label j = 0; j < i; ++j)
+ {
+ EValsIm_[j] /= h;
+ f += EValsIm_[j]*EValsRe_[j];
+ }
+
+ cmptType hh = f/(2.0*h);
+
+ for (label j = 0; j < i; ++j)
+ {
+ EValsIm_[j] -= hh*EValsRe_[j];
+ }
+
+ for (label j = 0; j < i; ++j)
+ {
+ f = EValsRe_[j];
+ g = EValsIm_[j];
+
+ for (label k = j; k <= i - 1; ++k)
+ {
+ EVecs_[k][j] -=
+ (f*EValsIm_[k] + g*EValsRe_[k]);
+ }
+
+ EValsRe_[j] = EVecs_[i - 1][j];
+ EVecs_[i][j] = Zero;
+ }
+ }
+
+ EValsRe_[i] = h;
+ }
+
+ // Accumulate transformations
+ for (label i = 0; i < n_ - 1; ++i)
+ {
+ EVecs_[n_ - 1][i] = EVecs_[i][i];
+ EVecs_[i][i] = cmptType(1);
+ cmptType h = EValsRe_[i + 1];
+
+ if (h != 0.0)
+ {
+ for (label k = 0; k <= i; ++k)
+ {
+ EValsRe_[k] = EVecs_[k][i + 1]/h;
+ }
+
+ for (label j = 0; j <= i; ++j)
+ {
+ cmptType g = Zero;
+
+ for (label k = 0; k <= i; ++k)
+ {
+ g += EVecs_[k][i + 1]*EVecs_[k][j];
+ }
+
+ for (label k = 0; k <= i; ++k)
+ {
+ EVecs_[k][j] -= g*EValsRe_[k];
+ }
+ }
+ }
+
+ for (label k = 0; k <= i; ++k)
+ {
+ EVecs_[k][i + 1] = Zero;
+ }
+ }
+
+ for (label j = 0; j < n_; ++j)
+ {
+ EValsRe_[j] = EVecs_[n_ - 1][j];
+ EVecs_[n_ - 1][j] = Zero;
+ }
+
+ EVecs_[n_ - 1][n_ - 1] = cmptType(1);
+ EValsIm_[0] = Zero;
+}
+
+
+template
+void Foam::EigenMatrix::symmTridiagonalQL()
+{
+ for (label i = 1; i < n_; ++i)
+ {
+ EValsIm_[i - 1] = EValsIm_[i];
+ }
+
+ EValsIm_[n_-1] = Zero;
+
+ cmptType f = Zero;
+ cmptType tst1 = Zero;
+ cmptType eps = pow(2.0, -52.0);
+
+ for (label l = 0; l < n_; l++)
+ {
+ // Find SMALL subdiagonal element
+ tst1 = max(tst1, mag(EValsRe_[l]) + mag(EValsIm_[l]));
+ label m = l;
+
+ // Original while-loop from Java code
+ while (m < n_)
+ {
+ if (mag(EValsIm_[m]) <= eps*tst1)
+ {
+ break;
+ }
+
+ ++m;
+ }
+
+ // If m == l, EValsRe_[l] is an eigenvalue, otherwise, iterate
+ if (m > l)
+ {
+ label iter = 0;
+
+ do
+ {
+ iter += 1;
+
+ // Compute implicit shift
+ cmptType g = EValsRe_[l];
+ cmptType p = (EValsRe_[l + 1] - g)/(2.0*EValsIm_[l]);
+ cmptType r = Foam::hypot(p, cmptType(1));
+
+ if (p < 0)
+ {
+ r = -r;
+ }
+
+ EValsRe_[l] = EValsIm_[l]/(p + r);
+ EValsRe_[l + 1] = EValsIm_[l]*(p + r);
+ cmptType dl1 = EValsRe_[l + 1];
+ cmptType h = g - EValsRe_[l];
+
+ for (label i = l + 2; i < n_; ++i)
+ {
+ EValsRe_[i] -= h;
+ }
+
+ f += h;
+
+ // Implicit QL transformation.
+ p = EValsRe_[m];
+ cmptType c = cmptType(1);
+ cmptType c2 = c;
+ cmptType c3 = c;
+ cmptType el1 = EValsIm_[l + 1];
+ cmptType s = Zero;
+ cmptType s2 = Zero;
+
+ for (label i = m - 1; i >= l; --i)
+ {
+ c3 = c2;
+ c2 = c;
+ s2 = s;
+ g = c*EValsIm_[i];
+ h = c*p;
+ r = std::hypot(p, EValsIm_[i]);
+ EValsIm_[i + 1] = s*r;
+ s = EValsIm_[i]/r;
+ c = p/r;
+ p = c*EValsRe_[i] - s*g;
+ EValsRe_[i + 1] = h + s*(c*g + s*EValsRe_[i]);
+
+ // Accumulate transformation
+ for (label k = 0; k < n_; ++k)
+ {
+ h = EVecs_[k][i + 1];
+ EVecs_[k][i + 1] = s*EVecs_[k][i] + c*h;
+ EVecs_[k][i] = c*EVecs_[k][i] - s*h;
+ }
+ }
+
+ p = -s*s2*c3*el1*EValsIm_[l]/dl1;
+ EValsIm_[l] = s*p;
+ EValsRe_[l] = c*p;
+ }
+ while (mag(EValsIm_[l]) > eps*tst1); // Convergence check
+ }
+
+ EValsRe_[l] = EValsRe_[l] + f;
+ EValsIm_[l] = Zero;
+ }
+
+ // Sorting eigenvalues and corresponding vectors
+ for (label i = 0; i < n_ - 1; ++i)
+ {
+ label k = i;
+ cmptType p = EValsRe_[i];
+
+ for (label j = i + 1; j < n_; ++j)
+ {
+ if (EValsRe_[j] < p)
+ {
+ k = j;
+ p = EValsRe_[j];
+ }
+ }
+
+ if (k != i)
+ {
+ EValsRe_[k] = EValsRe_[i];
+ EValsRe_[i] = p;
+
+ for (label j = 0; j < n_; ++j)
+ {
+ p = EVecs_[j][i];
+ EVecs_[j][i] = EVecs_[j][k];
+ EVecs_[j][k] = p;
+ }
+ }
+ }
+}
+
+
+template
+void Foam::EigenMatrix::Hessenberg()
+{
+ DiagonalMatrix orth_(n_);
+
+ label low = 0;
+ label high = n_ - 1;
+
+ for (label m = low + 1; m <= high - 1; ++m)
+ {
+ // Scale column
+ cmptType scale = Zero;
+
+ for (label i = m; i <= high; ++i)
+ {
+ scale = scale + mag(H_[i][m - 1]);
+ }
+
+ if (scale != 0.0)
+ {
+ // Compute Householder transformation
+ cmptType h = Zero;
+
+ for (label i = high; i >= m; --i)
+ {
+ orth_[i] = H_[i][m - 1]/scale;
+ h += orth_[i]*orth_[i];
+ }
+
+ cmptType g = Foam::sqrt(h);
+
+ if (orth_[m] > 0)
+ {
+ g = -g;
+ }
+
+ h -= orth_[m]*g;
+ orth_[m] -= g;
+
+ // Apply Householder similarity transformation
+ // H = (I - u*u'/h)*H*(I - u*u')/h)
+ for (label j = m; j < n_; ++j)
+ {
+ cmptType f = Zero;
+
+ for (label i = high; i >= m; --i)
+ {
+ f += orth_[i]*H_[i][j];
+ }
+
+ f /= h;
+
+ for (label i = m; i <= high; ++i)
+ {
+ H_[i][j] -= f*orth_[i];
+ }
+ }
+
+ for (label i = 0; i <= high; ++i)
+ {
+ cmptType f = Zero;
+
+ for (label j = high; j >= m; --j)
+ {
+ f += orth_[j]*H_[i][j];
+ }
+
+ f /= h;
+
+ for (label j = m; j <= high; ++j)
+ {
+ H_[i][j] -= f*orth_[j];
+ }
+ }
+
+ orth_[m] = scale*orth_[m];
+ H_[m][m-1] = scale*g;
+ }
+ }
+
+ // Accumulate transformations
+ for (label i = 0; i < n_; ++i)
+ {
+ for (label j = 0; j < n_; ++j)
+ {
+ EVecs_[i][j] = (i == j ? 1.0 : 0.0);
+ }
+ }
+
+ for (label m = high - 1; m >= low + 1; --m)
+ {
+ if (H_[m][m - 1] != 0.0)
+ {
+ for (label i = m + 1; i <= high; ++i)
+ {
+ orth_[i] = H_[i][m - 1];
+ }
+
+ for (label j = m; j <= high; ++j)
+ {
+ cmptType g = Zero;
+
+ for (label i = m; i <= high; ++i)
+ {
+ g += orth_[i]*EVecs_[i][j];
+ }
+
+ // Double division avoids possible underflow
+ g = (g/orth_[m])/H_[m][m - 1];
+
+ for (label i = m; i <= high; ++i)
+ {
+ EVecs_[i][j] += g*orth_[i];
+ }
+ }
+ }
+ }
+}
+
+
+template
+void Foam::EigenMatrix::realSchur()
+{
+ // Initialize
+ label nn = n_;
+ label n = nn - 1;
+ label low = 0;
+ label high = nn - 1;
+ cmptType eps = pow(2.0, -52.0);
+ cmptType exshift = Zero;
+ cmptType p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
+
+ // Store roots isolated by balance and compute matrix norm
+ cmptType norm = Zero;
+
+ for (label i = 0; i < nn; ++i)
+ {
+ if ((i < low) || (i > high))
+ {
+ EValsRe_[i] = H_[i][i];
+ EValsIm_[i] = Zero;
+ }
+
+ for (label j = max(i - 1, 0); j < nn; ++j)
+ {
+ norm += mag(H_[i][j]);
+ }
+ }
+
+ label iter = 0;
+
+ while (n >= low)
+ {
+ // Look for single SMALL sub-diagonal element
+ label l = n;
+
+ while (l > low)
+ {
+ s = mag(H_[l - 1][l - 1]) + mag(H_[l][l]);
+
+ if (s == 0.0)
+ {
+ s = norm;
+ }
+
+ if (mag(H_[l][l - 1]) < eps*s)
+ {
+ break;
+ }
+
+ l--;
+ }
+
+ // Check for convergence
+ if (l == n) // One root found
+ {
+ H_[n][n] = H_[n][n] + exshift;
+ EValsRe_[n] = H_[n][n];
+ EValsIm_[n] = Zero;
+ --n;
+ iter = 0;
+ }
+ else if (l == n - 1) // Two roots found
+ {
+ w = H_[n][n - 1]*H_[n - 1][n];
+ p = (H_[n - 1][n - 1] - H_[n][n])/2.0;
+ q = p*p + w;
+ z = Foam::sqrt(mag(q));
+ H_[n][n] += exshift;
+ H_[n - 1][n - 1] += exshift;
+ x = H_[n][n];
+
+ // Scalar pair
+ if (q >= 0)
+ {
+ if (p >= 0)
+ {
+ z = p + z;
+ }
+ else
+ {
+ z = p - z;
+ }
+
+ EValsRe_[n - 1] = x + z;
+ EValsRe_[n] = EValsRe_[n - 1];
+
+ if (z != 0.0)
+ {
+ EValsRe_[n] = x - w/z;
+ }
+
+ EValsIm_[n - 1] = Zero;
+ EValsIm_[n] = Zero;
+ x = H_[n][n - 1];
+ s = mag(x) + mag(z);
+ p = x/s;
+ q = z/s;
+ r = Foam::sqrt(p*p + q*q);
+ p /= r;
+ q /= r;
+
+ // Row modification
+ for (label j = n - 1; j < nn; ++j)
+ {
+ z = H_[n - 1][j];
+ H_[n - 1][j] = q*z + p*H_[n][j];
+ H_[n][j] = q*H_[n][j] - p*z;
+ }
+
+ // Column modification
+ for (label i = 0; i <= n; ++i)
+ {
+ z = H_[i][n - 1];
+ H_[i][n - 1] = q*z + p*H_[i][n];
+ H_[i][n] = q*H_[i][n] - p*z;
+ }
+
+ // Accumulate transformations
+ for (label i = low; i <= high; ++i)
+ {
+ z = EVecs_[i][n - 1];
+ EVecs_[i][n - 1] = q*z + p*EVecs_[i][n];
+ EVecs_[i][n] = q*EVecs_[i][n] - p*z;
+ }
+ }
+ else // Complex pair of roots
+ {
+ EValsRe_[n - 1] = x + p;
+ EValsRe_[n] = x + p;
+ EValsIm_[n - 1] = z;
+ EValsIm_[n] = -z;
+ }
+
+ n -= 2;
+ iter = 0;
+ }
+ else // No convergence yet
+ {
+ // Form shift
+ x = H_[n][n];
+ y = Zero;
+ w = Zero;
+
+ if (l < n)
+ {
+ y = H_[n - 1][n - 1];
+ w = H_[n][n - 1]*H_[n - 1][n];
+ }
+
+ // Wilkinson's original ad-hoc shift
+ if (iter == 10)
+ {
+ exshift += x;
+ for (label i = low; i <= n; ++i)
+ {
+ H_[i][i] -= x;
+ }
+
+ s = mag(H_[n][n - 1]) + mag(H_[n - 1][n - 2]);
+ x = 0.75*s;
+ y = 0.75*s;
+ w = -0.4375*s*s;
+ }
+
+ // New ad hoc shift
+ if (iter == 30)
+ {
+ s = (y - x)/2.0;
+ s = s*s + w;
+
+ if (s > 0)
+ {
+ s = Foam::sqrt(s);
+
+ if (y < x)
+ {
+ s = -s;
+ }
+
+ s = x - w/((y - x)/2.0 + s);
+
+ for (label i = low; i <= n; ++i)
+ {
+ H_[i][i] -= s;
+ }
+
+ exshift += s;
+ x = 0.964;
+ y = 0.964;
+ w = 0.964;
+ }
+ }
+
+ iter += 1;
+
+ // Look for two consecutive SMALL sub-diagonal elements
+ label m = n - 2;
+
+ while (m >= l)
+ {
+ z = H_[m][m];
+ r = x - z;
+ s = y - z;
+ p = (r*s - w)/H_[m + 1][m] + H_[m][m + 1];
+ q = H_[m + 1][m + 1] - z - r - s;
+ r = H_[m + 2][m + 1];
+ s = mag(p) + mag(q) + mag(r);
+ p /= s;
+ q /= s;
+ r /= s;
+
+ if (m == l)
+ {
+ break;
+ }
+
+ if
+ (
+ mag(H_[m][m - 1])*(mag(q) + mag(r))
+ < eps*(mag(p)*(mag(H_[m - 1][m - 1])
+ + mag(z) + mag(H_[m + 1][m + 1])))
+ )
+ {
+ break;
+ }
+
+ --m;
+ }
+
+ for (label i = m + 2; i <= n; ++i)
+ {
+ H_[i][i - 2] = Zero;
+
+ if (i > m + 2)
+ {
+ H_[i][i - 3] = Zero;
+ }
+ }
+
+ // Double QR step involving rows l:n and columns m:n
+ for (label k = m; k <= n - 1; ++k)
+ {
+ label notlast = (k != n - 1);
+
+ if (k != m)
+ {
+ p = H_[k][k - 1];
+ q = H_[k + 1][k - 1];
+ r = (notlast ? H_[k + 2][k - 1] : 0.0);
+ x = mag(p) + mag(q) + mag(r);
+
+ if (x != 0.0)
+ {
+ p /= x;
+ q /= x;
+ r /= x;
+ }
+ }
+
+ if (x == 0.0)
+ {
+ break;
+ }
+
+ s = Foam::sqrt(p*p + q*q + r*r);
+
+ if (p < 0)
+ {
+ s = -s;
+ }
+
+ if (s != 0)
+ {
+ if (k != m)
+ {
+ H_[k][k - 1] = -s*x;
+ }
+ else if (l != m)
+ {
+ H_[k][k - 1] = -H_[k][k - 1];
+ }
+
+ p = p + s;
+ x = p/s;
+ y = q/s;
+ z = r/s;
+ q /= p;
+ r /= p;
+
+ // Row modification
+ for (label j = k; j < nn; ++j)
+ {
+ p = H_[k][j] + q*H_[k + 1][j];
+
+ if (notlast)
+ {
+ p += r*H_[k + 2][j];
+ H_[k + 2][j] -= p*z;
+ }
+
+ H_[k][j] -= p*x;
+ H_[k + 1][j] -= p*y;
+ }
+
+ // Column modification
+ for (label i = 0; i <= min(n, k + 3); ++i)
+ {
+ p = x*H_[i][k] + y*H_[i][k + 1];
+
+ if (notlast)
+ {
+ p += z*H_[i][k + 2];
+ H_[i][k + 2] -= p*r;
+ }
+
+ H_[i][k] -= p;
+ H_[i][k + 1] -= p*q;
+ }
+
+ // Accumulate transformations
+ for (label i = low; i <= high; ++i)
+ {
+ p = x*EVecs_[i][k] + y*EVecs_[i][k + 1];
+
+ if (notlast)
+ {
+ p += z*EVecs_[i][k + 2];
+ EVecs_[i][k + 2] -= p*r;
+ }
+
+ EVecs_[i][k] -= p;
+ EVecs_[i][k + 1] -= p*q;
+ }
+ }
+ }
+ }
+ }
+
+ // Backsubstitute to find vectors of upper triangular form
+ if (norm == 0.0)
+ {
+ return;
+ }
+
+ for (n = nn-1; n >= 0; --n)
+ {
+ p = EValsRe_[n];
+ q = EValsIm_[n];
+
+ // scalar vector
+ if (q == 0)
+ {
+ label l = n;
+ H_[n][n] = cmptType(1);
+
+ for (label i = n-1; i >= 0; --i)
+ {
+ w = H_[i][i] - p;
+ r = Zero;
+
+ for (label j = l; j <= n; ++j)
+ {
+ r += H_[i][j]*H_[j][n];
+ }
+
+ if (EValsIm_[i] < 0.0)
+ {
+ z = w;
+ s = r;
+ }
+ else
+ {
+ l = i;
+
+ if (EValsIm_[i] == 0.0)
+ {
+ if (w != 0.0)
+ {
+ H_[i][n] = -r/w;
+ }
+ else
+ {
+ H_[i][n] = -r/(eps*norm);
+ }
+ }
+ else // Solve real equations
+ {
+ x = H_[i][i + 1];
+ y = H_[i + 1][i];
+ q = (EValsRe_[i] - p)*(EValsRe_[i] - p)
+ + EValsIm_[i]*EValsIm_[i];
+
+ t = (x*s - z*r)/q;
+ H_[i][n] = t;
+
+ if (mag(x) > mag(z))
+ {
+ H_[i + 1][n] = (-r - w*t)/x;
+ }
+ else
+ {
+ H_[i + 1][n] = (-s - y*t)/z;
+ }
+ }
+
+ // Overflow control
+ t = mag(H_[i][n]);
+
+ if ((eps*t)*t > 1)
+ {
+ for (label j = i; j <= n; ++j)
+ {
+ H_[j][n] /= t;
+ }
+ }
+ }
+ }
+ }
+ else if (q < 0) // Complex vector
+ {
+ label l = n - 1;
+
+ // Last vector component imaginary so matrix is triangular
+ if (mag(H_[n][n - 1]) > mag(H_[n - 1][n]))
+ {
+ H_[n - 1][n - 1] = q/H_[n][n - 1];
+ H_[n - 1][n] = -(H_[n][n] - p)/H_[n][n - 1];
+ }
+ else
+ {
+ complex cDiv = complex(0.0, -H_[n - 1][n])
+ /complex(H_[n - 1][n-1]-p, q);
+
+ H_[n - 1][n - 1] = cDiv.Re();
+ H_[n - 1][n] = cDiv.Im();
+ }
+
+ H_[n][n - 1] = Zero;
+ H_[n][n] = cmptType(1);
+
+ for (label i = n - 2; i >= 0; --i)
+ {
+ cmptType ra, sa, vr, vi;
+ ra = Zero;
+ sa = Zero;
+
+ for (label j = l; j <= n; ++j)
+ {
+ ra += H_[i][j]*H_[j][n-1];
+ sa += H_[i][j]*H_[j][n];
+ }
+
+ w = H_[i][i] - p;
+
+ if (EValsIm_[i] < 0.0)
+ {
+ z = w;
+ r = ra;
+ s = sa;
+ }
+ else
+ {
+ l = i;
+
+ if (EValsIm_[i] == 0)
+ {
+ complex cDiv = complex(-ra, -sa)/complex(w, q);
+ H_[i][n - 1] = cDiv.Re();
+ H_[i][n] = cDiv.Im();
+ }
+ else
+ {
+ // Solve complex equations
+ x = H_[i][i + 1];
+ y = H_[i + 1][i];
+ vr = (EValsRe_[i] - p)*(EValsRe_[i] - p)
+ + EValsIm_[i]*EValsIm_[i] - q*q;
+
+ vi = 2.0*(EValsRe_[i] - p)*q;
+
+ if ((vr == 0.0) && (vi == 0.0))
+ {
+ vr = eps*norm*(mag(w) + mag(q) + mag(x) + mag(y)
+ + mag(z));
+ }
+
+ complex cDiv =
+ complex(x*r - z*ra + q*sa, x*s - z*sa - q*ra)
+ /complex(vr, vi);
+
+ H_[i][n - 1] = cDiv.Re();
+ H_[i][n] = cDiv.Im();
+
+ if (mag(x) > (mag(z) + mag(q)))
+ {
+ H_[i + 1][n - 1] = (-ra - w*H_[i][n - 1]
+ + q*H_[i][n])/x;
+
+ H_[i + 1][n] = (-sa - w*H_[i][n]
+ - q*H_[i][n - 1])/x;
+ }
+ else
+ {
+ complex cDiv = complex(-r - y*H_[i][n - 1], -s
+ - y*H_[i][n])/complex(z, q);
+
+ H_[i + 1][n - 1] = cDiv.Re();
+ H_[i + 1][n] = cDiv.Im();
+ }
+ }
+
+ // Overflow control
+ t = max(mag(H_[i][n - 1]), mag(H_[i][n]));
+
+ if ((eps*t)*t > 1)
+ {
+ for (label j = i; j <= n; ++j)
+ {
+ H_[j][n - 1] /= t;
+ H_[j][n] /= t;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Vectors of isolated roots
+ for (label i = 0; i < nn; ++i)
+ {
+ if (i < low || i > high)
+ {
+ for (label j = i; j < nn; ++j)
+ {
+ EVecs_[i][j] = H_[i][j];
+ }
+ }
+ }
+
+ // Back transformation to get eigenvectors of original matrix
+ for (label j = nn - 1; j >= low; --j)
+ {
+ for (label i = low; i <= high; ++i)
+ {
+ z = Zero;
+
+ for (label k = low; k <= min(j, high); ++k)
+ {
+ z = z + EVecs_[i][k]*H_[k][j];
+ }
+
+ EVecs_[i][j] = z;
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template
+Foam::EigenMatrix::EigenMatrix(const SquareMatrix& A)
+:
+ n_(A.n()),
+ EValsRe_(n_, Zero),
+ EValsIm_(n_, Zero),
+ EVecs_(n_, Zero),
+ H_()
+{
+ if (n_ <= 0)
+ {
+ FatalErrorInFunction
+ << "Input matrix has zero size."
+ << abort(FatalError);
+ }
+
+ if (A.symmetric())
+ {
+ EVecs_ = A;
+ tridiagonaliseSymmMatrix();
+ symmTridiagonalQL();
+ }
+ else
+ {
+ H_ = A;
+ Hessenberg();
+ realSchur();
+ }
+}
+
+
+template
+Foam::EigenMatrix::EigenMatrix
+(
+ const SquareMatrix& A,
+ bool symmetric
+)
+:
+ n_(A.n()),
+ EValsRe_(n_, Zero),
+ EValsIm_(n_, Zero),
+ EVecs_(n_, Zero),
+ H_()
+{
+ if (n_ <= 0)
+ {
+ FatalErrorInFunction
+ << "Input matrix has zero size."
+ << abort(FatalError);
+ }
+
+ if (symmetric)
+ {
+ EVecs_ = A;
+ tridiagonaliseSymmMatrix();
+ symmTridiagonalQL();
+ }
+ else
+ {
+ H_ = A;
+ Hessenberg();
+ realSchur();
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template
+const Foam::SquareMatrix
+Foam::EigenMatrix::complexEVecs() const
+{
+ SquareMatrix EVecs(EVecs_.n());
+
+ std::transform
+ (
+ EVecs_.cbegin(),
+ EVecs_.cend(),
+ EVecs.begin(),
+ [&](const cmptType& x){ return complex(x); }
+ );
+
+ for (label i = 0; i < EValsIm_.size(); ++i)
+ {
+ if (mag(EValsIm_[i]) > VSMALL)
+ {
+ for (label j = 0; j < EVecs.m(); ++j)
+ {
+ EVecs(j, i) = complex(EVecs_(j, i), EVecs_(j, i+1));
+ EVecs(j, i+1) = EVecs(j, i).conjugate();
+ }
+ ++i;
+ }
+ }
+
+ return EVecs;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/EigenMatrix/EigenMatrix.H b/src/OpenFOAM/matrices/EigenMatrix/EigenMatrix.H
new file mode 100644
index 0000000000..9a40843d4e
--- /dev/null
+++ b/src/OpenFOAM/matrices/EigenMatrix/EigenMatrix.H
@@ -0,0 +1,267 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Code created 2014-2018 by Alberto Passalacqua
+ Contributed 2018-07-31 to the OpenFOAM Foundation
+ Copyright (C) 2018 OpenFOAM Foundation
+ Copyright (C) 2019-2020 Alberto Passalacqua
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Class
+ Foam::EigenMatrix
+
+Description
+ EigenMatrix (i.e. eigendecomposition or spectral decomposition) decomposes
+ a diagonalisable nonsymmetric real square matrix into its canonical form,
+ whereby the matrix is represented in terms of its eigenvalues and
+ eigenvectors.
+
+ The eigenvalue equation (i.e. eigenvalue problem) is:
+
+ \f[
+ A v = \lambda v
+ \f]
+
+ where
+ \vartable
+ A | a diagonalisable square matrix of dimension m-by-m
+ v | a (non-zero) vector of dimension m (right eigenvector)
+ \lambda | a scalar corresponding to v (eigenvalue)
+ \endvartable
+
+
+ If \c A is symmetric, the following relation is satisfied:
+
+ \f[
+ A = v*D*v^T
+ \f]
+
+ where
+ \vartable
+ D | diagonal real eigenvalue matrix
+ v | orthogonal eigenvector matrix
+ \endvartable
+
+ If \c A is not symmetric, \c D becomes a block diagonal matrix wherein
+ the real eigenvalues are present on the diagonal within 1-by-1 blocks, and
+ complex eigenvalues within 2-by-2 blocks, i.e. \f$\lambda + i \mu\f$ with
+ \f$[\lambda, \mu; -\mu, \lambda]\f$.
+
+ The columns of \c v represent eigenvectors corresponding to eigenvalues,
+ satisfying the eigenvalue equation. Even though eigenvalues of a matrix
+ are unique, eigenvectors of the matrix are not. For the same eigenvalue,
+ the corresponding eigenvector can be real or complex with non-unique
+ entries. In addition, the validity of the equation \f$A = v*D*v^T\f$
+ depends on the condition number of \c v, which can be ill-conditioned,
+ or singular for invalidated equations.
+
+ References:
+ \verbatim
+ OpenFOAM-compatible implementation:
+ Passalacqua, A., Heylmun, J., Icardi, M.,
+ Madadi, E., Bachant, P., & Hu, X. (2019).
+ OpenQBMM 5.0.1 for OpenFOAM 7, Zenodo.
+ DOI:10.5281/zenodo.3471804
+
+ Implementations for the functions:
+ 'tridiagonaliseSymmMatrix', 'symmTridiagonalQL',
+ 'Hessenberg' and 'realSchur' (based on ALGOL-procedure:tred2):
+ Wilkinson, J. H., & Reinsch, C. (1971).
+ In Bauer, F. L. & Householder A. S. (Eds.),
+ Handbook for Automatic Computation: Volume II: Linear Algebra.
+ (Vol. 186), Springer-Verlag Berlin Heidelberg.
+ DOI: 10.1007/978-3-642-86940-2
+
+ Explanations on how real eigenvectors
+ can be unpacked into complex domain:
+ Moler, C. (1998).
+ Re: Eigenvectors.
+ Retrieved from https://bit.ly/3ao4Wat
+
+ TNT/JAMA implementation:
+ Pozo, R. (1997).
+ Template Numerical Toolkit for linear algebra:
+ High performance programming with C++
+ and the Standard Template Library.
+ The International Journal of Supercomputer Applications
+ and High Performance Computing, 11(3), 251-263.
+ DOI:10.1177/109434209701100307
+
+ (No particular order) Hicklin, J., Moler, C., Webb, P.,
+ Boisvert, R. F., Miller, B., Pozo, R., & Remington, K. (2012).
+ JAMA: A Java Matrix Package 1.0.3.
+ Retrived from https://math.nist.gov/javanumerics/jama/
+ \endverbatim
+
+Note
+ - This implementation is an integration of the \c OpenQBMM \c eigenSolver
+ class (2019) without any changes to its internal mechanisms. Therefore, no
+ differences between EigenMatrix and \c eigenSolver (2019) classes should be
+ expected in terms of input-process-output operations.
+ - The \c OpenQBMM \c eigenSolver class derives almost completely from the
+ \c TNT/JAMA implementation, a public-domain library developed by \c NIST
+ and \c MathWorks from 1998 to 2012, available at
+ http://math.nist.gov/tnt/index.html (Retrieved June 6, 2020). Their
+ implementation was based upon \c EISPACK.
+ - The \c tridiagonaliseSymmMatrix, \c symmTridiagonalQL, \c Hessenberg and
+ \c realSchur methods are based on the \c Algol procedures \c tred2 by
+ Bowdler, Martin, Reinsch, and Wilkinson, Handbook for Auto. Comp.,
+ Vol. II-Linear Algebra, and the corresponding \c FORTRAN subroutine
+ in \c EISPACK.
+
+See also
+ Test-EigenMatrix.C
+
+SourceFiles
+ EigenMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef EigenMatrix_H
+#define EigenMatrix_H
+
+#include "scalarMatrices.H"
+#include "DiagonalMatrix.H"
+#include "SquareMatrix.H"
+#include "complex.H"
+#include
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class EigenMatrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template
+class EigenMatrix
+{
+ static_assert
+ (
+ std::is_floating_point::value,
+ "EigenMatrix operates only with scalar base type."
+ );
+
+ // Private Data
+
+ //- Number of rows and columns in input matrix
+ const label n_;
+
+ //- Real eigenvalues or real part of complex eigenvalues
+ DiagonalMatrix EValsRe_;
+
+ //- Zero-matrix for real eigenvalues
+ //- or imaginary part of complex eigenvalues
+ DiagonalMatrix EValsIm_;
+
+ //- Right eigenvectors matrix where each column is
+ //- a right eigenvector that corresponds to an eigenvalue
+ SquareMatrix EVecs_;
+
+ //- Copy of nonsymmetric input matrix evolving to eigenvectors matrix
+ SquareMatrix H_;
+
+
+ // Private Member Functions
+
+ //- Householder transform of a symmetric matrix to tri-diagonal form
+ void tridiagonaliseSymmMatrix();
+
+ //- Symmetric tri-diagonal QL algorithm
+ void symmTridiagonalQL();
+
+ //- Reduce non-symmetric matrix to upper-Hessenberg form
+ void Hessenberg();
+
+ //- Reduce matrix from Hessenberg to real Schur form
+ void realSchur();
+
+
+public:
+
+ // Generated Methods
+
+ //- No default construct
+ EigenMatrix() = delete;
+
+ //- No copy construct
+ EigenMatrix(const EigenMatrix&) = delete;
+
+ //- No copy assignment
+ EigenMatrix& operator=(const EigenMatrix&) = delete;
+
+
+ // Constructors
+
+ //- Construct from a SquareMatrix
+ explicit EigenMatrix(const SquareMatrix& A);
+
+ //- Construct from a SquareMatrix and symmetry flag
+ // Does not perform symmetric check
+ EigenMatrix(const SquareMatrix& A, bool symmetric);
+
+
+ // Access
+
+ //- Return real eigenvalues or real part of complex eigenvalues
+ const DiagonalMatrix& EValsRe() const
+ {
+ return EValsRe_;
+ }
+
+ //- Return zero-matrix for real eigenvalues
+ //- or imaginary part of complex eigenvalues
+ const DiagonalMatrix& EValsIm() const
+ {
+ return EValsIm_;
+ }
+
+ //- Return right eigenvectors matrix where each column is
+ //- a right eigenvector that corresponds to an eigenvalue
+ const SquareMatrix& EVecs() const
+ {
+ return EVecs_;
+ }
+
+ //- Return right eigenvectors in unpacked form
+ const SquareMatrix complexEVecs() const;
+
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "EigenMatrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //