Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
1/78
Organization (S): EDF/MTI/MN
Handbook of Référence
R5.01 booklet: Modal analysis
Document: R5.01.01
Algorithm of resolution for the problem
generalized
Summary
In this document, we present the algorithms of resolution for the generalized modal problems which
are established in Code_Aster via operators MODE_ITER_INV and MODE_ITER_SIMULT:
· method of the iterations opposite,
· method of Lanczos,
· method IRA (known as of “Sorensen”),
· method of Bathe & Wilson.
One gives to the reader the properties and the limitations, theoretical and practical, of the modal methods approached
while connecting these considerations, which can sometimes appear a little “éthérées”, to a precise parameter setting of
operators. For each method, one recapitulates in the form of tables the aforementioned parameter setting with his values by
defect and of the references to the paragraphs of the document.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
2/78
Contents
1 Introduction - Description of the document ................................................................................................. 4
2 Context ................................................................................................................................................. 6
2.1 Problems ................................................................................................................................. 6
2.2 Taking into account of the conditions limit ........................................................................................... 7
2.3 Properties of the matrices ................................................................................................................... 8
2.4 Properties of the clean modes ........................................................................................................ 9
2.5 Estimate of the real spectrum ............................................................................................................. 10
2.6 Establishment of the test of Sturm ....................................................................................................... 14
2.7 Spectral transformation ............................................................................................................... 16
2.8 Modal calculation .................................................................................................................................. 17
2.9 Establishment in Code_Aster .................................................................................................. 19
3 Method of the powers opposite (MODE_ITER_INV) ....................................................................... 22
3.1 Introduction .................................................................................................................................... 22
3.2 Localization and separation of the eigenvalues ............................................................................. 23
3.2.1 Method of bisection ......................................................................................................... 23
3.2.2 Method of the secant ......................................................................................................... 24
3.3 Method of the powers opposite ................................................................................................ 25
3.3.1 Principle ................................................................................................................................ 25
3.3.2 Method of iteration of the quotient of Rayleigh ........................................................................ 27
3.3.3 Establishment in Code_Aster ......................................................................................... 28
3.3.4 Display in the file message ....................................................................................... 29
3.3.5 Summary of the parameter setting ............................................................................................... 30
4 Method of subspace (MODE_ITER_SIMULT) .............................................................................. 31
4.1 Introduction .................................................................................................................................... 31
4.2 Analyze of Rayleigh-Ritz ............................................................................................................... 31
4.3 Choice of the space of projection ..................................................................................................... 33
4.4 Choice of the spectral shift ........................................................................................................... 35
5 Method of Lanczos (METHODE = “TRI_DIAG”) ............................................................................. 36
5.1 Introduction .................................................................................................................................... 36
5.2 Theoretical algorithm of Lanczos .................................................................................................. 36
5.2.1 Principle ................................................................................................................................ 36
5.2.2 Estimates of errors and convergences ........................................................................... 38
5.3 Algorithm of Lanczos practices .................................................................................................... 40
5.3.1 Problem of orthogonality ..................................................................................................... 40
5.3.2 Capture multiplicities ...................................................................................................... 41
5.3.3 Phenomenon of Lanczos ....................................................................................................... 41
5.4 Complementary processing ....................................................................................................... 42
5.4.1 Detection of spaces invariants ............................................................................................. 42
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
3/78
5.4.2 Strategies of restartings ................................................................................................. 43
5.5 Establishment in Code_Aster ................................................................................................... 44
5.5.1 Alternative of Newmann & Pipano ........................................................................................... 44
5.5.2 Parameter setting ......................................................................................................................... 46
5.5.3 Warning on the quality of the modes ............................................................................... 47
5.5.4 Perimeter of use ............................................................................................................ 47
5.5.5 Display in the file message ........................................................................................ 48
5.5.6 Summary of the parameter setting ................................................................................................ 49
6 Algorithm WILL GO (METHOD = “SORENSEN”) ....................................................................................... 50
6.1 Introduction .................................................................................................................................... 50
6.2 Algorithm of Arnoldi ....................................................................................................................... 50
6.2.1 Principle ................................................................................................................................. 50
6.2.2 Estimates of errors and convergence ............................................................................. 52
6.3 The stakes ...................................................................................................................................... 53
6.4 Algorithm “Implicit Restarted Arnoldi” (WILL GO) .................................................................................. 54
6.5 Establishment in Code_Aster ................................................................................................... 56
6.5.1 ARPACK ............................................................................................................................... 56
6.5.2 Adaptations of the algorithm of Sorensen ............................................................................. 56
6.5.3 Parameter setting ......................................................................................................................... 57
6.5.4 Display in the file message ........................................................................................ 58
6.5.5 Summary of the parameter setting ................................................................................................ 59
7 Method of Bathe and Wilson (METHODE = `JACOBI) ........................................................................ 60
7.1 Principle .......................................................................................................................................... 60
7.2 Tests of convergence .................................................................................................................... 60
7.3 Establishment in Code_Aster ................................................................................................... 60
7.3.1 Dimension of the subspace .................................................................................................. 60
7.3.2 Choice of the initial vectors ................................................................................................... 61
7.3.3 Parameters in Code_Aster ........................................................................................... 62
8 Conclusion - Synthesis .......................................................................................................................... 63
9 Bibliography ......................................................................................................................................... 65
Appendix 1 Généralités on algorithm QR .............................................................................................. 67
Appendix 2 Orthogonalization of Gram-Schmidt ....................................................................................... 72
Appendix 3 Méthode de Jacobi .................................................................................................................. 75
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
4/78
1
Introduction - Description of the document
A majority of study concerning the dynamic behavior of structures are carried out in
carrying out a transitory analysis on modal basis. To exhume these modes of vibrations, one
string of algorithms were developed since about fifty years. In order to face
continual increase in the size of the problems and with the degradation of conditionings of
the discretized operators, only most effective and most robust, in practice, were built-in
in the two modal operators of Code_Aster.
The optimal perimeters of use of these operators can be dissociated. When it is about
to determine some eigenvalues (typically a half-dozen) or to refine some
estimates, operator MODE_ITER_INV is completely indicated. He gathers algorithms
heuristic and those of type powers (cf [§3]).
On the other hand, to capture a significant part of the spectrum, one A resorts to MODE_ITER_SIMULT.
This last federates the methods known as of “subspace” (Lanczos [§4], [§5], IRAM [§6], Bathe & Wilson
[§7]) which projects the operator of work in order to obtain an approximated spectrum of more reduced size (treated
then by a total method of type QR or Jacobi).
Until now, these algorithms stumbled regularly on the same shelves: detection
correct of multiple modes, modes of rigid body and a general way, processing of
packed spectrum. All this led to the appearance of “phantom” modes sometimes badly easy to detect
(modes corresponding to multiplicities missed and being able to generate correct residues with the direction
of Aster, and this, more especially as the criteria of the post-modal checks were sometimes permissive
(residue in 102 instead of the current 106), decontaminated even insufficient (test of Sturm limited to the values
clean positive) which causes distortions of results downstream from calculation, during projections on
base modal.
To be been free from the recurring problems to this type of approach, one thus proposed to enrich
MODE_ITER_SIMULT (starting from V5) of algorithm IRA (“Implicit Restarted Arnoldi” [§6]). This
alternative of Arnoldi, initiated by D.C. Sorensen in 1992, makes real great strides for the resolution of large
modal systems on parallel supercomputers.
It tries to bring an elegant remedy for the recurring numerical problems raised by the others
approaches. In short, IRAM gets a better total robustness Dans V5, it made it possible to balance
all of software anomalies related to the modal problems generalized (in V5, it allowed
to balance all software anomalies related to the modal problems generalized) while improving
complexities calculation and memory for a fixed precision, and, it allows a real control on
quality of the modes via a suitable parameter setting. Its use (method by defect in
MODE_ITER_SIMULT) is thus with advising in all the cases of figures.
This document is articulated around the following parts:
· initially, one recalls the context of modal calculation in Code_Aster to
through matrices and limiting conditions used, particular properties of modes
clean exhumed and of the difficulties of estimate of the spectrum. One recapitulates also what it is necessary
to know with minima about the spectral transformations, the families of modal algorithms and
on the broad outline of the course of a modal calculation in the code,
· in the second time, one describes the method of the powers opposite and his phases
preliminaries of localizations of eigenvalues (method of bisection and the secant)
installations in MODE_ITER_INV,
· the third part treats the general framework of the methods of subspace installation in
MODE_ITER_SIMULT. One details in particular the implications of the choice of the shift and the type of
subspace of projection in the parameter setting of the operator,
· the three following chapters show, in the order, the three methods of this operator:
Lanczos, IRAM and Bathe & Wilson,
· in the appendices one approaches more in details of the processes “of second level” to which make
call three preceding methods. It is about the algorithms QR and Jacobi which are known as
general practitioners because they allow to capture all the spectrum of an operator. The first,
algorithm QR, is fundamental because it intervenes in the majority of the methods. It is
the algorithm of reference which offers a very great robustness but calculation complexities and
memory prohibitory. One approaches also the various algorithms of orthonormalisation put in
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
5/78
place in the code. Indeed, one will cease hammering, throughout this document, only
their quality and their robustness are crucial for the good unfolding of the algorithms.
This last version of the document was almost entirely rewritten while taking as a starting point the the indices
written precedents, respectively by D. Seligmann [R5.01.01] index A and B. Quinnez [R5.01.01]
index B, this in order to try to approach the code while clarifying more in details them
characteristics of the methods and the subjacent numerical phenomena. A particular effort was
brought to put in prospect the choices led in Code_Aster compared to search,
passed and current, like clarifying the general philosophy of a modal calculation.
One gives to the reader the properties and the limitations, theoretical and practical, of the modal methods
approached while connecting these considerations, which can sometimes appear a little “éthérées”, to one
precise parameter setting of the operators. For each method, one recapitulates in the form of tables the aforementioned
parameter setting with its default values and of the references to the paragraphs of the document.
At the time of the first passages, one strongly engages the user has to modify only the parameters
principal noted in fat in these tables and with reading the traces of the file message. Others,
concerning more the mysteries of the algorithms, were initialized empirically with values
standards.
One tried constantly to bind different the items approached and to limit to the bare minimum the recourse
with long mathematical demonstrations. In any event, the many references which enamel
the text must make it possible to seek precise information.
The object of this document is not to detail all the aspects approached, of the complete works having
already fulfilled this mission. One will quote in particular F. Chatelin [bib3], G.H. Golub [bib6], P. Lascaux
[bib11], B.N. Parlett [bib18] and Y. Saad [bib32], and one recommend the synthesis particularly
brought up to date and exhaustive made by J.L. Vaudescal [bib23].
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
6/78
2 Context
2.1 Problems
We consider the problem generalized with the eigenvalues:
To find (, U) such as A U = B U, U 0, where A and B are symmetrical matrices with coefficients
realities. This type of problem corresponds, in mechanics, in particular with:
· The study of the free vibrations of a not deadened and nonrevolving structure. For this
structure, one seeks the smallest eigenvalues or those which are in one
interval given to know if an exiting force can create a resonance. In this case,
stamp A is the matrix of rigidity, noted K, (possibly increased matrix of
geometrical rigidity noted kg, if the structure is prestressed) and B is the matrix of mass
or of noted inertia Mr. Les eigenvalues obtained are the squares of the associated pulsations
at the sought frequencies.
The system to be solved can be written: (K + kg) U = 2 M U where =
2 F is the pulsation,
F the Eigen frequency and U the vector of associated clean displacement.
· The search for linear mode of buckling. Within the framework of the linearized theory, in
supposing a priori that the phenomena of stability are suitably described by
system of equations obtained by supposing the linear dependence of displacement by report/ratio
at the level of critical load, the search of the mode of buckling U associated on this level with
critical load, is brought back to a problem generalized to the eigenvalues of the form:
(K+ kg) U = 0 with K stamps rigidity and kg stamps geometrical rigidity.
Note:
· this type of generalized modal problem is treated in Code_Aster by two operators:
MODE_ITER_INV and MODE_ITER_SIMULT. Each one having its perimeter of application, its
functionalities and its limitations,
· in V5, the user can specify the class of membership of his calculation by initializing it
key word TYPE_RESU with “DYNAMIQUE” (default value) or with “MODE_FLAMB”. Display of
results will then be formatted by taking account of this specificity. In the first case one
will speak about frequencies whereas in the second, one will speak about critical load,
· in the presence of depreciation and gyroscopic effects, the study of dynamic stability
of a structure leads to the resolution of a modal problem of a nature higher, known as quadratic:
(K +i C-2 M) U = 0. It is solved by the two modal operators and is the object
of a specific note [R05.01.02].
Now that links between the mechanics of the structures and the resolution of modal problems
generalized were recalled, we will be interested in the processing of the conditions limit in
code and with their incidences on the matrices of mass and rigidity.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
7/78
2.2
Taking into account of the limiting conditions
There are two ways, during the construction of the matrices of rigidity and mass, to take into account
the boundary conditions (this description in dynamic term of problem is extrapolated easily with
buckling):
·
double dualisation, by using ddls of Lagrange [R3.03.01], makes it possible to check
C U = 0 (CLL for Condition Limite Linéaire),
with C stamps real of size p X N (K and M are of command N). Matrices of rigidity and of
mass dualized have the form then
K
CT
CT
M 0 0
~
~
K = C - Id
Id
M = 0 0 0
C Id - Id
0 0 0
with and strictly positive realities which will be initialized to 1 without that involving a loss
of general information.
The dimension of the problem was increased by 2p, because to N ddls known as “physics”, one has
added ddls of Lagrange. There are two ddls of Lagrange by linear relation assigned to p
limiting conditions.
· The setting has zero of p lines and columns of the matrices of rigidity and mass. This is not
valid that for blockings of ddls. One cannot take into linear account of relation
and one will speak about kinematic blocking (CLB for Condition Limite de Blocage). Matrices
of rigidity and mass become:
~
K 0
~
M 0
K =
M =
0 Id
0 0
The dimension of the problem remains unchanged but it is however necessary to withdraw the participations of
ddls blocked with the components of the initial matrices (K is obtained starting from K in
eliminating the lines and the columns from the ddls which are blocked; idem for M).
When limiting conditions are imposed, the number of eigenvalues (with all their multiplicities)
really implied in physics (with the flat of modeling close (grid, frequency of
rupture…) phenomenon is thus lower than size N of the transformed problem:
3p'
·
N
= N -
with p'= 2 p
ddl - active
(double dualisation),
2
·
N
= N - p
ddl-credits
(kinematic blocking).
Framed below the display dedicated to these parameters in the file message shows.
------------------------------------------------------------------------
THE NUMBER OF DDL
TOTAL EAST: 220 N
LAGRANGE EAST: 58 p'= 2p
THE NUMBER OF ACTIVE DDL EAST: 133 nddl_actifs
-----------------------------------------
Example 1: Management of the ddls
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
8/78
In addition, in the modal calculation algorithms, one must ensure oneself of the membership of the solutions with
acceptable space. One brings back oneself there via auxiliary processing. Thus when blockings are used
kinematics (CLB), it is necessary in the various algorithms and for each iteration, to use a vector “of
positioning " ubloq, defined by:
· if
ième ddl is not blocked ubloq (I) = 1,
· if not
ubloq (I) = 0,
u1 (I) = u0 (I). U
(I) (I =… N
bloq
) u1
1
If one uses the method of double dualisation, one needs a vector of positioning of the ddls for
lagrange ulagr definite like ubloq. It is only used at the time of the choice of the random initial vector.
So that this vector u0 checks the conditions limit (CLL) one operates in the following way:
u1 (I) = u0 (I). U
(I) (I = 1… N
lagr
) u2
~
K u2 = u1
Thereafter, to simplify the notations, we will not make the distinction between the initial matrices and theirs
hanging dualized (noted with a tilde) that if necessary. Very often, they will be indicated by A
and B in order to approach the usual modal notation without being attached to such or such class of
problems.
2.3
Properties of the matrices
As we wrote previously the matrices considered are symmetrical and with
real coefficients. According to the cases of figure indexed in the table below, they can be
defined positive (noted > 0), semi-definite positive (0), indefinite (0 or 0) even singular
(S).
Structure
Lagranges *
Buckling
Fluid
free
structure
K
0 and S
< 0 or > 0
> 0
0 and S
M
> 0
0 and S
0 or 0
> 0
(resp. Kg)
Table 2.3-a: Propriétés of the matrices of the generalized problem
* This column relates to of course the properties of the matrices dualized made up to leave
initial matrices [R3.03.01]
Columns of this mutually being excluded table, in practice, a problem of using buckling
Lagranges doubles to model some of its limiting conditions, sees its dualized matrices
~
~
(K and K G) to become potentially indefinite.
Note:
This range of properties must be taken into account at the time of the choice of the couple operator of work -
(pseudo) produced scalar. This framework can thus reinforce, with effectiveness and transparency,
robustness and the perimeter of the modal calculation algorithm in all the encountered cases of figure
by Code_Aster.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
9/78
The following paragraphs go enable us to measure the incidence of these properties on the spectrum
of generalized problem.
2.4
Properties of the clean modes
Let us recall first of all that if the matrix of modal problem standard A U = U is real
symmetrical, then its own elements are real; The clean elements of a matrix are its
clean values and its vectors. In addition, A being normal, its own vectors are orthogonal.
In the case of the problem generalized A U = B U, this condition is not sufficient. Thus,
let us consider the following generalized problem:
1
1 U
1
0
U
1
1
=
1
0 U
0
- 1
U
2
2
1
1
- ±
its own modes are ± = (1± I 3) and u± =
2
1 + 2 1.
±
If one adds the assumption “one of matrices A or B is definite positive”, then the problem
generalized has its real solutions. There is even the characterization (sufficient condition) more precise
following.
Theorem 1
Are A and B two symmetrical matrices real. If there exists and such as A + B
that is to say definite positive, then the generalized problem has its real own elements.
Proof:
This result is obtained immediately by multiplying the problem generalized by and by carrying out one
spectral shift. The problem then is obtained (A +)
B U = (+) B U. Like
(A +)
B is definite positive, it admits a single decomposition of Cholesky in the form
CCT with C stamps regular.
1
The problem is written C B then
1
C-T Z = µ Z with Z = CT U and µ = +, which allows
to conclude, because matrix C-1B C-T is symmetrical.
Notice
This characterization is not necessary, thus the generalized problem associated the matrices
With = diag (1, - 2, -)
1 and B = diag (- 2 1
,)
1 admits a spectrum real all while not answering
condition of definite positivity.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
10/78
Proposal 2
If the matrices A and B are real and symmetrical, the clean vectors of the generalized problem are
With and B - orthogonal, which means that they check the relations
C
B U
I
J = I
J has J
C A U
I
J = J I
J has J
where J has is a scalar depend on the standard of the jième clean vector, ij is the symbol of Kronecker
and U J is the clean vector associated the eigenvalue J.
Proof:
Immediate for distinct eigenvalues, by writing A and B - product scalar between two
couples (I, J) and (J, I), then by using the symmetry of the matrices (cf [bib9]).
Note:
· it is shown that A and B orthogonalities of the clean vectors are a consequence of
hermiticity of the matrices. They are clearly a generalization of the properties of the problem
square standard (even normal),
· orthogonality compared to the matrices does not mean especially that the clean vectors are
orthogonal for the traditional euclidian norm. This one can be only the fruit of
particular symmetries (cf TP n°1 [bib25]),
· this property simplifies calculations of modal recombinations (DYNA_TRAN_MODAL
[R5.06.01]), when one handles matrices of generalized rigidity and mass which are
diagonals. The kj quantities = J aj and mj = aj are called, respectively, modal rigidity
and modal jième mode masses.
Knowing that the modes are real, we now will worry we about their estimate.
2.5
Estimate of the real spectrum
Because of membership of the spectrum to the real axis, the problems of counting of eigenvalues go
to be largely simplified. One does not have to treat régionnement complex plan and one can
to be based on the corollary of the law of inertia of following Sylvester.
Corollary 3
Are A and B two real matrices symmetrical, B being of more definite positive. The number of
eigenvalues, strictly lower than, of the problem generalized A U = B U is then equal to
a number of strictly negative diagonal coefficients of the matrix D such as (A - B) = LDLT.
Proof:
Cf paragraph n°1 of the article of Y. Haugazeau [bib7].
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
11/78
Remarks
· this corollary extends to the square matrices and results from the properties of the series Sturm
(N)
and of the principle of inclusion (one notes p
()
(N)
(N)
= die (your - B) the characteristic polynomial
(N) (N)
(N) (N) (N)
problem generalized A
U
= B U obtained by removing N last
()
lines and columns of matrices A and B. the continuation of polynomials (p N) constitutes a continuation of
N
Sturm because the roots of the ième polynomial frame those of (i+1) the ème. This property
of interlacing of the spectrum of symmetrical real submatrices “principle is called
of inclusion ") [bib7], [bib9],
· the possible multiple eigenvalues are taken into account their multiplicity,
· thereafter, one will call modal position of and one will note it pm (), this number of
strictly negative diagonal coefficients.
This corollary thus makes it possible to easily determine the number of eigenvalues contained in one
interval [, µ] and the modal position of these eigenvalues in the spectrum. It is enough to carry out two
decompositions LDLT, that of (A - B) and that of (A - µ B) and to enter the difference
number of strictly negative terms between the two diagonal matrices. In the jargon of the code,
one indicates (improperly besides!) this test under the term of “test of Sturm”.
It must however be extended to the quite particular shapes of the matrices met in
Code_Aster, and in particular, it is necessary to be able to take into account buckling with or without Lagrange.
With this intention, one widened the criterion with a matrix B unspecified and one generalized it while taking in
count the dualized matrices.
Let us recall that one defines usually the signature of a matrix (and that of the quadratic form
associated) like the triplet of natural entireties (R, S, T) where R indicates the number of eigenvalues > 0, S it
a number of null eigenvalues and T concealment be < 0. The latter entirety is thus what one notes pm ()
in our case of figure.
Property 4
Are A and B the two real matrices symmetrical (of command N) related to the generalized modal problem
~
~
(S): With U = B U. Let us note A and B, their matrices associated resulting from the double dualisation with p
Lagranges making it possible to check C U = 0 (CLL), with C stamps real of size p X N.
~
~
Then, the signature of (A - B) is written, by noting card [,
B] the number of values has
clean of the generalized problem included in the interval [has, B]:
if B is indefinite and A is definite positive
then card {}
0 = 0 and
R
= card] -, [+ card [, 0+ [
if < 0: S = p + card {}
éq 2.5-1
T
= (
pm) =
card],]
0 + p
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
12/78
R
= card] -, + [
if = 0: S = p
éq 2.5-2
T
= (
pm) = p
R
= card] -,]
0 + card], + [
if > 0: S = card {} + p
éq 2.5-3
T
= (
pm) =
card [,
0 [+ p
if B is definite positive and if A is semi-definite positive
then card]
,
- [
0 = 0 and
R
= card], 0+ [
if = 0: S = p + card {}
0
éq 2.5-4
T
= (
pm) = p
R
= card], + [
if > 0: S = p + card {}
éq 2.5-5
T
= (
pm) =
card [,
0 [+ p
Proof:
To impose boundary conditions linear, a technique of double dualisation is used
[R3.03.01] which leads to the transformed generalized system
WITH CT
CT
B 0 0 U
(~
~
-) ~ = C - Id
Id - 0 0 0 v = 0
(~
With
B U
S)
C
Id
- Id
0 0 0 W
J
One will note thereafter the VI null vector column of size p (i= 1. .p; j= 2,3) except with index I, for
which it is worth 1, and 0n, the null vector column of size N. Now let us consider the three families of
vectors following:
u1
I
~
·
N - p V1 vectors = v1
I
I which is the clean vectors of the system (S). Clearly,
v1i ~ ~
according to proposal 2, they are A and B - orthogonal (according to the properties of A and B) and
one notes their eigenvalues I.
0
N
·
p independent vectors V2 = v2
I
I,
v2
I
0
N
·
p independent vectors V3 = v3
I
I.
- v3i
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
13/78
U
These N + p vectors form a base B of space E = ~u = v/Cu = 0 of the solutions
W
acceptable of the dualized problem (as a free family of a space of finished size).~
~
Let us consider, for a given real number, the quadratic form associated matrix (A - B)
(~)
~
~ ~
= (-), ~, ~
U
With
B U U
U.E.
While breaking up on the basis B generated by the preceding vectors
~
U =
a1 V
1
2
2
3
3
I
I
+
V has
I I + has V
I I,
i=,
1 N p
i=,
1 p
i=,
1 p
one obtains
~
T
1 2
1
1
3 2
(U) = (A) (-) (U B U) + (- 4) (has
I
I
I
I
I)
i=,
1 N p
i=,
1 p
J
(one used in particular the properties of orthogonality of the families of vectors VI and the relation
(u1, CTv2
1
2
I) = (C U, v
I
I
I) = 0). From where, by noting Li the linear form which associates ~
U its ième
coordinated in the base B
~
~
T
2
1
1
2
~
2
~
(U) = Li (U) (I -) (U B U
I
I) + 0 Ln-p+i (U) + (- 4) Ln+i (U)
i=,
1 N p
i=,
1 p
i=,
1 p
Clearly them (Li)
are linearly independent from where an immediate reading of the terms of
i=,
1 n+ p
the signature according to the sign of the factors. Moreover, it is noticed that this decomposition comprises
obligatorily p zeros and p minus signs. The relations of the property result all then
naturally, by using the principle of inertia of Sylvester (who ensures us that this decomposition
is invariant), the relations of A and B - orthogonality of property 2, the properties of the table
[Table 2.3-a] and by noticing that:
~T ~ ~
U With U
C A U
~T ~ ~
I
J
I
J
U B U = C B U
I
J
I
J =
=
for J
0.
J
J
Restriction 0 of the case of figure (2) comes owing to the fact that if B is definite positive then the spectrum
problem (S) is positive and thus the shift with a strictly negative value does not have any interest (one
find within the framework of application of corollary 3).
According to table [Table 2.3-a] this property applies to the matrices handled by
Code_Aster and one can build the following corollary.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
14/78
Corollary 5 (theoretical wide Sturm)
In the matric configurations of Code_Aster, the number of clean modes of the problem
generalized (S) whose clean vector checks the linear conditions limiting (CLL) and of which the value
clean is contained in the interval], µ [is
If µ 0
card
], [= p (
m
) - p (
m
µ
µ
),
if not
card], [= p (
m
) + p (
m
µ
µ
) - 2p.
If no dualisation of Lagranges is required, p= 0 is posed.
Proof:
It is immediate by combining the results of the table with those of property 4, and, while noticing
that the generalized modal configurations of Code_Aster can be only of two types:
· buckling
A is definite positive and B is indefinite, the spectrum can be negative and one uses
relations of the type (1),
· dynamics
A is semi-definite positive and B is definite positive, the spectrum is positive and one
can use the relations of the type indifferently (1) or (2).
Of course, if no dualisation of Lagranges is required, the same reasoning applies in
~
~
posing p = 0,
=,
~
With
WITH B = B and U = U.
To finish, the principle of inertia of Sylvester ensures us that the signature of factorization
(~
~
-) ~~ ~
With
B = LDLT is identical to that of the shiftée matrix. In spite of this transformation, them
modal positions thus remain many perennial information.
However this corollary is tributary of arithmetic exact, in practice it is necessary to adapt it and control
its application.
2.6
Establishment of the test of Sturm
This test is confronted into arithmetic finished with two concomitant problems:
· the factorization of the shiftée matrix when is very close to an eigenvalue,
~
· the calculation of the strictly negative terms of D to evaluate the modal position
(
pm).
The first point requires, au préalable, the determination of a criterion of membership of the shift with
spectrum of the problem. In Code_Aster that Ci is founded on the loss of decimals at the time of
~
~
~ ~ ~
factorization of matrix (A - B) = LDLT. More precisely, is regarded as being one
eigenvalue, if during factorization one loses more decimal NPREC_SOLVEUR (in any rigor it
is only one test of bad numerical conditioning). It is then necessary to modify the value of in
shifting this shift of PREC_SHIFT % following the algorithm:
For I = 1, NMAX_ ITER_ SHIFT
If loss of more than decimal NPREC_ SOLVEUR
then (1+ PREC_ SHIFT),
If not
exit;
Fine loops.
Algorithm 1: Shift of the shift
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
15/78
If at the end of tentative NMAX_ITER_SHIFT, the matrix is still not numerically invertible,
one continues all the same calculation with the last shiftée value. Framed below shows
trace of such a shift in the file message.
ONE MODIFIES THE VALUE OF SHIFT OF: 5.00000E+00POURCENT
THE VALUE OF SHIFT BECOMES: 5.96840E+04
ONE MODIFIES THE VALUE OF SHIFT OF: 5.00000E+00POURCENT
THE VALUE OF SHIFT BECOMES: 5.66998E+04
VALEUR_MIN IN FREQUENCY EAST: 2.00000E-01
VALEUR_MAX IN FREQUENCY EAST: 5.78813E+01
THE VALUE OF SHIFT IN FREQUENCY EAST: 3.78975E+01
Example 2: Shift of the shift
In algorithm 1 if F
then = - F
corig
corig. This parameter F corig corresponds to a value
threshold in lower part of which it is considered that one has a numerically null eigenvalue (in
traditional mechanics, that corresponds to a mode of rigid body, cf [§5.5.4]). The imposition
= - fcorig thus makes it possible to dissociate these modes of the remainder of the spectrum and to avoid instabilities
numerical during the test of Sturm. This threshold is skeletal with key word SEUIL_FREQ.
Note:
· the preceding key words are accessible starting from the key word factor CALC_FREQ from both
modal operators,
· Y. Haugazeau [bib7] proposes more generally to deal with these problems of pivots
numerically “small” by building a matrix, unitairement similar (i.e.
similar via a unit matrix of passage (orthogonal in reality) to the shiftée matrix, of which
factorization would have less instability.
By taking into account these elements and knowing that, numerically, the calculation of the pivots
strictly negative of the diagonal matrix includes in fact also the elements (theoretically) null of
the signature, one can rewrite the preceding corollary.
Corollary 5bis (numerical wide Sturm)
According to the assumptions of corollary 5, there is the numerical criterion of accountancy of the modes according to:
If µ > 0
card [,] = p (
m
) - p (
m
µ
µ
),
if µ < 0
card [,] = p (
m
) + p (
m
µ
µ
) - 4p
Proof (heuristic):
The fact is applied that numerically the operator of factorization provides the “modal position
numerical "
(
pm) = S + T with property 4 and corollary 5. In addition, establishment of the criterion
neither the nullity of the product allows, nor the estimate of card {}.
Note:
This criterion is used in preprocessings in MODE_ITER_INV (options “AJUSTE” or “SEPARE”) and
in MODE_ITER_SIMULT (option “BANDE”), in postprocessings to check the validity of the number
modes (MODE_ITER_SIMULT) and in auxiliary controls (IMPR_STURM).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
16/78
Now that we are able to enter the spectrum of the generalized problem, it remains with
to determine! The generic algorithms being intended for the standard problems, it is necessary to transform
our initial problem.
2.7 Transformation
spectral
These techniques make it possible to answer objective triple:
· to exhume a standard modal problem,
· to direct the search of the spectrum,
· to separate the eigenvalues.
The spectral calculation algorithms converging of as much better than the spectrum (of work) which they treat
is separated, these techniques can be regarded as prepacking of the problem of
departure. They make it possible to return the separation of certain modes much more important than
those of other modes, and to improve their convergence thus.
Most widespread of these transformations is the technique known as of “shift and invert” which consists with
to work with operator A such as:
µ
-
With U = B U (A -)
B 1
1
B U =
U
!# “
#
$
# #
-
With
!“$
µ
Be reproduced 2.7-a: Effet of the “shift and invert” on the separation of the eigenvalues.
The figure [Figure 2.7-a] shows that this separation and this orientation of the spectrum of work in µ are
had with the particular properties of the hyperbolic function. In addition, it is observed that only them
eigenvalues are affected by the transformation. At the end of the modal process it is thus enough to
to pass by again in the plan of by a suitable change of variable.
Note:
· the variable is usually indicated by the term of “shift” or spectral shift,
· the matrix of work A must of course be invertible, that can become one besides of
motivations of this shift (cf [§5.5.1]).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
17/78
For memory, let us note that with a shift complexes several scenarios are possible:
· to work completely into arithmetic complex,
· into arithmetic real,
· mixer two steps by isolating the real and imaginary contributions from A
~
~
1 1
1
With =
(
Re With
) µ =
+
,
2
- -
~
~
~
~
1 1
1
With =
(
Im A
) µ =
-
.
2i
- -
Each one of its approaches has its advantages and its disadvantages. For the quadratic problems,
it is for the moment the second solution which was adopted in the code.
Remarks
· this choice of the operator of work is indissociable among that of (pseudo) - scalar product. It
allows to direct itself towards such or such algorithm and can thus influence the robustness of calculation,
· another class of spectral transformation, with double shifts, allows to select them
eigenvalues located on the right vertical axis. It is about the rational transformation of
- 1
-
Cayley:
(A - B
2
1
) (A - B
2
) U =
U
.
!# # # “
# # # # $
#
- 1
With
!“$
µ
We now will enumerate the various families of methods allowing to solve it
standard modal problem.
2.8 Calculation
modal
The methods of calculation modal can gather in (at least) four families:
· The algorithms of the type QR (cf [Annexe 1]) which were presented per H. Rutishauser (1958)
and formalized jointly by J.C. Francis and V.N. Kublanovskaya (1961). It is one
fundamental algorithm often implied in the other methods.
Perimeter of use: Calculation of all the spectrum.
Advantages: Good convergence, robustness, calculation of the form of Schur (any matrix
complex A admits the decomposition of Schur Q * A Q = T, with Q and T of the matrices
respectively, unit and triangular higher. In the real case, T is only diagonal by
block 1x1 or 2x2 (real form of Schur). Columns of the matrix unit, known as of Schur,
vectors of Schur are called).
Disadvantages: Prohibitory memory complexities and calculation, sensitivities to the “balancing (them
techniques of balancing (English balancing) consist in transforming reversibly
the operator of work so that its matric standard decreases. Thus its handling will be
less sensitive to the effects of round (cf [§10.3]))“.
Alternatives: With implicit shift or clarifies, simple or double…
· The algorithms of the powers type which were historically developed the first for
to solve generic modal problems. These are basic algorithms of which them
others are an improvement. They are implied in operator MODE_ITER_INV (cf [§3]).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
18/78
Perimeter of use: Calculation of the extreme values of the spectrum.
Advantages: Simplicity, very good estimate of the clean vector in some iterations.
Disadvantages: Convergence can become problematic, bad capture of
multiplicities, of the clusters… (together of eigenvalues “very close”).
Alternatives: Powers opposite, (Bi) iteration of the quotient of Rayleigh (cf [§3.3.2])…
·
methods of subspace which consist in projecting the operator of work on a space
H such as the spectrum of the projected operator is a good approximation of the part of
initial spectrum which one seeks. These algorithms are the hard core of the operator
MODE_ITER_SIMULT (cf [§4]).
Perimeter of use: Calculation of part of the spectrum.
Advantages: Reduction of the size of the problem and memory complexities and calculation,
require only the calculation of a product matrix-vector and not the knowledge of
all the matrix.
Disadvantages: Use the many pre one and postprocessings, convergence can become
problems, more or less easily capture the multiplicities and the clusters according to
alternatives.
Alternatives: Iterations of subspace, Bathe and Wilson (1971), Lanczos (1950), Arnoldi
(1951), Davidson (1975), Sorensen (1992)…
·
other approaches more or less empirical and are specialized. They are often
connected to other problems: search for roots of polynomials, functions
unspecified… One can thus quote the method of bisection used in preprocessings in
MODE_ITER_INV, but also that of Jacobi, Laguerre etc…
Note:
Many parallels can be led between these families (method QR is not thus
that a method of iterations of subspaces applied to entire space), but they
also lead to processes similar to those developed for other problems:
· optimization; The method of the quotient of Rayleigh is with the method of the powers opposite,
what the method of Newton is for a method of descent traditional,
· resolution of linear systems; The method of the combined gradient is a method of
subspace for the positive definite symmetrical systems,
· search for roots of polynomials: the method of the powers is a method of
Bernoulli applied to the matrix “companion” of the polynomial (cf [bib11] pp502/503).
The following paragraph will synthesize what precedes in the total flow chart by resolution of one
generalized modal problem of Code_Aster.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
19/78
2.9
Establishment in Code_Aster
This one can break up into four phases:
1)
The first operation consists in determining the shift like certain parameters of the problem.
That is carried out in a more or less transparent way according to the option of calculation chosen by
the user:
· MODE_ITER_INV + option: “SEPARE” or “AJUSTE” the shift is determined by the first
phase of the algorithm and the number of clean modes sought by frequential tapes
(provided by the criterion of Sturm) is limited by NMAX_FREQ,
· MODE_ITER_INV + option: “PROCHE” the shift is fixed by the user and the number of modes
clean is equal to the number of shifts,
· MODE_ITER_SIMULT + option: “PLUS_PETITE” the shift is null and the number of modes is
parameterized by NMAX_FREQ,
· MODE_ITER_SIMULT + option: “BANDE” the shift is equal in the middle of the tape fixed by
the user and the number of modes are determined by the criterion of Sturm,
· MODE_ITER_SIMULT + option: “CENTER” the shift is fixed by the user and the number of
clean modes is parameterized by NMAX_FREQ.
2)
In one second phase of preprocessings, one factorizes partially the matrix of work
WITH = (A - B)-1B
while being interested only in its reversed part. The action of this operator on an unspecified vector U,
noted u2, will be built thus more effectively
u1 = B U
(
u2
With - B) u2 = u1
This factorization undergoes the same risks besides as the criterion of Sturm when the shift is
near to an eigenvalue. One carries out the same shifts then according to algorithm 1.
Notice
In V5 when NMAX_ITER_SHIFT is reached calculation stops in fatal error
(only for this preprocessing), this problem of factorization which can give place to
numerical instabilities.
3)
Modal calculation itself is carried out: the standard problem is solved, then one returns to
initial problem. In MODE_ITER_INV two alternatives are available: method of the iterations
opposite (option: “DIRECT”) and its acceleration by quotient of Rayleigh (“RAYLEIGH”). For it
who is MODE_ITER_SIMULT, it allows the use of three distinct methods: method of
Bathe and Wilson (“JACOBI”), that of Lanczos (“TRI_DIAG”) and finally, that of Sorensen
(“SORENSEN”).
Each one of these methods have internal tests of stops. Without counting that methods
of projection employ auxiliary modal methods:Jacobi (cf [Annexe 3]) for the first
and QR/QL (cf [Annexe 1]) for the others. They require also tests of stops.
The user often has access to these parameters, although it is warmly recommended, at least
initially, to preserve their default values.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
20/78
4)
This last part gathers the tests of total stops which check the good unfolding of
calculation. They are of two types:
1
2 2
· (One normalizes recalls that U
U
= max I, A
max
With
=
ij, U =
U
2
I and
I
I
J
I
1
With = max I
(AA)
* 2
= max I
(A) if A square
2
residue of the initial problem
I
I
U
U U
If >
SEUIL_
FREQ then
With U - B U 2? < SEUIL,
With U 2
If not
With U - B U
? < SEUIL,
2
End if.
Algorithm 2: Test of the standard of the residue
This sequence is parameterized by key words SEUIL and SEUIL_FREQ, pertaining
respectively with the key word factor VERI_MODE and CALC_FREQ. In addition, it
postprocessing is activated by initialization with “OUI” (default value) of STOP_ERREUR in
the key word factor VERI_MODE. When this rule is activated and non-observance, calculation
arrète, if not the error is just announced by an alarm. One would know of course only too
to recommend not to decontaminate this parameter preferential treatment!
· Counting of the eigenvalues
This postprocessing is set up in MODE_ITER_SIMULT and it makes it possible to check that it
a many eigenvalues contained in a tape test [1,2] are equal to the number
detected by the algorithm. The inclusion of the initial tape (in option “BANDE” it acts of
values indicated by the user, if not they are the extreme values of the calculated spectrum)
[I, F] in this tape test is led in order to detect possible problems of
clusters or of multiplicities at the initial boundaries (cf TP n°1 [bib25]).
+
1
I
2
-
I
F
F
Initial tape
Bandage tested
Appear 2.9-a: Comptage of the eigenvalues
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
21/78
While noting +
-
I and F, respectively, largest and the smallest eigenvalue not requested by
the user and including the initial tape (cf [Figure 2.9-a]), one has the algorithm of construction of
bandage following test:
+
-
I - I
F - F
If + exists
-
I
(resp.f) and if
< PREC_ SHIFT resp.
I
F
then
+
-
+
+
= I
I
resp.
1
= F
F,
2
2
2
If not
=
1
I (1 - S (
ign I) PREC_ SHIFT) (resp. =
2
F (1+ S (
ign F) PREC_ SHIFT),
End if.
Algorithm 3: Construction of the tape test
This sequence is parameterized by key word PREC_SHIFT of the key word factor VERI_MODE.
Framed below the trace of the error messages shows when these post-checkings
start.
CHECKING A POSTERIORI OF THE MODES
<E> <MODE_ITER_SIMULT> <VERIFICATION OF THE MODES> FOR THE CONCEPT >MOD_4< IT
MODE NUMBER 2 OF CRITICAL LOAD 8.7243E+07 A a STANDARD Of ERROR
OF 6.3919E-01
<E> <MODE_ITER_SIMULT> <VERIFICATION OF THE MODES> FOR THE CONCEPT >MOD_4< IN
The INTERVAL (- 1.0608E+08,1.8130E+08) IT THERE A THEORETICALLY 3 LOAD (S)
CRITICAL (S) AND ONE OF A CALCULEE (S) 4.
Example 3: Post-checkings
The activation of the second postprocessing is subordinated to the initialization of STURM with “OUI” in
key word factor VERI_MODE. When this rule is activated and not respected, the continuation of the events
is controlled by STOP_ERREUR (if “OUI” calculation stops, if not the error is just announced by one
alarm). It would not be known of course that too much to recommend not to decontaminate these parameters 'passes
droit'!
Now that the context of the generalized modal problems of Code_Aster was brushed, us
let us be interested more particularly in the methods of the power type and their establishment
in operator MODE_ITER_INV.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
22/78
3
Method of the powers opposite (MODE_ITER_INV)
3.1 Introduction
To calculate several eigenvalues of the generalized problem, the general method is not used
iterations opposite just as it is.
One can, for example, to combine it with a technique deflation (technique of consistent filtering with
to transform the operator of work so that it has the same eigenvalues except some
preset modes which see themselves affecting a zero value) so as to filter automatically
spectral information updated more to find it with the following iteration. With deflation by
restriction (of other types of deflation exist, such as the method of Ducan-Collar which uses for
to filter spectral information the first line of the matrix and the clean vector. These techniques
vectorial spread of course per blocks) of Wielandt one builds the operator repeatedly of
work (in the symmetrical case), by noting (µk, the U.K.) the mode to be filtered,
With
With
1
- µ
U C
K + =
K
K
K
K.
This strategy, which does not apprehend the multiplicities, must be supplemented by a criterion of Sturm.
In addition, the fact of working into arithmetic finished and of not building the operator indeed
With K with each iteration, constrained with mâtiner this step of process of orthogonalization
powerful.
All these complications resulted in choosing another way, which breaks up into two parts:
1) localization of the eigenvalues (determination of an approximate value of each value
clean contained in an interval given by a technique of bisection, refined or not,
by a method of the secant),
2) improvement of these estimates and the calculation of their own vectors associated by one
method of iterations opposite.
The search for a value approached for each eigenvalue considered is selected in
key word factor CALC_FREQ via OPTION:
· if
OPTION = “SEPARE”, in each interval of frequencies defined by key word FREQ,
an approximate value of each eigenvalue contained in this interval is calculated in
using the method of dichotomy (cf [§3.2.1]),
· if
OPTION = “AJUSTE”, one carries out the same operations first of all as previously and
then, on the basis of these approximations, one refines the result by the method of the secant
(cf [§3.2.2]).
For these two options, one calculates at the same time the modal position of each value
clean what makes it possible to detect the multiple modes. Either only the NMAX_FREQ are retained
the lowest frequencies contained in the maximum interval specified by the user, is
one calculates all the values of this interval (if NMAX_FREQ = 0).
· if
OPTION = “PROCHE”, the frequencies given by key word FREQ, are considered
like the approximate values of the sought eigenvalues.
Note:
· of course, as one already specified (cf [§2.8]), this operator is to be used only for
to determine or refine some eigenvalues. For a wider search it is necessary
to use operator MODE_ITER_SIMULT,
· with option “PROCHE” one cannot calculate multiple modes.
· it is an expensive algorithm because it calls much upon the test of Sturm and thus with its
associated factorizations.
· the terminals of the intervals of search are provided by FREQ or following CHAR_CRIT
the initialization of TYPE_RESU.
We now will detail the various algorithms (and their parameter settings) which are put in
place in the first part of the process.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
23/78
3.2
Localization and separation of the eigenvalues
3.2.1 Method of bisection
As one already saw previously, the corollary 5bis of the law of Inertie de Sylvester (cf [§2.5])
allows to determine the number of eigenvalues contained in an interval given while carrying out
two decompositions LDLT. This criterion can thus result in refining the interval until not including
that an eigenvalue. This control being set up, one passes from one iteration to the other by using it
principle of the dichotomy.
On a starting interval [has, B], one thus operates in the following way, knowing
(
pm has) and
(
pm b):
1
Calculation of * = (+).
2
has
B
Calculation of
(
pm *).
If
(
pm *) =
(
pm has) (resp. p (
m b) then
to start again on [*, B] (resp.[, has *]),
If not
to start again on [, has *]
and on [*, B],
End if.
Algorithm 4: Method of Bissection
The process is stopped if one cut out more NMAX_ITER_SEPARE time the interval of departure, or if
- B has
for a given interval, one has
SEPARATE PREC
*
_
(in this case one does not refine any more
seek in this interval).
One obtains finally a list of frequencies. In each interval defined by the arguments of this
list, is an eigenvalue having a certain multiplicity. Like approximation of this value
clean, the medium of the interval is taken.
Note:
· one could have used as criterion the change of sign of the characteristic polynomial, but
in addition to the fact that it is very expensive to evaluate, it does not make it possible, just as it is, to detect them
multiplicities,
· the initialization of the process can be carried out in an empirical way according to the needs for
the user. To include part of the spectrum one can also use the régionnements
plan complexes theorems of Gerschgörin-Hadamard (on A, A T…). Accordingly,
the method of bisection can prove more effective than a QR in the presence of cluster. Its
convergence, although linear, is indeed raised by ½ whereas that of QR can tend
towards 1 [bib11].
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
24/78
has
*
B Intervalle of will tra
1 era division
p ()
has
has
*
B
*
B
the 2nd division
p ()
has
B
B
has
B
has
*
*
*
the 3rd division
p ()
has
B
has
B
has
B
has
*
*
*
*
B
the 4th division
Appear 3.2.1-a: Méthode of bisection
These approximated eigenvalues can be improved by a unidimensional minimization
seeking to cancel the characteristic polynomial (
p) = dét (A - B). But only the calculation of one
value of this determinant requires NR! operations.
3.2.2 Method of the secant
The method of the secant is a simplification of the method of Newton-Raphson. With the stage
unspecified K, knowing a value K and by approximating the nonlinear function (
p) by its
tangent in this point, one determines the next value k+1 as being the intersection of this line
with the axis of, and so on, according to the iterative diagram
p ()
K -
K - 1
K +1 = K - (
p K) (
p
K) - (
p k-1)
k-1
k+1
K
Appear 3.2.2-a: Méthode of the secant
The tangent being approximated by a difference finished in order not to have to calculate of derived from
(
p), only the estimate of the polynomial is necessary.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
25/78
k+1 - K
It is considered that one reached convergence when
<
PREC_ AJUSTE and, in addition, one
K
limits itself to NMAX_ITER_AJUSTE iterations if this criterion is not reached.
Note:
This method has a quadratic convergence when it is close to the solution, in the case
opposite, it can diverge. From where interest to combine the method of bisection with this
approach.
We now will detail the algorithm of the powers opposite (coupled to an acceleration of
Rayleigh) constituting the second part of the process.
3.3
Method of the powers opposite
3.3.1 Principle
To determine the eigenvalue of the problem generalized A U = B U nearest in module to
, one applies the method of the powers to operator (A -
-
1
B) B. In fact, one only builds
the factorized matrix (this notation should not be confused that symbolizing the “shift and invert”
WITH = (A -
)
B - 1B
) A = (A -)
B and that amount dealing with the generalized problem
(- 1
1
With) Drunk =
U
-
. Method of the powers opposite converging proritairement towards
µ
!“$
eigenvalues of stronger module (in µ), one will thus capture the closest to the shift.
The principle is as follows, knowing an estimate of the sought eigenvalue and on the basis of one
initial vector standardized y0, one builds a clean vector series approximate (yk) by
K
recurring formula
For F
K = 1,… anger
~
With yk = B yk-1,
~
µk = yk,
~
yk
yk =
,
µk
Fine loops.
Algorithm 5: Method of the powers opposite
With 1, the 2… first eigenvalues (of clean vector ui) closest in module to,
it is shown that one has a linear convergence of y
U
K
1 and one quadratic convergence of
1
~
yk (I)
µ
1
K
and
I = 1. .n (the factor of convergence of these continuations is of
1 -
yk-1 (I)
(
)
1 -
1 -
the command of
).
min I -
i1
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
26/78
Note:
· this result is not acquired (with the factor of convergence near) only when the eigenvalue
dominant (here that nearest in module to the shift) is single. In the contrary case,
results remain however possible, even in the complex framework, but the analysis
rigorous of the convergence of this algorithm is still incomplete [bib6], [bib23],
· in theory, these results require that the initial vector is not orthogonal with
clean subspace on the left required. In practice, the round-off errors avoid it
problem (cf [bib11] pp500-509),
· even if the estimate of the eigenvalue is coarse, the algorithm provides one quickly very
good estimate of the clean vector,
· the major disadvantage of this method is that it is necessary to carry out a factorization of A for
each eigenvalue to calculate.
In general one works in euclidian norm or infinite standard, but to facilitate calculations
post-modal one seeks here B-to standardize the clean vectors (when B is indefinite, one works
with the associated pseudo-standard). The basic algorithm can be rewritten while posing Z = By
0
0
For K = 1… to make
With yk = zk-1,
z~k = B yk, T.
(
y
Z
y
K
K 1
K) =
-
yT. z~,
K
K
T
~
K =
(
sign Y.Z
K
K),
z~
Z
K
K = K
,
1
yT. z~ 2
K
K
Fine loops.
Algorithm 6: Method of the powers opposite with B-pseudo-standard
1
Let us note that (yk)
. This presentation avoids the products matrix-vector by the matrix B
1 -
during the calculation of the scalar products and already the coefficient of Rayleigh of the following paragraph precedes.
It is observed that the factor of convergence is all the more small as the spectral shift is close
sought eigenvalue and thus which A - B is close to the singularity. That is not in fact not
prejudicial with the process because the error made by solving the system is “mainly” in
direction generated by the clean vector which is the sought direction. That means that at the time of
resolution of A y = Z
K
K - 1, one does not find the solution exact y K but only the round-off errors
lead to solution close to the form ~
y = y + W
K
K
. This one is proportional to the solution
exact, but as standardization is arbitrary, all happens correctly [bib18].
Note:
This bad conditioning, far from having an unfavourable effect, improves even convergence of
the algorithm.
This algorithm is thus used to improve the clean vector associated with the value approximated with
phase 1. To refine this estimate of the eigenvalue, a quotient of Rayleigh is introduced.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
27/78
3.3.2 Method of iteration of the quotient of Rayleigh
Let us recall that the quotient of Rayleigh applied to the generalized problem is defined by the number
R (X), with X a vector not no one of N, such as:
xT A X
R (X) = xTB X
This quotient has the remarkable property of stationnarity in the vicinity of any clean vector and
to reach a extremum (local) which is the corresponding eigenvalue: for each X fixed, R (X)
minimize (A - B) X.
2
What we can translate by “if X is an approximation of a clean vector of the system
With X = B X, then R (X) is an approximation of the eigenvalue associated vector X ", and
reciprocally, we saw that if one had a good estimate of an eigenvalue,
method of the iterations opposite made it possible to obtain a good estimate of the clean vector
corresponding.
From where the idea to combine these two properties by considering the algorithm of iteration reverses with
spectral shift for which one revalues, with each iteration, the eigenvalue via the quotient of
Rayleigh. One then obtains the algorithm known as of iterations of the quotient of Rayleigh (in B-pseudo-standard)
For K = 1… to make
(A - R (yk-1)) Byk = zk-, 1
z~k = B y,
K T.
R (
y
Z
y
K
K 1
K) =
- +
T
~
R (yk-1),
y. Z
K
K
T
K =
(
sign Y. z~
K
K),
z~
Z
K
K =
,
K
1
yT. z~ 2
K
K
Fine loops.
Algorithm 7: Method of iterations of the quotient of Rayleigh with B-pseudo-standard
For the standard modal problem, one can show [bib18] that the convergence of this algorithm is
cubic if the operator of work is normal (a matrix A is known as normal if AA * =A * A.
It is thus the case of the operators square, antihermitiens or unit) (a fortiori in the case
square) and at best quadratic in the other cases.
If one used this method without modifying it, it would be necessary for each iteration of the process of improvement
of each eigenvalue, to carry out a factorization LDLT, which would be very expensive. From where the idea of
not to even carry out (this strategy of the coupling of methods to the complementary characteristics
antagonists is often used in numerical analysis. For example, in optimization, the method of
Levenberg-Marquadt couples one steepest-descent and Newton) this shift of Rayleigh, that if one
is in a vicinity (arbitrary concept to define) solution.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
28/78
Notice
Within a more general framework, certain authors introduced an algorithm known as “of Bi-iterations of
T
y With X
quotient of Rayleigh ". Based on the stationnarity of the Bi-quotient R (X,
B
y) = T
in the vicinity
y B X
clean vectors on the right and on the left, it provides (in nonsquare) the two types of vectors
clean. Its cost is however crippling because it requires twice more factorizations
(cf B.N. Parlett, 1969 [bib18]).
3.3.3 Establishment in Code_Aster
This spectral shift is activated only if key word OPTION of the key word factor CALC_MODE is initialized
in “RAYLEIGH”. By defect, there are “DIRECT” and the shift is then traditional (total correction rather
that progressive of the eigenvalue). The algorithm set up in the code cuts out as follows
(in B-pseudo standard):
· Initialization of the eigenvalue starting from the estimate of the first phase: .
· Calculation of an initial vector y0 random checking the limiting conditions.
· B - orthonormalisation of y0 compared to the modes previously calculated (if it is a mode
multiple according to the first phase) by Gram-Schmidt Modifié (since the development
of this operator, more effective and more robust methods of orthogonalization are
spread, such as Gram-Schmidt Modifiés Itératifs (IGSM) used in
MODE_ITER_SIMULT (cf [Annexe 3] and B.N.Parlett [bib18])) (GSM).
·
T
Calculation of 0 = sign (y0 By0).
· For
k=1,
NMAX_ITER to make
To solve (A -
B) yk =
B y
K - 1
K - 1.
B-orthonormalisation (possible) of yk.
T
y B y
Calculation of the correction of the eigenvalue C
K
K
=
- 1.
T
yk B y K
If yT B y
K
K
1 PREC
- 1 -
then
=
+ C,
exit;
If not
If OPTION = “RAYLEIGH” and if C
01
. then
=
+ C;
End if.
End if.
Fine loops.
Algorithm 8: MODE_ITER_INV
The standard of maximum error acceptable PREC and the maximum number of authorized iterations NMAX_ITER
are arguments of the key word factor CALC_MODE.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
29/78
Note:
· the clean vector B-being normalized, one considers to have reached convergence when the value
absolute of the unsuitable scalar product is close to the unit,
· to avoid taking an B-orthogonal initial vector with the sought eigenvalue one uses
a method of hard copy random for the components of this vector,
· in addition to be able to determine multiple or close modes, one is used
B-orthogonalization with the already calculated modes,
· the option “of acceleration” of the algorithm by quotient of Rayleigh being expensive, it is not
used with each iteration that if one is in the vicinity of the required eigenvalue.
3.3.4 Display in the file message
In the file message, the results are displayed in the form of table
Position
Frequency Deadened
Dichotomy
Method
Secant
Method
Opposite
modal
sow
Numbers
Numbers
Precision Numbers
Precision
iterations
iterations
iterations
*
*
0.
*
*
*
*
*
Table 3.3.4-a: Trace of MODE_ITER_INV in the file message
------------------------------------------------------------------------
THE NUMBER OF DDL
TOTAL EAST: 192
LAGRANGE EAST: 84
THE NUMBER OF ACTIVE DDL EAST: 66
------------------------------------------------------------------------
CHECKING OF THE FREQUENCY SPECTRUM (CONTINUATION OF STURM)
THE NUMBER OF FREQUENCIES IN THE TAPE (1.00000E-02, 6.00000E-02) IS 4
------------------------------------------------------------------------
MODAL CALCULATION: METHOD Of ITERATION OPPOSITE
OPPOSITE SECANT DICHOTOMY
NUMBER FREQUENCY (HZ) AMORT NB_ITER/NB_ITER/PRECISION/NB_ITER/PRECISION
4 1.97346E-02 0.00000E+00 4 6 2.97494E-07 4 1.22676E-07
5 2.40228E-02 0.00000E+00 4 5 4.21560E-05 3 4.49567E-09
6 4.40920E-02 0.00000E+00 3 5 2.19970E-05 3 2.62910E-09
7 5.23415E-02 0.00000E+00 3 5 2.34907E-07 5 1.32212E-07
------------------------------------------------------------------------
CHECKING A POSTERIORI OF THE MODES
------------------------------------------------------------------------
Example 4: MODE_ITER_INV
With option “PROCHE”, the columns “Dichotomie” and “Sécante” do not appear, while with
option “SEPARE”, only the column “Sécante” disappears. The last column precision gathers
intermediate data and does not represent, as for the other modal operator
MODE_ITER_SIMULT, the standard of error of the residue. It is an artefact which will be brought to disappear.
Let us recapitulate now the parameter setting of operator MODE_ITER_INV.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
30/78
3.3.5 Summary of the parameter setting
Key word factor
Key word
Default value
References
“DYNAMIC” TYPE_RESU
“DYNAMIQUE”
[§2.1]
“MODE_FLAMB”
[§2.1]
CALC_FREQ
FREQ
[§3.1]
CHAR_CRIT
[§3.1]
OPTION “SEPARATE”
“AJUSTE”
[§2.1]
“AJUSTE”
[§3.1]
“PROCHE”
[§3.1]
NMAX_FREQ
0
[§3.1]
NMAX_ITER_SEPARE
30
[§3.2.1]
PREC_SEPARE
1.E-04
[§3.2.1]
NMAX_ITER_AJUSTE
15
[§3.2.2]
PREC_AJUSTE
1.E-04
[§3.2.2]
NPREC_SOLVEUR
8
[§2.6]
NMAX_ITER_SHIFT
5
[§2.6]
PREC_SHIFT
0.05
[§2.6]
SEUIL_FREQ
1.E-02
[§2.9]
CALC_MODE
OPTION “DIRECT”
“DIRECT”
[§3.3.3]
“RAYLEIGH”
[§3.3.3]
PREC
1.E-05
[§3.3.3]
NMAX_ITER
30
[§3.3.3]
VERI_MODE
STOP_ERREUR “YES”
“OUI”
[§2.9]
“NON”
SEUIL
1.E-02
[§2.9]
Table 3.3.5-a: Récapitulatif of the parameter setting of MODE_ITER_INV
Note:
· in V5, the user can specify the class of membership of his calculation by initializing it
key word TYPE_RESU. According to this value, one informs vector FREQ or CHAR_CRIT,
· one finds all the “tripaille” of parameters related to the preprocessings of the test of Sturm
(NPREC_SOLVEUR, NMAX_ITER_SHIFT, PREC_SHIFT) and with postprocessings of checking
(SEUIL_FREQ, VERI_MODE),
· at the time of the first passages, it is strongly advised to modify only them
principal parameters noted in fat. The others relate to more the mysteries of
the algorithm and they were initialized empirically with values standards.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
31/78
4
Method of subspace (MODE_ITER_SIMULT)
4.1 Introduction
If one only wishes to calculate the p clean elements of a generalized problem of command N where
N << p (for example, p smaller eigenvalues or all eigenvalues included/understood in
a given interval), one saw that it was preferable to have resorts to methods of subspace. They
are based on the analysis of Rayleigh-Ritz which consists in effectively projecting the problem considered
on under space of dimension m (p < m < N) and to seek certain clean elements of it
problem projected (much easier to treat) using robust algorithms (QR or QL for
Lanczos and IRAM, Jacobi for the method of Bathe and Wilson).
The criteria of effectiveness of the aforesaid projection are:
· small size of the space of projection (directly related to calculation complexities and memory),
· facility of its construction,
· the robustness of orthogonal projection (in the nonsquare case, of oblique projections
were developed by F. Chatelin [bib2],
· the setting in the canonical shape of the projected matrix,
· the good approximation of the part of the initial spectrum sought by that of the projected operator,
· and, of course, the minimization of calculation complexities and memory and that of the effects of round-offs
(those are especially related to the problems of orthogonality raised by the 3rd point).
One will first of all present the analysis of Rayleigh-Ritz before detailing three methods resulting from this
analyze: method of Lanczos, that known as of Sorensen (IRA) and that of Bathe and Wilson.
4.2
Analyze of Rayleigh-Ritz
Let us consider the standard modal problem of command N, A U = U, and the subspace H of N
generated by a orthonormée base (Q, Q,…, Q
1
2
m). The latter constitutes the matrix
orthogonal Q
T
m allowing to define the operator of projection P
= Q Q
m
m
Mr. method of
Galerkin used by this analysis consists in solving the following problem
m
= m
×
(
~
~
~
U Q X,
,
X
, ~u)
With
to find (
)
such as
To find
× H such as
*
~
Q WITH
Q X
m
m
= X
P
m (~ ~~
With - U) = 0
! “
# $
#
B
m
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
32/78
~
The search of the elements of Ritz (, ~
U) which thus interest cost to us to find the modes
clean of the matrix of Rayleigh B = QT AQ
m
m
Mr. Les eigenvalues remain unchanged in
two formulations, on the other hand knowing a clean vector X of Bm one must go up with space all
entirety via the transformation ~
U = Q X
m
. The step can thus be summarized in the form:
· Choice of a space H and a base (h1, H2,…, hm) represented per H.
· Orthonormalisation of the base
R
0
· Stamp of Rayleigh
N
X
H
=
Qm
m
~
· Clean elements of Bm: (, X)
~
· Elements of Ritz: (, ~
U = Q X
m),
~
· Test of error with the residue: R = Au~ - u~
.
Algorithm 9: Procedure of Rayleigh-Ritz
Note:
· the elements of Ritz are in fact the clean modes of the matrix (of command NR) of
the approximation of Galerkin Cm = Pm A Pm,
· calculation complexity of this process, the command at worst of O (m2 (4n+m)) (much less
with a good space, cf Lanczos and IRAM), are without common measurement with that of a good
QR (O (N2)) or that of an iteration of the quotient of Rayleigh encapsulated in a process
heuristics (the cost of factorizations is prevalent on that of the products matrix-vector
algorithm itself) such as that developed in MODE_ITER_INV (>> pn3). But
it is quite clear that these methods are more complementary than concurrent.
To consider the error made by using the clean elements of the Bm matrix there is the result
according to.
Theorem 6
That is to say (, U) a clean element of the matrix diagonalisable (it is the case of the normal matrices and of
nondefective matrices (matrices of which all the eigenvalues have even arithmetic multiplicity and
~
geometrical) A and, Bm the matrix of Rayleigh associated, then this one admits an eigenvalue
checking
, U
(-
~
I P
(I - P
m) U
m) U
-
2
P U 2
Pmu
m
2
H
~ ~
where is a number dependant on, A and P
, U
Mr.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
33/78
In the square case the numerator of this increase is squared. The associated clean vector, ~
U,
check as for him
(I - P) U
(
sin, ~
U U)
m
2.
P U
m
2
Proof:
Cf [bib11] pp531-537.
It is shown that for m rather large, can be raised by a number which does not depend that on A and on
min I -. The estimate of the second member (I - Pm) U depends on the choice of the subspace.
J
2
I J
Note:
The number min I - (and its alternatives) is omnipresent in the analyzes of errors, of
J
I J
convergence or of conditioning spectral. One sees, once more, the importance of
separation of spectrum of work on the good numerical behavior of the algorithms.
Of course, to take account of spectral information already obtained, even to improve it or filter it,
one reiterates this procedure of Rayleigh-Ritz by modifying the subspace of projection. One will connect
then the construction of this new space (recurring of the precedent) and this process of projection. Two
types of spaces (each one attached to distinct methods) are most often used.
4.3
Choice of the space of projection
Two choices of space of projection H are most often set up:
· The first, that of the method of Bathe and Wilson (cf [§7] METHODE = “JACOBI”),
consist in choosing under space H of dimension m then to build successively:
H = A H
1
H = A H = A2 H
2
1
…
H
I
I = A H I - = A H
1
This method, which holds at the same time of generalization in form block of the method of
powers and of the truncation of the algorithm QR, conduit with an impoverishment of space
of work: dim H dim H
I
i-1. Not to always find the same clean mode
dominating it is necessary to insert in the process a reorthogonalisation (with all the problems
of calculation complexity and round-offs that that implies).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
34/78
· The second, that of the method of Lanczos (cf [§5], METHODE = “TRI_DIAG”) and of IRAM
(cf [§6], METHODE = “SORENSEN”) starting from an initial vector y consists, to build
continuation of spaces of Krylov (Hi) I = 1, m where Hi is the space generated by the family of vectors
(y, Ay,…, Have-1y). This last is called the space of Krylov of A and command I initiate by y:
H = K (A, y) Vect
i-1
I
I
(y, Ay,…, A y).
They check the following property: dim H dim H
I
I +1. Contrary to the preceding choice, one
thus a certain enrichment of the workspace has. In addition, it will be seen that they
allow to project in an optimal way within the meaning of the criteria defined previously.
Historically, the first type of space was very much used in mechanics of the structures. But for
to decrease calculation complexity related on the size of the subspaces and the orthogonalizations, one prefers to them
from now on methods of the type Lanczos/Arnoldi.
Note:
One meets an impressive variety of algorithms using a space of the first type. They
are called “'iterations of subspace”, “iterations orthogonal” or “iterations
simultaneous ". Beyond the various terms, it should especially be retained that these adaptations
the strategies of restartings (technique concern consisting in starting again an algorithm
modal with an initial vector comprising spectral information already. Typically, one chooses
a linear combination of the clean vectors or already exhumed vectors of Schur
(cf [§5.4.2]), techniques of acceleration (technical consisting in replacing the matrix of work
With by a matric polynomial P (A) which with the characteristic to be of great amplitude in the areas
spectral of interest (cf [§6.4]) and of factorization implemented to enrich the workspace.
There are even versions based on the powers opposite making it possible to calculate p more
small modes (cf [bib11] pp538-45, [bib23] pp49-54).
For these methods of projection, by supposing that initial space H is not too poor according to
the p dominant directions, the factor of convergence of the ième mode Lorsqu' they are arranged
classically, i.e. by order descending of module (when they are classically arranged,
i.e. by order descending of module) at the end of K iterations is written:
K
O
m+
1
.
I
It expresses two phenomena clearly:
· the priority convergence of the dominant modes,
· improvement of these convergences (and thus of their standards of error for a number
iterations fixed) when the size of the subspace increases.
To transform the modal problem generalized into a standard problem one A resorts to
spectral transformations. They also make it possible to direct the search of the spectrum and to improve
convergence (cf [§2.7]).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
35/78
4.4
Choice of the spectral shift
To calculate the smallest eigenvalues close to a given frequency or all the values
clean included/understood in a given frequency band, one carries out a spectral shift of
generalized problem. That is to say the value of shift, one transforms initial problem A U = B U into
a shifted standard problem.
With B U
1 = µ U with A = A - B and µ = -
This transformation spectral, known as of “simple shift”, is used in the method of Bathe and Wilson.
As the process detects the smallest eigenvalues in µ, those gradually are captured
who are closest to (in module).
One carrying out the reverse transformation (“shift and invert”) it occurs
1
With U = U
µ with A = (A -
- 1
B) B
and µ = -
It is the problem dealt with with Lanczos and IRAM which makes it possible to calculate the greatest eigenvalues
µ thus the eigenvalues closest to the shift (in module).
In command MODE_ITER_SIMULT, there are three ways of choosing this shift (for buckling,
transposition at the levels of critical loads is immediate):
·
= 0, one seek the smallest eigenvalues of the starting problem. This corresponds to
OPTION = “PLUS_PETITE” under the key word factor CALC_FREQ.
2
·
= 0 with 0 = (
2 f0), one seeks the frequencies close to the frequency FREQ = f0.
(OPTION = “CENTER”).
0 +
2
2
·
=
1 with 0 = (
2 f0) and 1 = (
2 f1), one seeks all the frequencies
2
included/understood in the tape [F, F
0
1] (OPTION = “BANDE” with FREQ = {f0, f1}).
The number of frequencies to be calculated is given in general by the user using NMAX_FREQ under
the key word factor CALC_FREQ. But for option “BANDE”, it is automatically given in
using the corollary 5bis (cf [§2.6]).
Note:
It is pointed out that the factorization of the matrix of work, just like the tests of Sturm of the option
“BANDE” are dependant on numerical instabilities when the shift is close to an eigenvalue.
Palliative, skeletal processing by the user, were established (cf [§2.6], [§2.9]).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
36/78
5
Method of Lanczos (METHODE = “TRI_DIAG”)
5.1 Introduction
This algorithm, developped at the point by Lanczos [bib10] in 1950, was hardly used until the medium of
the Seventies. First of all very simple and effective, it is the seat nevertheless of great instabilities
numerical being able to cause the capture of phantom modes! Those are dependant mainly on
problems of maintenance of orthogonality enters the vectors generating the space of Krylov.
comprehension of this type of behavior vis-a-vis arithmetic finished computers was long with
to exhume.
Since, many palliative solutions were born and an abundant literature covers the subject.
Let us quote for example the work of J.K. Cullum [bib3] and Al which is entirely devoted to him, that of
B.N. Parlett [bib18] and the brought up to date and exhaustive synthesis of J.L. Vaudescal [bib23] (pp55-78).
In the continuation of this paragraph, we first of all will delay we on the theoretical framework of
operation of the algorithm. Then, before detailing the alternative installation in Code_Aster,
we will stick to the realistic framework of arithmetic finished.
5.2
Theoretical algorithm of Lanczos
5.2.1 Principle
Its perimeter of application covers the couples operator with work (pseudo) produced scalar ensuring
the hermiticity of A. It makes it possible to repeatedly build a matrix of Rayleigh Bm of size
flexible, tridiagonale and square (with a true scalar product, if not it loses this
last property). This particular form largely facilitates the calculation of the modes of Ritz: with QR
one loses a command busy magnitude then, when one seeks p clean modes while projecting on one
subspace of size m, O (pm2) with O (20pm) [bib6].
The algorithm consists in gradually building a family of vectors of Lanczos q1, q2,…, qm
while projecting orthogonally, with the iteration K, vector A qk on the two preceding vectors qk and
qk-1. The new vector becomes qk+1 and thus, gradually, one ensures structurally
the orthogonality of this family of vectors. The iterative process is summarized thus.
Q = 0
0
, 0 = 0.
Calculation of Q/Q
1
1 = 1.
For K = 1, m to make
Z = A Q
-
K
Q
,
k-1
K - 1
=
K
(Z, qk),
v = Z - Q,
K
K
= v,
K
If K 0
then
v
Q
=
,
K +1
K
If not
Déflatio;
N
End if.
Fine loops.
Algorithm 10: Theoretical Lanczos
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
37/78
While noting, EM the mième vector of the canonical base, the vector residue of the factorization of
Lanczos is written R
= Q
and
m
m
m+1 Mr.
B
m
m
N
With
Qm
=
Qm
+
0
m
N
m
R
T
m=m qm+1em
Appear 5.2.1-a: Factorization de Lanczos
The matrix of Rayleigh takes the form then (in complex with a square scalar product)
1
0
0
1
1
…
0
B
2
m =
.
0
…
…
m
- 1
0
0
m
- 1 m
By construction, this algorithm ensures us of the following results (into arithmetic exact):
·
(qk) k=1, m constitute an orthonormal family,
· they generate the space of Krylov of command m initiated by q1,
K (A, Q
- 1
1) = Vect
m
m
(Q, A Q,…, A Q
1
1
1),
· they make it possible to build a matrix Bm tridiagonale, square and of flexible size,
· the spectrum of Bm does not admit that the m dominant simple modes on which q1 has one
component.
Note:
· the algorithm thus does not allow, in theory, the capture of multiple modes. That can
to be explained by noticing that any matrix of irreducible Hessenberg (a matrix A is
said of Hessenberg higher (resp. lower) if Have, j= 0 for i>j+1 (resp. j>i+1). It is
of more irreducible if Have 1, I
0 I
+
) (thus a fortiori a matrix tridiagonale) does not admit that
simple modes,
· the cost of an iteration is weak, about that of a product matrix-vector A qk, that is to say
for the m first iterations a calculation complexity of O (Nm (c+5)) (with C the number
means of nonnull terms on the lines of the matrix of work). In addition, complexity
necessary memory is weak bus one does not need to build A in entirety, one has right need
to know its action on a vector. This characteristic is very interesting for
to solve problems of big sizes [bib23],
· in general the sphere of activity directs the choice of a vector of initial Lanczos, this one must by
example to belong to a space of acceptable solutions (checking constraints,
conditions limit…) like with the image unit of the operator. This last point is important
because it makes it possible more quickly to enrich modal search while not being limited to the core,
· the algorithm of Lanczos is only one means of approaching the subspaces of Krylov. Its
generalization with the nonsymmetrical configurations is called algorithm of Arnoldi (cf [§6.2]).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
38/78
5.2.2 Estimates of errors and convergences
Because of particular form of Bm, the extraction of the required elements of Ritz is simplified, and in
more, the following result allows us quickly (via a product of two realities!) to estimate theirs
qualities.
Property 7
~
The euclidian norm of the residue of the element of Ritz (, ~
U = Q X
m) is equal to
~
R
= (A
~
T
- I) U =
E X
2
2
m
m
Proof:
Commonplace by taking the euclidian norm of the factorization of Lanczos
WITH Q X
Q B X
Q
and
=
+
X
m
m
m
m m+1 m where qm+1 is normalized with the unit.
Note:
A more general result, independent of any method of resolution, us ensures that for
~
each clean mode (, X) of Bm, there exists a clean mode (, U) of A such as:
R 2
- ~
2
with = min
T
I - ~
~
~
U With U
(
R
sin U, ~u)
I
2
Therefore, more than minimization of the residue, the criterion of stop of the methods of the Lanczos type/Arnoldi
m < makes it possible to approach the original spectrum as well as possible. These increases are not accessible
that in the square case.
In the general case (and in particular nonnormal) they are more difficult even impossible with
to build without additional information (level of nonnormality, “defectivity”.).
The estimate of the residue is not then any more the good criterion to estimate the quality of the approximated modes.
Let us interest being maintained in the quality of convergence of the algorithm. Let us particularize theorem 6
for this algorithm. By limiting the standard of error of (I - Pm) U one obtains increases
2
following:
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
39/78
Theorem 8
That is to say (I, ui) the ième clean mode dominating of A diagonalisable, there exists a mode of Ritz then
(~, ~
I ui) such as:
2
~
2
J -
N
- 2
I - I
m
T - I (I)
j<i J - I
(
J - N
sin
, ~u
- 1
ui I)
m
T - I (I)
j<i J - I
-
I
i+1
with I = 1 + 2
, T
-
the semi (X) polynomial of Tchebycheff of semi degree, the constant of
1+i
N
theorem 6 and =
your (
N Q, U
1
I).
P U
m I 2
Proof:
Cf [bib11] pp554-557.
These increases indicate (cf [bib23] pp59-63) that:
· if the initial vector does not have any contribution along the clean vectors, one cannot capture
the aforementioned modes (+),
· one
has
firstly convergence of the extreme modes and as much better than the spectrum
is separate (these methods are less sensitive to the clusters than the methods of the types
reiterated powers), without packing of eigenvalues (famous “the clusters”),
· the error decrease when m increases (like the reverse of the polynomial Tm-I (X), therefore like
the reverse of exponential),
· the estimate on the eigenvalues is better than that their associated own vectors.
Note:
The convergence of the method was studied by P. Kaniel then improved (and corrected) by
S.C. Paige; One will find a synthesis of these studies in the paper of Y. Saad [bib20].
Let us look at now how behaves the algorithm into arithmetic finished. We will see that it
is the seat of surprising phenomena.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
40/78
5.3
Practical algorithm of Lanczos
5.3.1 Problem
of orthogonality
At the time of implementation the effective of this algorithm, one realizes that very quickly the vectors of
Lanczos lose their orthogonalities. The matrix of Rayleigh is not then any more the matrix projected of
the initial operator and the captured spectrum are sullied more and more with errors. A long time this phenomenon
inescapable was allotted wrongly to the effects of rounded the arithmetic one finished. In 1980, D.C. Paige showed
that the loss of orthogonality was especially ascribable with the convergence of a clean mode. Thus, as of
that one captures a mode dominating, one disturbs the layout of the preceding vectors of Lanczos.
following result expresses this paradox clearly.
Property 9 (Analyze de Paige)
By taking again the notations of algorithm 10 and by noting the precision machine, with the iteration k+1
the orthogonality of the new vector of Lanczos is expressed in the form:
vTq + A
I
qT Q
2
K 1 I
(I K
+
=
)
K
Proof:
Cf [bib26].
The problem occurs in fact during the standardization of the new vector of lanczos qk+1. When it
mode converges, according to property 7, that comes from the smallness of the coefficient K and thus in spite of
theoretical orthogonality (v T IQ = 0 (I K)), effective orthogonality is far from being acquired
(qT Q
K +1 I >> (I K)).
The digital processing of these problems has been the search object many for thirty years
and of many palliative strategies were developed. The choice of such or such method depends
type of required spectral information and average data processing available, the synthesis of
Besides JL.Vaudescal proposes a very good overflight of these alternatives (cf [bib23] pp66-78):
· Algorithm of Lanczos without reorthogonalisation, developed by J.K. Cullum and
R.A. Willoughby [bib3] who consists with expurger the calculated spectrum of his phantom modes of
studying interlacings of the eigenvalues of the matrix of work and one of its
submatrices. This alternative admits a weak overcost calculation and memory to determine them
eigenvalues, but it complexes (sometimes largely) the search of the vectors
clean.
· Algorithms of Lanczos with total reorthogonalisation (B.N. Parlett) or selective
(J. Scott) which consists with each step with réorthogonaliser the last vector of Lanczos obtained
compared to all the already calculated vectors or simply compared to the vectors of Ritz
converged (this alternative thus makes it possible to control the loss of orthogonality dynamically
acceptable). These alternatives are much more expensive in complexity calculation (them
automatic réorthogonalisations) and memories (storage of the preceding vectors) but
they prove more effective and more robust.
In the alternative of Newmann-Pipano [bib14] (METHOD = “TRI_DIAG”) of Code_Aster, it is
strategy of reorthogonalisation supplements which was selected: the vector qk+1 is thus not completely
calculated as in the theoretical algorithm.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
41/78
Note:
· these problems of loss of orthogonality occur with all the more of acuity that the size of
subspace increases. It is an additional argument for the installation of one
iterative process limiting this size (cf [§5.4.2]),
· these strategies are declined vectoriellement or per blocks, they are established in
simple or iterative algorithms. The explicit latter being able to profit from restarts or
implicit, controlled or automatic. Alternative IRAM profits thus from an implicit restart
calculated automatically,
· the central point of these algorithms is consisted the method of orthogonalization put in
place. All the process is dependant on its success and its robustness. With
orthogonal of type Housholder or Givens, very expensive transformations but very
robust, one prefers from now on the algorithms of the type Gram-Schmidt Itératifs (IGSM) which are
a better compromise between effectiveness and complexity calculation (cf [Annexe 2]). It is it besides
choice which was made for the alternatives of Lanczos/Arnoldi installation in the code,
· the alternative of selective reorthogonalisation compared to the converged modes returns to
to carry out an implicit deflation by restriction on the operator of work not to have with
to recompute (cf [§3.1]).
5.3.2 Capture multiplicities
It was seen that in theory the algorithm of Lanczos did not allow, some is his strategy of
reorthogonalisation, the theoretical capture of multiple modes. In practice, for once, effects
from round-off come to the rescue and “powder” with small components along almost
all clean vectors. One can thus capture multiplicities, however they can be
erroneous and require a complementary postprocessing of checking.
Note:
In fact, only a version per blocks (G. Golub & R. Underwood, 1979) can allow us one
correct detection of the multiplicities, at least as long as the size of the blocks is sufficient. This
version was wide (Mr. Sadkane [bib22], 1993) with the algorithm of Arnoldi but it would be a pity
to deprive itself of the alternative of Sorensen (IRAM) who is more effective and more robust.
5.3.3 Phenomenon of Lanczos
The success of this algorithm rested at the beginning on what is called the “phenomenon of Lanczos”.
This conjecture predicts that for a sufficiently large size of subspace (m >> N) one is
able to detect all the spectrum (with load to thereafter distinguish the “good grain from the ryegrass”) of
the operator of work. Taking into account weak the pre-necessary report of “basic” Lanczos
(grosso-modo a matrix tridiagonale and some vectors), this is particularly interesting for
to treat hollow systems of very big sizes (other algorithms of reiterated the powers type or
QR require them it knowledge of all the matrix of work).
We now will detail (a little) some elements whose interest largely exceeds the framework
of this algorithm.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
42/78
5.4 Processing
complementary
5.4.1 Detection of spaces invariants
In algorithm 10, the nullity of the coefficient l-1 prevents the standardization of the new vector and
an additional processing called of deflation requires. The subspace of calculated Krylov is then under
space invariant and its spectrum is thus exact. In other words, the l-1 first eigenvalues
exhumed by the algorithm of support (QR or QL) are those of A.
To continue, it is then necessary to choose a new random vector observing the limiting conditions,
the orthogonaliser with all the already calculated vectors of Lanczos and to build a second family of
vectors of Lanczos that one orthogonalise compared to the first family constituting space:
H = K (A, Q
1
1) = Vect
L -
L
L
(Q, A Q,…, A Q
1
1
1).
The matrix tridiagonale obtained with the following form:
1
1
…
0
1
2
…
…
0
Bm =
0 L
L
L L
+1
…
0
…
…
M-1
M-1 m
The intermediate factorization of Lanczos (with the command L) is written:
WITH Q = Q B
L
L
L
The values (K, K) k=1, L are obtained starting from the first random vector and the values
(K, K) k=l+1, m are obtained starting from a second random vector. To detect the nullity of the term
extra-diagonal the numerical criterion is used
L
PREC_ LANCZOS
- 1
L then l-1 = 0,
where PREC_LANCZOS is initialized under the key word factor CALC_FREQ.
Note:
· a more robust criterion retained by large mathematical bookstores EISPACK [bib26],
LAPACK [bib26] and ARPACK [bib27] use the preceding diagonal term, a parameter
flexible C and the precision machine,
L 1 < C
-
(l-1 + L)
it is this criterion which at summer appointed for IRAM (but the user cannot modify the parameter C,
this one is fixed at a standard value by the code), various messages preventing besides
the user of the detection of a space invariant.
· Obtaining such a space is completely sympathetic nerve since it ensures us of very
good quality of the first modes. Moreover, one reduces the size of the problem by dividing it into
two parts: one solved, the other to solve. Many techniques seek with
to reproduce artificially this phenomenon (cf [§3.1], [§5.3.1]).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
43/78
5.4.2 Strategies of restartings
To limit the endemic problems of orthogonalization, to limit to them pre-necessary report and
to take dynamically into account spectral information already obtained, of the strategies of
restartings were coupled with the algorithm of Lanczos. It loses its “simple” character and becomes
“iterative”. One reiterates the “restarts” until satisfying the desired criteria of convergence. The algorithm
does not undergo the convergence of the modes imposed by theorem 8 and it can support intéractivement
that of certain modes.
In theory however, plus the size of the subspace is large, better is convergence. One
compromise is thus to find between the size of the subspace and the frequency of the restartings.
Various vectors of restartings can be used, being generally written like a sum
balanced p required clean modes (Y. Saad 1980)
p
Q = Q X
1
I
m I
.
I =1
Indeed, if one starts again with an initial vector pertaining to the subspace invariant generated by
sought modes, one is then sure to obtain (with a good algorithm of orthogonalization) a standard
residual about the precision machine and of the almost exact clean modes. In addition to the decision
to start the restarting, all the problems lie in the search for these weights I. Us
will see that with IRAM, the restarts set up via implicit polynomial filters solve
elegantly and automatically this question.
Note:
· this philosophy of the restarts was initiated by W. Karush (1951),
· restartings based on the vectors of Schur associated with the spectral decomposition
sought seem more stable (J. Scott 1993),
· criteria of effectiveness were required to decide restart appropriateness (B. Vital
1990),
· instead of a simple linear combination of clean vectors, one can determine a vector
of restarting via polynomials (Tchebycheff in reality and Faber in complex) allowing
to center modal search in such or such zone. One speaks then about polynomial accelerations
explicit [bib22].
The following paragraph will show some of the concepts presented hitherto to clarify
alternative installation in the code.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
44/78
5.5
Establishment in Code_Aster
5.5.1 Alternative of Newmann & Pipano
This alternative developed by Mr. NEWMAN & A. PIPANO [bib3], [bib14], in 1977, is a method of
Simple Lanczos, into arithmetic real, with total reorthogonalisation (via a GSM). It uses it
crossbred “shift and invert” traditional of the scalar pseudo-product introduced by shifted matrix (A-B).
(A)
B -
1
1
B U =
U
!# “
#
$
# #
µ -
With
!“$
(X, y) = yt (A)
B X,
with =
(
sign X, X), 1
X = (X, X) 2.
This choice returns the couple operator of work - symmetrical scalar product and it are adapted to the matrices
particular of the mechanics of the structures. One can thus deal with problems:
· of free structure and fluid structure (A can be singular),
· of buckling (B is indefinite).
The scalar pseudo-product introduced by the shifted matrix is thus regular (answers this
prerogative cf [§2.6], [§2.9]) and it makes it possible to seek vectors of Lanczos (A - B)-orthogonal.
The total gonalisation is carried out only if it proves to be necessary and this criterion is skeletal.
Note:
· B-orthogonality and a fortiori that with the Euclidean direction cannot be but the fruit any more of
particular configurations. This does not modify the properties of orthogonalities of the modes
clean (cf proposal 2),
· as one already announced [§1.7] there is a whole zoology of couples operator - product
scalar, this one being only one possibility among others. Thus, the libraries [bib29]
traditional propose to deal with the problems of buckling in “buckling mode” via the same one
“shift and invert” and the scalar pseudo-product introduced by A. Mais because of introduction
quasi-systematic of Lagranges, this matrix becomes indefinite even singular, which
disturb the process largely. The same causes produce the same effects when
for a calculation of dynamics, one uses the scalar B-product.
The price to be paid for this gain in robustness is the loss of possible symmetry of the matrix of Rayleigh
1 1
0
0
1
…
0
=
B
2
I
I I
m =
with
.
0
…
…
m 1
-
I = sign (Q, Q
I +1
I +1)
0
0
M-1 m
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
45/78
After balancebeing balanced and putbeing put in the form of Hessenberg higher, Bm is diagonalisée by one
implicit method QL if it remains in spite of very symmetrical or by a method QR if not
(cf [Annexe 1]). The difference in cost calculation between these solveurs remains negligible vis-a-vis the costs of
réorthogonalisations. The diagram of established Lanczos becomes then:
Random hard copy of qinit.
q~
- 1
1 = (A -
)
B
(qinit (I) .ulagr (I)) .i
Q % = A
1
(q~1 (I) .ubloq (I)) .i
1 =
(
sign qT
% 1 (A -)
B Q % 1)
1
-
Q = Q
2
1
1% qT
1% 1 (A -
)
B Q % 1
Q = 0, =,
0 =.
0
0
0
0
For K = 1, m to make
Z = A Q -
Q
,
K
k-1
k-1
k-1
= T
K
qk (A -)
B Z,
v = Z - Q,
K
K
K
Réothogonalisation from
v
report/ratio with Q
(I)
(IGSM),
I =1, K
T
K = v (A -)
B v,
T
K =
(
sign v (A -)
B v),
If K
then
0 Kv
qk+1 =
,
1
K 2
If not
Deflation;
End if.
Fine loops.
Algorithm 11: Alternative of Newman-Pipano
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
46/78
The reorthogonalisation is carried out thanks to an alternative “house” of Méthode of Gram-Schmidt
Iterative (IGSM) according to the following process:
For I = 1, K
has
T
= qk+1 (A - B) Q,
I
If (qT 1 (A)
B Q
K
I
PREC_ORTHO
+
-
) then
For J = 1, NMAX_ITER_ORTHO
X = Q
T
K +1 - (qk +1 (A -)
B IQ) IQ,
B = xT (A -)
B IQ
If (B PREC_ORTHO) then
Q
X
K +1 =
,
I + 1 I
, exit;
If not
If (B has) then
Q
X
K +1 =
,
= B has;
If not
Failure, emission of a message of alarm,
I + 1 I
, exit;
End if.
End if.
Fine loops in J.
End if.
Fine loops out of I.
Algorithm 12: Procedure of reorthogonalisation of “TRI_DIAG”
Note:
After some tests, it seems that this alternative specific to the code is less effective than
the IGSM of the type Kahan-Parlett [bib18] selected for the IRAM (cf [Annexe 2]).
5.5.2 Parameter setting
To be able to activate this method, it is necessary to initialize key word METHODE with “TRI_DIAG”. In addition,
the size of the subspace of projection is determined, either by the user, or empirically to leave
formula:
m = min (my (
X 4 p, p + 7), nddl-credits)
where
p
is the number of eigenvalues to calculate,
nddl-credits
is the number of active degrees of freedom of the system (cf [§2.2])
The user can always impose to him even dimension by indicating it with the key word
DIM_SOUS_ESPACE of the key word factor CALC_FREQ.
Parameters of the total reorthogonalisation, PREC_ORTHO and NMAX_ITER_ORTHO, the number
maximum of iterations QR, NMAX_ITER_QR, and the criterion of deflation (cf [§5.4.1]), PREC_LANCZOS, are
accessible by the user under the key word factor CALC_FREQ. When this phenomenon of deflation
product, a specific message specifies its row L.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
47/78
5.5.3 Warning on the quality of the modes
Let us recall that this alternative is not iterative. From the number of desired frequencies, one
estimate the dimension of the subspace of calculation and one hopes that the clean modes of the problem
projected standard will be a good approximation of those of the generalized problem. One checks has
posteriori the validity of the results (cf [§2.9]) but one does not have any active control on this precision.
If they are not satisfactory, one has as only only resorts to carry out a new test in
increasing the dimension of the subspace.
Note:
· method IRA whom we will approach proposes a flexible and dynamic control of
precision of the results, it is an iterative method,
· in V5, in parallel of the modal position and eigenvalue (in frequency if one is in
dynamics), one displays from now on the standard of the residue of the initial problem calculated in
postprocessing (after having standardized the clean vectors ad infinitum (in order to exhume
clean vectors “neutral” with respect to the string of standardizations which they can undergo at the time
calculations of the modal parameters of the structure (factors of participation…)).
5.5.4 Perimeter
of use
Attention, this method is usable one to the only deal with quadratic problem. In the event of error
at the time of the initialization of key word METHODE, calculation stops and a trace is left in the files of
output.
In V5, an algebraic method calculating the null eigenvalues of the modal problem
generalized was introduced. Physically those correspond to the movements with energy of
null deformation of a free structure (cf [bib25] TP n°2). Except the numerical difficulties of
handling of quasi-null quantities, because of their multiplicities, their capture correct was until
often problematic present for Lanczos established in Code_Aster: phantom modes
appeared corresponding to missed multiplicities!
This algorithm of detection of the modes of rigid body, which one activates by initializing OPTION with
“MODE_RIGIDE” (default value `SANS `), intervenes in preprocessing of modal calculation properly
known as. It is based on the analysis of the matrix of rigidity and breaks up into three phases:
· detection of the null pivots of this matrix,
· blocking of these pivots,
· resolution of a linear system whose are solutions the associated clean vectors.
During the process of Lanczos it is enough consequently to orthogonaliser, progressively with their
determination, basic vectors with the latter.
However, the introduction of method IRA reduces the interest (if it is not as comparison) of one
such option (which is not free in calculation complexity since it requires inversions of systems).
What good is it to deprive itself of such an algorithm which is by no means affected by the presence of these modes
a little particular multiples!
This method remains nevertheless useful as solution of help in the event of failure of IRAM (that
can occur in the presence of clusters in the vicinity of the terminals or for operators
pathological (badly conditioned, nonnormal, defective…) who should however be rather rare
(that often reveals a badly posed problem). One can also plan to encapsulate it in one
auxiliary operator (like IMPR_STURM).
Note:
The quadratic processing of problem is incompatible with options “MODE_FLAMB” and
“MODE_RIGIDE”.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
48/78
5.5.5 Display in the file message
The example below resulting from the list of case tests of the code (sdll112a) recapitulates the whole of the traces
managed by the algorithm. The iteration count effective QR (or QL) can be only identical for
all eigenvalues. It is an artefact of information which will be brought to disappear.
------------------------------------------------------------------------
THE NUMBER OF DDL
TOTAL EAST: 86
LAGRANGE EAST: 20
THE NUMBER OF ACTIVE DDL EAST: 56
------------------------------------------------------------------------
The SELECTED OPTION EAST: PLUS_PETITE
THE VALUE OF SHIFT IN FREQUENCY EAST: 0.00000E+00
------------------------------------------------------------------------
INFORMATION ON CALCULATION REQUIRES:
A NUMBER OF REQUESTS MODES: 10
THE DIMENSION OF REDUCED SPACE EAST: 0
IT IS LOWER THAN THE NUMBER OF MODES, ONE TAKES IT EQUALIZES A 40
------------------------------------------------------------------------
FREQUENCIES CALCULEES INF. AND SUP. SONT:
FREQ_INF: 1.54569E+01
FREQ_SUP: 1.01614E+02
THE FIRST HIGHER FREQUENCY NOT SELECTED EAST: 1.29375E+02
------------------------------------------------------------------------
MODAL CALCULATION: METHOD Of SIMULTANEOUS ITERATION
METHOD OF LANCZOS
NUMBER FREQUENCY (HZ) STANDARD Of ERROR ITER_QR
1 1.54569E+01 1.29798E-12 4
2 1.54569E+01 7.15318E-13 4
3 3.35823E+01 3.98618E-13 4
4 3.35823E+01 4.25490E-12 4
5 4.73076E+01 4.26463E-12 4
6 4.73076E+01 1.43391E-12 4
7 5.45850E+01 9.52006E-12 4
8 8.80156E+01 3.45489E-13 4
9 1.01614E+02 6.13949E-12 4
10 1.01614E+02 3.18663E-12 4
------------------------------------------------------------------------
CHECKING A POSTERIORI OF THE MODES
IN the INTERVAL (1.54182E+01, 1.16326E+02)
IT THERE A WELL 10 FREQUENCY (S)
------------------------------------------------------------------------
Example 5: MODE_ITER_SIMULT with “TRI_DIAG”
Now let us recapitulate the parameter setting available of operator MODE_ITER_SIMULT with this
option METHOD = “TRI_DIAG”.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
49/78
5.5.6 Summary of the parameter setting
Key word factor
Key word
Default value
References
“DYNAMIC” TYPE_RESU
“DYNAMIQUE”
[§2.1]
“MODE_FLAMB”
[§2.1]
METHOD “TRI_DIAG”
“SORENSEN”
[§5.5.3]
OPTION “MODE_RIGIDE”
“SANS”
[§5.5.4]
“SANS”
[§5.5.4]
CALC_FREQ
FREQ
[§4.4]
CHAR_CRIT
[§4. 4]
OPTION “PLUS_PETITE”
`PLUS_PETITE ''
[§4.4]
“BANDE”
[§4.4]
“CENTER”
[§4.4]
NMAX_FREQ
10
[§4.4]
DIM_SOUS_ESPACE
Calculated
[§5.5.2]
PREC_ORTHO
1.E-12
[§5.5.1], [§5.5.2]
NMAX_ITER_ORTHO
1.E-04
[§5.5.1], [§5.5.2]
PREC_LANCZOS
1.E-04
[§5.5.1], [§5.5.2]
NMAX_ITER_QR
15
[§5.5.1], [§5.5.2]
NPREC_SOLVEUR
8
[§2.6]
NMAX_ITER_SHIFT
5
[§2.6]
PREC_SHIFT
0.05
[§2.6]
SEUIL_FREQ
1.E-02
[§2.9]
VERI_MODE
STOP_ERREUR “YES”
“OUI”
[§2.9]
“NON”
[§2.9]
PREC_SHIFT
5.E-03
[§2.9]
SEUIL
1.E-06
[§2.9]
STURM “YES”
“OUI”
[§2.9]
“NON”
[§2.9]
Table 5.5.6-a: Récapitulatif of the parameter setting of MODE_ITER_SIMULT with “TRI_DIAG”
Note:
· in V5, the user can specify the class of membership of his calculation by initializing it
key word TYPE_RESU. According to this value, one informs vector FREQ or CHAR_CRIT,
· one finds all the “tripaille” of parameters related to the preprocessings of the test of Sturm
(NPREC_SOLVEUR, NMAX_ITER_SHIFT, PREC_SHIFT) and with postprocessings of checking
(SEUIL_FREQ, VERI_MODE),
· at the time of the first passages, it is strongly advised to modify only them
principal parameters noted in fat. The others relate to more the mysteries of
the algorithm and they were initialized empirically with values standards,
· in particular, to improve quality of a mode, the only flexible parameter is
dimension of the subspace, DIM_SOUS_ESPACE.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
50/78
6
Algorithm WILL GO (METHOD = “SORENSEN”)
6.1 Introduction
We saw that one of the crucial problems of the method of Lanczos is the loss of orthogonality
inescapable of its basic vectors (cf [§5.3.1]). A generalization of this algorithm to the case not
square imagined by W.E. Arnoldi [bib31] in 1951 makes it possible to solve partially this
problems. It was given to the style of the day by Y. Saad [bib32] in 1980 and one abundant
literature covers the subject. Put aside the article founder and papers of Y. Saad and Mr. Sadkane [bib22],
one recommends the brought up to date and exhaustive synthesis of J.L. Vaudescal [bib23] (pp79-112).
This method being at the base of algorithm IRAM known as “of Sorensen” (IRAM for `Implicit
Restarted Arnoldi Method'), we first of all will detail his operation, its
behaviors and its limitations. Thereafter, we will determine the stakes to which IRAM must answer
(and it does it in the majority of the standard cases!) and we will consider its theoretical mysteries
and numerical. We will conclude by the summary from its effective parameter setting in Code_Aster and
for an example of file message.
6.2 Algorithm
of Arnoldi
6.2.1 Principle
Its perimeter of application covers all the couples operator with work (pseudo) produced
scalar. During this opening is the filling of the matrix of Rayleigh which becomes
form Hessenberg higher. It is not very prejudicial, because one can thus apply to him
directly algorithm QR (the first stage of a good QR, except balancing, consists in reducing
the matrix of work in the form of Hessenberg. That makes it possible to gain a command magnitude at the time of
the resolution itself (cf [Annexe 1]), from where a gain about O (10m3/3).
The algorithm is very similar to that of Lanczos, it consists in building a family gradually
of vector of Arnoldi q1, q2,…, qm while projecting orthogonally, with the iteration K, vector A qk on
K vectors precedent. The new vector becomes qk+1 and thus, gradually, one ensures
the orthogonality of this family of vectors.
With the difference of Lanczos, the orthogonality of the new vector compared to all the precedents is
thus explicitly and not assured implicitly. This one is managed by the algorithm of Gram-Schmidt
Modified (GSM) (cf appendix 2) which proves sufficiently robust in the majority of the cases.
While noting, EM the mième vector of the canonical base, the vector residue of the factorization of Arnoldi
R is written
= B
Q
and
m
m+1, m
m+1 Mr.
B
m
m
N
With
Qm
=
Qm
+
0
m
R
T
N
m
m=bm+1, m qm+1em
Appear 6.2.1-a: Factorization d' Arnoldi
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
51/78
The iterative process is summarized as follows:
Calculation of Q/Q
1
1 = 1.
For K = 1, m to make
Z = A Q,
K
For L = 1, K to make (GSM)
B =
lk
(Z, ql),
Z = Z - blk ql,
Fine loops;
B +
= Z,
K 1, K
If B +
0 then
K 1, K
v
Q
=
K +1
,
bk+1, K
If not
Deflation;
End if.
Fine loops.
Algorithm 13: Theoretical Arnoldi
The matrix of Rayleigh is written then:
B
B
…
B
11
12
m
1
B
B
…
B
B
21
22
2m
m =
.
0
…
…
bm-1, m
0
0
B
m, m
B
- 1
mm
Except the shape of this matrix and a less acuity with the problems of orthogonality, this algorithm
us ensures the same theoretical and numerical properties that Lanczos. However if one
let grow indefinitely the size of the subspace until convergence of the eigenvalues
wished, the effects of round-off despite everything will take again the top and one goes to the front of large troubles.
From where need, as for Lanczos, to make this process iterative via restartings.
Note:
· this method can be seen like an implicit alternative of the algorithm of Lanczos with
total reorthogonalisation (cf [§5.3.1]),
· the implicit orthogonalization of the algorithm can be led by more expensive algorithms
but more robust such as QR or IGSM. This is strongly necessary when the operator of
work presents a too strong defect of normality,
· because of structure of Bm, memory complexity is requested than with Lanczos, by
against complexity calculation remains of the same order of magnitude O (Nm (c+3+m)) (with C it
numbers average nonnull terms on the lines of the matrix of work),
· to improve this first point, Y. Saad [bib32] showed that the structure of Hessenberg
higher can see to cancel its extreme on-diagonals if one carries out only partially
reorthogonalisation,
· the vectorial algorithms having a tendency natural to miss by the multiplicities, one prefers
often to use a version blocks (Mr. Sadkane [bib22], 1993). But the size of those influence
on the quality of the results, for this reason one prefers the vectorial versions to them or
blocks of IRAM,
· the choice of the vector of initial Arnoldi is carried out same manner as for Lanczos.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
52/78
6.2.2 Estimates of errors and convergence
With regard to the evaluation of the quality of approximation of the clean modes obtained, there is one
as simple and effective criterion as for Lanczos.
Property 10
~
The euclidian norm of the residue of the element of Ritz (, ~
U = Q X
m) is equal to
R = (A
~
T
- ~ I) U = B
E X
2
m +
2
1, m
m
Proof:
Commonplace by taking the euclidian norm of the factorization of Arnoldi
WITH Q X = Q B X +
Q
and
B
X
m
m
m
m +1, m
m +1 m and as qm+1 is normalized with the unit.
Always by focusing us on the standard of the additional projector (I - Pm) U one can
2
to generalize the theorem of convergence 8 with the nonsquare case.
Theorem 11
That is to say (1, u1) the first (traditional arrangement, by order descending of module) clean mode
N
dominating of A diagonalisable and is Q =
U
1
K K the initial vector of Arnoldi broken up on
K =1
~
base clean vectors, it exists a mode of Ritz then (, ~
1 u1) such as:
~
2
m
1 - 1
1
1
N
with =
, the constant of theorem 6, =
K
and
P U
1
m I 2
K =1, K 1
1
1
m +
-
1
m +1
m
K -
1
1 =
-
. This result is declined in the same way on the other modes.
J = K =, K J
K
2
2
J
Proof:
By taking again to the demonstration of Y. Saad ([bib33] pp209-210) and the result of theorem 6.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
53/78
These increases very different from those obtained with Lanczos guide the same ones however
phenomena:
· if the initial vector does not have any contribution along the sought clean vectors, one cannot
to capture (I +),
· one
has
firstly convergence of the peripheral modes of the spectrum, and this, of as much
better than it is separate (property of m
I),
· the decrease of the error is proportional to the increase in m (property of m
I).
Note:
· when an eigenvalue is badly conditioned I + then it should be increased m so that
m
I decreases,
· similar results were exhumed in the case of a defective operator (cf Zia 94 [bib23]).
Extremely from this lesson, we now will recapitulate the stakes to which must answer
IRAM.
6.3
stakes
Algorithm IRA tries to bring an elegant remedy for the recurring problems raised by the others
approaches:
· minimization of the space of projection: it proposes with minima m > p+1 instead of the m= 4p
of Lanczos,
· optimal management of the overcosts of orthogonalization establishing a compromise enters the size of
subspace and the frequency of the restartings,
· transparent, dynamic and effective management of these restarts,
· taking into automatic account of spectral information,
· fixing
pre-necessary report and of the quality of the results.
There is thus more question to be posed concerning the strategy of reorthogonalisation, the frequency of
restarts, their establishments, criteria of detection of possible phantom modes… “super-IRAM”
charge of all!
In short, it gets:
· one
better
total robustness,
·
complexities calculations O (4nm (m-p)) and memories O (2np+p2) improved (especially by
report/ratio in simple Lanczos such that of Newmann & Pipano) for a fixed precision,
· a more rigorous capture of the multiplicities, clusters and rigid modes of body
(thus less parasitic modes!).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
54/78
Note:
· on this last point, only a version per blocks of IRAM or a version incorporating of the “purging
and lock " (techniques of capture and filtering, cf D. Sorensen & R.B. Lehoucq, 1997) can
us to guarantee a correct detection of the spectrum of a standard operator (i.e not too badly
conditioned),
· it seems, that in practice in Code_Aster, the report/ratio in calculation complexity between IRAM and
Lanczos/Bathe & Wilson are with minima of command 2 in favor of the first. With sizes of
reasonable problems (a few tens of thousands of ddls and p= O (100)) this one can
to go up up to 10 (without the encapsulation of MACRO_MODE_MECA). In certain cases
semi-industrial, it made it possible to unroll a search for spectrum which had failed with
Jacobi and who were inaccessible with Lanczos (taking into account the time limits),
· a class of algorithm known as of “Jacobi-Davidson” (cf R.B. Morgan 1990) seems even more
promising to treat pathogenic cases. It uses an algorithm of the Lanczos type that it
packaged via a method of Rayleigh.
In the following paragraph we will clarify the operation of IRAM.
6.4
Algorithm “Implicit Restarted Arnoldi” (WILL GO)
This algorithm was initiated by D.C. Sorensen [bib30] in 1992 and makes real great strides for the resolution
great modal systems on parallel supercomputers. Its framework of application is very with
general fact. It deals with as well the real problems as complex, square or not. It is summarized in
a succession of factorizations of Arnoldi whose results control automatically
static and implicit restartings, via polynomial filters modelled by implicit QR.
First of all it carries out a factorization of Arnoldi of command m= p+q (in theory Q = 2 is enough, in
practical q= p is preferable. It is this last default value besides which was retained
(cf [§6.5.3]) of the matrix of work. Then once this preprocessing carried out it reiterates a process of
filtering of the part of the undesirable spectrum (numerically and méthodologiquement, it is easier
to exclude to incorporate).
It starts by determining the spectrum of the matrix of Rayleigh (via indétronable QR) and it in
dissociate the nondesired part (while referring to the criteria fixed by the user) which it uses then
like shift to set up a series of Q QR implicit with simple shift (cf [Annexe 1]).
factorization is written then:
In Q+
Q+ B+
=
+ R Q
m
m
m
m
where Q+ = Q Q
+
T
m
m
, B = Q B Q
m
m
and Q = Q Q
Q
1 2.
Q the unit matrix associated the QR. Afterwards
to have updated the matrices, one truncates them until the command p
In Q+ = Q+ B+ + R+
p
p
p
p
and thus, at the price of Q news iterations of Arnoldi, one can find a factorization of Arnoldi of command
m which is viable. All the subtlety of the process rests on this last sequence. Let us note:
Q
(
~
WITH) = (A -
Id
p+i),
i=1
the matric polynomial of command Q generated by the operator of work. In fact, the implicit QR have acts
on the p first lines of the matrix of Arnoldi, Rayleigh and residue, so that
factorization complementary to command Q produces the same effect as a factorization of command m
initiated by the vector
~
Q = (
With) Q
1
1.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
55/78
The initial vector was implicitly modified (it does not have construction and effective application there of
polynomial) so that it generates the desired modes preferentially and this, by withdrawing them
components considered to be “unsuitable”. As one already pointed out (cf [§5.4.2]) this type of restart
allows to decrease the residue by reducing the components of the initial vector following the modes
undesirable. Into arithmetic exact, there would be immediately Rm = 0 and the process would stop there!
Iteration after iteration, one thus improves quality of the modes sought while resting on
auxiliary modes. Of course this one is estimated at each stage and conditions the opportunity of
to simulate another restart or not. The process can be summarized (very macroscopically!) like
follows:
Factorization of Arnoldi of command m: AQ = Q B + R.
m
m m
m
For K = 1, NMAX_ITER_SOREN to make
~ ~
~
~
~
To calculate
1
, 2
… p
, p+1… p+q,
! “
#
$
#
! “
# # $
#
Preserved for Utilisées like
améliorat ion
shifts
QR with implicit shifts,
Update Q, B and R
m
m
m,
Truncation of these matrices has the command p,
AQ = Q B + R AQ = Q B + R,
p
p p
p
m
m m
m
Quality estimate of the p modes.
Fine loops.
Algorithm 14: Method IRA (known as of Sorensen)
The maximum number of iterations is controlled by key word NMAX_ITER_SOREN of the key word factor
CALC_FREQ.
It should be noticed that the Arnoldi algorithm is prudently supplemented by a reorthogonalisation
total (started that if that proves to be necessary). This overcost is all the more acceptable as it is
established via algorithm IGSM of Kahan-Parlett (cf [Annexe 2]) which is particularly effective.
All this makes it possible to ensure us of the good behavior of projection with respect to the initial spectrum.
In addition, the evaluation of the quality of the modes is not carried out simply by building p
residues R
= (A
~
- ~ I) U
I
I
I
via property 10. One already mentioned that in the case not
2
2
square, they could not suffir with this task, in particular in the event of strong defect of normality. For of
to discharge, without having recourse to other information (to obtain a criterion exact it would be necessary to be able
to consider conditionings spectral of spaces invariants and the angles which they make between them. It
who is sometimes difficult to obtain, even a posteriori!) a priori, the preceding property is used
supplemented by a criterion due to Z.Bai et al. [bib34]
~
B
T
+1,
E X < my (
X B, PREC_SOREN
m
m
m
m
)
where is the precision machine and PREC_SOREN is a key word initialized under CALC_FREQ. The user has
thus a quality control (partial) of the modes, that of which it did not lay out with the other methods
established in the code. Taking into account the various standards used, this error is different from
that resulting from the total postprocessing (cf [§2.9]) which is displayed in the columns results.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
56/78
Note:
· the technique of polynomial acceleration used is more effective than that of Tchebycheff,
since the latter is explicit and requires m produced matrix-vector,
· to avoid the deterioration of the vectors of Ritz (and thus of the approximate clean vectors) by
eigenvalues of very large modules (partners with the core of B in “shift and invert”)
a filtering of the type Ericsson & Ruhe [bib35] was established. More robust techniques
exist but they require information a priori concerning in particular the blocks of
Jordan associated with the core with the operator (cf Meerbergen & Spence, 1996),
· this algorithm can be seen like the truncated shape of implicit algorithm QR of
J.C.F. Francis (cf [Annexe 1]).
The following paragraph will clarify the choices which led to the alternative installation in the code.
6.5
Establishment in Code_Aster
6.5.1 ARPACK
The package of origin [bib29] (ARPACK for ARnoldi PACKage) coding the method is available in
freeware on Internet. These designers, D. Sorensen, R. Lehoucq and C.Yang of Rice University of
Houston, wanted it at the same time:
· simple of access (FORTRAN77, “reverse communication”),
· flexible (it is based on libraries LINPACK, LAPACK [27] and BLAS [bib36] (BLAS
is the acronym for Basic Linear Algebra Subprograms)),
· and rich in functionalities (decomposition of Schur, shifts flexible, many
spectral transformations).
Its effectiveness is multiplied by ten by the use of very optimized BLAS from level 2 and 3 and by the setting in
place “reverse communication”. The user is thus a Master of his structures of data and
of its procedures of processing concerning the operator and the scalar product of work. It is him which
provides to routines ARPACK this type of information. That thus makes it possible to manage as well as possible, with
tools and the procedures ASTER, the products matrix-vector, factorizations, resolutions of
system…
6.5.2 Adaptations of the algorithm of Sorensen
To deal with generalized modal problems, this package proposes a whole series of
spectral transformations, in reality or complex, square or not. Into square, the algorithm of
Sorensen is based on the couple Lanczos-QL (IRLM for Implicit Restarted Lanczos Method).
two approaches (square or not) are not designed besides to treat pseudo-products
scalars related to indefinite matrices.
Precisely, for at the same time circumscribing numerical problems involved in their properties enough
heterogeneous in the code and, in addition, to ensure itself of a better total robustness, we have
chosen to work in nonsymmetrical (IRAM with Arnoldi and QR), on the couple operator of
following scalar work-product:
(A -
)
B -
1
1
B U =
U
!# “
#
$
# #
µ -
With
!“$
(X, y) = yTx
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
57/78
One could have dealt with the problems of buckling in “buckling mode” via the same “shift and invert” and
the scalar pseudo-product introduced by A. Mais because of quasi-systematic introduction of
Lagranges, this matrix becomes indefinite even singular, which disturbs the process largely.
The same causes produce the same effects when, for a calculation of dynamics, one A resorts to
scalar B-product. Rather than to modify all the package while introducing a scalar pseudo-product,
we thus chose a simple scalar product “Euclidean” more robust and much less
expensive.
It was thus necessary to modify the procedures of “reverse-communication” of the package, because it did not envisage
this option (with matrices standards one classically prefers to enrich the components with one
matric scalar product, even in nonsymmetrical). Contrary to the alternative of Newmann &
Pipano (cf [§5.5.1]) introduced for Lanczos, we deliberately placed ourselves in one
nonsymmetrical configuration. But in order to avoid as much as doing it can the problems of orthogonality
recurring, even into symmetrical, we would have chosen the version of IRAM using Arnoldi.
The disadvantage of this step is that it is necessary, in preliminary postprocessing of IRAM,
B-orthonormaliser approximate clean vectors to find property 2 numerically
exploited by the modal recombinations. This stage does not disturb the base of clean modes
exhumed and it is very effectively carried out via the IGSM of Kahan-Parlett.
The taking into account of the limiting conditions and, in particular of the double dualisations, was led
as for Lanczos according to the procedure describes in [§2.2]. In particular, once the vector
initial in acceptable space one is applies the operator of work to him. This traditional process
allows to purge the vectors of Lanczos (and thus the vectors of Ritz) of the components of the core.
In addition, in certain configurations for which the number of frequencies requested p
is equal to the number of ddls active (cf [§2.2]), one owed bluffer the algorithm which stopped in error
fatal! Indeed, it detected generally well space invariant awaited (from size p), but because of
particular structure of the clean vectors associated Lagranges (cf proof of property 4, [§2.5])
it had much evil to generate an initial vector which is proportional for them.
One would have needed particular processing taking account of the classification of these Lagranges, which
would have been all the more expensive as they are not fundamentally necessary to solve it
problem requested! All the spectral information being already present at deepest of the algorithm,
it is not thus the sorrow to complete the two remaining iterations (when the user asks
p = nddl-credits, one imposes m= p+2 automatically). It is enough to short-circuit the natural wire of
the algorithm, to withdraw the interesting modes of Ritz and post-to treat them to return in space
initial.
Note:
This type of case of figure in which one seeks a number of clean modes very near to
ddls numbers leaves the “ideal” perimeter of use of this type of algorithm (cf [§2.8]). A good
QR would be without any doubt more effective, but it is a good means of testing the algorithm.
6.5.3 Parameter setting
To be able to activate this method, it is necessary to initialize key word METHODE with “SORENSEN”. Size of
subspace of projection is determined, either by the user, or empirically from
formulate:
m = min (my (
X 2p, p + 2), nddl-credits)
where
p
is the number of eigenvalues to calculate,
nddl-credits
is the number of active degrees of freedom of the system (cf [§2.2])
The user can always impose to him even dimension by indicating it with the key word
DIM_SOUS_ESPACE of the key word factor CALC_FREQ.
The parameter of the IGSM of Kahan-Parlett (cf [Annexe 2]) PARA_ORTHO_SOREN, the maximum number
iterations of the total process, NMAX_ITER_SOREN, and the criterion of quality control of the modes,
PREC_SOREN, are accessible by the user under the key word factor CALC_FREQ. When this last
key word is null, the algorithm initializes it with the precision machine (= 2.22.E-16 on SGI).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
58/78
6.5.4 Display in the file message
The example below resulting from the list of case tests of the code (ssll103b) recapitulates the whole of the traces
managed by the algorithm. One finds in particular, for each critical load (or frequency),
the estimate of its quality via the standard of error.
Here, IRAM reiterated only only once and used 30 IGSM (in its first phase). The resolution
total consumed 91 products (makes some, less than that because of “implicit” introduction of the product
Euclidean scalar) matrix-vector and 31 inversions of system (just increase because the operator of
work is already factorized).
------------------------------------------------------------
THE NUMBER OF DDL
TOTAL EAST: 68
LAGRANGE EAST: 14
THE NUMBER OF ACTIVE DDL EAST: 47
------------------------------------------------------------
The SELECTED OPTION EAST: PLUS_PETITE
THE VALUE OF SHIFT CRITICAL LOAD EAST: 0.00000E+00
------------------------------------------------------------
INFORMATION ON CALCULATION REQUIRES:
A NUMBER OF REQUESTS MODES: 10
THE DIMENSION OF REDUCED SPACE EAST: 30
=============================================
= METHOD OF SORENSEN (CODE ARPACK) =
= VERSION: 2.4 =
= DATE: 07/31/96 =
=============================================
A NUMBER OF RESTARTINGS = 1
A NUMBER OF PRODUCTS COP * X = 31
A NUMBER OF PRODUCTS B * X = 91
A NUMBER OF REORTHOGONALISATIONS (STAGE 1) = 30
A NUMBER OF REORTHOGONALISATIONS (STAGE 2) = 0
A NUMBER OF RESTARTINGS OF A NULL V0 = 0
------------------------------------------------------------
CRITICAL LOADS CALCULEES INF. AND SUP. SONT:
CHARGE_CRITIQUE_INF: - 9.96796E+06
CHARGE_CRITIQUE_SUP: - 6.80007E+05
------------------------------------------------------------
MODAL CALCULATION: METHOD Of SIMULTANEOUS ITERATION
METHOD OF SORENSEN
NUMBER CRITICAL LOAD NORMALIZES ERROR
1 - 6.80007E+05 5.88791E-12
2 - 7.04572E+05 1.53647E-12
3 - 7.09004E+05 1.16735E-12
4 - 7.10527E+05 1.72306E-12
5 - 7.11205E+05 2.41783E-12
6 - 7.11542E+05 7.88981E-13
7 - 7.11703E+05 5.71621E-13
8 - 1.50492E+06 1.17776E-11
9 - 6.02258E+06 2.42221E-11
10 - 9.96796E+06 3.55014E-12
------------------------------------------------------------
CHECKING A POSTERIORI OF THE MODES
IN the INTERVAL (- 1.00178E+07, - 6.76607E+05)
IT THERE A WELL 10 LOAD (S) CRITICAL (S)
------------------------------------------------------------
Example 6: MODE_ITER_SIMULT with “SORENSEN”
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
59/78
Note:
The introduction of this method made it possible to balance many cards of anomalies related to
multiplicities, of the clusters or search of eigenvalues of orders of magnitude very
different on which Lanczos and Bathe & Wilson stumbled.
Now let us recapitulate the parameter setting available of operator MODE_ITER_SIMULT with this
option METHOD = “SORENSEN”.
6.5.5 Summary of the parameter setting
Key word factor
Key word
Default value
References
“DYNAMIC” TYPE_RESU
“DYNAMIQUE”
[§2.1]
“MODE_FLAMB”
[§2.1]
METHOD “SORENSEN”
“SORENSEN”
[§6.4]
CALC_FREQ
FREQ
[§4.4]
CHAR_CRIT
[§4. 4]
OPTION “PLUS_PETITE”
`PLUS_PETITE ''
[§4.4]
“BANDE”
[§4.4]
“CENTER”
[§4.4]
NMAX_FREQ
10
[§4.4]
DIM_SOUS_ESPACE
Calculated
[§6.5.3]
PREC_SOREN
0.
[§6.4], [§6.5.3]
NMAX_ITER_SOREN
20
[§6.4], [§6.5.3]
PARA_ORTHO_SOREN
0.717
[§6.5.3], [Appendix 2]
NPREC_SOLVEUR
8
[§2.6]
NMAX_ITER_SHIFT
5
[§2.6]
PREC_SHIFT
0.05
[§2.6]
SEUIL_FREQ
1.E-02
[§2.9]
VERI_MODE
STOP_ERREUR “YES”
“OUI”
[§2.9]
“NON”
[§2.9]
PREC_SHIFT
5.E-03
[§2.9]
SEUIL
1.E-06
[§2.9]
STURM “YES”
“OUI”
[§2.9]
“NON”
[§2.9]
Table 6.5.5-a: Récapitulatif of the parameter setting of MODE_ITER_SIMULT with “SORENSEN”
Note:
· in V5, the user can specify the class of membership of his calculation by initializing it
key word TYPE_RESU. According to this value, one informs vector FREQ or CHAR_CRIT,
· one finds all the “tripaille” of parameters related to the preprocessings of the test of Sturm
(NPREC_SOLVEUR, NMAX_ITER_SHIFT, PREC_SHIFT) and with postprocessings of checking
(SEUIL_FREQ, VERI_MODE),
· at the time of the first passages, it is strongly advised to modify only them
principal parameters noted in fat. The others relate to more the mysteries of
the algorithm and they were initialized empirically with values standards,
· in particular, to improve quality of a mode, the fundamental parameter is
dimension of subspace DIM_SOUS_ESPACE.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
60/78
7
Method of Bathe and Wilson (METHODE = `JACOBI)
7.1 Principle
The method of Bathe and Wilson is a method of simultaneous iterations which consists in extending
the algorithm of the iterations opposite. One works starting from the problem shifted A X = (-) B X.
One distinguishes four parts in the algorithm [bib1], [bib4]:
· to choose
p initial vectors independent X, X,…, X
1
2
p and to build matrix X that they
generate,
· to calculate the clean elements in the subspace of Ritz while solving
(A - B
T
T
I
) ui = 0 where A = Q AQ and B = Q BQ. One returns then to space
initial (for the clean vectors) via the transformation X = qu where U = {
[ui}],
· to test the convergence of the clean modes I.
7.2
Tests of convergence
The method of Bathe and Wilson converges towards p smaller eigenvalues provided that the p
initial vectors are not B-orthogonal with the one of the clean vectors. In addition, matrices
With and B tend towards diagonal matrices. For this reason and as the matrices are full,
one uses the method of Jacobi (cf [Annexe 3]) to find the elements clean of the subspace of
Ritz.
To test the convergence of the eigenvalues, one classifies them after each iteration by command
growing in absolute value and one looks at if, for each eigenvalue, the following test is checked
k+1 K
K
PREC BATHE
+
-
1
_
where the exhibitor K indicates the iteration count. If after NMAX_ITER_BATHE iterations, one does not have
converged for all the eigenvalues, a message of alarm is transmitted in the file message.
7.3
Establishment in Code_Aster
7.3.1 Dimension of the subspace
If one wishes to calculate p eigenvalues, it is recommended to use under space of dimension Q
higher. One will check convergence only for the R smaller eigenvalues where p R Q. It
seem that r= p is not sufficient: one can find the good eigenvalues but the vectors
clean are not correct (convergence is slower for the clean vectors than for
eigenvalues). R = (p + Q)/2 seems a good choice. For Q one usually takes [bib1]
Q =
(
min p + 8, 2 p).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
61/78
7.3.2 Choice of the initial vectors
To choose Q initial vectors, one operates in the following way:
With
· first vector such as X
II,
I
1 = Bii
· for the other vectors of 2 with (Q -)
1
0
0
1
- - line 2
I
X = 1
- - line I
X
2
1 Q - 1 = 0
,…
0
0
With
where I are the indices the corresponding to smallest successive values of
II,
Bii
· the last
vector
xq, random vector.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
62/78
7.3.3 Parameters in Code_Aster
To be able to use the method of Bathe and Wilson, command MODE_ITER_SIMULT should be chosen
and to select METHODE = “JACOBI”. Two parameters relating to convergence directly
method are accessible under the key word factor CALC_FREQ using the key words
PREC_BATHE and NMAX_ITER_BATHE. One finds there also those managing the internal method of resolution
modal, PREC_JACOBI and NMAX_ITER_JACOBI.
Key word factor
Key word
Default value
References
“DYNAMIC” TYPE_RESU
“DYNAMIQUE”
[§2.1]
“MODE_FLAMB”
[§2.1]
METHOD “JACOBI”
“SORENSEN”
[Appendix 3]
CALC_FREQ
FREQ
[§4.4]
CHAR_CRIT
[§4. 4]
OPTION “PLUS_PETITE”
`PLUS_PETITE ''
[§4.4]
“BANDE”
[§4.4]
“CENTER”
[§4.4]
NMAX_FREQ
10
[§4.4]
DIM_SOUS_ESPACE
Calculated
[§7.3.1]
PREC_BATHE
1.E-10
[§7.2]
NMAX_ITER_BATHE
40
[§7.2]
PREC_JACOBI
1.E-12
[Appendix 3]
NMAX_ITER_JACOBI
12
[Appendix 3]
NPREC_SOLVEUR
8
[§2.6]
NMAX_ITER_SHIFT
5
[§2.6]
PREC_SHIFT
0.05
[§2.6]
SEUIL_FREQ
1.E-02
[§2.9]
VERI_MODE
STOP_ERREUR “YES”
“OUI”
[§2.9]
“NON”
[§2.9]
PREC_SHIFT
5.E-03
[§2.9]
SEUIL
1.E-06
[§2.9]
STURM “YES”
“OUI”
[§2.9]
“NON”
[§2.9]
Table 7.3.3-a: Récapitulatif of the parameter setting of MODE_ITER_SIMULT with “JACOBI”
Note:
· in V5, the user can specify the class of membership of his calculation by initializing it
key word TYPE_RESU. According to this value, one informs vector FREQ or CHAR_CRIT,
· one finds all the “tripaille” of parameters related to the preprocessings of the test of Sturm
(NPREC_SOLVEUR, NMAX_ITER_SHIFT, PREC_SHIFT) and with postprocessings of checking
(SEUIL_FREQ, VERI_MODE),
· at the time of the first passages, it is strongly advised to modify only them
principal parameters noted in fat. The others relate to more the mysteries of
the algorithm and they were initialized empirically with values standards.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
63/78
8
Conclusion - Synthesis
The optimal perimeters of use of the modal operators of Code_Aster can be dissociated.
When it is a question of determining some eigenvalues (typically a half-dozen) or
to refine some estimates, operator MODE_ITER_INV is completely indicated. He gathers
heuristic algorithms and those of type powers (cf [§3]), which were historically developed
first to solve generic modal problems. He couples a phase of localization of
eigenvalues (via a technique of bisection and a method of the secant) with an improvement
of these estimates and a calculation of the clean vectors associated by a method of iterations opposite
(crossbred, or not, of quotient of Rayleigh).
On the other hand, to capture a significant part of the spectrum, one A resorts to MODE_ITER_SIMULT.
This last federates the methods known as of “subspace” (Lanczos [§4], [§5], IRAM [§6], Bathe & Wilson
[§7]) which projects the operator of work in order to obtain an approximated spectrum of more reduced size (treated
then by a total method of type QR or Jacobi, cf [Annexe 1] and [Annexe 3]). In addition to theirs
numerical qualities (reduced calculation complexities and memory, facilitated coupling of techniques of
deflation, of restarting and acceleration…) and mathematics (good convergence, controls
quality of the spectrum of the projected operator…), it should be noted that they necessarily do not require
knowledge of the operator of work but that of its action on a vector (this characteristic is
very useful to treat large systems but it is not used in Code_Aster where one assembles
all matrices before treating them).
Until now, these algorithms stumbled regularly on the same shelves: detection
correct of multiple modes, modes of rigid body and a general way, processing of
packed spectrum. All this led to the appearance of “phantom” modes sometimes badly easy to detect
(modes corresponding to multiplicities missed and being able to generate correct residues with the direction
of Aster, and this, more especially as the criteria of the post-modal checks were sometimes permissive
(residue in 10-2 instead of the 10-6 current), decontaminated even insufficient (test of Sturm limited to the values
clean positive)) who cause distortions of results downstream from calculation, during projections on
base modal. To be been free from the recurring problems to this type of approach, one thus proposed
to enrich MODE_ITER_SIMULT (starting from V5) by algorithm IRA (“Implicit Restarted Arnoldi”
[§6]).
This alternative of Arnoldi, initiated by D.C. Sorensen in 1992, makes real great strides for the resolution of
great modal systems on parallel supercomputers. Its framework of application is completely
General, it deals with as well the real problems as complex, square or not. He summarizes himself in one
succession of factorizations of Arnoldi whose results control restartings automatically
statics and implicit, via polynomial filters modelled by implicit iterations QR.
Algorithm IRA tries to bring an elegant remedy for the recurring problems raised by the others
approaches. There is not any more a question to be posed concerning the strategy of reorthogonalisation,
frequency of the restarts, their establishments, the criteria of detection of possible phantom modes…
In short, numerically, IRAM gets a better total robustness while improving them
calculation complexities and memory (especially compared to simple Lanczos such as that of Newman &
Pipano) for a fixed precision.
From a functional point of view, the establishment of this method allowed a gain in calculation complexity
(observed) with minima of command 2. From now on, the user has a real control on the quality of the modes via one
suitable parameter setting. One, moreover, did not reinforce the severity of the post-modal checks and extended them
applicability.
The use of IRAM (by defect in MODE_ITER_SIMULT) is thus with advising in all them
case of figures, including for the search for some eigenvalues (one can even draw part
excellent properties of the method of the powers opposite concerning the calculation of vectors
clean and, to refine the result, by starting again MODE_ITER_INV with for estimate values
clean exhumed by MODE_ITER_SIMULT). Beyond their numerical specificities and
functional calculuses which are included in this document, one can synthesize the modal methods of
Code_Aster in the shape of table below (the default values are materialized in
fat).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
64/78
Operator/
Algorithm
Key word
Advantages
Disadvantages
Perimeter
of application
MODE_ITER_INV
1ère phase
(heuristics)
Calculation of some
Bisection
“SEPARE”
modes
Calculation of some
Bisection +
“AJUSTE”
Better precision
Cost calculation
modes
Secant
Improvement of
Initialization by
“PROCHE”
Resumption of values
No the capture of
some estimates
the user
clean estimated
multiplicity
by another
process.
Cost calculation of this
phase quasi-no one
2nd phase
(method of
powers properly
said)
Basic method
Powers
“DIRECT”
Very good
Not very robust
opposite
construction of
clean vectors
Option of acceleration
Quotient of
“RAYLEIGH”
Improve
Cost calculation
Rayleigh
convergence
MODE_ITER_SIMULT
Calculation of part of
Bathe & Wilson
“JACOBI”
Not very robust
spectrum
Lanczos
“TRI_DIAG”
Not very robust
(Newman- Pipano)
IRAM (Sorensen)
“SORENSEN”
Increased robustness.
Better
calculation complexities and
memory.
Control
quality of the modes.
Table 8-a: Récapitulatif of the modal methods of Code_Aster
IRAM made it possible to balance all the software anomalies related to the generalized modal problems,
but it seems that alone its variation per blocks or a version incorporating of the “purging and
lock " [bib23] can guarantee to us “an quasi-infallible” detection of the spectrum of a standard operator
(i.e too badly not conditioned).
A class of algorithm known as of “Jacobi-Davidson” [bib37] seems even more promising to treat
pathogenic cases. It uses an algorithm of the Lanczos type that it packaged via a method
of Rayleigh.
Let us conclude by noting that algorithm IRA was not extended yet to the quadratic case in
modal operators of the code [R5.01.02].
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
65/78
9 Bibliography
[1]
K.J. BATHE: Broad solution methods for generalized eigenvalue problems in
structural engineering. California university Berkeley, 1971.
[2]
F. CHATELIN: Eigenvalues of matrices. Masson, 1988.
[3]
J.K. CULLUM & R.A. WILLOUGHBY: Broad Lanczos algorithms for symetric
eigenvalue computations. Flight 1 Theory, Birkhäuser 1985.
[4]
DHATT & TOUZOT: A presentation of the finite element method. Maloine
S.A. Edition, 1984.
[5]
A. DUBRULLE, R.S. MARTIN & J.H. WILKINSON: The Implicit QL Algorithm
pp241-248, in Handbook for automatic computation. Flight 2, Linear will algebra,
Springer-Verlag, 1971.
[6]
G.H. GOLUB & c.f. Van LOAN: Matrix computations, The Johns Hopkins
university near, 1989.
[7]
Y. HAUGAZEAU: Application of the theorem of Sylvester to the localization of
eigenvalues of Ax = Bx in the symmetrical case. Numerical RAIRO Analyze,
vol.14, n°1 pp25-41, 1980.
[8]
D. HO: Tchebychev technical acceleration for broad scale nonsymmetric matrices.
Numerische mathematik. Flight 56, pp731-734.
[9]
J.F. IMBERT: Analyze structures by finite elements. Sup Aéro, CEPADUES
Editions.
[10]
C. LANCZOS: Year iteration method for the solution off the eigenvalue problem off
linear differential and integral operators. Newspaper. off research. off the national office
off standard, flight 45, n°4 pp255-282, Oct. 1950.
[11]
P. LASCAUX & R. THEODOR: Matric numerical analysis applied to Art of
the engineer. Masson, 1986.
[12]
R.S. MARTIN, G. PETERS & J.H. WILKINSON: The QR Algorithm for Real
Hessenberg Matrices pp358-371, in Handbook for automatic computation. Flight.
2, Linear will algebra, Springer-Verlag, 1971.
[13]
R.S. MARTIN & J.H. WILKINSON: Similarity Reduction off has General Matrix to
Hessenberg From pp339-358, in Handbook for automatic computation. Vol. 2
Linear will algebra, Springer-Verlag, 1971.
[14]
Mr. NEWMANN & A. PIPANO: Modal Fast extraction in Nastran via FEER
computer programs, Nastran, to use manual, NASA Langley Research Center
pp485-506, 1977.
[15]
B. NOUR-OMID, B.N. PARLETT & R.L. TAYLOR: Lanczos versus subspace
iteration for solution off eigenvalue problems, Int. J. for numerical methods in
engineering. Vol. 19, pp859-871, 1983.
[16]
INTER-F OJALVO & Mr. NEWMANN: Vibrations modes off broad structures by year
automatic matrix-reduction method. AIAA, vol.8 n°7, 1970, pp1234-1239.
[17]
E.E. OSBORNE: One pre-conditioning off matrices. J. Assoc. Comput. Mach., flight 7,
pp338-345, 1960.
[18]
B.N. PARLETT: The symetric eigenvalue problem. Prentice hall, Englewoods
Cliffs, 1980.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
66/78
[19]
B.N. PARLETT & C. REINSCH: Balancing has matrix for calculation off eigenvalues
and eigenvectors. pp315-326, in Handbook for automatic computation. Vol. 2,
Linear will vlgebra, Springer-Verlag, 1971.
[20]
Y. SAAD: One the spleens off convergences off the Lanczos and the Block-Lanczos
methods. SIAM J. Numer. Anal., flight 17 N° 5, pp687-706, october 1980.
[21]
H.D. SIMON: The Lanczos algorithm with partial reorthogonalisation, Mathematics
off computation. Flight 42, n° 165, pp115-142, 1984.
[22]
Mr. SADKANE: With block Arnoldi-Tchebychev method for computing the leading
ergenpairs off broad sparse unsymmetric matrices. Numerische mathematik, flight 64,
pp181-193, 1993.
[23]
J.L. VAUDESCAL: Introduction to the methods of resolution of problems with
eigenvalues of big size. Note HI-72/00/01A, 2000.
[24]
J.H. WILKINSON & C. REINSCH: Handbook for automatic computation, flight 2,
Linear Algebra, Springer-Verlag, New York, 1971.
[25]
O. BOITEAU & B. QUINNEZ: Method of spectral analysis in Code_Aster.
Booklet of Cours of dynamics, March 2000.
[26]
D.C. PAIGE: Accuracy and effectivness off the Lanczos algorithm for the symmetric
eigenproblem. Linear will algebra and its applications, flight 34, 1980.
[27]
J.J. DONGARRA & Al: LAPACK' S users guide. SIAM, 1992.
[28]
J.J. DONGARRA & Al: EISPACK' S users guide. Springer verlag, flight 6, 1976.
[29]
R.B. LEHOUCQ, D.C. SORENSEN & C. YANG: ARPACK' S users guide. 1996.
[30]
D.C. SORENSEN: Implicit applications off polynomial filters in A k-step Arnoldi
method. SIAM J. Anal Matrix. Appl., flight 13, pp357-385, 1992.
[31]
W.E. ARNOLDI: The principle off minimized iterations in the solution off the matrix
eigenvalue problem. Quarter. Appl. Maths, flight 9, pp17-19, 1951.
[32]
Y. SAAD: Variation one Arnoldi' S method for computing eigenelements off broad
unsymetric matrices. Linear will algebra and its applications, flight 34, pp269-2395, 1980.
[33]
Y. SAAD: Broad Numerical meyhods for eigenvalue problems. Manchester
university, close series in algorithms and architectures for advanced scientific
computing, 1991.
[34]
Z. BAI, J. DEMMEL & A.M. KENNEY: One computing condition numbers for the
nonsymmetric eigenproblem. ACM transactions one mathematical software, flight 19,
pp202-223, 1993.
[35]
T. ERICSSON & A. RUHE: Spectral The transformation Lanczos method for the
numerical solution off broad sparse generalized symmetric eigenvalue problems.
Mathematics off computations, flight 35, pp1251-1268, 1980.
[36]
J.J. DONGARRA & Al: Year extended set off FORTRAN BASIC linear will algebra
subprograms. ACM transactions one mathematical software, flight 14, pp1-17, 1998.
[37]
R.B. MORGAN & D.S. SCOTT: Preconditionning the Lanczos algorithm for sparse
symmetric eigenvalue problems. SIAM J. Sci. Comp, flight 14, 1993.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
67/78
Appendix 1 Généralités on algorithm QR
A1.1 Principe
The algorithms of the type QR (cf [§2.8]) were had a presentiment of per H. Rutishauser (1958) and were formalized
jointly by J.C. Francis and V.N. Kublanovskaya (1961). This fundamental method is
often implied in the other approaches adapted better to deal with the problems of large
sizes (in particular methods of projection).
For K = 1,… to make
H = Q R (factorization QR),
K
K
K
H
= R Q;
K +1
K
K
Fine loops.
Algorithm 1.1: Theoretical QR
The process leads repeatedly towards a matrix triangular HK higher (or triangular
by blocks) whose diagonal terms are the eigenvalues of the initial operator H= H1. The notation
H is not in fact not innocent, because there is any interest with transforming beforehand orthogonally (
manner to modify only the shape of the shifté operator and not his spectrum) the operator of work A
in the form of higher Hessenberg, that is to say into arithmetic real
H = QT AQ
1
0
0
This can be carried out via various orthogonal transformations (Householder, Givens,
Gram-Schmidt…) and their cost (about O (10n3/3)) negligible is compared with the gain that they
allow to realize with each iteration of the total process: O (N2) (with Householder or
Fast-Givens one more precisely has O (2n2) against O (4n2) with simple Givens) against O (n3). It
gain of a command magnitude can be even improved when the operator is tridiagonal symmetrical
(it is the case of Lanczos with a true scalar product): O (20n).
Convergence towards a simple triangular matrix is carried out only if all the eigenvalues
are distinct modules and that if the initial matrix is not “pathologically” too poor in
clean directions. The convergence of the ième mode (arranged classically by order descending of
modulate) is carried out then in:
J
max
,
J I
I
what can prove very slow if no complementary process is implemented.
Note:
· determination of the spectrum of the matrix of Rayleigh with the alternative of Newmann & Pipano
(cf [§5.5.1]) is carried out via a QR (or a QL into symmetrical) simple of this type (with, with
precondition, balancing). The only parameter accessible by the user is the maximum number
acceptable iterations NMAX_ITER_QR,
· for IRAM (cf [§6.4]), this calculation is carried out via a method QR with double explicit shift
whereas the polynomial filters managing the restarts use a QR with double implicit shift.
By prudence, no parameter is accessible for the user in Code_Aster!
· one should not confuse the method, algorithm QR, and one of its conceptual tools,
factorization QR,
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
68/78
· this class of algorithm is very much used to determine the complete spectrum of an operator,
because it is very robust (it is the reference in this field). However it is very
greedy memory from there what places makes its use crippling on great systems,
· the perimeter of application of algorithm QR is much more general than that of Jacobi
(it is the second standard algorithm providing all the spectrum of an operator for can
that one is ready to entirely store it) who is limited to the square matrices.
To accelerate the convergence of the simple algorithm which can be very slow (in the presence of clusters
for example) a multitude of alternatives, based on the choice of shifts answering certain criteria,
were born.
A1.2 strategy of the shift
This strategy consists in causing artificially a phenomenon of deflation (cf [§3.1], [§5.4.1])
within the matrix of work. That offers triple favors:
· of
to be able to isolate a real eigenvalue even two combined complex eigenvalues,
· all
in
reducing the size of the problem to be treated,
· and
in
accelerating convergence.
In its version with simple explicit shift, the method is rewritten then in the following form:
For K = 1,… to make
To choose the shift µ,
Sk = HK - µ I,
S = Q R (factorization QR),
K
K
K
H
= R Q + µ;
K +1
K
K
I
Fine loops.
Algorithm 1.2: QR with simple explicit shift
Note:
This process spreads intuitively with several shifts. One builds then, for each
total iteration K, as many auxiliary matrices Sik of shift µi.
The convergence of the process is largely improved in the direction where under-diagonal terms
cancel themselves asymptotically in:
J - µ
max
.
J I
I
- µ
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
69/78
In theory, if this shift µ is eigenvalue of the problem, then deflation is exact. In practice, them
effects of round-off disturb this phenomenon, it is what is called the property of direct instability of
the algorithm. The principal difficulty lies in the choice of (or of) the shifts. In addition, one does not keep
not the same shift for all the iterations. One must change some when it is associated a value
clean converged. Indeed, it will have numerically caused its ousting of the spectrum of work in
causing a deflation with the preceding iteration.
Since the Sixties a whole zoology of shifts developed. While simplest
use the last diagonal term Hn, N, that of J.H. Wilkinson [bib24] consists in determining
analytically the eigenvalue µ of the diagonal block
H
H
n-1, n-1
n-1, N
H
H
N, n-1
N, N
nearest to this term. This technique makes it possible to even obtain a quadratic convergence
cubic (in the symmetrical case) and it proves particularly effective to capture modes
double or of the distinct eigenvalues of the same module. However this strategy can appear
ineffective in the presence of combined complex clean modes.
The same author then proposed an alternative including a double shift corresponding to both
complex eigenvalues µ1 and µ2 (beforehand given) of the accused block. But in addition to
numerical difficulties with réfreiner appearance of invading complex components (which in theory
do not take place to be), one has in addition much evil to preserve the character of “Hessenberg
higher " of the matrix of work. As well as possible, that slowed down the convergence of algorithm QR, with
worse, that distorts its operation completely.
In order to bearing with these underhand numerical disadvantages, J.C. Francis developed a version with
double implicit shift. It is very well clarified in [bib6] (pp377-81) and one will be satisfied just here
to summarize philosophy of it.
To minimize the effects of round-offs, it would be necessary to constitute the auxiliary matrix resulting Sk directly
simultaneous application of the two shifts µ1 and µ2
S = H2 - (1
µ + µ2) H + (µ1µ2) I
K
K
K
before factorizing it in form QR and building the new one reiterated Hk+1
S = Q R
K
K
K,
H + = QTH Q
K
K
K
K.
1
But only the cost of the initial assembly (about O (n3)) the tactics make inoperative. This
alternative, based on the Q-implicit theorem, consists in applying to the matrix of work of
particular transformations of Householder allowing to find a matrix of Hessenberg
“primarily” equalizes in Hk+1, i.e. of the type:
~
H
1
1 = E H
E
E
K +
K +1
with = dia (
G ± 1. ±
. )
1.
The matrix of work thus preserves, at lower cost, its particular structure and its spectrum.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
70/78
All these alternatives are in fact very sensitive to the techniques of balancing established for
préconditionner the initial operator. The following paragraph will summarize this technique of “setting with
the scale " (“balancing” for the Anglo-Saxons) of the terms of the matrix of work which is very
employee in modal calculation.
A1.3 Equilibrage
It is a question of mitigating the effects of round by attaching the perimeter of expansion of the terms of
the operator of work, i.e. to prevent that they do not become too small or, on the contrary, too large.
On this subject, E.E. Osborne [bib17] (1960) noticed that generally the error on the calculation of
clean elements of A is about A 2 where is the precision machine. It proposed then
to transform the initial matrix into a matrix
~
With = D A
1 D (with D a diagonal matrix),
such as:
~
With << A.
2
2
In fact one calculates a succession of matrices repeatedly
WITH = D-1 A
D
K
K - 1
K - 1
K - 1
such as:
With
With
K
K
F
checking have
= have with have and has
2
2
I, respectively, the ième column and line of Af.
Note:
· this technique was generalized to any matric standard induced by the discrete standard of lq
Q
with lq the whole of complex continuations (one) N such as one <,
N
· its employment is very widespread in scientific computation and in particular among the direct solveurs of
system of linear equations.
The base of calculation of the computer, noted, intervenes in the determination of the terms of matrices Dk.
In order to minimize the round-off errors, the elements of Dk are chosen so that they are powers
of this base.
Are Rk and Ck the p standards (in practice p= 2 is often taken), respectively line and
column I of the Ak-1 matrix (I is selected in such way that I - 1 K - 1 [N]). By supposing that
R C
K
K 0, one shows whereas there is a single entirety signed such as:
2 - 1 R
< K 2 +1.
Ck
Either F =, one defines matrix Dk such as:
Id + (-)
1 E and
p
p
p
p
F
I I
if (Ck F) + (R K/F) < C + R
D
(K K)
K = Id
if not
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
71/78
where 0 < 1 is a constant and I.E.(internal excitation) the ième canonical vector. One builds then repeatedly:
D
= D D
K
K
K - 1,
,
With = D-1A
D
K
K
K - 1 K.
by initializing D0= Id. the process stops as soon as D
Id
K =
.
Note:
· before carrying out balancing, a search of the isolated eigenvalues is carried out in
detecting the presence of lines and quasi-null columns (all the terms are null, except
that placed on the diagonal). When there are some one can, by carrying out permutations
suitable, to put the matrix of work in the following form blocks:
*
*
*
Y
Z
0
*
0
X
T
*
*
0
0
*
0
*
It is then necessary to carry out balancing only on the central block X because the terms
diagonal of the two higher triangular matrices are eigenvalues of the matrix
of work.
·
p = 1
= 0 9
. 5
in Code_Aster, one chose
and
(cf [bib19]),
· before applying method QR, one starts by balancing the matrix of work and then,
one transforms it in the form of Hessenberg higher. Once calculation QR carried out, one
go up vectors of Ritz to the clean vectors via the reverse of the transformations
orthogonal due to the setting in Hessenberg form and balancing. It is this process which
is set up as well in Lanczos in IRAM. However in addition to the fact that it
require the storage of these orthogonal matrices, it is also and especially extremely
sensitive to the effects of round-offs. Thus, it would be largely preferable to estimate the vectors
clean via a method of the powers opposite initialized by the eigenvalues
exhumed.
A1.4 method QL
Factorization QL of a matrix A consists with orthonormaliser its vectors columns of
starting with the last (contrary to algorithm QR which begin with the first) thus giving
a matrix L triangular lower regular. It is shown besides that algorithm QL applied to
invertible operator A is equivalent to algorithm QR applied to A * (transposed matrix
combined reverse of A).
The method installation in Code_Aster is identical to the method of simple shift of
J.H. Wilkinson seen previously by adapting the search for this shift to the miner 2x2 highest.
Note:
Contrary to QR which captures the eigenvalues by order ascending of module, one obtains here
preferentially dominant modes, then others, by order descending of module.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
72/78
Appendix 2 Orthogonalization of Gram-Schmidt
A2.1 Introduction
One saw on several occasions in this document that the quality of orthogonalization of a family of
vectors is crucial for the good unfolding of the algorithms and the quality of the modes obtained.
This task is besides to permanently realize from where the importance of a fast algorithm.
The simple algorithms of orthonormalisation are deduced from the traditional process of Gram-Schmidt
but they are often “conditionally stable” (the quality of their work depends on
matric conditioning of the family with orthonormaliser). This defect of robustness can prove
problems to treat situations particularly badly conditioned. One prefers to them then
more expensive but robust algorithms, containing projections or of rotations: transformations
spectral of Householder and Givens.
In practice, for the establishment of method IRA in Code_Aster, we retained one
iterative version of the process of Gram-Schmidt (the IGSM of Kahan-Parlett [bib18]) Elle carries out one
good compromise between the robustness and calculation complexity since it is unconditionally
stable and orthogonaliser allows except for the precision machine, in, to the maximum, twice more
time that traditional Gram-Schmidt (GS).
In the following paragraphs we will detail the operation of the basic algorithm, thus
that that of these two principal alternatives. First was installation in MODE_ITER_INV and
second is used in MODE_ITER_SIMULT. Comparative very percussion of these methods is
declined in [bib23] (pp33-36) starting from a very simple example.
A2.2 Algorithme of Gram-Schmidt (GS)
Being given K vectors independent of N, (xi)
one wishes to obtain K orthonormal vectors
I =1, K
(compared to an unspecified scalar product) (yi)
space which they generate. In others
I =1, K
terms, one wishes to obtain another orthonormal family generating same space. The process
of traditional orthonormalisation of Gram-Schmidt is as follows:
For I = 1, K to make
For J = 1, I - 1 to make
Calculation of rji = (y, X
J
I);
Fine loops.
Calculation of y % I = xi - R y
ji
J,
J =,
1 i-1
Calculation of II
R = y % I,
y %
Calculation of y
I
I =
;
II
R
Fine loops.
Algorithm 2.1: Algorithm of Gram-Schmidt (GS)
This process is simple but very unstable because of the round-off errors, which causes to produce
nonorthogonal vectors. In particular when the initial vectors are almost dependant
that creates important variations magnitude in the second stage of the process
% y = X - R y
I
I
ij
J
J =,
1 i-1
From a numerical point of view, the management of these variations is very difficult.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
73/78
Note:
· the problem is completely similar to that met by the inversions of systems put in
place in the code (cf [§2.6]),
· by noting Q the matrix generated by (yi)
, one explicitly was thus built
I =1, K
factorization QR of the initial matrix X related to (xi)
. It is in fact the goal of very proceeded
I =1, K
of orthogonalization.
A2.3 Algorithme of Gram-Schmidt Modifié (GSM)
In order to évincer these instabilities numerical one reorganizes the preceding algorithm. Mathematically
equivalent with the preceding process, this one is then much more robust because it avoids the variations of
magnitude important between the vectors handled in the algorithm.
In the initial process, the orthogonal vectors yi are obtained without taking account of the i-1
preceding orthogonalizations. With Gram-Schmidt Modifié, one orthogonalise more
gradually by taking account of preceding deteriorations according to the process below.
For I = 1, K to make
% yi = X, I
For J = 1, I - 1 to make
Calculation of rji = (y, J % yi),
Calculation of % yi = % yi - R y;
ji
J
Fine loops.
Calculation of II
R = % y,
I
% y
Calculation of y
I
I =
;
II
R
Fine loops.
Algorithm 2.2: Algorithm of Gram-Schmidt Modifié (GSM)
The orthonormality of the basic vectors is much better with this process and it can be even obtained with
the precision machine and except for a constant (depend on the conditioning of X). However, for
to treat situations particularly badly conditioned, this “conditional” stability can prove
quickly problematic, from where the recourse to the following iterative algorithm.
Note:
GSM is twice more effective than a method of Householder to obtain a factorization
QR of the initial matrix X. It requires only O (2nk2) operations (with N the number of lines
matrix).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
74/78
A2.4 Algorithme of Gram-Schmidt Itératif (IGSM)
To ensure itself nevertheless of orthogonality except for the precision machine, one recommends to realize
one second orthogonalization. And if, at the end of the latter, orthogonality is not assured it
is not any more the sorrow to start again, the handled quantities are then certainly very close and
their variations oscillate around zero. This approach is based on a theoretical analysis due to
W. Kahan and taken again by B.N. Parlett [bib18] (cf pp105-110).
For I = 1, K to make
% yi = X, I
For J = 1, I - 1 to make
Calculation of rji = (y, J % yi),
Calculation of ~yi = % yi - R y,
ji
J
If
~yi % y then
I
% yi ~y, I
Exit loops in J;
If not
Calculation of r~ji = (y, ~y
J
I),
Calculation of ~~yi = ~yi - ~r y,
ji
J
If
~
~
yi ~
y then
I
% yi ~~y, I
Exit loops in J;
If not
% y 0,
I
Passage to I following;
End if.
Fine loops.
Calculation of II
R = % y,
I
y
Calculation of y
I
I = %;
II
R
Fine loops.
Algorithm 2.3: Algorithm of Gram-Schmidt Itératif of the type Kahan-Parlett (IGSM)
During the use of IRAM in Code_Aster, the criterion is parameterized by the key word
PARA_ORTHO_SOREN (cf [§6.5]). It is shown that its interval of validity is [1.2e, 0.83-] with
precision machine and one generally allot value 0.717 to him (by defect).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
75/78
Note:
· the larger the value of the parameter is, the less the reorthogonalisation starts, but
that affects the quality of the process,
· contrary to the version “house” developed for Lanczos (cf [§5.5.1]) here criteria
stops concentrate on the standards of the vectors orthogonalized rather than on the products
scalars inter-vectors which tend more to reflect the effects of round-off. That, joined to
suppression of the iterations higher than two, therefore useless, can explain the addition
of effectiveness of the version of Kahan-Parlett,
· according to D.C.Sorensen the paternity of this method is rather to allot to J. Daniel and Al
(paper of 1976 subjected to Mathematics off Computation, vol.30, pp772-795).
Appendix 3 Méthode de Jacobi
A3.1 Principe
Method of Jacobi [bib4], [bib6], [bib11] allows to calculate all the eigenvalues of one
generalized problem whose matrices are definite positive and symmetrical (the matrices obtained with
each iteration by the method of Bathe &Wilson check these properties; cf [§7]). It consists with
to transform matrices A and B of problem A U = B U into diagonal matrices, while using
successively orthogonal similar transformations (matrices of rotation of Givens).
process can be schematized in the following way:
A1 = A
B1 = B
A2 = QT A Q
1
B2 = QTB Q
1
1
1
1
1
…
Ak = QT Ak- Q
1
Bk
T
K -
K -
K -
= Q
B
Q
1
1
1
K - 1
K - 1
stamp
stamp
Ak
AD
Bk
Data base
K
diagonal
K
diagonal
Algorithm 3.1: Process of Jacobi
- 1
AD
The eigenvalues are given by = AD (data base) that is to say = II and stamps it clean vectors
Bdii
check:
1/AD
11
D
1 2
K
1/A
X = Q Q… Q…
22
D
1/Year
Each Qk matrix is selected so that a term (I, J) nondiagonal and not no one of Ak or of
Bk is null after the transformation (for the choice of this term, cf [§0]).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
76/78
A3.2 Quelques choice
In this algorithm, one realizes that the important points are as follows:
· How to choose the terms to be cancelled?
· How to measure the diagonal character of the matrices when K tends towards the infinite one?
· How to measure convergence?
A3.2.1 Termes nondiagonal to cancel
For the choice of the terms to be cancelled, there are several methods:
· the first consists with each stage K, to choose the largest element modulates some not
diagonal of the Ak matrix or Bk and to cancel it by carrying out a rotation (cf [§0]). This choice
ensure the convergence of the method of Jacobi but is relatively expensive (search for
the maximum element),
· the second solution consists in successively cancelling all the nondiagonal elements of
these matrices while following the natural command Ak,…, Ak, Ak
K
13
1n
23,… . When one arrives at Year-1, N, one
start again the cycle. This method converges slowly.
An alternative of this method, consists very of following the natural command of the terms, to cancel
only those which are higher than a given precision K. At the end of a cycle, one decreases
value of this criterion and one starts again. It is this strategy which is used in Code_Aster.
A3.2.2 Test of convergence
To test the convergence and the diagonal character of the matrices, one operates thus. It is checked that all
coupling coefficients defined by:
Aij
Bij
F
=
F
=
I J
Aij
Bij
AiiA jj
BiiB jj
are lower than a given precision (diagonal character of the matrices). One also controls
convergence of the eigenvalues via the indicator
K
K 1
I
- -
I
F = Max
I
K - 1
I
so that there remains lower than a given precision jaco.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
77/78
A3.2.3 Algorithme in Code_Aster
The algorithm implemented in Code_Aster is summarized with:
Initialization of the matrix of the vectors specific to the matrix identity.
Initialization of the eigenvalues.
To define the precision of necessary dynamic convergence.
To define the precision total glob.
For each cycle k=1, n_ max jaco
K
To define the dynamic tolerance: K = (dyn),
L = 0,
For each line i= 1, N
For each column j= 1, N
K, L
K, L
Aij
Bij
Calculation of the coupling coefficients F K L, =
F K L, =
I J,
Aij
K, L
K, L
Bij
K, L K, L
Aii A jj
Bii B jj
If F K L
With, or F K L
K
ij
B, then
ij
Calculation of the coefficients of the rotation of Givens,
Transformation of the Ak matrices, L and Bk, L,
Transformation of the clean vectors,
L = L + 1
End if.
Fine loops.
Fine loops.
K, L
With
Calculation of the eigenvalues K
II
I =
K, L.
Bii
K
K 1
I
- -
I
Calculation of F = Max
-
I
K 1
.
I
K, L
K, L
Aij
Bij
Calculation of the total coupling coefficients F = Max
F = Max
With
B
I, J
K, L
K, L
I, J
K, L K, L
WITH A
B B
I J
II
jj
I J
II
jj
If F
and
F
and
F
With
glob
B
glob
glob then
Correction of the clean modes (let us divide by Bkii),
Exit;
End if.
Fine loops.
Algorithm 3.2: Method of Jacobi established in Code_Aster
One notes n_ max jaco the maximum number of allowed iterations.
Note:
The matrices A and B being symmetrical, only the half of the terms is stored in form
of a vector.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Code_Aster ®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C Page:
78/78
A3.2.4 Matrice of rotation
One seeks with each stage to cancel the terms being in position I and J (I J) of the matrices
Ak, L
,
and Bk L by multiplying them by a matrix of rotation which with the following form:
G = 1 L = 1, N; G = has; G
L
ij
ji = B
+1
+
other terms being null. One wishes to have Ak, L
= Bk, L 1
ij
ij
= 0 what leads to:
has K, L
K, L
K, L
With
II + (1 + b) A has
+ B
ij
With jj = 0
has K, L
K, L
K, L
Bii +
(1+ A b) B + B
ij
B jj = 0
because Ak, L +1
WP Ak, L G
Bk, l+
=
1
and
= GTBk, L G. If the two equations are proportional one a:
K, L
With
has = 0 and B
ij
= - K, L
With jj
if not:
C
C
= 2 B has = - 1
D
D
with:
C
K, L
K, L
K, L
K, L
= A B - B A
C
K, L
K, L
K, L
K, L
1
II
ij
II
ij
2 = A jj Bij - B jj Aij
2
K, L
K, L
K, L
K, L
C
C
C = A
B
- B A
D
3
3
3
jj
II
jj
=
+
(
sign c3)
C C
II
1 2
2
2 +
Note:
· if D = 0 then one is in the case proportional,
C 2
· if the matrix B is definite positive then
3
C C
1 2
2 +
is positive [bib11] what gives one well
feel with parameter D.
A3.2.5 Récapitulatif of the parameter setting
To reach the parameters which intervene in the method of Jacobi, the user must select
following key words:
·
dyn key word PREC_JACOBI under the key word factor CALC_FREQ,
·
n_ max jaco key word NMAX_ITER_JACOBI under the key word factor CALC_FREQ.
The total precision glob is equal to the precision of convergence required in the method of
Bathe and Wilson. There is thus glob = bath (cf [§7]).
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75/01/001/A
Outline document