OpenFOAM: API Guide
v2112
The open source CFD toolbox
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
-------------------------------------------------------------------------------
14
License
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
30
Class
31
Foam::EigenMatrix
32
33
Description
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
116
Note
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
132
See also
133
Test-EigenMatrix.C
134
135
SourceFiles
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
151
namespace
Foam
152
{
153
154
/*---------------------------------------------------------------------------*\
155
Class EigenMatrix Declaration
156
\*---------------------------------------------------------------------------*/
157
158
template
<
class
cmptType>
159
class
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_;
178
179
//- Right eigenvectors matrix where each column is
180
//- a right eigenvector that corresponds to an eigenvalue
181
SquareMatrix<cmptType>
EVecs_;
182
183
//- Copy of nonsymmetric input matrix evolving to eigenvectors matrix
184
SquareMatrix<cmptType>
H_;
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
202
public
:
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
229
const
DiagonalMatrix<cmptType>
&
EValsRe
()
const
230
{
231
return
EValsRe_;
232
}
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
249
const
SquareMatrix<complex>
complexEVecs
()
const
;
250
251
};
252
253
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
254
255
}
// End namespace Foam
256
257
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
258
259
#ifdef NoRepository
260
#include "
EigenMatrix.C
"
261
#endif
262
263
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
264
265
#endif
266
267
// ************************************************************************* //
Foam::EigenMatrix::EValsRe
const DiagonalMatrix< cmptType > & EValsRe() const
Return real eigenvalues or real part of complex eigenvalues.
Definition:
EigenMatrix.H:248
Foam::EigenMatrix::operator=
EigenMatrix & operator=(const EigenMatrix &)=delete
No copy assignment.
Foam::EigenMatrix::EigenMatrix
EigenMatrix()=delete
No default construct.
DiagonalMatrix.H
A
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
complex.H
Foam::DiagonalMatrix< cmptType >
EigenMatrix.C
Foam::EigenMatrix::EValsIm
const DiagonalMatrix< cmptType > & EValsIm() const
Definition:
EigenMatrix.H:255
Foam::EigenMatrix::complexEVecs
const SquareMatrix< complex > complexEVecs() const
Return right eigenvectors in unpacked form.
Definition:
EigenMatrix.C:1076
Foam
Namespace for OpenFOAM.
Definition:
atmBoundaryLayer.C:33
Foam::EigenMatrix::EVecs
const SquareMatrix< cmptType > & EVecs() const
Definition:
EigenMatrix.H:262
Foam::SquareMatrix< cmptType >
Foam::EigenMatrix
EigenMatrix (i.e. eigendecomposition or spectral decomposition) decomposes a diagonalisable nonsymmet...
Definition:
EigenMatrix.H:178
scalarMatrices.H
SquareMatrix.H
src
OpenFOAM
matrices
EigenMatrix
EigenMatrix.H
Generated by
1.8.17
OPENFOAM® is a registered
trademark
of OpenCFD Ltd.