EigenMatrix.H
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Code created 2014-2018 by Alberto Passalacqua
9 Contributed 2018-07-31 to the OpenFOAM Foundation
10 Copyright (C) 2018 OpenFOAM Foundation
11 Copyright (C) 2019-2020 Alberto Passalacqua
12 Copyright (C) 2020 OpenCFD Ltd.
13-------------------------------------------------------------------------------
14License
15 This file is part of OpenFOAM.
16
17 OpenFOAM is free software: you can redistribute it and/or modify it
18 under the terms of the GNU General Public License as published by
19 the Free Software Foundation, either version 3 of the License, or
20 (at your option) any later version.
21
22 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
23 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
24 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
25 for more details.
26
27 You should have received a copy of the GNU General Public License
28 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
29
30Class
31 Foam::EigenMatrix
32
33Description
34 EigenMatrix (i.e. eigendecomposition or spectral decomposition) decomposes
35 a diagonalisable nonsymmetric real square matrix into its canonical form,
36 whereby the matrix is represented in terms of its eigenvalues and
37 eigenvectors.
38
39 The eigenvalue equation (i.e. eigenvalue problem) is:
40
41 \f[
42 A v = \lambda v
43 \f]
44
45 where
46 \vartable
47 A | a diagonalisable square matrix of dimension m-by-m
48 v | a (non-zero) vector of dimension m (right eigenvector)
49 \lambda | a scalar corresponding to v (eigenvalue)
50 \endvartable
51
52
53 If \c A is symmetric, the following relation is satisfied:
54
55 \f[
56 A = v*D*v^T
57 \f]
58
59 where
60 \vartable
61 D | diagonal real eigenvalue matrix
62 v | orthogonal eigenvector matrix
63 \endvartable
64
65 If \c A is not symmetric, \c D becomes a block diagonal matrix wherein
66 the real eigenvalues are present on the diagonal within 1-by-1 blocks, and
67 complex eigenvalues within 2-by-2 blocks, i.e. \f$\lambda + i \mu\f$ with
68 \f$[\lambda, \mu; -\mu, \lambda]\f$.
69
70 The columns of \c v represent eigenvectors corresponding to eigenvalues,
71 satisfying the eigenvalue equation. Even though eigenvalues of a matrix
72 are unique, eigenvectors of the matrix are not. For the same eigenvalue,
73 the corresponding eigenvector can be real or complex with non-unique
74 entries. In addition, the validity of the equation \f$A = v*D*v^T\f$
75 depends on the condition number of \c v, which can be ill-conditioned,
76 or singular for invalidated equations.
77
78 References:
79 \verbatim
80 OpenFOAM-compatible implementation:
81 Passalacqua, A., Heylmun, J., Icardi, M.,
82 Madadi, E., Bachant, P., & Hu, X. (2019).
83 OpenQBMM 5.0.1 for OpenFOAM 7, Zenodo.
84 DOI:10.5281/zenodo.3471804
85
86 Implementations for the functions:
87 'tridiagonaliseSymmMatrix', 'symmTridiagonalQL',
88 'Hessenberg' and 'realSchur' (based on ALGOL-procedure:tred2):
89 Wilkinson, J. H., & Reinsch, C. (1971).
90 In Bauer, F. L. & Householder A. S. (Eds.),
91 Handbook for Automatic Computation: Volume II: Linear Algebra.
92 (Vol. 186), Springer-Verlag Berlin Heidelberg.
93 DOI: 10.1007/978-3-642-86940-2
94
95 Explanations on how real eigenvectors
96 can be unpacked into complex domain:
97 Moler, C. (1998).
98 Re: Eigenvectors.
99 Retrieved from https://bit.ly/3ao4Wat
100
101 TNT/JAMA implementation:
102 Pozo, R. (1997).
103 Template Numerical Toolkit for linear algebra:
104 High performance programming with C++
105 and the Standard Template Library.
106 The International Journal of Supercomputer Applications
107 and High Performance Computing, 11(3), 251-263.
108 DOI:10.1177/109434209701100307
109
110 (No particular order) Hicklin, J., Moler, C., Webb, P.,
111 Boisvert, R. F., Miller, B., Pozo, R., & Remington, K. (2012).
112 JAMA: A Java Matrix Package 1.0.3.
113 Retrived from https://math.nist.gov/javanumerics/jama/
114 \endverbatim
115
116Note
117 - This implementation is an integration of the \c OpenQBMM \c eigenSolver
118 class (2019) without any changes to its internal mechanisms. Therefore, no
119 differences between EigenMatrix and \c eigenSolver (2019) classes should be
120 expected in terms of input-process-output operations.
121 - The \c OpenQBMM \c eigenSolver class derives almost completely from the
122 \c TNT/JAMA implementation, a public-domain library developed by \c NIST
123 and \c MathWorks from 1998 to 2012, available at
124 http://math.nist.gov/tnt/index.html (Retrieved June 6, 2020). Their
125 implementation was based upon \c EISPACK.
126 - The \c tridiagonaliseSymmMatrix, \c symmTridiagonalQL, \c Hessenberg and
127 \c realSchur methods are based on the \c Algol procedures \c tred2 by
128 Bowdler, Martin, Reinsch, and Wilkinson, Handbook for Auto. Comp.,
129 Vol. II-Linear Algebra, and the corresponding \c FORTRAN subroutine
130 in \c EISPACK.
131
132See also
133 Test-EigenMatrix.C
134
135SourceFiles
136 EigenMatrix.C
137
138\*---------------------------------------------------------------------------*/
139
140#ifndef EigenMatrix_H
141#define EigenMatrix_H
142
143#include "scalarMatrices.H"
144#include "DiagonalMatrix.H"
145#include "SquareMatrix.H"
146#include "complex.H"
147#include <algorithm>
148
149// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
150
151namespace Foam
152{
153
154/*---------------------------------------------------------------------------*\
155 Class EigenMatrix Declaration
156\*---------------------------------------------------------------------------*/
157
158template<class cmptType>
159class EigenMatrix
160{
161 static_assert
162 (
163 std::is_floating_point<cmptType>::value,
164 "EigenMatrix operates only with scalar base type."
165 );
166
167 // Private Data
168
169 //- Number of rows and columns in input matrix
170 const label n_;
171
172 //- Real eigenvalues or real part of complex eigenvalues
173 DiagonalMatrix<cmptType> EValsRe_;
174
175 //- Zero-matrix for real eigenvalues
176 //- or imaginary part of complex eigenvalues
177 DiagonalMatrix<cmptType> EValsIm_;
179 //- Right eigenvectors matrix where each column is
180 //- a right eigenvector that corresponds to an eigenvalue
182
183 //- Copy of nonsymmetric input matrix evolving to eigenvectors matrix
185
186
187 // Private Member Functions
188
189 //- Householder transform of a symmetric matrix to tri-diagonal form
190 void tridiagonaliseSymmMatrix();
191
192 //- Symmetric tri-diagonal QL algorithm
193 void symmTridiagonalQL();
194
195 //- Reduce non-symmetric matrix to upper-Hessenberg form
196 void Hessenberg();
197
198 //- Reduce matrix from Hessenberg to real Schur form
199 void realSchur();
200
201
202public:
203
204 // Generated Methods
205
206 //- No default construct
207 EigenMatrix() = delete;
208
209 //- No copy construct
210 EigenMatrix(const EigenMatrix&) = delete;
211
212 //- No copy assignment
213 EigenMatrix& operator=(const EigenMatrix&) = delete;
214
215
216 // Constructors
217
218 //- Construct from a SquareMatrix<cmptType>
219 explicit EigenMatrix(const SquareMatrix<cmptType>& A);
220
221 //- Construct from a SquareMatrix<cmptType> and symmetry flag
222 // Does not perform symmetric check
223 EigenMatrix(const SquareMatrix<cmptType>& A, bool symmetric);
224
225
226 // Access
227
228 //- Return real eigenvalues or real part of complex eigenvalues
230 {
231 return EValsRe_;
233
234 //- Return zero-matrix for real eigenvalues
235 //- or imaginary part of complex eigenvalues
236 const DiagonalMatrix<cmptType>& EValsIm() const
237 {
238 return EValsIm_;
239 }
240
241 //- Return right eigenvectors matrix where each column is
242 //- a right eigenvector that corresponds to an eigenvalue
243 const SquareMatrix<cmptType>& EVecs() const
244 {
245 return EVecs_;
246 }
247
248 //- Return right eigenvectors in unpacked form
250
251};
252
253// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
254
255} // End namespace Foam
256
257// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
258
259#ifdef NoRepository
260 #include "EigenMatrix.C"
261#endif
263// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
264
265#endif
266
267// ************************************************************************* //
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
A templated (N x N) diagonal matrix of objects of <Type>, effectively containing N elements,...
EigenMatrix (i.e. eigendecomposition or spectral decomposition) decomposes a diagonalisable nonsymmet...
Definition: EigenMatrix.H:179
const SquareMatrix< cmptType > & EVecs() const
Definition: EigenMatrix.H:262
EigenMatrix(const EigenMatrix &)=delete
No copy construct.
const DiagonalMatrix< cmptType > & EValsIm() const
Definition: EigenMatrix.H:255
const DiagonalMatrix< cmptType > & EValsRe() const
Return real eigenvalues or real part of complex eigenvalues.
Definition: EigenMatrix.H:248
EigenMatrix()=delete
No default construct.
EigenMatrix & operator=(const EigenMatrix &)=delete
No copy assignment.
const SquareMatrix< complex > complexEVecs() const
Return right eigenvectors in unpacked form.
Definition: EigenMatrix.C:1076
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
Definition: SquareMatrix.H:66
Namespace for OpenFOAM.