13.1. R5.01.01 : Modal solvers and solving generalized eigenvalue problems (GEP)

Whether it is to study the vibrations of a structure or to research its buckling modes, the mechanician must often solve a modal problem: either generalized (GEP), or quadratic (QEP) [R5.01.02]. With this intention, Code_Aster provides several methods via the operator CALC_MODES : inverse powers and Rayleigh coefficient, Lanczos, IRA, Bathe & Wilson, and QZ. They each have their scope of use, their advantages, their disadvantages and their development history.

To deal effectively with large modal problems (in mesh size and / or in the number of modes sought), we advise the use of macro-command CALC_MODES with the option ‘BANDE’ divided into several sub-intervals. This process breaks up the modal computation of a standard GEP (symmetric and real), into a succession of independent sub-calculations, which are less expensive, more robust and more precise. The gains can be significant even in seuqential calculations: factors of 2 to 5 in time, 2 or 3 in peak RAM, and 10 to 10,000 on the average mode error.

In addition, its multi-level parallelism can provide additional gains of the order of 20 in time and 2 in peak RAM, if we use about sixty processors.

In the first part of the document, we summarize the general problem of solving a modal problem, the different classes of methods and their variations in public domain libraries. In the second part, we discuss things to keep in mind and the general architecture of a modal calculation in Code_Aster. Then we detail the computational and functional aspects of each of the approaches available in the code.

A specific chapter details the implementation of parallelism and intensive computation within the framework of modal computations of the standard GEP type.

13.1.1. Introduction

A majority of studies concerning the dynamic behavior of structures are performed with a transient analysis on a modal basis. To calculate these vibration modes, many algorithms have been developed for about sixty years. In order to cope with the continual increase in the size of problems and the degradation of the conditioning of the discretized operators, only the most efficient and the more robust in practice were incorporated into the modal operator of Code_Aster : CALC_MODES.

The optimal regimes of use of the search criteria for this operator can be distinct. When the user wants to determine a few eigenvalues (typically half a dozen) or to refine a few estimates, the criterion OPTION = ‘PROCHE’ (NEAR) or ‘SEPARE’ (SEPARATED) or ‘AJUSTE’ (ADJUSTED) is quite useful. It brings together heuristic algorithms and those of power type (see Section 13.1.4).

On the other hand, to capture a significant part of the spectrum, one resorts to the criterion OPTION = ‘BANDE’ (INTERVAL) or ‘CENTRE’ or ‘PLUS_PETITE’ (SMALLEST) or ‘PLUS_GRANDE’ (LARGEST). The latter use “subspace” methods (Lanczos Section 13.1.5 / Section 13.1.6, IRAM Section 13.1.7, Bathe & Wilson Section 13.1.8) which project the operator in order to obtain an approximated modal problem of smaller size (which is then treated by a global method such as QR or Jacobi).

This operator also makes it possible to calculate in a robust way the whole spectrum of the problem thanks to the criterion OPTION = ‘TOUT’ (ALL). With this intention, we use a global reference method (QZ method Section 13.1.9) which calculates all modes exhaustively. Given its cost, however, it should be reserved for special uses: small problem size (less than 10,000 degrees of freedom) or algorithm benchmarks.

The different search criteria can also be complemented because the methods implemented with OPTION = ‘PROCHE’/’SEPARE’/’AJUSTE’ are very efficient in optimizing eigen modes that are already almost converged. In one or two iterations, they can improve the calculated eigenvectors in the first attempt via OPTION = ‘BANDE’/’CENTRE’/’PLUS_PETITE’/’PLUS_GRANDE’/’TOUT’. The projection on modal basis will only be better!

Table 13.1.1 Summary of the modal methods of Code_Aster

Operator/ Application

Algorithm

Keyword

Benefits

Disadvantages

Method of inverse powers

OPTION =

Only real+symmetric (GEP and QEP)

OPTION = ‘SEPARE’,’AJUSTE’ or ‘PROCHE’

First phase (heuristic)

Calculation of a few modes

Bisection (not applicable in QEP)

‘SEPARE’

Calculation of a few modes

Bisection + Secant (GEP) or Müller-Traub (QEP)

‘AJUSTE’

Better precision

Calculation cost

Improvement of some estimates

Initialization by the user

‘PROCHE’

Recovery of value estimates by another process. Cost of calculating this phase is almost zero

No capture of multipvalidy of modes

Second phase (method of inverse powers)

Onlysymmetrical real (GEP and QEP)

Basic method

Inverse powers

‘DIRECT’

Very good construction of eigenvectors

Not very robust

Acceleration option

Rayleigh quotient (not applicable in QEP)

‘RAYLEIGH’

Improves convergence

Calculation cost

Method ofsimultaneous iterations

SOLVEUR_MODAL =_F(METHODE=

OPTION = ‘PLUS_*’,’CENTRE’,’BANDE’ or ‘TOUT’

Calculation of part of spectrum

Bathe & Wilson

‘JACOBI’

Not very robust. Only real, symmetric (GEP)

Lanczos (Newman-Pipanoin GEP and Jennings in QEP)

‘TRI_DIAG’

Specific detection of rigid modes

Only real, symmetric (GEP and QEP)

IRAM (Sorensen)

‘SORENSEN’

Increased robustness. Best complexity of calculation and memory. Control of quality of modes. Default method.

No non-symmetric or complex symmetric.

Calculation of the full spectrum then filtering from one part

QZ

‘QZ’

Reference method in terms of robustness

Very expensive in CPU and in memory. Used in small cases (< 1000 degrees of freedom). Not for unsymmetric and complex symmetric

Options ‘BANDE’ ‘CENTRE’ ‘PLUS_PETITE’ ‘PLUS_GRANDE’ ‘TOUT’ + AMELIORATION = ‘OUI’

Subspace methods selected through SOLVEUR_MODAL then method of inverse powers

Addition implied by OPTION =’PROCHE’ on the modes calculated by the 1st step (subspace method)

Improves construction of eigenvectors

Calculation cost increases

Note

  • The effective establishment and the maintenance of the modal solvers in Code_Aster is the fruit of a work team: D. Seligmann, B. Quinnez, G. Devesa, O.Boiteau, O. Nicolas, E.Boyère, I.Nistor …

  • We have constantly tried to link the various items discussed and to limit long mathematical demonstrations to a strict minimum. In any case, the many references in the text should allow you to search for specific information.

  • The purpose of this document is not to detail all the aspects covered, complete works having already fulfilled this mission. We will mention in particular F. Chatelin 42, GHGolub 48, P.Lascaux 57, BNParlett 64, Y.Saad 67, DSWatkins 73 and the synthesis carried out by JLVaudescal 85.

13.1.3. Context

13.1.3.1. Problem

We consider the generalized modal problem (GEP)

(13.1.6)\[\begin{split}\text{Find} & \quad (\lambda, \mathbf{u}) \quad \text{such that} \\ & (\mathbf{A} - \lambda \mathbf{B}) \mathbf{u} = \mathbf{0} \quad\mathbf{u} \ne \mathbf{0}\end{split}\]

where \(\mathbf{A}\) and \(\mathbf{B}\) are matrices with real or complex coefficients, symmetrical or not (in structure and / or in values). This type of problem corresponds, in mechanics, in particular to:

  • The study of free vibrations of a non-damped and non-rotating structure. For this structure, we look for the smallest eigenvalues ​​or those which are in a given interval to know if an excitation force can create resonance. In this standard case, the matrix \(\mathbf{A}\) is the stiffness matrix, denoted \(\mathbf{K}\), real and symmetric (possibly increased by the matrix of geometric stiffness denoted \(\mathbf{K}_g\), if the structure is prestressed) and \(\mathbf{B}\) is the mass matrix or inertia matrix denoted \(\mathbf{M}\) (real symmetric). The eigenvalues ​​obtained are the squares of the angular frequency associated with the frequencies sought.

    The system to be solved can be written: \((\mathbf{K} + \mathbf{K}_g) \mathbf{u} = \omega^2 \mathbf{M}\mathbf{u}\) where \(\omega = 2\pi f\) is the angular frequency, \(f\) is the natural frequency and \(\mathbf{u}\) is the associated natural displacement vector. The eigen modes \((\omega, \mathbf{u})\) are real.

    In the presence of hysteretic damping, \(\mathbf{K}\) becomes complex symmetric. So the eigen modes are potentially complex and mismatched.

    On the other hand, if \(\mathbf{K}\) and / or \(\mathbf{M}\) remain real but possibly not symmetrical 6, the modes are either real, or in pairs \((\lambda, \tilde\lambda)\).

    This type of problem is activated by the keyword TYPE_RESU = ‘DYNAMIQUE’

    Table 13.1.2 Scope of use of operator CALC_MODES according to the search criterion and the analysis algorithm, as a function of properties of the matrices of the GEP.

    \(\mathbf{A}/\mathbf{B}\)

    Real symmetric

    Real unsymmetric

    Complex

    Real symmetric

    Most common case. No restrictions on methods. Real modes. 7

    OPTION = ‘BANDE’, ‘CENTRE’, ‘PLUS_PETITE’, ‘PLUS_GRANDE’, ‘TOUT’ with SOLVEUR_MODAL = _F(METHODE = ‘SORENSEN’ or ‘QZ’); Real or complex modes \((\lambda, \tilde\lambda)\)

    Not applicable

    Real unsymmetric

    OPTION = ‘BANDE’, ‘CENTRE’, ‘PLUS_PETITE’, ‘PLUS_GRANDE’, ‘TOUT’ with SOLVEUR_MODAL = _F(METHODE = ‘SORENSEN’ or’QZ’); Real or complex \((\lambda, \tilde\lambda)\) 8

    SIMULT with ‘SORENSEN’/’QZ’. Real or complex modes \((\lambda, \tilde\lambda)\)

    Not applicable

    Complex symmetric

    OPTION = ‘BANDE’, ‘CENTRE’, ‘PLUS_PETITE’, ‘PLUS_GRANDE’, ‘TOUT’ with SOLVEUR_MODAL = _F(METHODE = ‘SORENSEN’ or ‘QZ’), Real or complex modes, any of \((\lambda, \tilde\lambda)\)

    Not Applicable

    Not Applicable

    Other complex (hermitian, unsymmetric)

    Not applicable

    Not applicable

    Not applicable

  • Search for linear buckling modes. Under linearized theory, assuming a priori that the stability phenomena are suitably described by the system of equations obtained by assuming a linear dependence of displacement at the critical load level, the search for the buckling mode \(\mathbf{u}\) associated with this critical load level \(\lambda\), comes down to a generalized eigenvalue problem ​​of the form: \((\mathbf{K} + \lambda\mathbf{K}_g) \mathbf{u} = \mathbf{0}\) with \(\mathbf{K}\) the material stiffness matrix and \(\mathbf{K}_g\) the geometrical stiffness matrix. To transform into the “mold” of a standard GEP calculation, the code calculates, initially, the real eigen modes \((-\lambda, \mathbf{u})\) real. Then it converts them to the format of a buckling calculation: \((\lambda, \mathbf{u})\).

    This type of problem is activated by the keyword TYPE_RESU = ‘MODE_FLAMB’. This option should be reserved for real symmetric GEPs (if not, the code detects it and produces a fatal error).

Note

  • CALC_MODES allows automatic choice of the algorithm according to the search criterion, and of reducing the costs (CPU and memory) of searching for a large part of the spectrum (only in GEP with real modes) when one seeks the modes on a total interval divided into sub-intervals.

  • The user can specify the class of computation by using the keyword TYPE_RESU with ‘DYNAMIQUE’ (default value) or with ‘MODE_FLAMB’. The display of results will then be formatted taking this choice into account. In the first case we will speak of frequency (FREQ) while in the second, one will speak of critical load (CHAR_CRIT).

  • In the presence of damping and gyroscopic effects, the study of the dynamic stability of a structure leads to the solution of a modal problem of higher order, known as quadratic (QEP): \((\mathbf{K} + i \omega \mathbf{C} -\omega^2 \mathbf{M}) \mathbf{u} = \mathbf{0}\). This is solved by the two modal operators and is the subject of a specific note 79.

Now that the links between structural mechanics and the resolution of generalized modal problems have been recalled, we are going to focus on the processing of boundary conditions in the code and their impact on the mass and stiffness matrices.

6

This property is found in a Hermitian complex.

7

Vibration calculation of an undamped, non-rotating, structure. Buckling calculations.

8

Vibration calculation with hysteretic damping.

13.1.3.2. Consideration of boundary conditions

There are two ways, during the construction of the stiffness and mass matrices, to take into account the boundary conditions (this description in terms of dynamic problem is easily extrapolated to buckling):

  • Double dualisation, by using Lagrange degrees of freedom 83, makes it possible to check the “linear boundary condition” (LBC)

    \[\mathbf{C} \mathbf{u} = \mathbf{0}\]

    with \(\mathbf{C}\) a real matrix of size \(p \times n\) (\(\mathbf{K}\) and \(\mathbf{M}\) are of order \(n\)). It drives the manipulation of larger matrices (called “dualized”) because they incorporate these new unknowns. The dualized stiffness and mass matrixes then have the form

    \[\begin{split}\tilde{\mathbf{K}} = \begin{bmatrix} \mathbf{K} & \beta\mathbf{C}^T & \beta\mathbf{C}^T \\ \beta\mathbf{C} & -\alpha\mathbf{I} & \alpha\mathbf{I} \\ \beta\mathbf{C} & \alpha\mathbf{I} & -\alpha\mathbf{I} \end{bmatrix} \quad \tilde{\mathbf{M}} = \begin{bmatrix} \mathbf{M} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \end{bmatrix}\end{split}\]

    with \(\alpha\) and \(\beta\) strictly positive reals (which are used to balance the terms of the matrix). The dimension of the problem increases by \(2p\) because with \(n\) “physical” degrees of freedom, we added Lagrange dofs. There are two Lagrange multipliers per linear relation assigned to \(p\) boundary conditions.

  • Setting to zero \(p\) rows and columns of the stiffness and mass matrices. This is valid only for prescribing degrees of freedom (simple Dirichlet, no relation of proportionality between dofs). We cannot take into account a linear relation and we will speak of kinematic constraints, “constrained boundary condition” (CBC). The stiffness and mass matrices become:

    \[\begin{split}\mathbf{K} = \begin{bmatrix} \tilde{\mathbf{K}} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \quad \mathbf{M} = \begin{bmatrix} \tilde{\mathbf{M}} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix}\end{split}\]

    The dimension of the problem remains unchanged but it is necessary to withdraw the participation of the constrained dofs with the components of the initial matrices (\(\tilde{\mathbf{K}}\) is obtained starting from \(\mathbf{K}\) by eliminating the rows and columns of the dofs which are constrained; same for \(\tilde{\mathbf{M}}\)).

When boundary conditions are imposed, the number of eigenvalues ​​(with all their multiplicities) actually involved in the physics of the phenomenon is therefore less than the size \(n\) of the transformed problem:

  • \(n_{\text{dof - active}} = n - \frac{3 p'}{2}\) with \(p' = 2 p\) (double dualization),

  • \(n_{\text{dof - active}} = n - p\) (kinematic constraint).

The list below shows the display dedicated to these parameters in the message file.

----------------------------------------------
The method of counting:  STURM.
The total number of D.O.F. is :         11424
The number of D.O.F. of Lagrange is :   672
The number of activ D.O.F. is :         10416
----------------------------------------------

On the other hand, in the modal calculation algorithms, one must make sure that the solutions being to an allowable space. We retrieve the correct solution via auxiliary treatments. When using kinematic constraints (CBC), it is necessary, in the different algorithms and at each iteration, to use a “vector of positions” \(\mathbf{u}_{constraint}\), defined by

  • if the \(i\)-th degree of freedom is not constrained, \(\mathbf{u}_{constraint}(i) = 1\)

  • otherwise \(\mathbf{u}_{constraint}(i) = 0\)

By pre-multiplying each manipulated vector \(\mathbf{u}^1(i) = \mathbf{u}^0(i)\cdot\mathbf{u}_{constraint}(i) \quad (i = 1 \dots n) \rightarrow \mathbf{u}^1\) we recover the displacement.

This trick introduces the constraint in the algorithm and implicitly directs the search of solution in the admissible space. In the same way, if one uses the method of double dualization, one needs a vector of the Lagrange degrees of freedom \(\mathbf{u}_{lagr}\) defined in the same way as \(\mathbf{u}_{constraint}\). This is only used when choosing the initial random vector. For the vector \(\mathbf{u}^0\), the boundary conditions (LBC) are applied as follows:

\[\mathbf{u}^1(i) = \mathbf{u}^0(i) \cdot \mathbf{u}_{lagr}(i)\quad (i = 1 \dots n) ~,~~ \mathbf{K}\mathbf{u}^2 = \mathbf{u}^1\]

On the other hand, we often include the additional constraint that this initial vector belong to the image set of the working operator. That makes it possible to enrich modal computation more quickly by not being limited to the kernel. Thus, in the case of Lanczos and IRAM, one will take as initial vector, not the previous \(\mathbf{u}^2\), but \(\mathbf{u}^3\) such that

\[\mathbf{u}^3 = (\mathbf{K} -\sigma\mathbf{M})^{-1}\mathbf{M}\mathbf{u}^2\]

Subsequently, to simplify the notations, we will not distinguish between the initial matrices and their dualized versions (denoted with a tilde only if necessary). Very often, they will be designated by \(\mathbf{A}\) and \(\mathbf{B}\) in order to get closer to the usual modal notation without being attached to any particular class of problems.

13.1.3.3. Matrix properties

In the case (the most common) where the matrices considered are symmetric and with real coefficients, a list of possible scenarios is described in the table below. Matrices can be defined as positive (denoted \(> 0\) ), positive semi-definite (\(\ge 0\)), indefinite (\(\le 0 or \ge 0\)) or even singular (S).

Table 13.1.3 Properties of the matrices of the GEP.

Free structure

Lagranges 9

Buckling

Fluid-structure

\(\mathbf{A}(\mathbf{K})\)

\(\ge 0\) and S

\(< 0\) or \(> 0\)

\(> 0\)

\(\ge 0\) and S

\(\mathbf{B}(\mathbf{M}/\mathbf{K}_g)\)

\(> 0\)

\(\ge 0\) and S

\(\le 0\) or \(\ge 0\)

\(> 0\)

The columns of this table being mutually exclusive, in practice, a buckling problem using double Lagranges to model some of its boundary conditions, sees its dualized matrices (\(\mathbf{K}\) and \(\mathbf{K}_g\)) become potentially undefined.

This range of properties must be taken into account when choosing the work operator and the scalar product combination. This framework can thus reinforce, with efficiency and transparency, the robustness and the scope of the modal calculation algorithm in all the cases encountered by Code_Aster.

If the matrices are complex symmetric or real unsymmetric, one can no longer define the matrix dot product. In that case, only the methods QZ (Section 13.1.9) and IRAM (Section 13.1.7) are available. The first does not need this type of mechanism and, we “bluff” the second by providing it with a “false” matrix scalar product, in fact the usual Euclidean scalar product (Section 13.1.7.5). This last trick is allowed with IRAM, because as an Arnoldi variant, it can work with a unsymmetric combination (work operator, scalar product).

The following paragraphs will allow us to measure the impact of these properties on the spectrum of typical problems.

9

This column relates to the properties of the dualized matrices made up from the initial matrices.

13.1.3.4. Properties of the eigen modes

Let us recall first of all that if the matrix of SEP, \(\mathbf{A} \mathbf{u} = \lambda \mathbf{u}\) is real symmetric, then its elements are real. The eigen elements of a matrix are its eigenvalues ​​and vectors. On the other hand, \(\mathbf{A}\) being normal, its eigenvectors are orthogonal.

In the case of the GEP, \(\mathbf{A} \mathbf{u} = \lambda \mathbf{B}\mathbf{u}\),and this condition is not sufficient. So, consider the following problem:

\[\begin{split}\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = \lambda \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix}\end{split}\]

The eigensolution is

\[\begin{split}\lambda_{\pm} = \tfrac{1}{2}(1 \pm i \gamma) \quad \text{and} \quad \mathbf{u}_{\pm} = \frac{1}{\sqrt{1+\lambda_{\pm}^2}}\begin{bmatrix} -\lambda_{\pm}\\ 1 \end{bmatrix}\end{split}\]

If we add the hypothesis that one of the matrices \(\mathbf{A}\) or \(\mathbf{B}\) is positive definite, the generalized problem has real solutions. We even have the following more precise characterization (sufficient condition).

Theorem 1

Let \(\mathbf{A}\) and \(\mathbf{B}\) be two real symmetric matrices. If there exists \(\alpha \in \mathbb{R}\) and \(\beta \in \mathbb{R}\) such that \(\alpha\mathbf{A} +\beta\mathbf{B}\) is positive definite, then the generalized eigenvalue problem has its real eigen elements.

Proof: This result is obtained by multiplying the generalized problem by \(\alpha\) and by carrying out a spectral shift \(\beta\). We then obtain the problem

\[(\alpha \mathbf{A} +\beta \mathbf{B}) \mathbf{u} = (\lambda\alpha+\beta ) \mathbf{B}\mathbf{u}\]

As we defined \(\alpha \mathbf{A} +\beta \mathbf{B}\) as positive, it admits a unique Cholesky decomposition of the form \(\mathbf{C}\mathbf{C}^T\) with \(\mathbf{C}\) regular. The problem is then written

\[\mathbf{C}^{-1}\mathbf{B}\mathbf{C}^{-T} \mathbf{z} = \mu \mathbf{z} \quad\text{with}\quad \mathbf{z} = \mathbf{C}^T \mathbf{u} \quad\text{and}\quad \mu =\frac{1}{\alpha \lambda+\beta}\]

leading to the proof because the matrix \(\mathbf{C}^{-1}\mathbf{B}\mathbf{C}^{-T}\) is symmetric.

Note

  • This characterization is not necessary, the generalized problem associated with the matrices \(\mathbf{A} = \text{diag}(1, -2, -1)\) and \(\mathbf{B} = \text{diag}(-2,1,1)\) has a real spectrum though the matrices are not positive definite.

  • In the case of complex and Hermitian matrices, this theorem remains valid. On the other hand, in complex non Hermitian (e.g., hysteretic damping), the eigen modes can be complex: eigen vectors with complex components and any real or complex eigenvalues.

  • For real unsymmetric matrices, the eigen modes can be complex: eigenvectors with complex components and eigenvalues real or complex \((\lambda , \tilde{\lambda})\).

Proposition 2

If the matrices \(\mathbf{A}\) and \(\mathbf{B}\) are real and symmetric, the eigenvectors of the generalized problem are orthogonal to \(\mathbf{A}\) and \(\mathbf{B}\), which means that they satisfy the relations

\[\begin{split}\begin{cases} \mathbf{u}_i^T \mathbf{B} \mathbf{u}_j &= \delta_{ij} a_j \\ \mathbf{u}_i^T \mathbf{A} \mathbf{u}_j &= \lambda_j \delta_{ij} a_j \\ \end{cases}\end{split}\]

where has \(a_j\) is a scalar depending on the norm of the j-th eigenvector, \(\delta_{ij}\) is the Kronecker symbol and \(\mathbf{u}_j\) is the eigenvector associated with the eigenvalue \(\lambda_j\).

Proof: Immediate for distinct eigenvalues, writing the \(\mathbf{A}\) and \(\mathbf{B}\) - dot product between two pairs \(( i, j )\) and \(( j, i )\), then using the symmetry of the matrices (cf 52).

Note

  • We show that the orthogonality of the eigenvectors of \(\mathbf{A}\) and \(\mathbf{B}\) are a consequence of the hermiticity of the matrices. They are a generalization of the properties of the standard Hermitian problem : in the case of a matrix with complex and Hermitian coefficients, the scalar product to be considered is a Hermitian product.

  • The orthogonality relative to the matrices does not mean that the eigenvectors are orthogonal for the classical Euclidean norm. This can only result for specific symmetries (cf. TP no 1 76).

  • This property simplifies the calculations of modal recombinations (DYNA_TRAN_MODAL 82), when one handles generalized matrices of stiffness and mass which are diagonal. The quantities \(k_j = \lambda_j a_j\) and \(m_j = a_j\) are called, respectively, modal stiffness and modal mass for the j-th mode.

  • For non-Hermitian matrices, Theorem 1 is no longer satisfied.

  • Knowing that the modes are often real, we will now be concerned with their estimation.

13.1.3.5. Real spectrum estimate

Document [R5.01.04] is dedicated to this problem involving the modal operators: INFO_MODE and CALC_MODES.

Let us just recall that in the most frequent case of real modes (symmetric real GEP), the problem of counting of eigenvalues ​​is solved using the Sturm test (Section 13.1.2.2 / Section 13.1.3.2). The situation is much less favorable when the spectrum resides in the complex plane (complex or unsymmetric GEP and QEP). In this case, only the operator INFO_MODE has an implemented method: the APM method (Section 13.1.2.3 / Section 13.1.3.3). But because of its enormous calculation costs and its innovative character, it is advisable to reserve it, for the moment, to simple problems of small size (< 10000 degrees of freedom).

Now that we are able to count the spectrum of the GEP, it remains to be built! The generic algorithms being intended for SEP, we must transform our initial problem.

13.1.3.6. Spectral transformation

These techniques make it possible to meet a threefold objective:

  • Identify a SEP,

  • Orient the search for the spectrum,

  • Separate the eigenvalues.

The spectral calculation algorithms converge better if the (working) spectrum that they process is adequately separated. These techniques can therefore be considered as preconditioning of the starting problem and make it possible to make the separation of certain modes much greater than those of other modes, and thus improve their convergence.

The most widespread of these transformations is the technique known as “shift and invert” which consists in working with an operator \(\mathbf{A}_\sigma\) such that

\[\mathbf{A}\mathbf{u} = \lambda \mathbf{B}\mathbf{u} \implies \underbrace{(\mathbf{A}-\sigma\mathbf{B})^{-1}\mathbf{B}}_{\mathbf{A}_\sigma} \mathbf{u} = \frac{1}{\underbrace{\lambda - \sigma}_{\mu}}\mathbf{u}\]
../../_images/fig9_r5.01.01.png

Fig. 13.1.9 Effect of “shift and invert” on the separation of eigenvalues.

The figure shows that this separation the spectrum in \(\mu\) are due to particular properties of the hyperbolic function. On the other hand, we observe that only the eigenvalues are affected by the transformation. At the end of the modal process, it is therefore sufficient to go back to the plane of \(\lambda\) by an appropriate change of variable.

Note

  • The variable \(\sigma\) is usually designated as “shift” or “spectral shift”.

  • The matrix \(\mathbf{A}_\sigma\) must be invertible; this can also be one of the reasons for the shift (Section 13.1.6.5).

Note that with a complex shift several scenarios are possible:

  • Work completely in complex arithmetic,

  • In real arithmetic, by isolating the real and imaginary parts of \(\mathbf{A}_\sigma\), for example via operators

    \[\begin{split}\mathbf{A}^{+}_\sigma & = \text{Re} (\mathbf{A}_\sigma) \implies \mu^{+} = \tfrac{1}{2} \left( \frac{1}{\lambda-\sigma} + \frac{1}{\lambda - \tilde{\sigma}} \right)\\ \mathbf{A}^{-}_\sigma & = \text{Im} (\mathbf{A}_\sigma) \implies \mu^{-} = \tfrac{1}{2i} \left( \frac{1}{\lambda-\sigma} - \frac{1}{\lambda - \tilde{\sigma}} \right)\end{split}\]

Each approach has its advantages and disadvantages. For the QEP of Code_Aster 79, the first approach is used for Sorensen (METHODE = ‘SORENSEN’ + APPROCHE = ‘COMPLEX’). The second is reserved for the other approaches: METHODE = ‘SORENSEN’ or ‘TRI_DIAG’ + APPROCHE = ‘REEL’ / ‘IMAGINAIRE’).

Note

  • This choice of the operator is inseparable from that of the (pseudo)-scalar product. It allows you to choose an algorithm and can thus influence the robustness of computation.

  • Other classes of spectral transformations exist. For example that of Cayley, with a double shift \((\sigma_1, \sigma_2)\), used to select the eigenvalues ​​located to the right of a vertical axis

    \[\underbrace{(\mathbf{A}-\sigma_1\mathbf{B})^{-1}(\mathbf{A}-\sigma_2\mathbf{B})}_{\mathbf{A}_\sigma} \mathbf{u} = \underbrace{\frac{\lambda-\sigma_2}{\lambda-\sigma_1}}_{\mu} \mathbf{u}\]

The following section will summarize the above in the overall problem solving flow chart for generalized modal analysis in Code_Aster.

13.1.3.7. Implementation in Code_Aster

The course of a modal computation in Code_Aster can be broken down into four phases.

13.1.3.7.1. Determination of the shift

The first operation consists in determining the shift as well as certain parameters of the problem. This is carried out in a more or less transparent way according to the calculation option chosen by the user:

  • OPTION = ‘SEPARE’ or ‘AJUSTE’ => the shift is determined by the first phase of the algorithm and the number of eigen modes sought by frequency interval (provided by the Sturm criterion) is bounded by NMAX_FREQ.

  • OPTION = ‘PROCHE’ => the shift is fixed by the user and the number of eigen modes is equal to the number of shifts.

  • OPTION = ‘PLUS_PETITE’ => the shift is null and the number of modes is parameterized by NMAX_FREQ.

  • OPTION = ‘BANDE’ => the shift is equal to the middle of the interval fixed by the user and the number of modes is determined by the Sturm criterion.

  • OPTION = ‘CENTRE’ => the shift is fixed by the user and the number of eigen modes is parameterized by NMAX_FREQ.

Note

  • This phase does not concern the QZ algorithm. Indeed, the latter calculates the whole spectrum of the GEP. It is just after the modal calculation that one takes into account the wishes of the user (OPTION = ‘BANDE’ or ‘CENTRE’ or ‘PLUS_PETITE’ or ‘PLUS_GRANDE’ or ‘ALL’) to select the desired modes.

  • For a GEP with complex modes (\(\mathbf{K}\) complex symmetric or unsymmetric matrices), only the options ‘CENTRE’, ‘PLUS_PETITE’, ‘PLUS_GRANDE’, ‘TOUT’ are available (with Sorensen and QZ).

13.1.3.7.2. Factoring of dynamic matrices

In the second phase of preprocessing, one factorizes 77 the “dynamic” matrices of the type

\[\mathbf{Q}(\sigma) := \mathbf{A} - \sigma \mathbf{B}\]

That is, we decompose the dynamic matrix \(\mathbf{Q}(\sigma)\) into a product of particular matrices (triangular, diagonal) which are easier to manipulate to solve linear systems of the type \(\mathbf{Q}(\sigma) \mathbf{x} = \mathbf{y}\). For symmetric matrices, the decomposition is of the form \(\mathbf{Q}(\sigma) = \mathbf{L}\mathbf{D}\mathbf{L}^T\). For unsymmetric matrices, we have \(\mathbf{Q}(\sigma) = \mathbf{L}\mathbf{U}\), with \(\mathbf{U},\mathbf{L}\) and \(\mathbf{D}\), respectively, upper triangular, lower triangular, and diagonal.

Once this numerical factorization (expensive) has been carried out, the solution of other linear systems comprising the same matrix but with a second different vector is very fast.

This scenario is satisfied:

  • Each time one invokes the Sturm test. For certain scenarios of Phase 1 of selection of the shift (OPTION = ‘SEPARE’ / ‘AJUSTE’ / ‘BANDE’) and in anticipation of Phase 4 of post-checking during the counting of the real eigenvalues ​​(GEP with real and symmetric matrices only).

  • When one must handle a matrix comprising an inverse. For example in the case of Lanczos or IRAM (OPTION [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘ALL’] + SOLVEUR_MODAL = _F (METHODE = ‘TRI_DIAG’ / ‘SORENSEN’)), one is interested in \(\mathbf{A}_\sigma = (\mathbf{A} -\sigma \mathbf{B})^{-1}\mathbf{B}\). With the method of the inverse powers (OPTION = ‘SEPARE’ / ‘AJUSTE’ / ‘PROCHE’), it is \(\mathbf{A}_\sigma = (\mathbf{A}-\sigma \mathbf{B})^{-1}\). The other two approaches, Bathe & Wilson or QZ (OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’] + SOLVEUR_MODAL = _F (METHOD = ‘JACOBI’ / ‘QZ’)) are not concerned with this preliminary factorization of a matrix.

Note

  • This factorization also has the same hazards as the Sturm criterion when the shift is close to an eigenvalue. One then proceeds to the same shifts according to algorithm 1 of R5.01.04.

  • When this phase of preprocessing implements algorithm 5 of R5.01.04 and that NMAX_ITER_SHIFT is reached, the computation produces an alarm if it is about a factorization for the Sturm test (risk of numerical instabilities), a fatal error if one seeks to factorize the work matrix (risk of false results).

  • When one uses the QZ algorithm on an unsymmetric or complex GEP, one does not need factorize of dynamic matrix. This distinction makes it possible to save a little calculation complexity, the QZ algorithm being already quite expensive in itself!

  • By chaining operators INFO_MODE [U4.52.01] and CALC_MODES, one can better direct / calibrate the spectral computation even to limit the extra cost which the initial Sturm test of the option ‘BANDE’ involves (cf. keyword TABLE_FREQ / CHAR_CRIT [U4.52.02]).

13.1.3.7.4. Postprocessing checks

This last part discusses the postprocessing check on the progress of the computation. They are of two types:

  • Norm 10 of the residual of the initial problem

    Algorithm 1. Test of the norm of the residual

    \[\mathbf{u} \leftarrow \frac{\mathbf{u}}{\| \mathbf{u} \|_\infty}\]

    If \(|\lambda|\) > SEUIL_FREQ then

    \[\frac{\| \mathbf{A}\mathbf{u} - \lambda\mathbf{B}\mathbf{u}\|_2}{\mathbf{A}\mathbf{u}\|_2} ? < SEUIL\]

    Else

    \[\| \mathbf{A}\mathbf{u} - \lambda \mathbf{B}\mathbf{u} \|_ 2 ? < SEUIL\]

    End if

    This sequence is parameterized by the keywords SEUIL and SEUIL_FREQ, belonging respectively with the keyword factor VERI_MODE and CALC_FREQ. On the other hand, this post-processing is activated by initialization with OUI (default) of STOP_ERREUR in the keyword factor VERI_MODE. When this rule is activated and not respected, calculation stops, otherwise the error is just indicated by an alarm. Of course, we cannot recommend highly enough not to deactivate this pass-through parameter.

  • Counting eigenvalues

    This postprocessing step is set up only in the case of symmetric real matrices and it is not activated by default for OPTION = ‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’. In this context, it checks that the number of eigenvalues ​​contained in a test interval \([\sigma1, \sigma2]\) is equal to the number detected by the algorithm. This counting procedure is activated in two stages: determination of the test interval (cf. figure and algorithm below) then computation of Sturm proper. This Sturm test calculation is the object of a specific documentation (cf [R5.01.04]), one thus details here only the preliminary stage which is specific to postprocessing in CALC_MODES with OPTION = ‘BANDE’ / ‘CENTRE’ / ‘PLUS_*’ / ‘TOUT’.

    The inclusion of the initial interval 11 in \([\sigma_i, \sigma_f]\) is carried out in order to detect possible problems of clusters or multiplicities at the initial bounds.

    ../../_images/fig10_r5.01.01.png

    Fig. 13.1.10 Counting of the eigenvalues ​​by the Sturm test for symmetric real GEP.

    Noting \(\sigma_i^+\) and \(\sigma_f^-\), respectively, the largest and the smallest eigenvalues not requested by the user and including the initial interval (figure above), one has the following algorithm for construction of the test interval:

    Algorithm 2. Construction of the test interval for the Sturm test

    If \(\sigma_i^+\) exists (or \(\sigma_f^-\)) and

    If \(\frac{|\sigma_i^+ - \sigma_i|}{|\sigma_i|}\) < PREC_SHIFT (or \(\frac{|\sigma_f^- - \sigma_f|}{|\sigma_f|}\))

    Then \(\sigma_1 = \frac{\sigma_i^+ + \sigma_i}{2}\) (or \(\sigma_2 = \frac{\sigma_f^- + \sigma_f}{2})\)

    Else \(\sigma_1 = \sigma_i ( 1 - \text{sign}(\sigma_i)\) PREC_SHIFT ) (or \(\sigma_2 = \sigma_f ( 1 + \text{sign}(\sigma_f)\) PREC_SHIFT ))

    End if

This test interval construction sequence, via the calculated and unrequested modes, is only active for the methods of Lanczos and Bathe & Wilson of CALC_MODES with OPTION = ‘BANDE’ / ‘CENTRE’ / ‘PLUS_*’ / ‘TOUT’. With the methods of Sorensen and QZ, one selects only the requested modes and therefore the limits of the test interval are determined by the second part of Algorithm 2. This algorithm is parameterized by the keyword PREC_SHIFT of the keyword factor VERI_MODE.

The boxes below display the trace of error messages when these post-checks are triggered.

       A posteriori verification of the modes
<E> <ALGELINE5_15>
for the concept MODE_B_1, the mode number 5
frequency 84.821272
has an error norm of 0.000102 greater than the admissible threshold 0.000001.

       A posteriori verification of the modes

<E> <ALGELINE5_23>
for the MODET concept,
in the interval [28.987648, 5071.142099]
theoretically there are 6 natural frequency (s) ()
and we calculated 5.
This problem can appear when there are multiple modes (structure with symmetries) or a high modal density.

The counting procedure can be deactivated via the keyword VERI_MODE / STURM. Likewise, when at least one of the two post-verification steps (residual norm or counting procedure) is not respected, the continuation of the process is controlled by STOP_ERREUR (if ‘OUI’ the computation stops, if not the error is just indicated by an alarm). Of course, we cannot recommend highly enough not to disable these ‘pass-thru’ settings!

Note

  • When the improvement procedure (keyword AMELIORATION = ‘OUI’) is used with the options ‘BANDE’, ‘CENTRE’, ‘PLUS_*’ or ‘TOUT’, the Sturm test is not carried out at the end of the first step of solution (by a subspace method). This is unnecessary, since these modes are likely to be modified in the second stage (with a few iterations of the inverse powers method). It is only after this step that all the checks: residual, Sturm test and check of belonging to the initial interval (if OPTION = ‘BANDE’) are performed.

  • The Sturm test can also be activated with the options ‘PROCHE’, ‘SEPARE’ and ‘AJUSTE’. It does not fulfill quite the same function in those cases. Indeed, these options are used to refine initial estimates of eigenvalues ​​and to calculate the associated eigenvectors: the presence of ‘gaps’ in the calculated spectrum is then possible and allowed. On the other hand, the other options (‘BANDE’ etc.) must calculate all the specified modes, without missing any: no gap in the spectrum is then tolerated.

    With ‘PROCHE’, ‘SEPARE’ and ‘AJUSTE’, the emission of alarm or error messages therefore does not reveal necessarily a numerical problem but can be normal according to the setting of the data. So if we specify FREQ = (100.0,300.0), and the ‘PROCHE’ option refines these values ​​to (99.0, 299.0) but the spectrum of the mechanical problem normally also comprises the value 200.0, the Sturm test (if it is active) will stop computation with an error message though the absence of this value 200.0 is normal. The user did not ask to refine it! With the initial setting FREQ = (100.0, 200.0, 300.0) this alarm would no longer appear. In this context, the activation of the Sturm test is thus to be handled with caution. Functionally, it is used for the deferral of verification tests when the keyword AMELIORATION is activated at the end of a subspace method.

10

Remember that \(\| \mathbf{u} \|_\infty = \text{max}_i \| \mathbf{u}_i\ |\), \(\| \mathbf{A} \|_\infty = \text{max}_i \sum_j| A_{ij}|\), \(\| \mathbf{u} \|_2 = ( \sum_i | \mathbf{u}_i |^2 )^{1/2}\) and \(\| \mathbf{A} \|_2 = \text{max}_i | \lambda_i (\mathbf{A}\mathbf{A}^\star) |^{1/2}\) (\(= \text{max}_i | \lambda_i(\mathbf{A}) |\) if \(\mathbf{A}\) is Hermitian)

11

In option ‘BANDE’ these are the values ​​indicated by the user, otherwise they are the extreme values ​​of thecalculated spectrum.

13.1.3.7.5. Display in the message file

Information relating to the selected modes is mentioned in the message file. For example, in the most current case of a real symmetric GEP (real eigen modes), one specifies the modal solver used, the list frequencies \(f_i\) retained (FREQUENCY) and their error norm (ERROR NORM, cf. Algorithm 1). An example from the benchmark test forma11a is shown below.

------------------------------------------------------------------------
  Calcul modal : Méthode globale de type QR
                 Algorithme QZ_SIMPLE

numéro    fréquence (HZ)     norme d'erreur
   1       1.67638E+02        3.57918E-11
   2       1.67638E+02        3.84429E-11
   3       1.05060E+03        1.32150E-12
   4       1.05060E+03        1.15921E-12
   5       1.48054E+03        1.41808E-13
   6       2.59704E+03        5.70995E-14
   7       2.94237E+03        4.07506E-13
   8       2.94237E+03        3.61074E-13
   9       4.47822E+03        4.61687E-14
  10       5.76991E+03        6.96812E-14
  11       5.76991E+03        1.00502E-13

Norme d'erreur moyenne   :  7.08175E-12
------------------------------------------------------------------------

For a GEP with complex modes (real unsymmetrical or complex symmetric 12), contrary to QEP 79, Code_Aster does not filter the combined eigenvalues ​​of any complexes. It keeps all of them and displays them in ascending order of the real part. Thus, the columns FREQUENCY and AMORTIZATION group together, respectively, \(Re(f_i)\) and \(\frac{Im(f_i )}{2 Re(f_i)}\). An example from the benchmark test zzzz208a is shown below.

------------------------------------------------------------------------
   MODAL CALCULATION: SIMULTANEOUS ITERATION METHOD
   SORENSEN METHOD

   NUMBER FREQUENCY (HZ) DAMPING     ERROR NORM
   1      5.93735E+02    5.00000E-03 2.73758E-15
   2      9.45512E+02    5.00000E-03 1.64928E-15
   3      3.51464E+03    5.00000E-03 2.14913E-15
   .....

Now that the context of the GEPs in Code_Aster has been skimmed, we will be interested more particularly with the power methods and their implementation in the operator CALC_MODES with OPTION = [‘PROCHE’, ‘AJUSTE’, ‘SEPARATE’].

12

The cumulation of the two, complex unsymmetrical or any complex, is not currently treated.

13.1.4. Method of inverse powers (CALC_MODES with OPTION = [‘PROCHE’, ‘AJUSTE’, ‘SEPARATE’])

13.1.4.1. Introduction

This method is accessible in Code_Aster only for matrices with real coefficients.

To calculate several eigenvalues ​​of the generalized eigenvalue problem, one does not use the general method of inverse powers as it is. We can, for example, combine it with a “deflation” technique 13 so as to automatically filter the updated spectral information so as not to search for it in the next iteration. With “deflation with restriction” 14) of Wielandt, one builds the work operator iteratively (in the symmetric case), by using the mode \((\mu_k , u_k)\) to create a filtered matrix

\[\mathbf{A}_{k+1} = \mathbf{A}_k -\mu_k \mathbf{u}_k \mathbf{u}_k^T.\]

This strategy, which does not consider multiplicity of eigenvalues, must be supplemented by a Sturm criterion. On the other hand, finite precision arithmetic and the requirement of not having to rebuild the matrix \(\mathbf{A}_k\) at every iteration, requires efficient orthogonalization algorithms.

These complications led to the choice of another approach that consists of two parts:

  • The localization of the eigenvalues (determination of an approximate value of each eigenvalue contained in a given interval by a secant bisection technique).

  • The improvement of these estimates and the computation of their associated eigenvectors by an inverse iteration method.

The search for an approximate value for each eigenvalue is selected in the keyword factor CALC_FREQ via an OPTION:

  • If OPTION = ‘SEPARE’, in each interval of frequencies defined by the keyword FREQ, an approximate value of each eigenvalue contained in the interval is calculated using the bisection method (Section 13.1.4.2.1).

  • If OPTION = ‘AJUSTE’, we first perform the operations for SEPARE and then, starting from the approximations, we refine the result by the secant method.

  • If OPTION = ‘PROCHE’, the frequencies given by the keyword FREQ are considered as the approximate values ​​of the eigenvalues ​​sought.

For the first two options, we calculate the modal position of each eigenvalue simultaneously. This allows us to detect multiple identical modes. We retain either only the NMAX_FREQ frequencies ​​contained in the maximum interval specified by the user, or we calculate all the values in this interval (if NMAX_FREQ = 0).

Note

  • Of course, as we have already specified, these options are to be used only to determine or refine some eigenvalues. For a more extensive search, use the other search options.

  • With the option ‘PROCHE’, one cannot calculate multiple modes.

  • It is an expensive algorithm because it performs the Sturm test multiple times and thus has to perform multiple matrix factorizations.

  • The limits of the search intervals are provided by FREQ or CHAR_CRIT according to the initialization of TYPE_RESU.

We will now detail the different algorithms (and their settings) which are set up int he first part of the process.

13

Filtering technique consisting of transforming the work operator in such a way that it has the same eigenvalues ​​except at certain predefined null modes.

14

Other types of deflation exist, such as the Ducan-Collar method, which uses the first row of the matrix and the eigenvector to filter spectral information. These vector techniques generalize well in blocks.

13.1.4.2. Localization and separation of eigenvalues

13.1.4.2.1. Bisection method

A corollary Sylvester’s law of inertia makes it possible to determine the number of eigenvalues ​​contained in a given interval by performing two \(\mathbf{L}\mathbf{D}\mathbf{L}^T\) decompositions. This criterion can be used to refine the interval until it contains only one eigenvalue. After this control has been set up, we iterate using the bisection principle.

See [R5.01.04] for more details and for the definition of the “modal position” (called “pm” in Code_Aster)

Starting with the interval \([\lambda_a , \lambda_b]\), we operate as follows:

Algorithm 3: Bisection method

Calculate \(\lambda^\star = \tfrac{1}{2}(\lambda_a + \lambda_b)\)

Calculate the modal position \(\text{pm}(\lambda^\star)\)

If \(\text{pm}(\lambda^\star) = \text{pm}(\lambda_a)\) [or \(\text{pm}(\lambda_b)\) ] then

Restart with the interval \([\lambda^\star, \lambda_b]\) [or \([\lambda_a, \lambda^\star]\) ]

Else

Restart with the intervals \([\lambda_a, \lambda^\star]\) and \([\lambda^\star, \lambda_b]\)

End If

We stop the process if we subdivide the starting interval more than NMAX_ITER_SEPARE times, or if for a given interval, we have

\[\left(\frac{|\lambda_a - \lambda_b|}{\lambda^\star}\right) \le \text{PREC_SEPARE}\]

In the latter case, we no longer refine the search in this interval). We finally get a list of frequencies. In each interval defined by the arguments of this list, is an eigenvalue with a certain multiplicity. As an approximation of this eigenvalue, we take the middle of the interval.

Note

  • We could have used the change of sign of the characteristic polynomial as a criterion, but in addition to the fact that it is very expensive to evaluate, it cannot detect multiplicities

  • The initialization of the process can be carried out empirically according to the needs of the user. To encompass a part of the spectrum one can also use the regions of the complex plane using Gerschgörin-Hadamard theorems (on \(\mathbf{A}, \mathbf{A}^T\), …). In this perspective, the bisection method can be more efficient than a QR in the presence of a cluster. Its convergence, although linear, is in effect increased by 1/2 whereas that of QR can tend towards 1 57.

../../_images/fig11_r5.01.01.png

Fig. 13.1.11 Bisection method

13.1.4.2.2. Secant method

The secant method is a simplification of the Newton-Raphson method. At any stage \(k\), knowing a value \(\lambda_k\) and by approximating the nonlinear function \(p(\lambda)\) by its tangent at this point, we determine the next value \(\lambda_{k+1}\) as being the intersection of this line with the \(\lambda\) -axis, and so on, according to the iterative scheme:

\[\lambda_{k+1} = \lambda_k - p(\lambda_k) \frac{\lambda_k - \lambda_{k-1}}{p(\lambda_k) - p(\lambda_{k-1})}\]
../../_images/fig12_r5.01.01.png

Fig. 13.1.12 Secant method.

The tangent is approximated by a finite difference so as not to have to calculate a derivative of \(p(\lambda)\); only the estimation of the polynomial is required. We consider that we have reached convergence when

\[\frac{|\lambda_{k+1} - \lambda_k |}{\lambda_k} < \text{PREC_AJUSTE}\]

and, moreover, one is limited to NMAX_ITER_AJUSTE iterations if this criterion is not attained.

Note

  • This method has an almost quadratic convergence when it is close to the solution, otherwise it may diverge. Hence we combine the bisection method with this approach. This strategy of combining bisection and secant was initiated by Van Wijngaardeen and Dekker (1960).

  • If one does not restart with the calculation of the last point in the formula of the secant (\(\lambda_{k+1} = \lambda_0\)), we get the classic method of “false position”. 15 This converges less quickly (while requiring an iterative function evaluation), but it is robust and stable (like the bisection method). We are guaranteed to improve the initial assessment.

  • We show that the secant method achieves the best compromise “speed of convergence / number of function evaluations per iteration”.

  • The secant method is based on a linear interpolation between the last two iterations \((\lambda_{k+1}, p(\lambda_{k+1}))\) and \((\lambda_k, p(\lambda_k))\). By generalizing this process to a higher order of interpolation, i.e. the last three iterated, we get the Müller-Traub method used for the QEP 79.

We will now detail the inverse powers algorithm (coupled with a Rayleigh acceleration) constituting the second part of the process.

15

Method also called “regula falsi” (in Latin). Very old method of finding zero of functions, probably invented by the Hindus mathematicians of the Vth century, and re-discovered in Europe by L.Fibonnacci in the XII th century.

13.1.4.3. Inverse powers method

13.1.4.3.1. Principle

To determine the eigenvalue of the generalized problem \(\mathbf{A} \mathbf{u} = \lambda \mathbf{B} \mathbf{u}\) closest in modulus to the shift \(\sigma\), one applies the method of powers to the operator \((\mathbf{A} - \sigma \mathbf{B})^{-1}\mathbf{B}\). In fact, we only build the factored matrix 16 \(\mathbf{A}^\sigma = (\mathbf{A} - \sigma \mathbf{B})^{-1}\) and this amounts to dealing with the generalized problem

\[(\mathbf{A}^\sigma)^{-1} \mathbf{B} \mathbf{u} = \frac{1}{\underbrace{\lambda - \sigma}_{\mu}}\mathbf{u}\]

The power method converges towards eigenvalues ​​of higher modulus (in \(\mu\)). We can therefore capture the \(\lambda\) closest to the shift. The principle is the following: knowing an estimate \(\sigma\) of the required eigenvalue and starting from a normalized initial vector \(\mathbf{y}_0\), we build a sequence of approximate eigenvectors \((\mathbf{y}_k)\) with a recurrence relation.

Algorithm 4: Method of inverse powers

For \(k = 1, \dots\) do

\(\mathbf{A}^\sigma \tilde{\mathbf{y}}_k = \mathbf{B} \mathbf{y}_{k-1}\)

\(\mu_k = \| \tilde{\mathbf{y}}_k \|\)

\(\mathbf{y}_k = \cfrac{\tilde{\mathbf{y}}_k}{\mu_k}\)

End For

With \(\lambda 1, \lambda 2, \dots\) as the first eigenvalues ​​(with eigenvectors \(\mathbf{u}_i\)) the closest in modulus to \(\sigma\), we have a linear convergence of \(\mathbf{y}_k \rightarrow \mathbf{u}_1\) and a quadratic convergence of

\[\mu_k \rightarrow \frac{1}{|\lambda_1 -\sigma|} \quad \text{and} \quad \frac{\tilde{\mathbf{y}}_k(i)}{\tilde{\mathbf{y}}_{k-1}(i)} \rightarrow \frac{1}{|\lambda_1-\sigma|} \quad (i = 1 \dots n)\]

The convergence factor of these sequences is the order of

\[\frac{|\lambda_1 -\sigma|}{\text{min}_{i \ne 1} |\lambda i -\sigma|}\]

Note

  • This result is acquired (except for the convergence factor) when the dominant eigenvalue (here the one closest to the shift mod) is unique. Otherwise, results are possible even in the complex framework, but the rigorous analysis of the convergence of this algorithm is still incomplete 48 85.

  • In theory, these results require that the initial vector is not orthogonal to the proper left subspace desired. In practice, rounding errors avoid this problem (cf. 57 pp500-509).

  • Even if the estimate of the eigenvalue is rough, the algorithm quickly provides a very good estimate of the eigenvector.

  • The major disadvantage of this method is that it is necessary to carry out a factorization for each value to be calculated.

In general one works in Euclidean norm or in infinite norm, but to facilitate post-modal computations one seeks here to \(\mathbf{B}\) -normalize the eigenvectors (when \(\mathbf{B}\) is undefined, one works with the associated pseudo-norm). The basic algorithm can be rewritten by setting \(\mathbf{z}_0 = \mathbf{B} \mathbf{y}_0\):

Algorithm 5: Method of inverse powers with \(\mathbf{B}\) - pseudo-norms

For \(k = 1, \dots\) do

\(\mathbf{A}^\sigma \tilde{\mathbf{y}}_k = \mathbf{z}_{k-1}\)

\(\tilde{\mathbf{z}}_k = \mathbf{B} \mathbf{y}_k\)

\(\rho(\mathbf{y}_k) = \cfrac{{\mathbf{y}}_k^T \cdot \mathbf{z}_{k-1}}{{\mathbf{y}}_k^T \cdot \tilde{\mathbf{z}}_{k}}\)

\(\varepsilon_k = \text{sign}(\mathbf{y}_k^T \cdot \tilde{\mathbf{z}}_k)\)

\(\mathbf{z}_k = \varepsilon_k \cfrac{\tilde{\mathbf{z}}_k}{|\mathbf{y}_k^T \cdot \tilde{\mathbf{z}}_{k}|^{1/2}}\)

End For

Note that \(\rho(\mathbf{y}_k) \rightarrow \frac{1}{|\lambda 1 -\sigma|}\). This approach avoids matrix-vector products with matrix \(\mathbf{B}\) during the computation of the dot products and already foreshadows the Rayleigh coefficient of the following paragraph. We observe that the convergence factor is smaller if the spectral shift \(\sigma\) is close to the required eigenvalue and therefore that \(\mathbf{A} - \sigma\mathbf{B}\) is close to a singularity. This fact is not detrimental to the process because the error made in solving the system is “mostly” in the direction generated by the eigenvector which is the direction sought. Therefore, when solving \(\mathbf{A}^\sigma \mathbf{y}_k = \mathbf{z}_{k+1}\), we do not find the exact solution \(\mathbf{y}_k\) but that rounding errors lead to a close solution of the form \(\tilde{\mathbf{y}}_k = \mathbf{y}_k + \mathbf{w}\). This is proportional to the exact solution, but as the normalization is arbitrary, everything is fine 64.

Note

  • This poor conditioning, far from having an unfavorable effect, even improves the convergence of the algorithm. This algorithm is therefore used to improve the eigenvector associated with the approximate value of phase 1. To refine this estimate of the eigenvalue, we introduce a Rayleigh quotient.

16

This notation should not be confused with that symbolizing the “shift and invert” operation \(\mathbf{A}_\sigma = (\mathbf{A} - \sigma\mathbf{B})^{-1}\mathbf{B}\).

13.1.4.3.2. Rayleigh quotient iteration method

The Rayleigh quotient applied to the generalized problem is defined by the quantity \(R(\mathbf{x})\), where \(\mathbf{x} \in \mathbb{R}^n\) is a non-zero vector such that

\[R(\mathbf{x}) = \frac{\mathbf{x}^T \mathbf{A}\mathbf{x}}{\mathbf{x}^T \mathbf{B}\mathbf{x}}\]

This quotient has the remarkable property of stationarity in the neighborhood of any eigenvector and reaches an extremum (local) which is the corresponding eigenvalue: for each fixed \(\mathbf{x}\), \(R(\mathbf{x})\) minimizes \(\lambda\rightarrow\|( \mathbf{A} -\lambda \mathbf{B})\mathbf{x}\|_2\).

By the requirement “if \(\mathbf{x}\) is an approximation of an eigenvector of the system \(\mathbf{A}\mathbf{x} = \lambda \mathbf{B}\mathbf{x}\) we mean that “\(R(\mathbf{x})\) is an approximation of the eigenvalue associated with the vector \(\mathbf{x}\)”, and conversely, we have seen that if we had a good estimate of an eigenvalue, the method of inverse iterations makes it possible to obtain a good estimate of the corresponding eigenvector. Hence the we can combine these two properties by considering the inverse iteration algorithm with a spectral shift for which one reevaluates, at each iteration, the eigenvalue via the Rayleigh quotient. We then get the “Rayleigh quotient iteration algorithm” (in \(\mathbf{B}\) -pseudo-norm)

Algorithm 6. Method of iterations of the Rayleigh quotient with \(\mathbf{B}\) -pseudo-norm.

For \(k = 1, \dots\) do

\([\mathbf{A} - R(\mathbf{y}_{k-1})]\mathbf{y}_k = \mathbf{z}_{k-1}\)

\(\tilde{\mathbf{z}}_k = \mathbf{B} \mathbf{y}_k\)

\(R(\mathbf{y}_k) = \cfrac{{\mathbf{y}}_k^T \cdot \mathbf{z}_{k-1}}{{\mathbf{y}}_k^T \cdot \tilde{\mathbf{z}}_{k}} + R(\mathbf{y}_{k-1})\)

\(\varepsilon_k = \text{sign}(\mathbf{y}_k^T \cdot \tilde{\mathbf{z}}_k)\)

\(\mathbf{z}_k = \varepsilon_k \cfrac{\tilde{\mathbf{z}}_k}{|\mathbf{y}_k^T \cdot \tilde{\mathbf{z}}_{k}|^{1/2}}\)

End For

For the standard modal problem, we can show 64 that the convergence of this algorithm is cubic in the event that the work operator is normal 17 (a fortiori in the Hermitian case) and at best quadratic in other cases. If we use this method without modifying it, we would need at each iteration of the process an improvement of each eigenvalue, perform a \(\mathbf{L}\mathbf{D}\mathbf{L}^T\) factorization, which would be very expensive. Hence we carry out this Rayleigh shift only if we are in a neighborhood (to be defined) of the solution. 18

Note

  • In a more general framework, certain authors have introduced an algorithm known as “bi-iterations of the Rayleigh quotient”. Based on the stationarity of the bi-quotient

    \[R_b(\mathbf{x}, \mathbf{y}) = \frac{\mathbf{y}^T \mathbf{A} \mathbf{x}}{\mathbf{y}^T \mathbf{B} \mathbf{x}}\]

    in the vicinity of eigenvectors on the right and on the left, it provides (in non-Hermitian matrices) the two types of eigenvectors. Its cost is however prohibitive because it requires twice as many factorizations (cf. BNParlett 1969 64).

17

A matrix \(\mathbf{A}\) is called normal if \(\mathbf{A}\mathbf{A}^\star = \mathbf{A}^\star \mathbf{A}\). This is the case with Hermitian, anti-Hermitian, or unitary operators.

18

This strategy of coupling methods with complementary or even antagonistic characteristics is often used in numerical analysis. For example, in optimization, the Levenberg-Marquadt method couples a steepest-descent and a Newton.

13.1.4.3.3. Implementation in Code_Aster

This spectral shift is activated only if the keyword OPTION of the factor keyword CALC_MODE is initialized to ‘RAYLEIGH’. By default, one has ‘DIRECT’ and the shift is then traditional (global correction rather than progressive eigenvalue). The algorithm set up in the code is divided as follows (in \(\mathbf{B}\) -pseudo norm):

Algorithm 7. Functional diagram of the method of the inverse powers in CALC_MODES with OPTION among [‘PROCHE’, ‘AJUSTE’, ‘SEPARE’].

  • Initialization of the eigenvalue starting from the estimate of the first phase: \(\lambda^\star\).

  • Calculation of a random initial vector \(\mathbf{y}_0\) that satisfies the boundary conditions.

  • \(\mathbf{B}\) -orthonormalization of \(\mathbf{y}_0\) compared to the previously calculated modes (if it is a multiple mode according to the first phase) by a Modified Gram-Schmidt 19 (GSM).

  • Calculation of \(\delta_0 = \text{sign}(\mathbf{y}_0^T \mathbf{B} \mathbf{y}_0)\).

  • For k = 1 , NMAX_ITER do

    • Solve \((\mathbf{A} - \lambda^\star \mathbf{B}) \mathbf{y}_k = \delta_{k-1} \mathbf{B}\mathbf{y}_{k-1}.\) This is \(\mathbf{B}\) -orthonormalization (possible) of \(\mathbf{y}_k\)

    • Calculate the correction of the eigenvalue

      \[c =\frac{\mathbf{y}_k^T \mathbf{B} \mathbf{y}_{k-1}}{\mathbf{y}_k^T \mathbf{B}\mathbf{y}_k}\]

  • If \(||\mathbf{y}_k^T \mathbf{B}\mathbf{y}_{k-1} | - 1| \le\) PREC then

    • \(\lambda^\star = \lambda^\star + c\)

      exit

  • Else

    • If OPTION = ‘RAYLEIGH’ and if \(|c| \le 0.1 \lambda^\star\) then

      • \(\lambda^\star = \lambda^\star + c\)

      • lambda*= lambda*+ c ;End if.End if.End loop.

  • End If

  • End If

  • End For

The maximum admissible error norm PREC and the maximum number of iterations NMAX_ITER are arguments of the keyword factor CALC_MODE.

Note

  • The eigenvector being \(\mathbf{B}\) -normal, we assume that we have reached convergence when the absolute value of the dot product is close to unity.

  • To avoid taking a \(\mathbf{B}\) -orthogonal initial vector to the desired eigenvalue one uses a random selection method for the components of this vector.

  • On the other hand, to be able to determine multiple or close modes, we use a \(\mathbf{B}\) -orthogonalization with the modes that have already been calculated.

  • The “acceleration” of the algorithm by Rayleigh quotient being expensive, it is not used in each iteration but only if one is in the vicinity of the required eigenvalue.

19

Since the development of this operator, more efficient and more efficient and robust orthogonalization methods have spread, such as the Gram-Schmidt Modified Iteratives (IGSM) used in CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’] (see Section 13.1.15 and 64).

13.1.4.4. Scope of use

GEP with symmetric real matrices.

The user can specify the type of calculation (dynamic or buckling) with the keyword TYPE_RESU. Depending to this value, one provides either the vector FREQ or CHAR_CRIT.

13.1.4.5. Display in message file

In the message file, the results are displayed as a table:

------------------------------------------------------------------------
THE NUMBER OF DOF
TOTAL IS:                    192
FROM LAGRANGE IS:            84
THE NUMBER OF ACTIVE DDL IS: 66
------------------------------------------------------------------------
CHECKING THE FREQUENCY SPECTRUM (CONTINUED FROM STURM)
THE NUMBER OF FREQUENCIES IN THE BAND (1.00000E-02, 6.00000E-02) IS 4
------------------------------------------------------------------------
MODAL CALCULATION: INVERSE ITERATION METHOD
BISECTION SECANT INVERSE
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
------------------------------------------------------------------------
POSTERIORI CHECKING OF THE MODES
------------------------------------------------------------------------

With the option ‘PROCHE’, the columns “Bisection” and “Secant” do not appear, while with the option ‘SEPARE’, only the “Secant” column disappears. The last precision column groups intermediate data and does not represent, as for the case of CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’], the error norm of the residual. It is an artifact that will be deprecated.

13.1.4.6. Summary of settings

Let us now recapitulate the parameter setting of the operator CALC_MODES with OPTION = [‘PROCHE’, ‘AJUSTE’, ‘SEPARE’].

Table 13.1.4 Summary of the parameter setting of CALC_MODES with OPTION among[‘PROCHE’, ‘AJUSTE’, ‘SEPARE’] (GEP).

Operand

Keyword

Default value

References

TYPE_RESU = ‘DYNAMIQUE’

‘DYNAMIQUE’

Section 13.1.3.1

‘MODE_FLAMB’

Section 13.1.3.1

OPTION = ‘SEPARE’

‘AJUSTE’

Section 13.1.4.1

‘AJUSTE’

Section 13.1.4.1

‘PROCHE’

Section 13.1.4.1

CALC_FREQ or CALC_CHAR_CRIT

FREQ

Section 13.1.4.1

CHAR_CRIT

Section 13.1.4.1

NMAX_FREQ

0

Section 13.1.4.1

NMAX_ITER_SHIFT

3

81 Section 13.1.3.2

PREC_SHIFT

0.05

81 Section 13.1.3.2

SEUIL_FREQ

1.E-02

Section 13.1.3.7

SOLVEUR_MODAL

NMAX_ITER_SEPARE

30

Section 13.1.4.2

PREC_SEPARE

1.E-04

Section 13.1.4.2

NMAX_ITER_AJUSTE

15

Section 13.1.4.2

PREV_ADJUST

1.E-04

Section 13.1.4.2

OPTION_INV = ‘DIRECT’

‘DIRECT’

Section 13.1.4.3

‘RAYLEIGH’

Section 13.1.4.3

PREC_INV

1.E-05

Section 13.1.4.3

NMAX_ITER_INV

30

Section 13.1.4.3

VERI_MODE

STOP_ERREUR = ‘OUI’

‘OUI’

Section 13.1.3.7

‘NON’

SEUIL

1.E-02

Section 13.1.3.7

Note

  • All the parameters related to the Sturm pre-test (NMAX_ITER_SHIFT, PREC_SHIFT) and with post-checks (SEUIL_FREQ, VERI_MODE) are provided above.

  • Beginners should avoid modifying the default values of most of these parameters. We recommend that the user change only FREQ, CHAR_CRIT, and NMAX_FREQ and leave the rest at their default settings.

13.1.5. Subspace methods (CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’])

13.1.5.1. Introduction

Warning

GEP with complex modes (Section 13.1.3.1) can be handled in Code_Aster only by the IRAM and QZ methods. This section is therefore limited to matrices with real coefficients. Only general principles of subspace methods are discussed.

If we want to calculate only \(p \ll n\) eigen elements of a generalized problem of order \(n\) (e.g, the p smallest eigenvalues ​​or all the eigenvalues ​​in a given interval), it is preferable to use subspace methods. These are based on a Rayleigh-Ritz approach which projects the problem onto a subspace of dimension \(m\) (\(p < m < n\)). The search for eigen-decompositions is then limited to this reduced problem. Robust algorithms (QR or QL for Lanczos and IRAM, Jacobi for the method of Bathe and Wilson) are available for the task.

The projection process is effective if:

  • The projection space is small

  • The space is easy to construct

  • The orthogonal projection is robust 20

  • The projected matrix can be expressed in a canonical form

  • The projection operation gives a good approximation of the part of the initial spectrum of interest

  • Computation and memory costs are minimal and rounding effects are small

We will first present the Rayleigh-Ritz analysis. We will then explore three variants: the Lanczos method, the Sorensen approach (IRA), and that of Bathe and Wilson.

20

In the non-Hermitian case, oblique projections were developed by F. Chatelin 42.

13.1.5.2. Rayleigh-Ritz analysis

Consider the modal problem of order \(n\), \(\mathbf{A}\mathbf{u} = \lambda \mathbf{u}\), and the subspace \(H \subset \mathbb{R}^n\) generated by an orthonormal basis \((\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_m)\). The latter constitutes the orthogonal matrix \(\mathbf{Q}_m\) which can be used to define the projection operator \(\mathbf{P}_m = \mathbf{Q}_m \mathbf{Q}_m^T\).

The Galerkin method used in this analysis solves the problem:

Problem

Find \((\tilde{\lambda}, \tilde{\mathbf{u}}) \in \mathbb{R} \times H\) such that \(\mathbf{P}_m (\mathbf{A} \tilde{\mathbf{u}} - \tilde{\lambda} \tilde{\mathbf{u}}) = \mathbf{0}\)

Alternatively, we can write

Problem

With \(\tilde{\mathbf{u}} = \mathbf{Q}_m\mathbf{x}\), find \((\tilde{\lambda}, \mathbf{x}) \in \mathbb{R} \times \mathbb{R}^m\) such that \(\underbrace{\mathbf{Q}_m^\star \mathbf{A} \mathbf{Q}_m}_{\mathbf{B}_m} \mathbf{x} = \tilde{\lambda} \mathbf{x}\)

The search for the Ritz elements of Ritz \((\tilde{\lambda}, \tilde{\mathbf{u}})\) therefore amounts to finding the eigen modes of the Rayleigh matrix \(\mathbf{B}_m = \mathbf{Q}_m^T\mathbf{A}\mathbf{Q}_m\).

The eigenvalues ​​remain unchanged in the two formulations. But, knowing an eigenvector \(\mathbf{x}\) of \(\mathbf{B}_m\) we can recover the whole space via the transformation \(\tildes{\mathbf{u}} = \mathbf{Q}_m \mathbf{x}\). The approach is summarized below:

Algorithm 8. Rayliegh-Ritz procedure

  • Choose a space \(H\) and a basis \((\mathbf{h}_1, \mathbf{h}_2, \dots, \mathbf{h}_m)\) represented by \(\mathbf{H}\).

  • Orthonormalize the basis

  • Calculate the Rayleigh matrix \(\mathbf{B}_m := \mathbf{Q}_m^\star \mathbf{A}\mathbf{Q}_m\)

  • Calculate of the eigen-decomposition of \(\mathbf{B}_m : (\tilde{\lambda}, \mathbf{x})\)

  • Calculate the Ritz elements: \((\tilde{\lambda}, \tilde{\mathbf{u}} = \mathbf{Q}_m \mathbf{x})\)

  • Test the residual error: \(\mathbf{r} = \mathbf{A} \tilde{\mathbf{u}} -\tilde{\lambda} \tilde{\mathbf{u}}\)

Note

  • The Ritz elements are the eigen modes of the matrix (of order \(n\)) of the Galerkin approximation \(\mathbf{C}_m = \mathbf{P}-m \mathbf{A}\mathbf{P}_m\).

  • The computational complexity is at worst \(O(m^2(4n+ m))\) (much less with a good space, e.g., Lanczos and IRAM), and is not commensurate with that of a good QR \((O(n^2))\) or that of an iteration of the Rayleigh quotient encapsulated in a heuristic process such as the one used in CALC_MODES with OPTION = [‘PROCHE’, ‘AJUSTE’, ‘SEPARE’] (\(\gg pn^3\)). But these methods are more complementary than competitive.

To estimate the error due to eigen-elements of the matrix \(\mathbf{B}_m\), we have the following result.

Theorem 1

Let \((\lambda, \mathbf{u})\) be an eigen element of a diagonalizable matrix \(\mathbf{A}\) and let \(\mathbf{B}_m\) be the associated Rayleigh matrix. Then the Rayleigh matrix has an eigenvalue \(\tilde{\lambda}\) satisfying

\[|\lambda- \tilde{\lambda}|\le\beta\frac{\|(\mathbf{I} - \mathbf{P}_m )\mathbf{u}\|_2}{ \| \mathbf{P}_m \mathbf{u} \|_2^2}\]

where \(\beta = \beta(\lambda, \mathbf{A}, \mathbf{P}_m)\).

In the Hermitian case, the numerator is also squared.

The associated eigenvector, \(\tilde{\mathbf{u}}\) satisfies

\[\sin(\mathbf{u}, \tilde{\mathbf{u}}) \le \beta \frac{\|(\mathbf{I} - \mathbf{P}_m )- \mathbf{u} \|_2}{\|\mathbf{P}_m \mathbf{u}\|_2}\]
../../_images/fig13_r5.01.01.png

Proof See 57 pp531-537.

We can show that for \(m\) large enough, \(\beta\) can be increased by a number which depends only on \(\mathbf{A}\) and \(\text{min}_{\lambda_i \ne \lambda_j} |\lambda_i -\lambda_j |\). The estimate of the second member \(\|(\mathbf{I} - \mathbf{P}_m ) \mathbf{u} \|_2\) depends on the choice of the subspace.

Note

  • The quantity \(\text{min}_{\lambda_i \ne \lambda_j} |\lambda_i -\lambda_j |\) (and its variants) is omnipresent in error and convergence analysis, or spectral conditioning. This expression points to the importance of spectrum separation on good numerical performance.

Of course, to take into account known spectral information, or even to improve or filter it, we iterate the Rayleigh-Ritz procedure by modifying the projection subspace. We can then chain the construction of this new space (recurring from the previous one) and the process of projection. Two types of spaces (each associated with different methods) are most often used and discussed next.

13.1.5.3. Choice of projection space

Two choices of projection space \(H\) are most often used:

  • The method of Bathe and Wilson (Section 13.1.8 METHODE = ‘JACOBI’), uses a sub-space \(H\) of dimension \(m\) then to recursively builds:

    \[\begin{split}\mathbf{H}_1 &= \mathbf{A}\mathbf{H} \\ \mathbf{H}_2 &= \mathbf{A}\mathbf{H}_1 = \mathbf{A}^2 \mathbf{H} \\ & \dots \\ \mathbf{H}_i &= \mathbf{A}\mathbf{H}_{i-1} = \mathbf{A}^i \mathbf{H}\end{split}\]

    This method, which generalizes the method of inverse powers to block form and also truncates the QR algorithm, leads to a reduced problem: \(\text{dim}\mathbf{H}_i \le \text{dim}\mathbf{H}_{i-1}\). In order to avoid finding the same dominant eigenmode repeatedly, it is necessary to add a reorthogonalisation step in the process.

  • The Lanczos method (Section 13.1.6 METHODE = ‘TRI_DIAG’) and the IRAM method (Section 13.1.7 METHOD = ‘SORENSEN’) starts with an initial vector \(\mathbf{y}\), to build the sequence of Krylov spaces \((H_i), i = 1, m\) where \(H_i\) is the space generated by the family of vectors \((\mathbf{y}, \mathbf{A}\mathbf{y}, \dots, \mathbf{A}^{i-1}\mathbf{y})\). The latter is called the Krylov space of \(\mathbf{A}\) and of order \(i\), initiated by \(\mathbf{y}\):

    \[H_i = K_i(\mathbf{A}, \mathbf{y}) \equiv \text{Vect}(\mathbf{y}, \mathbf{A}\mathbf{y}, \dots, \mathbf{A}^{i-1}\mathbf{y})\]

    These spaces satisfy \(\text{dim} H_i \le \text{dim} H_{i+1}\). Unlike the previous choice, we therefore have an enrichment of the problem space. On the other hand, they make it possible to project optimally according to the criteria defined above.

Historically, the first type of space has been widely used in structural mechanics. But to decrease the computational complexity linked to the size of the subspaces and to the orthogonalizations, we now prefer the Lanczos / Arnoldi type methods.

Note

  • There are an impressive variety of algorithms that use a space of the first type. These are called “subspace iterations”, “orthogonal iterations” or even “simultaneous iterations”. Note that all these adaptations use the strategies of restart 21, acceleration 22 and factorization techniques implemented to enrich the space. There are even versions based on inverse powers that allow the calculation of the \(p\) smallest modes (cf. 57 pp538-45, 85 pp49-54).

For these projection methods, assuming that the initial space \(H\) is not too poor relative to the \(p\) dominant directions, the convergence factor of the i-th mode after k iterations is written (when they are ordered in decreasing order of modulus):

\[O\left(\left|\frac{\lambda_{m+1}}{\lambda_i}\right|^k\right)\]

This exhibits two features:

  • Prioritized convergence of dominant modes

  • Improvement of these convergence rates (and therefore of their error norms for a fixed number of iterations) as the size of the subspace increases.

To transform the generalized modal problem into a standard problem one has recourse to a spectral transformation. This also makes it possible to orient the search for the spectrum and to improve convergence (Section 13.1.3.6) .

21

Technique consisting of restarting a modal algorithm with an initial vector already containing spectral information. Typically, one chooses a linear combination of the eigenvectors or Schur vectors that have already been identified.

22

Technique consisting of replacing the working matrix \(\mathbf{A}\) by a matrix polynomial \(P(\mathbf{A})\) which has larger amplitude in the spectral regions of interest (cf. Section 13.1.7.4).

13.1.5.4. Choice of spectral shift

To calculate the smallest eigenvalues ​​close to a given frequency or all the eigenvalues in a given frequency interval, we perform out a spectral shift of the generalized problem.

Let \(\sigma\) be the value of the shift that transforms the initial problem \(\mathbf{A} \mathbf{u} = \lambda \mathbf{B} \mathbf{u}\) into a standard shifted problem, \(\mathbf{A}^\sigma \mathbf{B}^{-1} \mathbf{u} = \mu \mathbf{u}\) with \(\mathbf{A}^\sigma = \mathbf{A} -\sigma \mathbf{B}\) and \(\mu = \lambda - \sigma\).

This spectral transformation, called “simple shift”, is used in the method of Bathe and Wilson. As the process detects the smallest eigenvalues ​​in \(\mu\), one gradually captures those which are the closest to \(\sigma\) (in modulus if the shift is real and partly real and imaginary if it is complex).

When we perform the inverse transformation (“shift and invert”), we have \(\mathbf{A}_\sigma \mathbf{u} = \mu \mathbf{u}\) with \(\mathbf{A}_\sigma = ( \mathbf{A} -\sigma \mathbf{B} )^{-1} \mathbf{B}\) and \(\mu =\frac{1}{\lambda - \sigma}\)

This is the problem treated with Lanczos and IRAM; which makes it possible to calculate the largest eigenvalues \(\mu\) (thus the eigenvalues ​​closest to shift \(\sigma\) in modulus). In the case of the IRAM method for matrices with complex coefficients, the shift is necessarily complex.

In the command CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’], there are four ways of choosing this shift (for buckling, the critical load levels are analogous):

  • \(\sigma = 0\), one seeks the smallest eigenvalues ​​of the starting problem. This corresponds to OPTION = ‘PLUS_PETITE’.

  • \(\sigma = \sigma_0\) with \(\sigma_0 = (2\pi f_0 )^2\), one seeks the frequencies close to the frequency FREQ = \(f_0\) (OPTION = ‘CENTRE’).

  • \(\sigma = \tfrac{1}{2}(\sigma_0 + \sigma_1)\) with \(\sigma_0 = (2\pi f_0)^2\) and \(\sigma_1 = (2\pi f_1 )^2\), we are looking for all frequencies included in the interval \([ f_0, f_1 ]\) (OPTION = ‘BANDE’ with FREQ = \(\{ f_0 , f_1 \}\) ).

  • \(\sigma = \sigma_0 + j \sigma_1\) with \(\sigma_0 = (2\pi f_0 )^2\) and \(\sigma_1= \frac{\eta}{2} (2\pi f 0 )^2\) we are looking for frequencies close to the frequency FREQ = \(f_0\) and reduced damping AMOR_REDUIT = \(\eta\) (OPTION = ‘CENTRE’). This approach is only available in dynamics with a GEP for complex modes: matrices \(\mathbf{A}\) and / or \(\mathbf{B}\) symmetric or real/ complex unsymmetric (cf. table 1).

The number of frequencies to be calculated is given in general by the user using NMAX_FREQ under the keyword factor CALC_FREQ. But for the ‘BANDE’ option, it is determined automatically using the corollary 3bis of [R5.01.04].

Note

  • The factorization of the initial matrix, just like the Sturm tests of the option BANDE, may be numerically unstable when the shift is close to an eigenvalue. For options that can reduce the problem see [R5.01.04].

13.1.6. Lanczos method (SOLVEUR_MODAL =_F(METHOD = ‘TRI_DIAG’))

13.1.6.1. Introduction

This algorithm, developed by Lanczos 55 in 1950, was hardly used until the mid-1970s. At first glance, it appears simple and effective. Nevertheless it can exhibit numerical instability which can cause ghost modes to be extracted! These are mainly related to the problem of maintaining orthogonality between the vectors generating the Krylov space.

Understanding this type of behavior in the face of finite precision of computers took a long time to unearth. Since then, many solutions have emerged and an abundant literature covers the subject, e.g., JK Cullum 43 which is entirely devoted to it, BNParlett 64 and the exhaustive synthesis by JL Vaudescal 85 (pp55-78).

In the rest of this section, we will first focus on the theoretical framework. Then, before detailing the Code_Aster implementation, we will connect the theory to a realistic finite precision framework.

13.1.6.2. Theoretical Lanczos algorithm

13.1.6.2.1. Principle

The Lanczos algorithm can be applied to the pair, [work operator, (pseudo) scalar product], such that the Hermiticity of \(\mathbf{A}_\sigma\) is ensured. A Rayleigh matrix \(\mathbf{B}_m\) of modular size, tridiagonal and Hermitian (with a real scalar product, otherwise it loses the latter property) is built iteratively. This particular form greatly facilitates the computation of the Ritz modes: for QR from \(O(pm^2)\) to \(O(20pm)\), when one searches for \(p\) eigen modes by projecting on a subspace of size \(m\). The algorithm consists of gradually building a family of Lanczos vectors \(\mathbf{q}_1, \mathbf{q}_2, \dots \mathbf{q}_m\) by projecting orthogonally, at iteration \(k\), the vector \(\mathbf{A}_\sigma \mathbf{q}_k\) on the two previous vectors \(\mathbf{q}_k\) and \(\mathbf{q}_{k-1}\). The new vector becomes \(\mathbf{q}_{k+1}\) and so, step by step, the orthogonality of this family of vectors is structurally ensured.

The iterative process is summarized below.

Algorithm 9. Theoretical Lanczos.

  • \(\mathbf{q}_0 = 0, \beta_0= 0\)

  • Calculate \(\mathbf{q}_1\) with \(\| \mathbf{q}_1\| = 1\)

  • For \(k = 1, m\) do

    • \(\mathbf{z} = \mathbf{A}_\sigma\mathbf{q}_k - \beta_{k-1}\mathbf{q}_{k-1}\)

    • \(\alpha_k= ( \mathbf{z} , \mathbf{q}_k )\)

    • \(\mathbf{v} = \mathbf{z} - \alpha_k\mathbf{q}_k\)

    • \(\beta_k = \| \mathbf{v} \|\)

    • If \(\beta_k \ne 0\) then

      • \(\mathbf{q}_{k+1} = \frac{\mathbf{v}}{\beta_k}\)

    • Else

      • Deflation

    • End If

  • End For

If \(\mathbf{e}_m\) is the m-th vector of the canonical basis, the residual vector can be written \(\mathbf{R}_m = \beta_m \mathbf{q}_{m+1} \mathbf{e}_m^T\).

../../_images/fig14_r5.01.01.png

Fig. 13.1.13 Lanczos factorization.

The Rayleigh matrix then takes the form

\[\begin{split}\mathbf{B}_m = \begin{bmatrix} \alpha_1 & \bar{\beta_1} & 0 &0 \\ \beta_1 & \alpha_2 & \dots & 0 \\ 0 & \dots & \dots & \bar{\beta}_{m-1} \\ 0 & 0 & \beta_{m-1} & \alpha_m \end{bmatrix}\end{split}\]

By construction, this algorithm assures us of the following results (in exact arithmetic):

  • The vectors \((\mathbf{q}_k)_{k = 1, m}\) constitute an orthonormal family

  • They generate a Krylov space of order \(m\) beginning with \(\mathbf{q}_1\),

    \[\mathbf{K}_m(\mathbf{A}_\sigma, \mathbf{q}_1) = \text{Vect}(\mathbf{q}_1, \mathbf{A}_\sigma\mathbf{q}_1, \dots , \mathbf{A}_\sigma^{m-1} \mathbf{q}_1 )\]

  • They allow us to build a matrix tridiagonal and Hermitian matrix \(\mathbf{B}_m\)

  • The spectrum of \(\mathbf{B}_m\) admits only the \(m\) dominant simple modes of which \(\mathbf{q}_1\) is a component.

Note

  • The algorithm therefore does not, in theory, allow the capture of multiple modes. This can be explained by noting that any irreducible Hessenberg matrix 23 (therefore a tridiagonal matrix) only admits simple modes.

  • The cost of an iteration is low, of the order of that of a matrix-vector product \(\mathbf{A}_\sigma \mathbf{q}_k\), i.e. for the first m iterations a computational complexity of \(O( nm ( c +5))\) (with \(c\) the mean number of non-zero terms on the rows of the work matrix). On the other hand, the memory complexity required is weak because we don’t need to build \(\mathbf{A}_\sigma\) in full, we just need to know its action on a vector. This feature is very useful for solving problems of large sizes 85.

  • In general the field of application determines the choice of an initial Lanczos vector. This vector must belong to a space of admissible solutions (satisfying constraints, boundary conditions, etc.) as well as to the image set of the operator. This last point is important because it allows us to enrich modal search more quickly by not being limited to the kernel.

  • The Lanczos algorithm is only one way to approach Krylov subspaces. Its generalization to unsymmetrical configurations is called the Arnoldi algorithm (Section 13.1.7.2).

23

A matrix \(\mathbf{A}\) is said to be upper (or lower) Hessenberg if \(\mathbf{A}_{i,j} = 0\) for \(i > j +1\) (or \(j > i +1\)). It is also called irreducible if \(\mathbf{A}_{i+1, i} \ne 0 \forall i\).

13.1.6.2.2. Error estimates and convergence

Due to the shape of \(\mathbf{B}_m\), the extraction of the Ritz elements is simplified and the following result allows us to quickly ( via a product of two reals) estimate their qualities.

Property 2.

The Euclidean norm of the residual of the Ritz element \((\tilde{\lambda}, \tilde{\mathbf{u}} = \mathbf{Q} m \mathbf{x} )\) is equal to

\[\|\mathbf{r}\|_2 = \|(\sigma - \tilde{\lambda} \mathbf{I}) \tilde{\mathbf{u}}\|_2 = |\beta_m ||\mathbf{e}_m^T \mathbf{x}|\]

Proof: Take the Euclidean norm of the Lanczos factorization

\[\mathbf{A}_\sigma \mathbf{B}_m \mathbf{x} = \mathbf{Q}_m \mathbf{B}_m \mathbf{x}+ \beta_m \mathbf{q}_{m+1} \mathbf{e}_m^T \mathbf{x}\]

where \(\mathbf{q}_{m+1}\) is a normalized.

Note

  • A more general result, independent of any solution method, assures us that for each eigenmode \((\tilde{\lambda}, \mathbf{x})\) of \(\mathbf{B}_m\), there exists an eigenmode \((\tilde{\lambda}, \mathbf{u})\) of \(\mathbf{A}_\sigma\) such that:

    • \(|\lambda- \tilde{\lambda}|\le \frac{\| \mathbf{r} \|_2^2}{\gamma}\)

      where \(\gamma = \text{min}_{\lambda_i \ne \tilde{\lambda}}|\lambda_i - \mathbf{u}^T \mathbf{A}_\sigma \tilde{\mathbf{u}} |\)

    • \(|\sin(\mathbf{u}, \tilde{\mathbf{u}})| \le \frac{\| \mathbf{r} \|_2}{\gamma}\)

    • Therefore, more than the minimization of the residual, the stopping criterion of the Lanczos/Arnoldi methods, \(|\beta_m|<\varepsilon\), lead to the best approximation of the original spectrum. These improvements are only accessiblein the Hermitian case.

    • In the general case (and in particular non-normal matrices) they are more difficult or even impossible to construct without additional information (level of non-normality, of “defectivity”, etc.). An estimate of the residue is then no longer the good criterion to estimate the quality of the approximated modes.

Let us now look at the quality of convergence of the algorithm. Theorem 2 can be specialized for that purpose by limiting the error norm of \(\|(\mathbf{I} - \mathbf{P}_m) \mathbf{u}\|_2\) the following bounds are obtained:

Theorem 3

Let \((\lambda_i , \mathbf{u}_i)\) be the i-th dominant eigenmode of the diagonalizable matrix \(\mathbf{A}_\sigma\). Then there is a Ritz mode \((\tilde{\lambda}_i , \tilde{\mathbf{u}}_i )\) such that

\[\begin{split}|\lambda_i-\tilde{\lambda}_i| & \le \frac{\alpha^2}{\beta} \left(\prod_{j < i} \frac{\lambda_j - \lambda_n}{\lambda_j - \lambda_i}\right)^2 T_{m-i}^{-2}(\gamma_i) \\ |\sin(\mathbf{u}_i, \tilde{\mathbf{u}}_i)| & \le \alpha \left(\prod_{j < i} \frac{\lambda_j - \lambda_n}{\lambda_j - \lambda_i}\right) T_{m-i}^{-1}(\gamma_i)\end{split}\]

with

\[\gamma_i= 1 + 2 \frac{\lambda_i - \lambda_{i+1}}{\lambda_{i+1} - \lambda_n}\]

and \(T_{m-i}(\mathbf{x})\) is the Tchebycheff polynomial of degree \(m-i\), \(\beta\) is the constant of Theorem 1 and

\[\alpha = \frac{\beta}{\| \mathbf{P}_m\mathbf{u}_i\|_2} \tan(\mathbf{q}_1, \mathbf{u}_i)\]

Proof: See 57 pp554-557.

These observations indicate that:

  • If the initial vector has no contribution along the eigenvectors, we cannot capture the related modes \((\alpha\rightarrow +\infty)\).

  • We primarily have convergence of the extreme modes which is better because the spectrum is separated (these methods are less sensitive to clusters than the iterated methods), without settling of eigenvalues into “clusters”.

  • The error decreases when \(m\) increases (as the inverse of the polynomial \(\mathbf{T}_{m-i}(\mathbf{x})\), therefore as the inverse of an exponential).

  • The estimate on the eigenvalues ​​is better than that of their associated eigenvectors.

Note

The convergence of the method was studied by P. Kaniel then improved (and corrected) by SCPaige; A synthesis of these studies can be found in the paper by Y. Saad 67.

13.1.6.3. Practical Lanczos algorithm

13.1.6.3.1. Problem of orthogonality

During the implementation of this algorithm, one realizes that the Lanczos vectors lose orthogonality. The Rayleigh matrix is ​​then no longer the projected matrix of the initial operator and the computed spectrum is increasingly erroneous.

For a long time this phenomenon was attributed to rounding effects. In 1980, CCPaige showed that the loss of orthogonality could be traced to the convergence of an eigen mode. Thus, as soon as we capture a dominant mode, we perturb the preceding Lanczos vectors. The following result expresses this paradox.

Property 4 (Paige Analysis)

We will use the notation used in (Algorithm 9, Theoretical Laczos). If \(\varepsilon\) is the machine precision, at the iteration \(k+1\), the orthogonality of the new Lanczos vector is expressed in the form

\[|\mathbf{q}_{k+1}^T \mathbf{q}_i| = \frac{|\mathbf{v}^T \mathbf{q}_i| + \varepsilon\|\mathbf{A}_\sigma\|_2}{|\beta_k|} \quad (i \le k)\]

Proof:See 63.

The problem arises during the normalization of the new Lanczos vector \(\mathbf{q}_{k+1}\). When this mode converges according to Property 2, convergence is judged from the smallness of the coefficient \(\beta_k\). Therefore, despite the theoretical orthogonality \((|\mathbf{v}^T \mathbf{q}_i| = 0 ( i \le k ) )\), the effective orthogonality is far from being attained \((|\mathbf{q}_{k+1}^T\mathbf{q}_i| \gg \varepsilon ( i \le k ) )\).

The numerical treatment of these problems has been the subject of much research for thirty years and of many palliative strategies have been developed. The choice of method depends on the required spectral information and available computer resources, the synthesis by JL.Vaudescal gives a good overview of these variants:

  • Lanczos algorithm without reorthogonalisation, developed by JKCullum and RAWilloughby 43 which consists in deleting phantom modes from the calculated spectrum by studying the interlacing of the eigenvalues ​​of the work matrix and one of its sub-matrices. This variant requires a lost-cost additional calculation and memory to determine the eigenvalues, but it complicates (sometimes greatly) the search for eigenvectors.

  • Lanczos algorithms with total (BNParlett) or selective (J. Scott) reorthogonalization which involve reorthogonalizing the last Lanczos vector obtained relative to all the vectors already calculated or simply compared to the converged Ritz vectors (this variant makes it possible to dynamically control the loss of orthogonality).

These variants are much more expensive but they turn out to be more efficient and more robust. In the variant of Newmann-Pipano 61 (METHODE = ‘TRI_DIAG’ of Code_Aster), it is the strategy of complete reorthogonalisation 24 which is used: the vector \(\mathbf{q}_{k+1}\) is therefore not calculated as in the theoretical algorithm.

Note

  • These problems of loss of orthogonality arise all the more acutely as the size of the sub-space increases. This is an additional argument for setting up an iterative process limiting this size.

  • These strategies are available for vectors or for blocks, they are implemented in simple or iterative algorithms. The latter can benefit from explicit or implicit restarts, managed or automatic. The IRAM variant thus benefits from an implicit restart calculated automatically.

  • The central feature of these algorithms is the orthogonalization method used. The whole process is dependent on the success and robustness of the process of computing the orthogonal vectors. Since orthogonal transformations of the Housholder or Givens type are very expensive (though robust), we now prefer iterative Gram-Schmidt algorithms (IGSM) which are a better compromise between efficiency and computational complexity (see Section 13.1.14) This choice was made for the variants of Lanczos / Arnoldi implemented in the code.

  • Selective reorthogonalisation amounts to carrying out an implicit deflation by restricting the work operator so as not to have to recompute the converged modes (Section 13.1.4.1).

24

In GEP as in QEP, one uses this strategy of complete reorthogonalisation (for Lanczos and IRAM).

13.1.6.3.2. Capturing multiplicities

The Lanczos algorithm (in theory) does not allow, whatever its strategy of reorthogonalisation, the capture of multiple modes. In practice, rounding effects come to the rescue and “sprinkle” small components along almost all eigenvectors. We can therefore capture multiplicities, however they can be erroneous and require additional verification.

Note

  • In fact, only a block version (G. Golub & R. Underwood 1979) can allow us to detect correct multiplicities, at least as long as the size of the blocks is sufficient. This version was extended (M.Sadkane 69 1993) to the Arnoldi algorithm but the Sorensen variant (IRAM) is more efficient and robust.

13.1.6.3.3. Lanczos phenomenon

The success of this algorithm was initially based on what is called the “Lanczos phenomenon”. This conjecture predicts that for a sufficiently large size of subspace (\(m \gg n\)) we are able to detect the entire spectrum of the operator.

Given the low memory prerequisites of a “basic” Lanczos (roughly a tridiagonal matrix and some vectors), this is particularly interesting for treating sparse systems of very large size (the iterated power or QR type algorithms require knowledge of the entire matrix).

13.1.6.4. Complementary treatments

13.1.6.4.1. Detection of invariant spaces

In (Algorithm 8), if \(\beta_{l-1} = 0\), the new vector cannot be normalized and requires an ancillary treatment called deflation . The computed Krylov subspace is then an invariant subspace and its spectrum is exact. In other words, the first \(l-1\) eigenvalues extracted by the algorithm (QR or QL) are those of \(\mathbf{A}_\sigma\).

We therefore have to choose a new random vector satisfying the boundary conditions, orthogonalize it along with all the Lanczos vectors already calculated and build a second family of Lanczos vectors:

\[\mathbf{H}_l = \mathbf{K}_l(\mathbf{A}_\sigma, \mathbf{q}_1) = = \text{Vect}(\mathbf{q}_1, \mathbf{A}_\sigma\mathbf{q}_1, \dots, \mathbf{A}_\sigma^{l-1}\mathbf{q}_1)\]

The resulting tridiagonal matrix obtained has the form:

\[\begin{split}\mathbf{B}_m = \begin{bmatrix} \alpha_1 & \bar{\beta}_1 & 0 & \dots & 0 & 0\\ \beta_1 & \alpha2 & \dots & 0 & 0 & 0\\ \dots & \dots & \dots & \approx 0 & \dots & \dots \\ 0 & 0 & \approx 0 & \alpha_l & \bar{\beta}_l & \dots \\ 0 & 0 & \dots & \beta_l & \alpha_{l+1} & \dots \\ 0 & 0 & \dots & \dots & \dots & \bar{\beta}_{m-1} \\ 0 & 0 & \dots & \dots & \beta_{m-1} & \alpha_m \end{bmatrix}\end{split}\]

The intermediate Lanczos factorization (to order \(l\)) is written

\[\mathbf{A}_\sigma \mathbf{Q}_l = \mathbf{Q}_l \mathbf{B}_l\]

The values \((\alpha_k , \beta_k )_{ k = 1, l}\) are obtained from the first random vector and the following values \((\alpha_k , \beta_k )_{k = l+1, m}\) via a second random vector. To detect the nullity of the extra-diagonal term we uses the numerical criterion

\[|\beta_{l-1}| \le \text{PREC_LANZOS} |\alpha_l| \quad \text{then} \quad \beta_{l-1} = 0\]

where PREC_LANCZOS is initialized under the keyword CALC_FREQ.

Note

  • A more robust criterion used by the mathematical libraries EISPACK, LAPACK and ARPACK uses the previous diagonal term, a tunable parameter \(c\) and machine precision \(\varepsilon\):

    \[|\beta_{l-1}|< c \varepsilon ( | \alpha_{l-1}|+| \alpha_l| )\]

  • It is this criterion which was retained for IRAM (but the user cannot modify the parameter \(c\), it is set to a standard value by the code), various messages also warn the user of the detection of an invariant space.

  • Obtaining such a space is quite nice since it assures us of the 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 be solved. Many techniques seek to artificially reproduce this phenomenon (Section 13.1.4.1).

13.1.6.4.2. Restart strategies

To limit endemic orthogonalization problems, limit the memory prerequisites and take into account the spectral information already obtained, restart strategies have been coupled with the Lanczos algorithm. It loses its “simple” character and becomes “iterative”. We iterate the “restarts” until the desired convergence criteria are met. The algorithm does not undergo the convergence of the modes imposed by Theorem 8 and it can interactively favor that of certain modes.

In theory, however, the larger the size of the subspace, the better the convergence. A compromise is therefore to be found between the size of the subspace and the frequency of restarts. Different restart vectors can be used, usually written as a weighted sum of the \(p\) eigenmodes (Y. Saad 1980)

\[\mathbf{q}_1= \sum_i^p \chi_i \mathbf{Q}_m \mathbf{x}_i\]

Indeed, if one restarts with an initial vector belonging to the invariant subspace generated by the modes sought, one is then sure to obtain (with a good algorithm of orthogonalization) a residual norm of the order of machine precision and almost exact eigenmodes. In addition to the decision to trigger the restart, the whole problem lies in finding these weights \(\chi\). We will see that with IRAM, the restarts set up via implicit polynomial filters solve this problem elegantly and automatically.

Note

  • The restart philosophy was initiated by W. Karush (1951).

  • Restarts based on Schur vectors associated with the desired spectral decomposition appear to be more stable (J. Scott 1993).

  • Efficiency criteria have been sought to decide on the advisability of a restart (B. Vital 1990).

  • Instead of a simple linear combination of eigenvectors, we can determine a restart vector of via polynomials (Chebyshev in real and Faber in complex) allowing to center the modal search in a given area. We then speak of explicit polynomial acceleration [Sad93].

13.1.6.5. Implementation in Code_Aster

13.1.6.5.1. Variant of Newmann & Pipano

This variant was developed by M.Newmann & A.Pipano 43 61 in 1977. It is a simple Lanczos method, in real arithmetic, with total reorthogonalisation (via a GSM). It uses a classical “shift and invert” with the pseudo-scalar product introduced by the shifted matrix \((\mathbf{A} - \sigma \mathbf{B})\)

\[\underbrace{( \mathbf{A} - \sigma \mathbf{B} )^{-1}}_{\mathbf{A}_\sigma} \mathbf{B} \mathbf{u} = \frac{1}{\underbrace{\mu - \sigma}_\lambda} \mathbf{u}\]

with:

\[\begin{split}( \mathbf{x} , \mathbf{y} ) & = \mathbf{y}^T ( \mathbf{A} -\sigma \mathbf{B} ) \mathbf{x} \\ \delta & = \text{sign} ( \mathbf{x} , \mathbf{x} ) \\ \| \mathbf{x} \| &= \delta | ( \mathbf{x} , \mathbf{x} ) |^{1/2}\end{split}\]

This choice makes the pair “operator-scalar product” symmetric and it is adapted to structural mechanics matrices. We can thus deal with problems:

  • Of free structure and fluid-structure (\(\mathbf{A}\) can be singular).

  • Buckling (\(\mathbf{B}\) is not defined).

The pseudo-scalar product introduced by the shifted matrix is ​regular and it makes it possible to search for orthogonal Lanczos vectors (\(\mathbf{A} - \sigma \mathbf{B}\)). Total reorthogonalisation is carried out only if necessary and this criterion can be tuned. The price to pay for this gain in robustness is the possible loss of symmetry of the Rayleigh matrix:

\[\begin{split}\mathbf{B}_m = \begin{bmatrix} \alpha_1 & \gamma_1 & 0 & \dots & 0 & 0\\ \beta_1 & \alpha2 & \dots & \dots & 0 & 0\\ 0 & 0 & \dots & \dots & \dots & \gamma_{m-1} \\ 0 & 0 & \dots & \dots & \beta_{m-1} & \alpha_m \end{bmatrix}\end{split}\]

where \(|\gamma_i= \delta_i\beta_i\) and \(\delta_i= \text{sign}(\mathbf{q}_{i+1}, \mathbf{q}_{i+1})\).

After being balanced and shaped into upper Hessenberg form, \(\mathbf{B}_m\) is diagonalized by an implicit QL method if it remains symmetrical despite everything or by a QR method otherwise (Section 13.1.13.1)

The difference in cost between these solvers remains negligible compared to the costs of reorthogonalization. The modified Lanczos scheme then becomes:

Algorithm 10. Lanczos method. Newman-Pipano variant.

  • Random initialization of \(\mathbf{q}_{\text{init}}\)

  • \(\tilde{\mathbf{q}}_1 = (\mathbf{A} - \sigma\mathbf{B} )^{-1} (\mathbf{q}_{\text{init}}(i) \cdot \mathbf{u}_{\text{lagr}}(i))_i\)

  • \(\hat{\mathbf{q}}_1 = \mathbf{A}_\sigma (\tilde{\mathbf{q}}_1(i) \cdot \mathbf{u}_{\text{block}}(i))_i\)

  • \(\delta_1 = \text{sign}(\hat{\mathbf{q}}_1^T (\mathbf{A} - \sigma \mathbf{B}) \hat{\mathbf{q}}_1)\)

  • \(\mathbf{q}_1 = \delta_1 \hat{\mathbf{q}}_1 |\hat{\mathbf{q}}_1^T (\mathbf{A} - \sigma \mathbf{B}) \hat{\mathbf{q}}_1|^{-1/2}\)

  • \(\mathbf{q}_0 = 0, \beta_0= 0, d_0= 0\)

  • For \(k = 1, m\) do

    • \(\mathbf{z} = \mathbf{A}_\sigma \mathbf{q}_k - \delta_{k-1} \beta_{k-1} \mathbf{q}_{k-1}\)

    • \(\alpha_k = \mathbf{q}_k^T (\mathbf{A} - \sigma \mathbf{B}) \mathbf{z}\)

    • \(\mathbf{v} = \mathbf{z} -\delta_k\alpha_k\mathbf{q}_k\)

    • Reothogonalization of \(\mathbf{v}\) with respect to \((\mathbf{q}_i)_{i = 1,k}\) (IGSM )

    • \(\beta_k= \mathbf{v}^T (\mathbf{A} -\sigma \mathbf{B}) \mathbf{v}\)

    • \(\delta_k= \text{sign} (\mathbf{v}^T (\mathbf{A} -\sigma \mathbf{B}) \mathbf{v})\)

    • If \(\beta_k \ne 0\) then

      • \(\mathbf{q}_{k+1} = \cfrac{\delta_k\mathbf{v}}{|\beta_k|^{1/2}}\)

    • Else

      • Deflation

    • End If

  • End For.

Note

  • The \(\mathbf{B}\)-orthogonality, and even more so Euclidean orthogonality, is not the result of special configurations. This does not modify the orthogonality properties of the eigen modes (see TODO Proposition 2).

  • There exists a whole zoo of “operator-scalar product” pairs, this one used in this method being only one possibility among others. Thus, the classical 86 libraries treat buckling problems in “buckling mode” via the same “shift and invert” and the pseudo-scalar products introduced by \(\mathbf{A}\). But due to the almost systematic introduction of Lagrange multipliers, this matrix becomes undefined or even singular, which greatly disturbs the process. The same causes produce the same effects when, for a dynamic calculation, we use the \(\mathbf{B}\)-scalar product.

The reorthogonalisation is carried out with a “homebrew” variant of the Iterative Gram-Schmidt Method (IGSM) described below.

Algorithm 11. Procedure for reorthogonalisation of TRI_DIAG.

  • For \(i = 1, k\)

    • \(a = \mathbf{q}_{k+1}^T (\mathbf{A} - \sigma \mathbf{B}) \mathbf{q}_i\)

    • If \(\left|\left<\mathbf{q}_{k+1}^T (\mathbf{A} - \sigma \mathbf{B}) \mathbf{q}_i\right>\right| \ge \text{PREC_ORTHO}\) then

      • For \(j = 1 , \text{NMAX_ITER_ORTHO}\)

        • \(\mathbf{x} = \mathbf{q}_{k+1} - (\mathbf{q}_{k+1}^T (\mathbf{A} - \sigma \mathbf{B})\mathbf{q}_i) \mathbf{q}_i\)

        • \(b = \mathbf{x}^T(\mathbf{A} - \sigma \mathbf{B})\mathbf{q}_i\)

        • If \((|b| \le \text{PREC_ORTHO})\) then

          • \(\mathbf{q}_{k+1} = \mathbf{x}\)

          • \(i +1 \implies i\)

          • Exit

        • Else

          • If \(b \le a\) then

            • \(\mathbf{q}_{k+1} = \mathbf{x}\)

            • \(a = b\)

          • Else

            • Failure, emission of an alarm message

            • \(i +1 \implies i\)

            • Exit

          • End If

        • End If

      • End For in \(j\)

    • End If

  • End For in \(i\)

Note

  • After some testing, it seems that this code specific variant is less efficient than the IGSM of Kahan-Parlett type 64 chosen for IRAM (see Section 13.1.14).

13.1.6.5.2. Parameter setting

To be able to activate this method, it is necessary to use METHODE = ‘TRI_DIAG’. On the other hand, the sizeof the projection subspace is determined either by the user or empirically from the formula

\[m = \text{min} (\text{max} (4p, p+7) , n_{\text{dofs - active}}\]

where

  • \(p\) is the number of eigenvalues ​​to calculate

  • \(n_{\text{dof - active}}\) is the number of active degrees of freedom of the system (Section 13.1.3.2)

The user can always impose the dimension himself by indicating it with the keyword DIM_SOUS_ESPACE of the keyword factor CALC_FREQ. Instead of a default value, we may prefer a multiplicative coefficient related to the number of frequencies contained in the study interval. This type of functionality concerns the keyword COEF_DIM_ESPACE.

The parameters for total reorthogonalisation, PREC_ORTHO and NMAX_ITER_ORTHO, the maximum number of QR iterations, NMAX_ITER_QR, and the deflation criterion, PREC_LANCZOS, are accessible under the keyword factor CALC_FREQ. When the deflation phenomenon occurs, a message specifies its rank 1.

13.1.6.5.3. Warning on the quality of the modes

Remember that this variant is not iterative. From the number of desired frequencies, we estimate the dimension of the computation subspace and one hopes that the eigen modes of the projected problem will be a good approximation of those of the generalized problem. The validity of the results is checked a posteriori (Section 13.1.3.5) but one has no active control on this precision. If they are not satisfactory, the only recourse is to carry out a new test by increasing the dimension of the subspace.

Note

The IRA method that we are going to discuss offers a flexible and dynamic control of the precision results, it is an iterative method.

13.1.6.5.4. Detection of rigid modes

An algebraic method calculating the zero eigenvalues of the generalized modal problem was introduced. Physically these correspond to the movements at zero strain energy of a free structure (cf 76 TP no 2). Apart from the computational difficulties of handling almost zero quantities, due to their multiplicities, their correct capture was so far often problematic for the Lanczos implemented in Code_Aster: phantom modes appeared corresponding to failed multiplicities! This algorithm of detection of the modes of rigid body, which one activates by initializing OPTION = ‘MODE_RIGIDE’ (default value ‘SANS’), intervenes in preprocessing of modal computation (only with ‘TRI_DIAG’). It is based on the analysis of the stiffness matrix and breaks down into three phases:

  • Detection of null pivots of this matrix.

  • Locking of these pivots.

  • Solution of a linear system of which the associated eigenvectors are solutions.

During the Lanczos process it is therefore sufficient to orthogonalize the base vectors. However, the introduction of the IRA method reduces the interest (if not for comparison) of such an option (which is has a cost since it requires matrix inversions).

This method nevertheless remains useful as a back-up solution in the event of IRAM failure. We can also consider encapsulating it in an auxiliary operator (like INFO_MODE).

Note

The option MODE_RIGIDE can be activated only with Lanczos for GEP with real modes.

13.1.6.6. Scope of use

GEP with real symmetric matrices. The user can specify the type of calculation (dynamic or buckling) with TYPE_RESU and one of FREQ or CHAR_CRIT.

13.1.6.7. Display in message file

The example below from the list of code tests (sdll112a) summarizes all the traces from the algorithm. The number of effective QR iterations (or QL ) can only be the same for all eigen values. This is an information artefact which will be removed in future versions.

------------------------------------------------------------------------
THE NUMBER OF DOF TOTAL IS:   86
FROM LAGRANGE IS:             20
THE NUMBER OF ACTIVE DOL IS:  56
------------------------------------------------------------------------
THE OPTION CHOSEN IS: PLUS_PETITE
THE FREQUENCY OFFSET VALUE IS: 0.00000E + 00
------------------------------------------------------------------------
INFORMATION ON THE CALCULATION REQUIRED :
NUMBER OF MODES REQUESTED: 10
THE DIMENSION OF THE REDUCED SPACE IS: 0
IT IS LESS THAN THE NUMBER OF MODES, IT IS TAKEN EQUAL TO 40
------------------------------------------------------------------------
THE CALCULATED FREQUENCIES INF. AND SUP. ARE:
FREQ_INF: 1.54569E + 01
FREQ_SUP: 1.01614E + 02
THE FIRST HIGH FREQUENCY NOT RETAINED IS: 1.29375E + 02
------------------------------------------------------------------------
MODAL CALCULATION: SIMULTANEOUS ITERATION METHOD
LANCZOS METHOD
NUMBER FREQUENCY (HZ) ERROR NORM 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
...
10 1.01614E + 02 3.18663E-124
------------------------------------------------------------------------
POSTERIORI CHECKING OF THE MODES
IN THE INTERVAL (1.54182E + 01, 1.16326E + 02)
THERE ARE 10 FREQUENCY (S)
------------------------------------------------------------------------

13.1.6.8. Summary of settings

Let us recapitulate the parameter settings available for the operator CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’] with this algorithm.

Table 13.1.5 Summary of the parameter setting of CALC_MODES with OPTION among[‘BANDE’, ‘CENTRE’, ‘PLUS_*’] (GEP) with SOLVEUR_MODAL = _F (METHOD = ‘TRI_DIAG’).

Operand

Keyword

Default value

References

TYPE_RESU = ‘DYNAMIQUE’

‘DYNAMIQUE’

Section 13.1.3.1

‘MODE_FLAMB’

Section 13.1.3.1

OPTION = ‘PLUS_*’

‘PLUS_PETITE’

Section 13.1.5.4

‘BANDE’

Section 13.1.5.4

‘CENTRE’

Section 13.1.5.4

SOLVEUR_MODAL

METHOD = ‘TRI_DIAG’

‘SORENSEN’

Section 13.1.6.5

MODE_RIGIDE = ‘NON’

‘NON’

Section 13.1.6.5

‘OUI’

Section 13.1.6.5

PREC_ORTHO

1.E-12

Section 13.1.6.5

NMAX_ITER_ORTHO

1.E-04

Section 13.1.6.5

PREC_LANCZOS

1.E-04

Section 13.1.6.5

NMAX_ITER_QR

15

Section 13.1.6.5

DIM_SOUS_ESPACE

Calculated

Section 13.1.6.5

COEF_DIM_ESPACE

Calculated

Section 13.1.6.5

CALC_FREQ

FREQ or CHAR_CRIT (if ‘CENTRE’ )

Section 13.1.5.4

NMAX_FREQ (if ‘PLUS_*’)

10

Section 13.1.5.4

NMAX_CHAR_CRIT (if ‘PLUS_*’)

10

Section 13.1.5.4

NMAX_ITER_SHIFT

3

81 Section 13.1.3.2

PREC_SHIFT

0.05

81 Section 13.1.3.2

SEUIL_FREQ or SEUIL_CHAR_CRIT

1.E-02

Section 13.1.3.7

VERI_MODE

STOP_ERREUR = ‘OUI’

‘OUI’

Section 13.1.3.7

‘NON’

Section 13.1.3.7

PREC_SHIFT

5.E-03

Section 13.1.3.7

SEUIL

1.E-06

Section 13.1.3.7

STURM = ‘OUI’

‘OUI’

Section 13.1.3.7

‘NON’

Section 13.1.3.7

Note

  • We find all the parameters related to the Sturm test (NMAX_ITER_SHIFT, PREC_SHIFT) and with verification checks (SEUIL_FREQ /THRESHOLD_CHAR_CRIT, VERI_MODE).

  • Beginners are advised not to change most of the default parameters as these have been determined empirically. Only the ones in bold should ideally be modified. During the first passages, it is strongly recommended to modify only the main parametersnoted in bold. The others relate more to the mysteries of the algorithm and they have been initializedempirically at standard values. * In particular, to improve the quality of a mode, the only adjustable parameter is the dimension of thesubspace (DIM_SOUS_ESPACE / COEF_DIM_ESPACE

13.1.7. IRA Algorithm (SOLVEUR_MODAL = _F(METHOD = ‘SORENSEN’))

13.1.7.1. Introduction

One of the problems of the Lanczos method is the loss of orthogonality of its basis vectors. A generalization of this algorithm to the non-Hermitian case by WEArnoldi 36 in 1951 partially solves this problem. It has been brought up to date by Y.Saad 67 in 1980 and an abundant literature covers the subject. Apart from the seminal article and the papers of Y.Saad and M.Sadkane 69, we recommend the updated and exhaustive synthesis of JLVaudescal 85. This method is the basis of the IRAM algorithm known as the “Sorensen” method (IRAM for ‘Implicit Restarted Arnoldi Method’).

We will first detail its operation, its behaviors and its limitations. Subsequently, we will identify the issues to which IRAM must respond (and it does so in most cases) and we will examine its theoretical and numerical properties. We will conclude with a summary of the related Code_Aster parameters and with an example of a message file when this method is used.

13.1.7.2. Arnoldi algorithm

13.1.7.2.1. Principle

IRA algorithm is intended to be applicable to all “operator-scalar product” pairs. The first step of a good QR, except for balancing, consists in reducing the work matrix to Hessenberg form. This is not very detrimental, because we can then directly apply the QR algorithm to it. Also, this makes it possible to gain an order of magnitude during the actual solution (Section 13.1.13.1).

The algorithm is very similar to that of Lanczos, it consists in gradually building a family of Arnoldi vectors \(\mathbf{q}_1, \mathbf{q}_2, \dots , \mathbf{q}_m\) by orthogonal projection, at iteration \(k\), of the vector \(\mathbf{A}_\sigma \mathbf{q}_k\) on the \(k\) previous vectors. The new vector becomes \(\mathbf{q}_{k+1}\) and so, step by step, we assure the orthogonality of this family of vectors.

Unlike Lanczos, the orthogonality of the new vector to all the previous ones is explicit and not implicit. This is managed by the modified Gram-Schmidt algorithm (GSM) (see Section 13.1.14) which proves to be sufficiently robust in most cases.

If \(\mathbf{e}_m\) is the m-th vector of the canonical basis, the residual vector of the Arnoldi factorization is written \(\mathbf{R}_m = b_{m+1, m} \mathbf{q}_{m+1} \mathbf{e}_m^T\).

../../_images/fig15_r5.01.01.png

Fig. 13.1.14 Arnoldi factorization.

The iterative process is summarized as follows:

Algorithm 12. Theoretical Arnoldi.

  • Calculate \(\mathbf{q}_1\) with \(\| \mathbf{q}_1\| = 1\)

  • For \(k = 1, m\) do

    • \(\mathbf{z} = \mathbf{A}_\sigma\mathbf{q}_k\)

    • For \(l = 1, k\) do (GSM)

      • \(b_{lk} = (\mathbf{z}, \mathbf{q}_l)\)

      • \(\mathbf{z} = \mathbf{z} - b_{lk} \mathbf{q}_l\)

    • End For \(l\)

    • \(b_{k+1, k} = \| \mathbf{z} \|\)

    • If \(b_{k+1, k} \ne 0\) then

      • \(\mathbf{q}_{k+1} = \mathbf{v} b_{k+1, k}\)

    • Else

      • Deflation

    • End if

  • End For \(k\).

The Rayleigh matrix is then written:

\[\begin{split}\mathbf{B}_m = \begin{bmatrix} b_{11} & b_{12} & \dots & b_{1m} \\ b_{21} & b_{22} & \dots & b_{2m} \\ 0 & \dots & \dots & b_{m-1,m} \\ 0 & 0 & b_{m,m-1} & b_{mm} \end{bmatrix}\end{split}\]

Apart from the shape of this matrix and a lower sensitivity to orthogonality issues, this algorithm ensures the same theoretical and numerical properties as Lanczos. However if we let the size of the subspace grow indefinitely until convergence to the desired eigenvalues, rounding effects will nevertheless regain the upper hand and we can run into problems. Hence the need, as for Lanczos, to make this process iterative through restarts.

Note

  • This method can be seen as an implicit variant of the Lanczos algorithm with total reorthogonalisation.

  • The implicit orthogonalization of the algorithm can be carried out by more expensive, but more robust, algorithms such as QR or IGSM. This is required when the work operator lacks normality.

  • Due to the structure of \(\mathbf{B}_m\), the memory complexity is more than with Lanczos. On the other hand, the computational complexity remains of the same order of magnitude \(O(nm(c+3+m))\) (with \(c\) the mean number of non-zero terms on the rows of the work matrix).

  • To improve the first point, Y.Saad 67 showed that the upper Hessenberg structure can see its extreme sub-diagonals canceled out if we only partially perform the reorthogonalization.

  • Since vector algorithms have a natural tendency to miss multiplicities, we often prefer to use a block version (M.Sadkane 1993). But the size of these affects the quality of results, it is for this reason that we prefer the vector versions of IRAM.

  • The choice of the initial Arnoldi vector is carried out in the same way as for Lanczos.

13.1.7.2.2. Error estimates and convergence

To evaluate of the approximation quality of the eigen modes, there is a simple and efficient and efficient criterion for IRAM.

Property 5

The Euclidean norm of the residual of the Ritz element \((\tilde{\lambda}, \tilde{\mathbf{u}} = \mathbf{Q}_m\mathbf{x})\) is equal to

\[\|\mathbf{r}\|_2 = \|(\mathbf{A}_\sigma -\tilde{\lambda}\mathbf{I}) \tilde{\mathbf{u}}\|_2 = | b_{m+1, m}|| \mathbf{e}_m^T \mathbf{x}|\]

Proof: Take the Euclidean norm of the Arnoldi factorization

\[\mathbf{A}_\sigma \mathbf{Q}_m \mathbf{x} = \mathbf{Q}_m \mathbf{B}_m \mathbf{x} + b_{m+1, m} \mathbf{q}_{m+1} \mathbf{e}_m^T \mathbf{x}\]

and normalize it like \(\mathbf{q}_{m+1}\). The proof follows after some algebra.

Using the norm of the projection, \(\| (\mathbf{I} - \mathbf{P}_m ) \mathbf{u} \|_2\) we can generalize the convergence Theorem 3 to the non-Hermitian case.

Theorem 6

Let \((\lambda_1, \mathbf{u}_1)\) be the first (classical ordering, in decreasing order of magnitude) dominant eigenmode of \(\mathbf{A}_\sigma\) (diagonalizable) and let \(\mathbf{q}_1 = \sum_{k = 1}^n \alpha_k\mathbf{u}_k\) be the initial Arnoldi vector decomposed in the eigenvector basis.

Then there exists a Ritz mode \((\tilde{\lambda}_1, \tilde{\mathbf{u}}_1)\) such that:

\[| \lambda_1-\tilde{\lambda_1}|\le\frac{\alpha^2}{\beta} \xi_1\delta_1^m\]

with \(\alpha =\beta\| \mathbf{P}_m\mathbf{u}_i\|_2\), \(\beta\) being the constant of Theorem 1, and

\[\xi_1 = \sum_{k=1, k \ne 1}^n \left|\frac{\alpha_k}{\alpha_1}\right|\]

and

\[\delta_1^m = \left(\sum_{j = 2}^{m +1} \prod_{k = 2, k \ne j}^{m +1} \left| \frac{\lambda_k- \lambda_1}{\lambda_k- \lambda_j}\right| \right)^{-1}\]

This result is applied in the same way to the other modes.

Proof: By taking again the proof of Y.Saad (68 pp209-210) and the result of Theorem 1.

These bounds, very different from those obtained with Lanczos, however guide the same phenomena:

  • If the initial vector has no contribution along the desired eigenvectors, one cannot capture them \((\xi_i \rightarrow +\infty)\).

  • One has primarily convergence of the peripheral modes of the spectrum, and it is better that the spectrum is separated (property of \(\lambda_i^m\)).

  • The decrease in error is proportional to the increase in \(m\) (property of \(\lambda_i^m\)).

Note

  • When an eigenvalue is badly conditioned, \(\xi_i \rightarrow+\infty\) then it is necessary to increase \(m\) for \(\lambda_{im}\) that decreases.

  • Similar results have been seen in the case of a defective operator (see Zia 94).

13.1.7.3. Challenges

The IRA Algorithm tries to provide an elegant cure for recurring problems that are raised by other approaches:

  • Minimization of the projection space: It proposes at least \(m > p +1\) instead of \(m = 4 p\) for Lanczos.

  • Optimal management of additional orthogonalization costs establishing a compromise between the size of the subspace and frequency of restarts.

  • Transparent, dynamic and efficient management of these restarts.

  • Automatic taking into account of spectral information .

  • Determination of the memory prerequisites and the quality of the results.

One thus has more question to be asked concerning the strategy of reorthogonalisation, the frequency of restarts, their locations, the criteria for detecting any phantom modes etc. The “super-IRAM” method takes care of all these!

In short, it provides:

  • Better overall robustness.

  • The calculations complexity \(O(4nm( m - p ))\) and memory complexity \(O (2np+ p 2)\) are improved (especially compared to a simple Lanczos such as Newmann & Pipano) for fixed precision.

  • A more rigorous capture of multiplicities, clusters and rigid body modes (hence fewer parasitic modes!).

Note

  • On this last point, only a block version of IRAM or a version incorporating “purge and lock” (capture and filtering techniques, cf. D. Sorensen & RBLehoucq, 1997) can guarantee correct detection of the spectrum of a standard operator (ie not too badly conditioned).

  • It seems, that in practice in Code_Aster, the report in complexity computation between IRAM and Lanczos / Bathe & Wilson is at least order 2 in favor of the first. With reasonable problem sizes (a few tens of thousands of dof and \(p = O (100)\)) this can go up to 10 (without the encapsulation of CALC_MODES on several sub-intervals of search). In some semi-industrial cases, it made it possible to carry out a spectrum search which had failed with Jacobi and which was unaffordable with Lanczos (given the time limits).

  • A class of algorithm known as “Jacobi-Davidson” (cf. RBMorgan 1990) seems even more promising to treat pathogenic cases. It uses a Lanczos-type algorithm that is preconditioned via a Rayleigh method.

13.1.7.4. ‘Implicit Restarted Arnoldi’ (IRA) algorithm

This algorithm was initiated by DCSorensen 70 in 1992 and is experiencing a boom for the solution of large modal systems on parallel supercomputers. Its scope of application is quite general. It deals with both real and complex problems, Hermitian or not.

The algorithm involves a succession of Arnoldi factorizations whose results automatically control static and implicit restarts, via polynomial filters modeled by implicit QRs.

First, an Arnoldi factorization of order \(m = p + q\) (in theory \(q = 2\) is sufficient, in practice \(q = p\) is preferable) is performed. It is this default value which is used (Section 13.1.7.5) for the working matrix.

Once this preprocessing is done, an iterative filtering process for the unwanted part of the spectrum is performed (numerically and methodologically it is easier to exclude than to incorporate).

Th algorithm then determines the spectrum of the Rayleigh matrix (via QR) (and removes the unwanted part using criteria set by the user) which it then uses like a shift for a series of \(q\) implicit QR with simple shift (Section 13.1.13.1).

The factorization is then:

\[\mathbf{A}_\sigma \mathbf{Q}_m^+ = \mathbf{Q}_m^+ \mathbf{B}_m^+ + \mathbf{R}_m \mathbf{Q}\]

where \(\mathbf{Q}_m^+ = \mathbf{Q}_m \mathbf{Q}\), \(\mathbf{B}_m^+ = \mathbf{Q}^T \mathbf{B}_m \mathbf{Q}\) and \(\mathbf{Q} = \mathbf{Q}_1 \mathbf{Q}_2 \dots \mathbf{Q}_q\) is the unitary matrix associated with the QRs. After having updated the matrices, they are truncated up to the order \(p\)

\[\mathbf{A}_\sigma \mathbf{Q}_p^+ = \mathbf{Q}_p^+ \mathbf{B}_p^+ + \mathbf{R}_p^+\]

and thus, at the cost of \(q\) new Arnoldi iterations, we can find an Arnoldi factorization of order \(m\) that is viable. All the subtlety of the process is based on this last sequence.

Note

\[\Phi(\mathbf{A}_\sigma) = \prod_{i=1}^q ( \mathbf{A}_\sigma-\tilde{\lambda}_{p+i} \mathbf{I})\]

is the matrix polynomial of order \(q\) generated by the work operator. In fact, the implicit QRs act on the first \(p\) rows of the Arnoldi, Rayleigh and residual matrix, so that the complementary factorization of order \(q\) produces the same effect as a factorization of order \(m\) initiated by vector

\[\tilde{\mathbf{q}}_1 = \Phi(\mathbf{A}_\sigma) \mathbf{q}_1\]

The initial vector has been implicitly modified (there is no construction and effective application of the polynomial) in order that it generates preferentially the desired modes, by subtracting the components judged “unsuitable”.

As one already pointed out this type of restart makes it possible to decrease the residue by reducing the components of the initial vector by the unwanted modes. In exact arithmetic, we would immediately have \(\mathbf{R}_m = 0\) and the process would end there!

Iteration after iteration, one thus improves the quality of the modes sought by leaning on auxiliary modes. Of course, this is estimated at each step and determines whether another restart is needed or not.

The process can be summarized as follows:

Algorithm 13. IRA method (known as Sorensen).

  • Arnoldi factorization order \(m\) : \(\mathbf{A}\mathbf{Q}_m = \mathbf{Q}_m\mathbf{B}_m + \mathbf{R}_m\)

  • For \(k = 1\), NMAX_ITER_SOREN do

    • Calculate \(\underbrace{\tilde{\lambda}_1,\tilde{\lambda}_2 \dots \tilde{\lambda}_p}_{\text{Kept for improvement}}\), \(\underbrace{\tilde{\lambda}_{p+1} \dots \tilde{\lambda}_{p+q}}_{\text{Used as shifts}}\)

    • QR with implicit shifts:

    • Update \(\mathbf{Q}_m\), \(\mathbf{B}_m\) and \(\mathbf{R}_m\)

    • Truncate these matrices to order \(p\):

    • \(\mathbf{A}\mathbf{Q}_p = \mathbf{Q}_p \mathbf{B}_p + \mathbf{R}_p \implies \mathbf{A}\mathbf{Q}_m = \mathbf{Q}_m \mathbf{B}_m + \mathbf{R}_m\)

    • Quality estimation of the \(p\) modes

  • End loop.

The maximum number of iterations is controlled by the keyword NMAX_ITER_SOREN of the factor keyword CALC_FREQ. It should be noted that the Arnoldi algorithm is carefully supplemented by a total reorthogonization (triggered only if necessary). This additional cost is all the more acceptable as it is implemented viathe Kahan-Parlett IGSM algorithm (see Section 13.1.14) which is particularly efficient. All this allows us to ensure the good behavior of the projection vis-à-vis the initial spectrum.

On the other hand, the evaluation of the quality of the modes is not carried out simply by constructing the \(p\) residues \(\|\mathbf{r}_i \|_2 = \|(\mathbf{A}_\sigma -\tilde{\lambda} \mathbf{I}) \tilde{\mathbf{u}}_i\|_2\) via Property 10. We have already mentioned that in the non-Hermitian case, they cannot suffice for this task, especially in the event of a strong lack of normality. To discharge it, without having use of other information 25 a priori , we use the previous property completed by a criterion due to Z. Baiet al 40:

\[| b_{m+1, m}| | \mathbf{e}_m^T \mathbf{x} | < \text{max} ( \varepsilon \| \mathbf{B}_m\|, \text{PREC_SOREN} |\tilde{\lambda} |)\]

where \(\varepsilon\) is the machine precision and PREC_SOREN is a keyword initialized under CALC_FREQ. The user has therefore a (partial) control of the quality of the modes, which they did not have with the other methods implemented in the code. Given the different norms used, this error is different from that resulting from total postprocessing (Section 13.1.3.6) which is displayed in the columns of results.

Note

  • The polynomial acceleration technique used is more efficient than that of Chebyshev, since the latter is explicit and requires \(m\) matrix-vector products.

  • To avoid the deterioration of the Ritz vectors (and thus of the approximate eigenvectors) due to eigenvalues ​​of very large magnitude (associated with the kernel of \(\mathbf{B}\) in “shift and invert”) a filter devised by Ericsson & Ruhe 46 was implemented. More robust techniques exist but they require a priori information concerning in particular the Jordan blocks associated with the kernel of the operator (cf. Meerbergen & Spence, 1996).

  • This algorithm can be seen as a truncated form of JCFFrancis’ implicit QR algorithm (Section 13.1.13.1).

25

To obtain an exact criterion, it would be necessary to be able to estimate the spectral conditioning of the invariant spaces and the angles they make between them. That is difficult to obtain, even retrospectively!

13.1.7.5. Implementation in Code_Aster

13.1.7.5.1. ARPACK

The ARPACK package (ARPACK for ARnoldi PACKage) implements the IRA method. It is a public product downloadable from the internet 86. The designers, D. Sorensen, R. Lehoucq and C. Yang from Rice University inHouston, wanted:

  • Easy to access: Fortran / C / C ++ interface via the “reverse communication” model.

  • Modular: based on LINPACK, LAPACK and BLAS libraries.

  • Rich in functionality : Schur decomposition, modular shifts, numerous spectral transformations, configurable stopping criteria.

  • Covering large scope of use : single / double precision, real / complex, matrix GEP symmetrical / unsymmetrical etc.

  • Efficient: to save CPU time and memory usage, the package is available in PARPACK parallelized version via MPI and BLACS. On the other hand, even with sequential ARPACK, the solution of linear systems can be carried out by a parallel linear solver. The memory management of large objects required by ARPACK is left to the care of the host code.

Its effectiveness is greatly increased by the use of BLAS Level 2 and 3 highly optimized and implementation of “reverse communication”. The user is therefore master of his data structures and of his processing procedures concerning the operator and the working dot product. It is he who provides ARPACK routines this type of information. This therefore makes it possible to manage best, with the tools and procedures of Code_Aster, the matrix-vector products and the solutions of linear systems.

For the moment, the memory management of the large objects required by ARPACK is entrusted to JEVEUX. It therefore benefits from the “out-of-core” faculties of Code_Aster 26 . On the other hand, the modal operators of Aster having access to the same linear solver than those of statics, they can also benefit from improvements in time and memory from parallelism (SOLVEUR / METHODE = ‘MUMPS’ or ‘MULT-FRONT’) [U4.50.01] [U2.08.06]). But note that these gains are limited to the linear solver part of modal computation. Achievable speed-ups are therefore less promising since this part can represent, in some cases, only around 50% of the calculation time 27.

The ARPACK site offers documentation (theoretical and user), examples of use and a download page. The finalized version of ARPACK (v2.5 of 08/27/96) used in Code_Aster seems to be the last. Despite thousands of uses in the world (academic and industrial), a proven reference 51 and a certain recognition (the “ARPACK applications” tab of the site and the links with other libraries such as PETSc and TRILINOS) this software development project was stopped in 1997. It has lasted a few more years via the C++ version (ARPACK++) of D. Sorensen and F. Gomez.

26

“Out-Of-Core” (OOC) means that the memory management software package can offload a part of the memory to disk (objects contained in the RAM). Unlike traditional “In-Core” memory management (IC) which keeps everything in RAM, which more quickly limits the size of problems.

27

In this extreme case of only 50% of the elapsed time spent in the linear solver part, the parallelism, even on tens of processors, will bring at best only a gain (or acceleration or ‘speed-up’) by a factor of 2 (in time)! In RAM memory, the gains are also limited by the ratio between the own needs of ARPACK and those of the linear solver.

13.1.7.5.2. Adaptations of Sorensen’s algorithm

To deal with generalized modal problems, this package proposes a whole series of spectral transformations, both real and complex, and Hermitian/non-Hermitian. In Hermitian, Sorensen’s algorithm is based on the Lanczos-QL pair (IRLM for Implicit Restarted Lanczos Method). The two approaches are not intended to deal with pseudo-scalar products linked to undefined matrices. More precisely, to at the same time circumscribe the numerical problems related to their rather heterogeneous properties in the code and, on the other hand, to ensure better overall robustness, we have chosen to work with the unsymmetrical version (IRAM with Arnoldi and QR), on the pair (operator-scalar product) :

\[\underbrace{(\mathbf{A} - \sigma \mathbf{B})^{-1}\mathbf{B}}_{\mathbf{A}_\sigma} \mathbf{u} = \frac{1}{\underbrace{\mu - \sigma}_{\lambda}} \mathbf{u} \quad \text{and} \quad (\mathbf{x}, \mathbf{y}) = \mathbf{y}^T \mathbf{x}\]

One could have treated the problems of buckling in “buckling mode” via the same “shift and invert” and the pseudo-scalar product induced by \(\mathbf{A}\). But due to the almost systematic introduction of Lagrange multipliers,this matrix becomes indefinite or even singular, which greatly disturbs the process. The same causes produce the same effects when, for a dynamic computation, we have recourse to the \(\mathbf{B}\)-scalar product. Rather than modifying the whole package by introducing a pseudo-dot product, we therefore opted for a simple “Euclidean” dot product that is more robust and much less expensive.

It was therefore necessary to modify the “reverse-communication” procedures of the package, because it did not provide for this option (with standard matrices one classically prefers to enrich the components with a matrix scalar product, even for unsymmetrical matrices). Unlike the variant of Newmann & Pipano introduced for Lanczos, we have deliberately placed ourselves in a non-symmetrical configuration. But in order to avoid as much as possible the recurring problems of orthogonality, even in symmetry, we would have opted for the version of IRAM using Arnoldi.

The disadvantage of this approach is that it is necessary, in preliminary postprocessing of IRAM, to \(\mathbf{B}\) -orthonormalizethe approximate eigenvectors to numerically find the Property 2 exploited by the model recombinations. This step does not disturb the basis of eigen modes and it is very effectively carried out via the Kahan-Parlett IGSM. Boundary conditions, in particular double dualisations, were applied as for Lanczos following the operating mode described in Section 13.1.3.2. In particular, once the initial vector is found in the admissible space, the work operator is applied to it. This classic process allows you to purge the Lanczos vectors (and therefore Ritz vectors) of components of the kernel.

On the other hand, in certain configurations for which the number of frequencies requested (\(p\)) is equal to the number of active degrees of freedom, we have had to trick the algorithm which stopped with a fatal error! Indeed, it generally detected the expected invariant space (of size \(p\)), but due to the particular structure of the eigenvectors associated with the Lagrange multipliers (cf. Proof of Property 4, Section 13.1.3.5) it is difficult to generate an initial vector which is proportional to them. It would have required special treatments taking into account the numbering of these Lagrange multipliers, which would have been all the more expensive since they are not fundamentally necessary to solve the requested problem!

The spectral information is already present in the depths of the algorithm, so there is no point in completing the two remaining iterations (when the user asks for \(p = n_{\text{dof - actives}}\), we automatically impose \(m = p +2\)). Just short-circuit the natural thread of the algorithm, remove the interesting Ritz modes and to post-process them to return to the initial space.

Note

  • The case where we seek a number of eigen modes very close to the number of degrees of freedom is outside the “ideal” scope of use of this type of algorithm. A good QR would undoubtedly be more efficient, but it’s a good way to test the algorithm.

13.1.7.5.3. Parameter setting

To be able to activate this method, the keyword METHOD should be initialized with SORENSEN. The size of the projection space is determined, either by the user, or empirically from the formula: \(m = \text{min} (\text{max} (2 p, p +2), n_{\text{dof - active}})\) where

  • \(p\) is the number of eigenvalues ​​to calculate

  • \(n_{\text{dof - active}}\) is the number of active degrees of freedom of the system (Section 13.1.3.2)

The user can always impose the dimension himself by indicating it with the key word DIM_SOUS_ESPACE of the keyword factor CALC_FREQ. We can prefer a default value, a multiplicative coefficient compared to thenumber of frequencies contained in the study interval. This type of functionality concerns the keyword COEF_DIM_ESPACE. The parameter of the IGSM of Kahan-Parlett (see Section 13.1.14) PARA_ORTHO_SOREN, the maximum number of iterations of the total process, NMAX_ITER_SOREN, and the quality control criterion of the modes, PREC_SOREN, are accessible by the user under the keyword factor CALC_FREQ. When this last keyword is null, the algorithm initializes it to machine precision.

13.1.7.6. Scope of use

GEP with unspecified real matrices (symmetrical or not) or with complex matrix \(\mathbf{A}\) and real symmetric \(\mathbf{B}\). The user can specify the type of calculation (dynamic or buckling if matrices are symmetric and real) by initializing the keyword TYPE_RESU. According to this value, one assigns the vector FREQ or CHAR_CRIT.

13.1.7.7. Display in message file

The example below from the list of code tests (ssll103b) summarizes messages produced by the algorithm. We find in particular, for each critical load (or frequency), the estimate of its quality via the error norm. Here IRAM iterated only once and used 30 IGSMs (in its first phase). The overall solution has consumed 91 (in fact, less than that due to the “implicit” introduction of the Euclidean scalar product) matrix-vector products and 31 matrix inversions (because the operator of work is already factorized).

------------------------------------------------------------
THE NUMBER OF DDLTOTAL IS:68
FROM LAGRANGE EAST:14
THE NUMBER OF ACTIVE DDL IS:47
-------------------------------------------------- ----------
THE OPTION CHOSEN IS: PLUS_SMALL
THE CRITICAL LOAD OFFSET VALUE IS: 0.00000E + 00
-------------------------------------------------- ----------
INFORMATION ON THE CALCULATION REQUIRED :
NUMBER OF MODES REQUESTED: 10
THE DIMENSION OF THE REDUCED SPACE IS: 30
==============================================
= SORENSEN METHOD (ARPACK CODE) =
= VERSION: 2.4=
= DATE: 07/31/96=
==============================================
NUMBER OF RESTART= 1
NUMBER OF PRODUCTS OP * X= 31
NUMBER OF PRODUCTS B * X= 91
NUMBER OF REORTHOGONALIZATIONS (STEP 1) = 30
NUMBER OF REORTHOGONALIZATIONS (STEP 2) = 0
NUMBER OF RESTART DUE TO A NULL V0 = 0
-------------------------------------------------- ----------
THE CALCULATED CRITICAL LOADS INF. AND SUP. ARE:
CRITICAL_LOAD_INF: -9.96796E + 06
CRITICAL_LOAD_SUP: -6.80007E + 05
-------------------------------------------------- ----------
MODAL CALCULATION: SIMULTANEOUS ITERATION METHOD
SORENSEN METHOD
NUMBER CRITICAL LOAD  ERROR NORM
1 -6.80007E + 05 5.88791E-12
2 -7.04572E + 05 1.53647E-12
3 -7.09004E + 05 1.16735E-12
...
10 -9.96796E + 06 3.55014E-12
-------------------------------------------------- ----------
POSTERIORI CHECKING OF THE MODES
IN THE INTERVAL (-1.00178E + 07, -6.76607E + 05)
THERE ARE A LOT 10 CRITICAL LOAD (S)
-------------------------------------------------- ----------

Note

The introduction of this method made it possible to settle many anomaly files linked to multiplicities, clusters or searches for eigenvalues ​​of very different orders of magnitude on which Lanczos and Bathe & Wilson were stumbling.

13.1.7.8. Summary of settings

Let us recapitulate now the parameter setting available of the operator CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’] with this algorithm.

Table 13.1.6 Parameters for CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’] for SOLVEUR_MODAL = _F(METHODE = ‘SORENSEN’)

Keyword factor

Keyword

Default value

References

TYPE_RESU = ‘DYNAMIQUE’

‘DYNAMIQUE’

Section 13.1.3.1

‘MODE_FLAMB’

Section 13.1.3.1

OPTION = ‘PLUS_*’

‘PLUS_PETITE’

Section 13.1.5.4

‘BANDE’

Section 13.1.5.4

‘CENTRE’

Section 13.1.5.4

SOLVEUR_MODAL

METHODE = ‘SORENSEN’

‘SORENSEN’

Section 13.1.7.4

PREC_SOREN

Section 13.1.7.4, Section 13.1.7.5

NMAX_ITER_SOREN

20

Section 13.1.7.4, Section 13.1.7.5

PARA_ORTHO_SOREN

0.717

Section 13.1.7.5, Section 13.1.14

DIM_SOUS_ESPACE

calculated

Section 13.1.7.5

COEF_DIM_ESPACE

calculated

Section 13.1.7.5

CALC_FREQ or CALC_CHAR_CRIT

FREQ or CHAR_CRIT (if ‘CENTRE’)

Section 13.1.5.4

NMAX_FREQ or NMAX_CHAR_CRIT (if ‘PLUS_*’)

10

Section 13.1.5.4

AMOR_REDUIT (if ‘CENTRE’ and FREQ)

Section 13.1.5.4

NMAX_ITER_SHIFT

3

81 Section 13.1.3.2

PREC_SHIFT

0.05

81 Section 13.1.3.2

SEUIL_FREQ or SEUIL_CHAR_CRIT

1.E-02

Section 13.1.3.7

VERI_MODE

STOP_ERREUR = ‘YES’

‘YES’

Section 13.1.3.7

‘NO’

Section 13.1.3.7

PREC_SHIFT

5.E-03

Section 13.1.3.7

SEUIL

1.E-06

Section 13.1.3.7

STURM = ‘YES’

‘YES’

Section 13.1.3.7

‘NO’

Section 13.1.3.7

Note

  • We list all the parameters related to the Strum test (NMAX_ITER_SHIFT, PREC_SHIFT) and verification tests (SEUIL_FREQ / SEUIL_CHAR_CRIT, VERI_MODE).

  • As always, beginners should only modify the parameters in bold and void changing the other parameters.

  • In particular, to improve the quality of a mode, the fundamental parameter is the dimension of the sub-space (DIM_SOUS_ESPACE / COEF_DIM_ESPACE)

13.1.8. Method of Bathe and Wilson (SOLVEUR_MODAL = _F(METHOD = ‘JACOBI’))

13.1.8.1. Principle

The method of Bathe and Wilson uses simultaneous iterations which extend the inverse iteration algorithm. We work with the shifted problem \(\mathbf{A}^\sigma \mathbf{x} := (\lambda - \sigma) \mathbf{B} \mathbf{x}\) with \(\mathbf{A}^\sigma := \mathbf{A} - \sigma \mathbf{B}\) (Section 13.1.5.4). There are four parts in the algorithm 37:

  • Choose \(p\) initial independent vectors \(\mathbf{x}_1, \mathbf{x}_2, \dots \mathbf{x}_p\) and build the matrix \(\mathbf{X}\) that they generate.

  • Calculate the eigen dcomposition, via the Jacobi global method (Section 13.1.15), in the Ritz subspace by solving

    \[(\tilde{\mathbf{A}} - \tilde{\lambda}_i \tilde{\mathbf{B}}) \mathbf{u}_i = \mathbf{0}\]

    where \(\tilde{\mathbf{A}} = \mathbf{Q}^T \mathbf{A}^\sigma \mathbf{Q}\) and \(\tilde{\mathbf{B}} = \mathbf{Q}^T \mathbf{B} \mathbf{Q}\)

  • One then returns to the initial space (for the eigenvectors) via the transformation

    \[\mathbf{X} = \mathbf{Q}\mathbf{U}\]

    where \(\mathbf{U} = [\{ \mathbf{u}_i \}]\)

  • Test the convergence of the eigen modes \(\lambda_i\).

13.1.8.2. Convergence tests

The method of Bathe and Wilson converges towards the \(p\) smallest eigenvalues ​​provided that the \(p\) initial vectors are not \(\mathbf{B}\) -orthogonal to one of the eigenvectors. Moreover, the matrices \(\tilde{\mathbf{A}}\) and \(\tilde{\mathbf{B}}\) tend towards diagonal matrices. For this reason and as the matrices are dense, we use Jacobi’s method (see Section 13.1.15) to find the eigen decomposition of the Ritz subspace.

To test the convergence of the eigenvalues, one classifies them after each iteration in ascending order in absolute value and we see if, for each eigenvalue, the following test is satisfied:

\[|\lambda^{k+1} - \lambda^k| \le \text{PREC_BATHE} |\lambda^{k+1} |\]

where the exponent \(k\) indicates the number of iterations. If after NMAX_ITER_BATHE iterations, one does not converge for all the eigenvalues, an alarm message is emitted in the message file.

13.1.8.3. Implementation in Code_Aster

13.1.8.3.1. Dimension of the subspace

If one wishes to calculate \(p\) eigenvalues, it is recommended to use a subspace of dimension larger than \(q\). One will check the convergence only for the \(r\) smallest eigenvalues ​​where \(p \le r \le q\). It seems that \(r = p\) is not sufficient: we can find good eigenvalues ​​but the eigenvectors are not correct (convergence is slower for the eigenvectors than for the eigenvalues).

The value \(r = (p + q) / 2\) seems a good choice.

For \(q\) we usually take 37 \(q = \text{min}(p +8, 2 p)\).

13.1.8.3.2. Choice of the initial vectors

To choose the \(q\) initial vectors, we operate as follows:

  • Choose the first vector such that

    \[\mathbf{x}_{1i} = \frac{\mathbf{A}_{ii}^\sigma}{\mathbf{B}_{ii}}`\]

  • For the other vectors from \(2\) to \((q-1)\):

    \[\begin{split}\mathbf{x}_2= \begin{bmatrix} \mathbf{0} & \\ 1 & \quad (\text{row}~~ i_1) \\ \mathbf{0} & \end{bmatrix} \quad \mathbf{x}_{q-1} = \begin{bmatrix} \mathbf{0} & \\ 1 & \quad (\text{row}~~ i_q) \\ \mathbf{0} & \end{bmatrix} \dots\end{split}\]

    where the \(i\) are the indices corresponding to the smallest successive values ​​of \(\mathbf{A}_{ii}^\sigma / \mathbf{B}_{ii}\)

  • The last vector \(\mathbf{x}_q\) is a random vector.

13.1.8.4. Scope of use

GEP with symmetric real matrices.

The user can specify the type of computation with TYPE_RESU and the vectors FREQ or CHAR_CRIT.

13.1.8.5. Summary of settings

The parameter settings for the operator CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’] when using this algorithm are listed in the table below. To be able to use the method of Bathe and Wilson, it is necessary to select METHODE = ‘JACOBI’ under the keyword factor SOLVEUR_MODAL. Both parameters directly concerning the convergence of the method are accessible under this keyword factor using keywords PREC_BATHE and NMAX_ITER_BATHE. There are also those that manage the internal of modal solution process, PREC_JACOBI and NMAX_ITER_JACOBI.

Table 13.1.7 Summary of the parameters of CALC_MODES with OPTION=[‘BANDE’, ‘CENTRE’, ‘PLUS_*’] (GEP) with SOLVEUR_MODAL = _F(METHOD = ‘JACOBI’).

Keyword factor

Keyword

Default value

References

TYPE_RESU = ‘DYNAMIQUE’

‘DYNAMIQUE’

Section 13.1.3.1

‘MODE_FLAMB

‘OPTION = ‘PLUS_*’

‘PLUS_PETITE’

Section 13.1.5.4

‘BANDE’

‘CENTRE’

SOLVEUR_MODAL

METHODE = ‘JACOBI’

‘SORENSEN’

Section 13.1.15

PREC_BATHE

1.E-10

Section 13.1.8.2

NMAX_ITER_BATHE

40

Section 13.1.8.2

PREC_JACOBI

1.E-12

Section 13.1.15

NMAX_ITER_JACOBI

12

Section 13.1.15

DIM_SOUS_ESPACE

Calculated

Section 13.1.8.3

COEF_DIM_ESPACE

Calculated

Section 13.1.8.3

CALC_FREQ or CALC_CHAR_CRIT

FREQ or CHAR_CRIT (if ‘CENTRE’)

Section 13.1.5.4

NMAX_FREQ

10

Section 13.1.5.4

NMAX_ITER_SHIFT

3

81 Section 13.1.3.2

PREC_SHIFT

0.05

81 Section 13.1.3.2

SEUIL_FREQ or SEUIL_CHAR_CRIT

1.E-02

Section 13.1.3.7

VERI_MODE

STOP_ERREUR = ‘OUI’

‘OUI’

Section 13.1.3.7

‘NON’

PREC_SHIFT

5.E-03

Section 13.1.3.7

SEUIL

1.E-06

Section 13.1.3.7

STURM = ‘OUI’

‘OUI’

Section 13.1.3.7

‘NON’

Section 13.1.3.7

Note

  • The table contains all parameters related to the Sturm test (NMAX_ITER_SHIFT, PREC_SHIFT) and with verification checks (SEUIL_FREQ / SEUIL_CHAR_CRIT, VERI_MODE).

  • Beginners are encouraged to modify only the parameters in bold and to retain the rest of the empirically derived parameters values without modification.

13.1.9. Global QZ method (SOLVEUR_MODAL =_F(METHOD = ‘QZ’))

13.1.9.1. Introduction

Most modal methods use a spectral “shift and invert” transformation to facilitate the search for part of the spectrum and to transform the initial GEP into a classic SEP (Section 13.1.3.7). This approach is often adequate but it involves factorization of matrices of the type \(\mathbf{A}_\sigma = (\mathbf{A} -\sigma \mathbf{B})^{-1}\mathbf{B}\) (Section 13.1.3.6). The result of this factorization is only an intermediate step performed many times in modal solvers. However these factorizations can have numerical errors related to the conditioning of the matrices and to the matrix inversion error 28 of the direct solver (see [R6.02.03], 77, Section 13.1.2.3). Hence the need for a robust modal method to avoid these issues and only use:

  • Transformations of the Householder or Givens type, that produce / spread some rounding errors (called “backward stable”)

  • Easy to handle canonical matrices (triangular, diagonal or Hessenberg) and / or very well conditioned (orthogonal in real or unitary in complex space).

28

The inverse error measures the propensity of the algorithm to transmit / amplify rounding errors.

13.1.9.2. GEP and the QZ method

The QZ method directly treats a GEP rather than an SEP resulting from its modification via a spectral transformation.

GEP

Find \((\lambda, \mathbf{u})\) such that

(13.1.7)\[(\mathbf{A} - \lambda\mathbf{B}) \mathbf{u} = \mathbf{0}\]

This method seeks to construct the following generalized Schur decomposition (counterpart of the standard Schur decomposition, see (13.1.5)).

Theorem 12 (Generalized Schur Decomposition)

Let \(\mathbf{A}\) and \(\mathbf{B}\) be two complex square matrices. Then there exist two unit square matrices \(\mathbf{Q}\) and \(\mathbf{Z}\) such that \(\mathbf{Q} ^\star\mathbf{A}\mathbf{Z} = \mathbf{T}\) and \(\mathbf{Q}^\star\mathbf{B}\mathbf{Z} = \mathbf{S}\) are upper triangular matrices. The eigenvalues ​​of the GEP (13.1.7) are then calculated easily thanks to the diagonal terms (complex) of the \(\mathbf{T}\) and \(\mathbf{S}\) matrices :

(13.1.8)\[\begin{split}\lambda(\mathbf{A}, \mathbf{B}) := \left\{ \lambda_i := \begin{cases} \mathbf{T}_{ii}/ \mathbf{S}_{ii} & \quad \text{if} \quad \mathbf{S}_{ii} \ne 0 \\ \infty & \quad \text{if} \quad \mathbf{S}_{ii} = 0 \end{cases} \right\}\end{split}\]

Proof:See 48 pp396 and 59.

This theorem indicates that, once the Generalized Schur Decomposition has been found, solving the GEP (13.1.7) amounts to solving the equivalent GEP (which is much simpler because it triangular)

Equivalent GEP

Find \((\lambda, \mathbf{v})\) such that

(13.1.9)\[(\mathbf{T} - \lambda\mathbf{S}) \mathbf{v} = \mathbf{0}\]

They have the same spectrum, \(\lambda(\mathbf{A}, \mathbf{B}) = \lambda(\mathbf{T}, \mathbf{S})\), and we pass eigenvectors from one to those of the other, via the canonical transformation

(13.1.10)\[\mathbf{u} = \mathbf{Z} \mathbf{v}\]

On the other hand, the manipulation of pairs of the type \((\alpha_i, \beta_i) = (\mathbf{T}_{ii}, \mathbf{S}_{ii})\), rather than the direct manipulation of an eigenvalue \(\lambda_i := \frac{\alpha_i}{\beta_i}\), is a numerical trick to better understand the complexity induced by the direct solution of a GEP.

Considering, for example, the rank of the matrix \(\mathbf{B}\), the spectrum can be:

  • finite, \(\lambda(\mathbf{A}, \mathbf{B}) = \{\lambda_1 \dots \lambda_n\}\), with possibly very large eigenvalues \(\pm \infty\)

  • empty, \((\lambda(\mathbf{A}, \mathbf{B}) = \emptyset)\)

  • infinite \((\lambda(\mathbf{A}, \mathbf{B}) = \mathbb{R}\) or \(\mathbb{C})\).

In the first two cases, the system is said to be “regular” as opposed to the last situation where it is said to be “singular”. The examples below illustrate the three scenarios:

(13.1.11)\[\begin{split}\mathbf{A} & := \begin{bmatrix} 1 & 2 \\ 0 & 3 \end{bmatrix},~~ \mathbf{B} := \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix} \implies \lambda (\mathbf{A}, \mathbf{B}) = \{ 1 \} \\ \mathbf{A} & := \begin{bmatrix} 1 & 2 \\ 0 & 3 \end{bmatrix},~~ \mathbf{B} := \begin{bmatrix} 0 & 0 \\ 1 & 0 \end{bmatrix} \implies \lambda (\mathbf{A}, \mathbf{B}) = \emptyset \\ \mathbf{A} & := \begin{bmatrix} 1 & 2 \\ 0 & 0 \end{bmatrix},~~ \mathbf{B} := \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix} \implies \lambda (\mathbf{A}, \mathbf{B}) = \mathbb{C}\end{split}\]

The representation of the solutions as \((\alpha, \beta)\) then makes it possible to manage these cases more easily. For example, the value \(\beta \sim 0\) can be interpreted as an infinite eigenvalue \(\lambda \sim \infty\) (probably not to retain in the work spectrum!). If \(\alpha ≈ 0\), we have an indeterminacy in the corresponding eigenvalue and this indicates a singular or very poorly conditioned system. In the case of the GEPs solved by Code_Aster this case appears when one treats the eigenvalues associated with constrained dofs and those of Lagrange multipliers (Section 13.1.3.2). The decomposition \((\alpha, \beta)\) then gives a framework to filter out these numerical artifacts from the physical spectrum:

(13.1.12)\[\lambda_{\text{physical}} (\mathbf{A}, \mathbf{B}) := \{ -\infty \ll | \lambda_i | \ll +\infty,~~ i = 1 \dots n_{\text{dof active}} \}\]

Note

  • The QZ method is “backward stable” for a standard GEP. But it loses this property when it is used on a linearized GEP (Section 13.1.2.1) which is an intermediate computation for a QEP (72 Section 13.1.5.1).

  • In the case of real matrices (most frequent case in Code_Aster), there is a similar decomposition similar called “Generalized Real” with orthogonal square matrices \(\mathbf{Q}\) and \(\mathbf{Z}\) such that \(\mathbf{Q}^T\mathbf{A}\mathbf{Z} = \mathbf{T}\) is upper quasi-triangular (i.e., block diagonal with \(1 \times 1\) or \(2 \times 2\)) and \(\mathbf{Q}^T \mathbf{B}\mathbf{Z} = \mathbf{S}\) is upper triangular. The philosophy of calculation remains the same, by introducing pairs \((\alpha, \beta)\) as intermediaries, before the actual calculation of the eigenvalues \(\lambda\). A notable difference in the complex case comes from the fact that we can, even in the presence of complex modes, continue to work in real arithmetic. This strategy is exploited by certain routines of LAPACK and is used in the QZ method of CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’].

  • In certain cases (SPD matrices, interval structures, etc.), the approach can be applied more effectively. For example, if the GEP is symmetric with real matrices and \(\mathbf{B}\) is, in addition, positive definite (thus there exists a Cholesky factorization such that \(\mathbf{B} = \mathbf{L}\mathbf{D}\mathbf{L}^T\)), the GEP (13.1.7) is solved directly via the equivalent real symmetric SEP:

    GEP

    Find \((\lambda, \mathbf{u})\) such that \((\mathbf{A} - \lambda \mathbf{B}) \mathbf{u} = 0\)

    SEP

    Find \((\lambda, \mathbf{u})\) such that

    (13.1.13)\[(\mathbf{L}^{-1} \mathbf{A} \mathbf{L}^{T})(\mathbf{L}^T\mathbf{u}) = \lambda (\mathbf{D}\mathbf{L}^T\mathbf{u})\]

  • This canonical SEP is then solved in a robust and more efficient way via an adapted QR (see Section 13.1.13). This strategy exploited by certain routines of LAPACK and used in the QZ method of CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’] and (SOLVEUR_MODAL = _F(METHOD = ‘QZ’, TYPE_QZ = ‘QZ_QR’)).

13.1.9.3. The QZ method

In general, the QZ algorithm can be broken down into three stages:

  • Using either a Givens or a \(2 \times 2\) Householder transformation (Section 13.1.15), the two complex (or real) matrices \(\mathbf{A}\) and \(\mathbf{B}\) are broken down into a form known as “Generalized Upper Hessenberg”. That is, we determine two unitary square matrices (or orthogonal for real matrices) \(\mathbf{U}_1\) and \(\mathbf{V}_1\) such that \(\mathbf{A} = \mathbf{U}_1 \mathbf{H} \mathbf{V}_1^\star\) and \(\mathbf{B} = \mathbf{U}_1 \mathbf{R} \mathbf{V}_1^\star\) (for real matrices, \(\mathbf{A} = \mathbf{U}_1 \mathbf{H} \mathbf{V}_1^T\) and \(\mathbf{B} = \mathbf{U}_1 \mathbf{R} \mathbf{V}_1^T\)) with upper Hessenberg matrix \(\mathbf{H}\) and upper triangular \(\mathbf{R}\).

    For example, with

    \[\begin{split}\mathbf{A} := \begin{bmatrix} 10 & 1 & 2 \\ 1 & 2 & -1 \\ 1 & 1 & 2 \end{bmatrix} \quad \text{and} \quad \mathbf{B} := \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{bmatrix}\end{split}\]

    we find

    \[\begin{split}\mathbf{U}_1 & := \begin{bmatrix} -0.12 & -0.99 & 0.03 \\ -0.49 & -0.02 & -0.86 \\ -0.86 & 0.12 & 0.49 \end{bmatrix}~,~~\mathbf{V}_1 := \begin{bmatrix} 1.00 & 0.00 & 0.00 \\ 0.00 & -0.89 & -0.44 \\ 0.00 & 0.44 & -0.89 \end{bmatrix} \\ \mathbf{H} & := \begin{bmatrix} -2.58 & 1.54 & 2.42 \\ -9.76 & 0.08 & 1.92 \\ 0.00 & 2.72 & -0.76 \end{bmatrix}~,~~ \mathbf{R} := \begin{bmatrix} -8.12 & 3.63 & 14.20 \\ 0.00 & 0.00 & 1.87 \\ 0.00 & 0.00 & 0.00 \end{bmatrix}\end{split}\]

  • By generalizing Francis’ implicit double-shift algorithm (based on Householder transformations, Section 13.1.13.1), we finalize the transformation. We then generate two square unitary matrices (orthogonal for real matrices) \(\mathbf{U}_2\) and \(\mathbf{V}_2\) such that \(\mathbf{H} = \mathbf{U}_2 \mathbf{t} \mathbf{V}_2^\star\) and \(\mathbf{R} = \mathbf{U}_2 \mathbf{S} \mathbf{V}_2^\star\) (for real matrices, \(\mathbf{H} = \mathbf{U}_2 \mathbf{T} \mathbf{V}_2^T\) and \(\mathbf{R} = \mathbf{U}_2 \mathbf{S} \mathbf{V}_2^T\)) with upper triangular \(\mathbf{T}\) and \(\mathbf{S}\) (for real matrices, \(\mathbf{T}\) is almost upper triangular). By combining these two steps, we establish the desired result, that the product of unitary transformations remains unitary:

    (13.1.14)\[\underbrace{(\mathbf{U}_1 \mathbf{U}_2)}_{\mathbf{Q}^T} \mathbf{A} \underbrace{(\mathbf{V}_1 \mathbf{V}_2)}_{\mathbf{Z}} = \mathbf{T} \quad \text{and} \quad \underbrace{(\mathbf{U}_1 \mathbf{U}_2)}_{\mathbf{Q}^T} \mathbf{B} \underbrace{(\mathbf{V}_1 \mathbf{V}_2)}_{\mathbf{Z}} = \mathbf{S}\]

  • By using the diagonal terms of the triangular matrices \(\mathbf{T}\) and \(\mathbf{S}\), one deduces (after filtering and some tests of validity) the physical eigenvalues ​​of the model vibratory problem. The right and left generalized Schur vectors (column vectors of \(\mathbf{Q}\) and \(\mathbf{Z}\)) allow us to recover the eigenvectors of the equivalent GEP (13.1.9) and then the original eigenvectors (via transformation (13.1.10)).

The QZ method is recommended when one seeks, in a robust way, the full the spectrum of a GEP. However, it is costly in CPU time \(O(30 n^3)\) against \(O(nm(m-p))\) for IRAM with for a matrix of size \(n\), \(p\) desired eigen modes, and \(m\) the size of the projection space. Further, it is very greedy in memory because it requires the full storage of the matrices. Hence a memory complexity of \(O(n^2)\) against, for example, \(O(2np + p^2)\) for IRAM. In short, this algorithm is reserved for small, dense, GEPs (less than 1000 degrees of freedom) for which we want to obtain a reliable estimate of the spectrum. As for the QR algorithm in Lanczos / IRAM or for Jacobi with Bathe & Wilson, the QZ method is exploited as an elementary part of an algorithmic process which will have (drastically) reduced the size of the initial problem.

For further information on the method, one can consult GHGolub 48 or the initial paper by CBMoler and GWStewart 59 (1973).

13.1.9.4. Implementation in Code_Aster

13.1.9.4.1. LAPACK

The linear algebra library LAPACK 87 (Section 13.1.2.3) implements dedicated routines (and drivers encapsulating them) for the direct solution of GEPs via the QZ method for a large set of problems: SPD, real non-symmetrical, complex Hermitian or any, single / double precision, full or partial storage etc. The library also offers a range of routines to break down the solution, pre-and post-process the data or even estimate numerical errors. However, it uses full storage of the matrices and only works “In-Core” and does not exploit any parallelism 29.

On the other hand, it does not appear that there is any public domain package that implements 30 the QZ method more efficiently than LAPACK. Code_Aster therefore uses this library for its implementation of the QZ method in CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’]. Note that many other products have made the same choice of LAPACK and its QZ: the eigenvalues ​​module of Python / linear algebra, eig and qz functions in Matlab, spec for Scilab, MSC / Nastran.

29

LAPACK has a parallel version (with a different scope): ScaLAPACK [Sca].

30

That is to say with sparse storage, or even in parallel mode.

13.1.9.4.2. Integration and post-verification

In the implementation which was used in Code_Aster, the GEP with real modes (real symmetric matrices) are processed in real arithmetic via the routines DGGEV and DGGEVX. The first routine is used during a standard run (TYPE_QZ = ‘QZ_SIMPLE’, default value). The second (TYPE_QZ = ‘QZ_EQUI’) is more sophisticated and reserved for difficult resolutions. In particular, it balances and permutes the terms of the input matrices in order to improve the spectral conditioning and thus decreases the numerical error in the solutions. Sometimes this pre-treatment can be detrimental by disproportionately increasing the number of terms that are close to machine precision, thus distorting the results. For this reason, this option is not enabled by default.

In the very particular case of a real symmetric, positive definite 31, GEP with matrix :math`mathbf{B}`, the LAPACK DSYGV driver is more efficient (see (13.1.13)) than DGGEV / X. It must be activated explicitly via the keyword TYPE_QZ = ‘QZ_QR’. For certain situations (buckling, double Lagrange, complex or non-symmetric matrix, linearized QEP), the non-SPD character of \(\mathbf{B}\) can be detected a priori. A fatal error then stops computation if the value QZ_EQUI was chosen.

In the third and last case of GEP with complex modes (symmetrical complex matrices 32 and / or real unsymmetric), the solution is directly carried out in complex arithmetic via the LAPACK drivers ZGGEV and ZGGEVX. The first is the standard default mode (activated by TYPE_QZ = ‘QZ_SIMPLE’) and the second, the expert mode (TYPE_QZ = ‘QZ_EQUI’).

At the end of modal calculation, QZ returns to Code_Aster pairs :math:(alpha_i , beta_i)` and their associated eigenvector \(\mathbf{u}_i\), which are sorted and filtered as follows:

  • Mechanics checks: If \(| \beta_i|\) is close to machine precision, the mode is not retained (it probably corresponds to a constrained dof or a constraining Lagrange multiplier) and an alarm message is sent (if INFO = 2).

    Once all the sorting has been carried out, one makes sure that the number of physically acceptable modes (13.1.12) corresponds well to the number of active dofs \(n_{\text{dof - actives}}\). This is determined in the initialization phase and is independent of the modal solver (Section 13.1.3.2). If the GEP has real modes, the non-satisfaction of this criterion leads to a fatal error. In other cases, it is only signaled by a message. It is therefore not necessarily synonymous with malfunctions. It all depends on subsequent sorting of the conjugate complex eigenvalues!

  • Numerical checks: If \(| \alpha_i |\ge \| \mathbf{A} \|\) or \(| \beta_i |\ge\| \mathbf{B} \|\) , one does not retain the mode and an alarm message is issued.

    In the case of a real symmetric GEP (with TYPE_QZ = ‘QZ_SIMPLE’ / ‘QZ_EQUI’), one checks if the eigenvalues ​​have a very small imaginary part (less than the \(f_{\text{corrected}}\) value informed by SEUIL_FREQ, default value equal to 1.0E-2). Otherwise, an alarm is issued.

    In the case of a real positive definite symmetric GEP treated via ‘QZ_QR’, one excludes the eigenvalues less than \(-2 f_{\text{corrected}}\) and greater than 0.5E308 (0.5 \(\times\) maximum value representable in the machine 33). These arbitrary values ​​result from numerical experiments on modal benchmarks of the Aster base. When an eigenvalue is found outside this interval, an alarm is emitted if INFO = 2.

Once the modes resulting from QZ have been filtered, these will be sorted and reordered according to the desiderata of the user:

  • CALC_FREQ / OPTION = ‘BANDE’: lawful option only with a GEP with real modes. We select the eigenvalues ​​located in the frequency interval desired by the user (keyword FREQ / CHAR_CRIT).

  • CALC_FREQ / OPTION = ‘PLUS_PETITE’: one retains the smallest eigenvalue in magnitude.

  • CALC_FREQ / OPTION = ‘CENTRE’: one selects the closest eigenvalues ​​in magnitude of shift parameterized by the keyword FREQ / CHAR_CRIT. We retain NMAX_FREQ.

  • CALC_FREQ / OPTION = ‘TOUT’: one keeps all the legal eigenvalues ​​resulting from the previous filters.

To finish, as for all the modal solvers of the code (Section 13.1.3.6), the selected modes will undergo two checks:

  • Test on the norm of the residual of the mode \((\lambda, \mathbf{u})\).

  • In the standard case of a GEP with real modes, Sturm test of counting of eigenvalues ​​in frequency interval.

Note

This work of integration of a QZ solver in Code_Aster continues and extends the developments undertaken by M.Begon 75. It was carried out to answer 78, in particular regarding the new robustness needs (on small problems) of the dynamic functionalities of the CADYRO code 34.

31

For example a dynamic calculation without hysteretic damping with boundary condition constraints modeled only by elimination (without recourse to the double Lagrange).

32

This case applies for the Aster hysteretic damping model.

33

Returned by the Aster environment variable R8MIEM.

34

Software performing vibration analyzes on shaft lines of rotating machines.

13.1.9.5. Range of applicability

GEP with unspecified real matrices (symmetrical or not) or with complex matrix \(\mathbf{A}\) and real symmetrical \(\mathbf{B}\).

The user can specify the type of calculation (dynamic or buckling) with TYPE_RESU. Depending on this value, the use needs to provide either of the vectors FREQ or CHAR_CRIT.

13.1.9.6. Display in message file

The message file contains information relating to the chosen method (here QZ), to its variant (here ‘QZ_SIMPLE’) and with the selected modes. In the most common case of a GEP with real eigen modes, we specify the list of frequencies \(f_i\) retained (FREQUENCY) and their error norm (NORMERROR).

A trace of CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’] in the file message with SOLVEUR_MODAL = _F(METHODE = ‘QZ’) is shown below. This is an extract from the forma11a benchmark.

-------------------------------------------------
MODAL CALCULATION: GLOBAL QR TYPE METHOD
QZ_SIMPLE ALGORITHM

NUMBER FREQUENCY (HZ) ERROR NORM
1      1.67638E+02   2.47457E-11
2      1.67638E+02   1.48888E-11
3      1.05060E+03   2.00110E-12
4      1.05060E+03   1.55900E-12
.....

For a GEP with complex modes, unlike the QEP 79, Code_Aster does not filter the eigenvalue conjugates. It keeps them and displays them in ascending order of the real part. Thus, the columns FREQUENCY and AMORTIZATION group together, respectively, \(\text{Re}( f_i )\) and \(\frac{\text{Im}(f_i)}{2\text{Re}(f_i)}\).

A trace of CALC_MODES from the sdld313c benchmark test is shown below.

------------------------------------------------------------------------
MODAL CALCULATION: GLOBAL QR TYPE METHOD
ALGORITHM QZ_EQUI
NUMBER FREQUENCY (HZ) DAMPING    ERROR NORM
1        6.44568E+00 5.00000E-02 1.28280E-15
2        1.55613E+01 5.00000E-02 6.26512E-16
AVERAGE ERROR NORM: 0.95466E-15
.....

13.1.9.7. Summary of settings

Let us recapitulate now the parameter settings available for the operator CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’] with this QZ algorithm.

Table 13.1.8 Summary of the parameter settings of CALC_MODES with OPTION among [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘ALL’] (GEP) with SOLVEUR_MODAL = _F ( METHOD = ‘QZ’).

Keyword factor

Keyword

Default value

References

TYPE_RESU = ‘DYNAMIQUE’

‘DYNAMIQUE’

Section 13.1.3.1

‘MODE_FLAMB’

Section 13.1.3.1

OPTION = ‘PLUS_*’

‘PLUS_PETITE’

Section 13.1.9.4

‘BANDE’

Section 13.1.9.4

‘CENTRE’

Section 13.1.9.4

‘TOUT’

Section 13.1.9.4

SOLVEUR_MODAL

METHODE = ‘QZ’

‘SORENSEN’

Section 13.1.9.2

TYPE_QZ = ‘QZ_SIMPLE’

‘QZ_SIMPLE’

Section 13.1.9.4

‘QZ_EQUI’

Section 13.1.9.4

‘QZ_QR’

Section 13.1.9.4

CALC_FREQ or CALC_CHAR_CRIT

FREQ or CHAR_CRIT (if ‘CENTRE’)

Section 13.1.9.4

NMAX_FREQ or NMAX_CHAR_CRIT (if ‘PLUS_*’)

10

Section 13.1.9.4

AMOR_REDUIT (if ‘CENTRE’ and FREQ)

Section 13.1.5.4

NMAX_ITER_SHIFT

3

81 Section 13.1.3.2

PREC_SHIFT

0.05

81 Section 13.1.3.2

SEUIL_FREQ or SEUIL_CHAR_CRIT

1.E-02

Section 13.1.3.7

VERI_MODE

STOP_ERREUR = ‘OUI’

‘OUI’

Section 13.1.3.7

‘NON’

Section 13.1.3.7

PREC_SHIFT

5.E-03

Section 13.1.3.7

SEUIL

1.E-06

Section 13.1.3.7

STURM = ‘OUI’

‘OUI’

Section 13.1.3.7

‘NON’

Section 13.1.3.7

Note

  • The table contains all parameters related to the Sturm test (NMAX_ITER_SHIFT, PREC_SHIFT) and with verification checks (SEUIL_FREQ / SEUIL_CHAR_CRIT, VERI_MODE).

  • Beginners are encouraged to modify only the parameters in bold and to retain the rest of the empirically derived parameters values without modification.

  • The QZ algorithm is supposed to provide the “most reliable possible” values ​​of the modes sought. However, when it calculates the whole spectrum of a problem of non-negligible size (> 100 dofs), it can provide “approximate” values. Not being an algorithm with restart (cf. IRAM) or with a projection stage (cf. Lanczos / IRAM), the user has no digital “lever” to improve the convergence of modes. At most, one can change computation option (TYPE_QZ), in adding or removing the balancing stage (Section 13.1.9.4)

13.1.10. Parallelism and high performance computing

13.1.10.1. Integrated operators

Operators CALC_MODES and INFO_MODE can benefit from a first level of parallelism: that intrinsic to the linear solver MUMPS. However, the parallel efficiency of MUMPS in a modal calculation is more limited than for other types of analyses. In general, a parallel efficiency in time of the order of 0.2 to 0.3 is typical for a small number of processors: from 2 to 16. Beyond that we do not gain anything. This can be explained by:

  • the quasi-singularity of certain dynamic work matrices

  • a very unfavorable “number of descent-ascent / number of factorizations” ratio

  • the significant cost of the analysis phase compared to that of factorization,

  • a ratio “cost in time / memory of the linear solver” / “cost in time / memory of the modal solver” that is very unfavorable.

For more technical and functional information, one can consult documents [R6.02.03], [U4.50.01] and [U2.08.03 / 06].

In order to progress in terms of performance, we break down the initial computation into sub-computations that are more efficient and more precise. This is the object of one of the functionalities of CALC_MODES with OPTION = ‘BANDE’ and the list of \(n > 2\) values ​​given under CALC_FREQ = _F(FREQ), detailed in the next paragraph. Moreover, this algorithmic rewriting of the problem exhibits two levels of parallelism that are relevant and more efficient to “boost” the modal calculations of Code_Aster.

13.1.10.2. Operator CALC_MODES, option ‘BANDE’ divided into sub-intervals

To deal effectively with large modal problems (in mesh size and / or in the number of modes sought), we advises the use of CALC_MODES with the option ‘BANDE’ divided into sub-intervals. This breaks up the modal computation of a standard GEP (symmetric and real) into a succession of independent calculations which are less expensive, more robust, and more precise.

In sequential mode alone, the gains can be significant : factors 2 to 5 in time, 2 or 3 in peak RAM and 10 to 10,000 on the average error of the modes. In addition, with multi-level parallelism, by reserving about sixty processors, it can provide additional gains of the order of 20 in time and 2 in RAM peak (see table below). This is without loss of precision, no scope restriction and with the same numerical behavior.

Table 13.1.9 Benchmark perf016a results

Benchmark perf016a (N = 4M, 50 modes) division into 8 sub-intervals

Time elapsed

Memory peak (RAM)

1 processor

5524s

16.9GB

8 processors

1002s

19.5GB

32 processors

643s

13.4GB

division into 4 sub-intervals

1 processor

3569s

17.2GB

4 processors

1121s

19.5GB

16 processors

663s

12.9GB
Table 13.1.10 Sesimic study results

Seismic study (N = 0.7M, 450 modes) division into 20 sub-intervals

Time elapsed

Peak memory (RAM)

1 processor

5200s

10.5GB

20 processors

407s

12.1GB

80 processors

270s

9.4GB

division into 5 sub-intervals)

1 processor

4660s

8.2GB

5 processors

1097s

11.8GB

20 processors

925s

9.5GB
../../_images/fig16_r5.01.01.png

Fig. 13.1.15 Test cases for parallel CALC_MODES with default parameters(SOLVEUR = ‘MUMPS’ in IN_CORE and RENUM = ‘QAMD’). Code_Aster v11.3.11 on the IVANOE machine (1 or 2 MPI processes per node).

The principle of CALC_MODES U4.52.02, with the option ‘BANDE’ divided into sub-intervals, rests on the fact that the calculation and memory costs of the modal algorithms depend more than linearly on the number of modes sought. Therefore, as for the decomposition of fields [R6.01.03], one will decompose the search for hundreds of modes into more reasonably sized packages.

A packet of the order of forty modes seems to be an empirical optimum in sequential mode. In parallel, we can continue to improve performance down to 15. The example of Fig. 13.1.16 illustrates a total CALC_MODES computation in the interval \([\text{freq_{\text{min}}, \text{freq}_{\text{max}}]\) which is often advantageously replaced by ten CALC_MODES calculations targeted on contiguous sub-interval

\[[\text{freq}_1 = \text{freq}_{\text{min}}, \text{freq}_2], [\text{freq}_2, \text{freq}_3], \dots [\text{freq}_{10}, \text{freq}_{11} = \text{freq}_{\text{max}}]\]

On the other hand, this type of decomposition makes it possible to:

  • reduce robustness problems

  • to improve and homogenize the modal errors

In practice, this modal operator dedicated to HPC is broken down into four main stages:

  1. Modal precalibration (via INFO_MODE) of the sub-intervals parameterized by the user:

    Potentially, a loop of nb_freq independent calculations on each modal position of frequencies ([R5.01.04]).

  2. Effective modal calculation ( via CALC_MODES + OPTION = ‘BANDE’ + CALC_FREQ = _F (TABLE_FREQ)) of the modes contained in each non-empty sub-interval (by aggregating the modal calibrations of step 1):

    Potentially, a loop of nb_sbande_nonvide < nb_freq independent calculations.

  3. Post-verification with a Sturm test on the extreme limits of the calculated modes (via INFO_MODE):

    Potentially, a loop of 2 independent calculations.

  4. Post-processing of all modes obtained: Normalization (via NORM_MODE) and filtering (via EXTR_MODE).

../../_images/fig17_r5.01.01.png

Fig. 13.1.16 Principle of the decomposition of calculations of CALC_MODES with the option ‘BANDE’ divided into sub-intervals.

In parallel simulations, each of the stages of calculation can use at least one level of parallelism:

  • The first two by distributing the calculations of each sub-interval over the same number of processors.

  • The third, by distributing the modal positions of the limits of the interval of verification on two processor packages.

  • The fourth step, inexpensive, remains sequential.

If the number of processors and the parameter setting allow it (in particular, if one uses the linear solver MUMPS), one can exploit a second level of parallelism.

Fig. 13.1.17 illustrates a modal computation seeking to take advantage of 40 processors by breaking down the initial computation into ten search sub-intervals. Each benefits from the support of 4 MUMPS occurrences for the inversion of linear systems required by modal solvers. For an exhaustive presentation of this multi-level parallelism, its challenges and some technical and functional details, one can consult documents [R5.01.04], [U4.52.01], U4.52.02 and [U2.08.06].

../../_images/fig18_r5.01.01.png

Fig. 13.1.17 Example of two levels of parallelism in the INFO_MODE of preprocessing and in the loop of the sub-intervals of CALC_MODES. Distribution on nb_proc = 40 processors with a division into 10 sub-intervals (parallelism known as “10x4”). One uses here the linear solver MUMPS and the parameter setting of parallelism by default (‘COMPLET’).

Note

  • In MPI parallelism, the main concerns are the distribution of tasks and their communication. For CALC_MODES with the option ‘BANDE’ (divided in sub-intervals), task distribution is carried out in the Python as well as in FORTRAN. The two communicate by internal keywords: PARALLELISME_MACRO. But all MPI calls are limited to only the F77 / F90 layers.

  • The global communication of the eigenvalues ​​and vectors is carried out at the end of the modal computation on each sub-interval at an intermediate level, between the simple communication of linear algebra results (e.g., MUMPS / PETSc) and the communication of Aster data structures in Python after filtering (optimal in term of performance but much more complicated to implement).

  • The ideal would have been to be able to empirically balance the frequency sub-intervals to limit load unbalances related to the distribution of modes in sub-intervals and those related to modal computation. This approach would make it possible to benefit from the decomposition of the macro, even on few processors. This feature has not been validated because of the complexity involved in manipulating Code_Aster concepts.

13.1.11. Bibliography

13.1.11.1. Books / articles / proceedings / theses

35(1,2)

P.Arbenz, ULHetmaniuk, RBLehoucq & RSTuminaro. A comparison of eigensolvers for large-scale 3d modal analysis using AMG-Preconditioned Iterative Methods . Int. J. Numer. Meth. Engng, vol. 64, pp204-236 (2005).

36(1,2)

WEArnoldi. The principle of minimized iterations in the solution of the matrix eigenvalue problem . Quarter. Appl. Math, flight9, pp17-19 (1951).

37(1,2)

KJBathe Solution methods for large generalized eigenvalue problems in structural engineering. Ed. California universityBerkeley (1971).

38

Bathe, K. J., & Wilson, E. L. (1972). Large eigenvalue problems in dynamic analysis. Journal of the Engineering Mechanics Division, 98(6), 1471-1485.

39

O.Bertrand. Eigenvalue enumeration procedures . Thesis of the University of Rennes 1 (1999).

40

Z.Bai, J. Demmel & AMKenney. On computing condition numbers for the nonsymmetric eigenproblem . ACM transactions on mathematical software, vol 19, pp202-223 (1993).

41

Bergamaschi & M.Putti. Numerical comparison of iterative eigensolvers for large symmetric positive definite matrices .Comput. Methods App. Mech. Engrg. 191 5233-5247 (2002).

42(1,2)

F.Chatelin. Eigenvalues of matrices . Ed. Masson (1988).

43(1,2,3)

JKCullum & RAWilloughby. Lanczos algorithms for large symetric eigenvalue computations. Vol 1 Theory, EdBirkhäuser (1985).

44

ER,Davidson. “The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices.” Journal of Computational Physics 17 (1975): 87-94.

45

A.Dubrulle, RSMartin & JHWilkinson. The Implicit QL Algorithm . Handbook for automatic computation. Vol 2, Linearalgebra, Springer-Verlag (1971).

46

T.Ericsson & A.Ruhe. The spectral transformation Lanczos method for the numerical solution of large sparse generalised symmetric eigenvalue problems. Mathematics of computations, vol 35, pp1251-1268 (1980).

47

Francis, J. G. (1961). The QR transformation a unitary analogue to the LR transformation—Part 1. The Computer Journal, 4(3), 265-271.

48(1,2,3,4,5)

GHGolub & CFVan Loan. Matrix computations . The Johns Hopkins university press (1989).

49

G.Gambolati, F.Sartoretto & P.Florian. An orthogonal accelerated deflation technique for large symmetric eigenproblems .Comp. Meth. Appl. Mech. Engrg., 94 13-23 (1992).

50

Y.Haugazeau. Application of Sylsvester’s theorem to the localization of eigenvalues . RAIRO digital analysis,flight. 14, 1, pp25-41 (1980).

51(1,2,3,4)

V.Hernandez, JERoman, A.Tomas & V. Vidal. A survey of software for sparse eigenvalue problems . SLEPc report available at http://www.grycap.upv.es/slepc.

52

JFImbert. Structural analysis by finite elements. Ed. CEPADUES (1991).

53

AVKnyazev. Toward the optimal preconditioned eigensolver: locally optimal block preconditioned conjugate gradient method. SIAM J. Sci. Comput., 23 pp517-541 (2001).

54

KUBLANOVSKAYA, V. N. (1961). On some algorithms for the solution of the compute eigenvalue problem. Zh. Vych, Math., 1, 555-570.

55(1,2)

C.Lanczos. An iteration method for the solution of the eigenvalue problem of linear differential and integral operators .Newspaper. of research. of the national bureau of standards, vol 45, no 4 pp255-282 (1950).

56

RBLehoucq & JAScott. An evaluation of software for computing eigenvalues of sparse nonsymmetric matrices. Argonne National Laboratory Report MCS-P547-1195.

57(1,2,3,4,5,6,7)

P.Lascaux & R.Theodor. Matrix numerical analysis applied to the Art of the engineer . Masson (1986).

58

RSMartin, G.Peters & JHWilkinson. The QR Algorithm for Real Hessenberg Matrices. Handbook for automatic computation ,. Flight. 2, Linear algebra, Springer-Verlag (1971).

59(1,2)

CBMoler & GWStewart. An algorithm for GEP. SIAM J. Num. Anal., 10, 241-256 (1973).

60

RBMorgan & DSScott. Preconditionning the Lanczos algorithm for sparse symmetric eigenvalue problems . SIAM J.Sci. Comp, vol 14 (1993).

61(1,2)

M.Newmann & A.Pipano. Fast modal extraction in NASTRAN via FEER computer programs, NASTRAN, user manual, NASA Langley Research Center pp485-506 (1977).

62

EEOsborne. On pre-conditioning of matrices . J.Assoc.Comput.Mach., Vol 7, pp338-345 (1960).

63

CCPaige. Accuracy and effectivness of the Lanczos algorithm for the symmetric eigenproblem Linear algebra and its applications, vol 34 (1980).

64(1,2,3,4,5,6,7,8,9)

BNParlett. The symetric eigenvalue problem. Prentice Hall (1980).

65

BNParlett & C. Reinsch. Balancing a matrix for calculation of eigenvalues and eigenvectors. Handbook for automatic computation. Flight. 2, Linear algebra, Springer-Verlag (1971).

66

Rutishauser, H. (1958). Solution of eigenvalue problems with the LR-transformation. National Bureau of Standards, Applied Mathematics Series, 49, 47-81.

67(1,2,3,4)

Y.Saad. On the rates of convergences of the Lanczos and the Block-Lanczos methods. SIAM J. Numer. Anal., Vol 17 No 5, pp687-706 (1980). Y.Saad. Variation on Arnoldi’s method for computing eigenelements of large unsymetric matrices . Linear algebra and its applications, vol 34, pp269-2395 (1980).

68

Y.Saad. Numerical methods for large eigenvalue problems. Ed. Manchester university (1991).

69(1,2)

Mr. Sadkane. A block Arnoldi-Tchebychev method for computing the leading ergenpairs of large sparse unsymmetric matrices. Numerische mathematik, vol 64, pp181-193 (1993).

70(1,2)

DCSorensen. Implicit applications of polynomial filters in a k-step Arnoldi method. SIAM J. Matrix Anal. Appl., Vol 13, pp357-385 (1992).

71(1,2)

GLGSleijpen & H. Van der Vorst. A Jacobi-Davidson method for linear eigenvalue problems . SIAM J. Matrix Anal. Appl.,17 pp401-425 (1996).

72

F.Tisseur & K. Meerbergen. The quadratic eigenvalue problem. SIAM Review vol 43, no 2, pp235-286 (2001).

73

DSWatkins. The matrix eigenvalue problem . SIAM (2007).

74

JHWilkinson & C. Reinsch. Handbook for automatic computation , vol 2, Linear Algebra, Springer-Verlag (1971).

13.1.11.2. EDF reports / reports

75

Mr Begon. Establishment of modal solvers in Code_Aster . INSTN Master’s Report (2006).

76(1,2)

O.Boiteau & B. Quinnez. Method of spectral analysis in Code_Aster . Dynamics Course Booklet availableon the website http://www.code-aster.org (2000).

77(1,2)

O.Boxeau. General on the direct linear solvers and use of the product MUMPS . Reference DocumentationCode_Aster R6.02.03 (2008).

78

O.Boiteau. Integration of a modal solver of the type QZ in Code_Aster and Extension of the solvers QZ and Sorensen to Unsymmetrical modal problems of Code_Aster. EDF R&D internal reports CR-I23 / 2008/030 and CR-I23 / 2008/044 (2008).

79(1,2,3,4,5)

O.Boxeau. Resolution of the quadratic modal problem (QEP) . Reference Documentation Code_Aster R5.01.02(2009).

80

O.Boiteau. General on the conjugate gradient: GCPC Aster and use of the product PETSc . Documentation ofReference Code_Aster R6.01.02 (2009).

81(1,2,3,4,5,6,7,8,9,10)

O.Boxeau. Procedure for enumeration of eigenvalues . Reference Documentation Code_Aster R5.01.04 (2012).

82

E.Boyère. Operator DYNA_TRAN_MODAL . Documentation of Use Code_Aster U4.53.21 (2007).

83

J.Pellet. Dualization of boundary conditions . Reference Documentation Code_Aster R3.03.01 (2001).

84

C.Rose. Multifrontal method . Reference Documentation Code_Aster R6.02.02 (2007).

85(1,2,3,4,5,6,7)

JLVaudescal. Introduction to the methods of solving large eigenvalue problems. EDF rating R&D HI-72/00/01 (2000).

13.1.11.3. Internet resources

86(1,2,3)

ARPACK website: http://www.caam.rice.edu/software/ARPACK/

87

LAPACK website: http://www.netlib.org/lapack/.

88

MatrixMarket website: http://math.nist.gov/MatrixMarket/index.html.

13.1.12. Description of the versions of the document

Please refer to the original PDF documentation for updated information on authors and version.

13.1.13. Appendix 1. General information on the QR algorithm

13.1.13.1. Principle

QR-type algorithms were introduced by H. Rutishauser (1958) and concurrently formalized by JCFrancis and VNKublanovskaya (1961). This fundamental method is often used within other approaches better suited to dealing with large-scale problems (particularly projection methods.

Algorithm 1.1. Theoretical QR.

  • For \(k = 1, \dots\) do

    • \(\mathbf{H}_k = \mathbf{Q}_k \mathbf{R}_k\) (QR factorisation)

    • \(\mathbf{H}_{k+1} = \mathbf{R}_k \mathbf{Q}_k\)

  • End For

The process leads iteratively to an upper triangular matrix \(\mathbf{H}_k\) (or block triangular) whose diagonal terms are the eigenvalues ​​of the initial operator \(\mathbf{H} = \mathbf{H}_1\). The \(\mathbf{H}\) notation is not arbitrary, because we are interested in orthogonal transformation (so as to modify only the form of the shift operator and not its spectrum) of the work operator \(\mathbf{A}\) into an upper Hessenberg form, or in real arithmetic

\[\mathbf{H}_1 = \mathbf{Q}_0^T \mathbf{A} \mathbf{Q}_0`\]

This can be done via various orthogonal transformations (Householder, Givens, Gram-Schmidt etc.) and their cost (of the order of \(O(10n^3/ 3)\)) is negligible compared to the gain that they allow to achieve at each iteration of the global process: \(O(n^2)\) (with Householder or Fast-Givens) against \(O(n^3)\). This gain of an order of magnitude can even be improved when the operator is symmetric tridiagonal (this is the case of Lanczos with a real scalar product): \(O(20n)\).

Convergence towards a simple triangular matrix is ​​carried out only if all the eigenvalues ​​are of distinct magnitude and that if the initial matrix is ​​not “pathologically” too poor in the eigen-directions. The convergence of the i-th mode (classically arranged in decreasing order of magnitude) is then checked with

\[\underset{j \ne i}{\text{max}} \frac{|\lambda_j|}{|\lambda_i|}\]

which can be very slow if no additional process is implemented.

Note

  • The determination of the spectrum of the Rayleigh matrix with the variant of Newmann & Pipano is performed via a simple QR (or a symmetric QL ) of this type (with a preliminary balancing procedure). The only parameter accessible by the user is the maximum number of iterations NMAX_ITER_QR.

  • For IRAM, this calculation is carried out via a QR method with explicit double shift while the polynomial filters handling restarts use a QR with implicit double shift. As a precaution, no parameter is accessible for the user in Code_Aster!

  • We must not confuse the method, the QR algorithm, and one of its conceptual tools, QR factorization.

  • This class of algorithm is widely used to determine the full spectrum of an operator, because it is very robust (it is the reference in this field). However it is very greedy in memory which makes its use prohibitive for large systems.

  • The scope of application of the QR algorithm is much more general than that of Jacobi (which is the second standard algorithm providing the full spectrum of an operator so that we need to store it entirely) which is limited to Hermitian matrices. To accelerate the convergence of the simple algorithm which can be very slow (in the presence of eigenvalue clusters, for example) a multitude of variants, based on the choice of shifts meeting certain criteria, have emerged.

13.1.13.2. The shift strategy

This strategy consists of artificially creating a deflation phenomenon (Section 13.1.4.1, Section 13.1.6.4 ) in the work matrix. This offers the following advantages:

  • to be able to isolate a real eigenvalue or even two conjugate complex eigenvalues

  • while reducing the size of the problem to be treated

  • and accelerating convergence.

In the version with simple explicit shift, the method is then rewritten in the following form:

Algorithm 1.2. QR with simple explicit shift.

  • For \(k = 1, \dots\) do

    • Choose the shift \(\mu\)

    • \(\mathbf{S}_k = \mathbf{H}_k - \mu \mathbf{I}\)

    • \(\mathbf{S}_k = \mathbf{Q}_k\mathbf{R}_k\) (QR factorisation)

    • \(\mathbf{H}_{k+1} = \mathbf{R}_k\mathbf{Q}_k+\mu \mathbf{I}\)

  • End For

Note

  • This process generalizes to several shifts. In that case we build, for each iteration \(k\), as many auxiliary matrices \(\mathbf{S}_k\) as the number of shifts \(\mu_i\). The convergence of the process is improved in the sense that the subdiagonal terms cancel each other asymptotically:

    \[\underset{j \ne i}{\text{max}} \frac{|\lambda_j-\mu|}{|\lambda_i-\mu|}\]

    In theory, if the shift \(\mu\) is an eigenvalue of the problem, then the deflation is exact. In practice, rounding effects perturb this phenomenon. This is called the direct instability property of the algorithm. The main difficulty lies in the choice of the shift (or shifts). On the other hand, we do not keep the same shift for all iterations. One must change it when it is associated with a converged eigenvalue. Indeed, we will have numerically caused its crowding out of the work spectrum by inducing deflation in the previous iteration.

    Since the sixties a whole zoo of shifts has developed. While the simplest uses the last diagonal term \(\mathbf{H}_{n,n}\), that of JH Wilkinson 74 analytically determines the eigenvalue \(\mu\) of the block diagonal matrix

    \[\begin{split} \begin{bmatrix} \mathbf{H}_{n-1, n-1} & \mathbf{H}_{n-1, n} \\ \mathbf{H}_{n, n-1} & \mathbf{H}_{n,n} \end{bmatrix}\end{split}\]

    closest to this term. This technique makes it possible to obtain a quadratic or even cubic convergence (in the symmetrical case) and it is particularly effective for capturing double modes or distinct eigenvalues ​​of the same magnitude. However, this strategy may prove ineffective in the presence of conjugated complex eigen modes.

    The same author then proposed a variant including a double shift corresponding to the two complex eigenvalues \(\mu_1\) and \(\mu_2\) (determined beforehand) of the offending block. But besides the computational issues relating to the appearance of invasive complex components (which in theory do not take place), it is not easy to maintain the “upper Hessenberg” characteristic of the work matrix. At best, this slows down the convergence of the QR algorithm. At worst, it completely distorts its operation.

    In order to remedy these drawbacks, JCFrancis has developed a version with an implicit double shift. It is explained in 48 (pp377-81) and we will content ourselves here with summarizing it.

    To minimize the effects of rounding, it would be necessary to directly constitute the auxiliary matrix \(\mathbf{S}_k\) resulting from simultaneous application of two shifts \(\mu_1\) and \(\mu_2\):

    \[\mathbf{S}_k = \mathbf{H}_k^2- (\mu_1+\mu_2) \mathbf{H}_k+ (\mu1\mu2) \mathbf{I}\]

    before factoring it in QR form and building the new iterated \(\mathbf{H}_{k+1}\):

    \[\begin{split}\mathbf{S}_k & = \mathbf{Q}_k \mathbf{R}_k \\ \mathbf{H}_{k+1} &= \mathbf{Q}_k^T \mathbf{H}_k \mathbf{Q}_k\end{split}\]

    But the only cost of the initial assembly (of the order of \(O(n^3)\)) makes the tactic ineffective. This variant, based on the \(\mathbf{Q}\) -implicit theorem, consists in applying Householder transformations to the work matrix and this allowing us to find a Hessenberg matrix that is “essentially” equal to \(\mathbf{H}_{k+1}\), i.e. of the type:

    \[\tilde{\mathbf{H}}_{k+1} = \mathbf{E}^{-1} \mathbf{H}_{k+1} \mathbf{E} \quad \text{where} \quad \mathbf{E} = \text{diag}(\pm 1, \dots, \pm 1)\]

    The working matrix thus preserves, at the lowest cost, its particular structure and its spectrum.

    All these variants are in fact very sensitive to the balancing techniques implemented to precondition the initial operator. The following paragraph will summarize this technique of “balancing” terms of the work matrix which is widely used in modal computation.

13.1.13.3. Balancing

Balancing (also called equilibration) is used to mitigate rounding errors by restricting the scope of expansion of the terms of the work operator, i.e., by preventing them from becoming too small or too large. EEOsborne 62 (1960) noticed that the error in the calculation of the eigen decomposition of \(\mathbf{A}\) is of the order of \(\varepsilon\| \mathbf{A} \|_2\) where \(\varepsilon\) is the machine precision. He then proposed to transform the initial matrix into a matrix \(\tilde{\mathbf{A}} = \mathbf{D}^{-1} \mathbf{A}\mathbf{D}\) (where \(\mathbf{D}\) is a diagonal matrix), such that \(\|\tilde{\mathbf{A}}\|_2 \ll\| \mathbf{A} \|_2\)

In fact, we iteratively calculate a succession of matrices

\[\mathbf{A}_k = \mathbf{D}_{k-1}^{-1} \mathbf{A}_{k-1} \mathbf{D}_{k-1}\]

with

\[\mathbf{A}_k \underset{k \rightarrow \infty}{\rightarrow} \mathbf{A}_f\]

satisfying \(\| \mathbf{a}^i\|_2 = \| \mathbf{a}_i \|_2\) with \(\mathbf{a}^i\) and \(\mathbf{a}_i\), respectively, the i-th column and row of \(\mathbf{A}_f\).

Note

  • This technique has been generalized to any matrix norm induced by a discrete norm 89 \(l^q\).

  • Its use is widespread in scientific computation and in particular among the direct solvers of linear systems of equations.

The computer arithmetic base, denoted \(\beta\), intervenes in the determination of the terms of the matrices \(\mathbf{D}_k\). In order to minimize rounding errors, we choose the elements of \(\mathbf{D}_k\) so that they are powers of this base. Let \(R_k\) and \(C_k\) be the \(p\) norms (in practice we often take \(p = 2\)), respectively of the row and column \(i\) of the matrix \(\mathbf{A}_{k-1}\) (\(i\) is chosen such that \(i -1 \equiv k -1 [ n ]\)). Assuming that \(R_k C_k \ne 0\), we can show that there exists a unique signed integer \(\sigma\) such that

\[\beta^{2\sigma - 1} < \frac{R_k}{C_k} \le \beta^{2\sigma+1}\]

Let \(f = \beta^\sigma\). We define the matrix \(\tilde{\mathbf{D}}_k\) such that

\[\begin{split}\tilde{\mathbf{D}}_k = \begin{cases} \mathbf{I} + (f-1) \mathbf{e}_i \mathbf{e}_i^T & \quad \text{if} \quad (C_k f)^p + \left(\frac{R_k}{f}\right)^p < \gamma (C_k^p + R_k^p) \\ \mathbf{I} & \quad \text{otherwise} \end{cases}\end{split}\]

where \(0 < \gamma \le 1\) is a constant and \(\mathbf{e}_i\) is the i-th canonical vector. We then build the \(\mathbf{A}, \mathbf{D}\) matrices iteratively:

\[\begin{split}\mathbf{D}_k &= \tilde{\mathbf{D}}_k \mathbf{D}_{k-1} \\ \mathbf{A} k &= \tilde{\mathbf{D}}_{k-1} \mathbf{A}_{k-1} \tilde{\mathbf{D}}_k\end{split}\]

with the initialization \(\mathbf{D}_0 = \mathbf{I}\). The process stops when \(\tilde{\mathbf{D}}_k = \mathbf{I}\).

Note

  • Before carrying out the balancing, a search for the isolated eigenvalues ​​is carried out by detecting the presence of almost zero rows and columns (all the terms are zero, except the one placed on the diagonal). When there is one we can, by carrying out suitable permutations, put the work matrix in the following block form:

    \[\begin{split}\left[ \begin{array}{@{}c|c@{}|@{}c} \begin{matrix} \star & & \star \\ & \star & \\ 0 & & \star \end{matrix} & \mathbf{Y} & \mathbf{Z} \\ \hline \mathbf{0} & \mathbf{X} & \mathbf{T} \\ \hline \mathbf{0} & \mathbf{0} & \begin{matrix} \star & & \star \\ & \star & \\ 0 & & \star \end{matrix} \end{array} \right]\end{split}\]

  • It is then necessary to carry out the balancing only on the central block \(\mathbf{X}\) because the diagonal terms of the two upper triangular matrices are eigenvalues ​​of the working matrix.

  • In Code_Aster, we chose \(p = 1\) and \(\gamma = 0.95\) 65.

  • Before applying the QR method, we begin by balancing the work matrix and then we transform into upper Hessenberg form. Once the QR calculation has been carried out, we convert the Ritz vectors to the eigenvectors via the inverse of the orthogonal transformations. It is this process that is implemented both in Lanczos and in IRAM. However, beside the fact that it requires the storage of these orthogonal matrices, it is also extremely sensitive to the effects of rounding. So it would be greatly preferable to estimate the eigenvectors via an inverse powers method initialized by the extracted eigenvalues.

89

The function space \(l^q\) is the set of complex sequences \((u_n)_n\) such that \(\sum_n | u_n |^q <\infty\).

13.1.13.4. The QL method

The QL factorization of a matrix \(\mathbf{A}\) consists of orthonormalizing its column vectors starting with the last (unlike the QR algorithm which starts with the first) thus giving a matrix \(\mathbf{L}\) that is lower triangular and regular. The QL algorithm applied to an invertible operator \(\mathbf{A}\) is equivalent to the QR algorithm applied to \(\bar{\mathbf{A}}^\star\) (conjugate transpose matrix of the inverse of \(\mathbf{A}\)). The method implemented in Code_Aster is identical to the simple shift method of JHWilkinson discussed previously, by adapting the search for this shift to the highest 2 \(\times\) 2 minor.

Note

Unlike QR, which captures the eigenvalues ​​in increasing order of magnitude, QL produces the dominant modes, and then the others in decreasing order of magnitude.

13.1.14. Appendix 2. Gram-Schmidt orthogonalization

13.1.14.1. Introduction

We have seen that the quality of orthogonalization of a family of vectors is crucial for the good progress of the algorithms and the quality of the modes obtained. This task is moreover to be carried out constantly, hence the importance of a fast algorithm.

The simple orthonormalization algorithms are deduced from the classical Gram-Schmidt process but they are often “conditionally stable” (the quality of their work depends on the conditioning of the matrix family to be orthonormalized). This lack of robustness can prove to be problematic in dealing with situations with poorly conditioned matrices. We then prefer more expensive but robust algorithms, based on projections or rotations: the spectral transformations of Householder and Givens.

In practice, for the implementation of the IRA method in Code_Aster, we have retained an iterative version of the Gram-Schmidt process (the IGSM of Kahan-Parlett 64). It achieves a good compromise between robustness and computational complexity since it is unconditionally stable and makes it possible to orthogonalize to machine precision, and take at worst up to twice as long to compute as a conventional Gram-Schmidt (GS).

In the following paragraphs we will detail the operation of the basic algorithm, as well as that of these two main variants. The first is set up in CALC_MODES with OPTION = [‘PROCHE’, ‘AJUSTE’, ‘SEPARE’] and the second is used in CALC_MODES with OPTION = [‘BANDE’, ‘CENTRE’, ‘PLUS_*’, ‘TOUT’].

A comparison of these methods is provided in 85 (pp33-36) using a simple example.

13.1.14.2. Gram-Schmidt algorithm (GS)

Given \(k\) independent vectors of \(\mathbb{R}^n\), \((\mathbf{x}_i)_{i = 1,k}\) we want to get \(k\) orthonormal vectors (with respect to any scalar product) \((\mathbf{y}_i)_{i=1,k}\) of the space they generate. In other words, we want to obtain another orthonormal family generating the same space. The process of classical Gram-Schmidt orthonormalization is as follows:

Algorithm 2.1. Gram-Schmidt algorithm (GS).

  • For \(i = 1, k\) do

    • For \(j = 1, i-1\) do

      • Calculate \(r_{ji} = (\mathbf{y}_j, \mathbf{x}_i)\)

    • End For \(j\)

    • Calculate \(\hat{\mathbf{y}}_i = \mathbf{x}_i - \sum_{j=1,i-1} r_{ji} \mathbf{y}_j\)

    • Calculate \(r_{ii} = \|\hat{\mathbf{y}}_i \|\)

    • Calculate \(\mathbf{y}_i = \frac{\hat{\mathbf{y}}_i}{r_{ii}}\)

  • End For

This process is simple but very unstable due to rounding errors, which has the effect of producing non-orthogonal vectors. In particular, when the initial vectors are almost dependent this creates large magnitude deviations in the second step of the process

\[\hat{\mathbf{y}}_i = \mathbf{x}_i - \sum_{j=1,i-1} r_{ji} \mathbf{y}_j\]

From a computational point of view, managing these discrepancies is very difficult.

Note

If \(\mathbf{Q}\) is the matrix generated by the \((\mathbf{y}_i)_{i=1,k}\), we have explicitly constructed a QR factorization of the initial matrix \(\mathbf{X}\) related to \((\mathbf{x}_i)_{i=1,k}\). This is in fact the goal of any orthogonalization process.

13.1.14.3. Modified Gram-Schmidt Algorithm (GSM)

In order to remove these numerical instabilities, we reorganizes the preceding algorithm. Mathematically equivalent to the previous process, it is then much more robust because it avoids deviations in magnitude between the vectors handled in the algorithm.

In the initial process, the orthogonal vectors \(\mathbf{y}_i\) are obtained without taking into account the \(i-1\) previous orthogonalizations. With the Modified Gram-Schmidt, we orthogonalize more progressively by taking into account the previous alterations following the process below.

Algorithm 2.2. Modified Gram-Schmidt Algorithm (GSM).

  • For \(i = 1, k\) do

    • \(\hat{\mathbf{y}}_i = \mathbf{x}_i\)

    • For \(j = 1, i-1\) do

      • Calculate \(r_{ji} = (\mathbf{y}_j, \mathbf{x}_i)\)

      • Calculate \(\hat{\mathbf{y}}_i = \mathbf{x}_i - r_{ji} \mathbf{y}_j\)

    • End For \(j\)

    • Calculate \(r_{ii} = \|\hat{\mathbf{y}}_i \|\)

    • Calculate \(\mathbf{y}_i = \frac{\hat{\mathbf{y}}_i}{r_{ii}}\)

  • End For \(i\)

The orthonormality of the basis vectors is much better with this method and it can even be obtained at the machine precision and within a constant (depending on the conditioning of \(\mathbf{X}\)). However, for particularly poorly conditioned matrices, this “conditional” stability can quickly prove to be problematic, hence the use of the following iterative algorithm.

Note

GSM is twice as efficient as a Householder method to obtain a QR factorization of the initial matrix \(\mathbf{X}\). It only requires \(O(2nk^2)\) operations (with \(n\) the number of rows of the matrix).

13.1.14.4. Iterative Gram-Schmidt Algorithm (IGSM)

To ensure the orthogonality with close to machine precision, a second orthogonalization is performed. And if orthogonality is not achieved after the second step, it is no necessary to start again as the quantities are then certainly very close and their differences oscillate around zero. This approach is based on a theoretical analysis by W. Kahan and expanded by BNParlett 64 (cf. pp105-110).

Algorithm 2.3. Kahan-Parlett Iterative Gram-Schmidt Algorithm (IGSM).

  • For \(i = 1, k\) do

    • \(\hat{\mathbf{y}}_i = \mathbf{x}_i\)

    • For \(j = 1, i-1\) do

      • Calculate \(r_{ji} = (\mathbf{y}_j, \hat{\mathbf{y}}_i)\)

      • Calculate \(\tilde{\mathbf{y}}_i = \hat{\mathbf{y}}_i - r_{ji} \mathbf{y}_j\)

      • If \(\|\tilde{\mathbf{y}}_i\| \ge \alpha \|\hat{\mathbf{y}}_i\|\) then

        • \(\hat{\mathbf{y}}_i \leftarrow \tilde{\mathbf{y}}_i\)

        • Exit loop in \(j\)

      • Else

        • Calculate \(\tilde{r}_{ji} = (\mathbf{y}_j, \tilde{\mathbf{y}}_i)\)

        • Calculate \(\tilde{\tilde{\mathbf{y}}}_i = \tilde{\mathbf{y}}_i - \tilde{r}_{ji} \mathbf{y}_j\)

        • If \(\|\tilde{\tilde{\mathbf{y}}}_i\| \ge \alpha \|\tilde{\mathbf{y}}_i\|\) then

          • \(\hat{\mathbf{y}}_i \leftarrow \tilde{\tilde{\mathbf{y}}}_i\)

          • Exit loop in \(j\)

        • Else

          • \(\hat{\mathbf{y}}_i \leftarrow \mathbf{0}\)

          • Go to the next \(i\)

        • End if

      • End if

    • End For \(j\)

    • Calculate \(r_{ii} = \|\hat{\mathbf{y}}_i\|\)

    • Calculate \(\mathbf{y}_i = \frac{\hat{\mathbf{y}}}{r_{ii}}\)

  • End For \(i\)

During the use of IRAM in Code_Aster, the criterion \(\alpha\) is parameterized by the keyword PARA_ORTHO_SOREN (see Section 13.1.7.5). We show that its validity interval is \([1.2\varepsilon, 0.83 - \varepsilon]\) with \(\varepsilon\) being the machine precision and it is assigned the value 0.717 (by default).

Note

  • The greater the value of parameter \(\alpha\), the less reorthogonalization is triggered, but that affects the quality of the solution.

  • Unlike the “in-house” version developed for Lanczos, the stopping criteria are concentrated on the norms of orthogonalized vectors rather than on vector dot products which have more tendency to reflect the effects of rounding. This, together with the removal of iterations greater than two, may explain the increased efficiency of the Kahan-Parlett version.

  • According to DCSorensen, this method is to be attributed to J. Daniel (paper of 1976 submitted to Mathematics of Computation, vol.30, pp772-795).

13.1.15. Appendix 3. Jacobi method

13.1.15.1. Principle

The Jacobi method 57 makes it possible to calculate all the eigenvalues ​​of a generalized problem whose matrices are positive definite and symmetric (the matrices obtained with each iteration by the method of Bathe & Wilson satisfy these properties; Section 13.1.8). The matrices \(\mathbf{A}\) and \(\mathbf{B}\) of the GEP \(\mathbf{A}\mathbf{u} = \lambda \mathbf{B}\mathbf{u}\) are transformed into diagonal matrices, using similar successively similar orthogonal transformations (Givens rotation matrices). The process can be schematized as follows:

Algorithm 3.1. Jacobi process

\[\begin{split}\mathbf{A}^1 & = \mathbf{A} \quad \mathbf{B}^1 &= \mathbf{B} \\ \mathbf{A}^2 & = \mathbf{Q}_1^T \mathbf{A}^1 \mathbf{Q}_1 \quad \mathbf{B}^2 & = \mathbf{Q}_1^T\mathbf{B}^1\mathbf{Q}_1 \\ \dots & &\\ \mathbf{A}^k & = \mathbf{Q}_{k-1}^T\mathbf{A}^{k-1}\mathbf{Q}_{k-1} \quad \mathbf{B}^k & = \mathbf{Q}_{k-1}^T\mathbf{B}^{k-1}\mathbf{Q}_{k-1} \\ \mathbf{A}^k & \underset{k \rightarrow \infty}{---\rightarrow}\underbrace{\mathbf{A}^d}_{\text{diagonal matrix}} \quad \mathbf{B}^k & \underset{k \rightarrow \infty}{---\rightarrow}\underbrace{\mathbf{B}^d}_{\text{diagonal matrix}}\end{split}\]

The eigenvalues ​​are given by \(\lambda = \mathbf{A}^d(\mathbf{B}^d)^{-1}\) or \(\lambda = \frac{\mathbf{A}_{ii}^d}{\mathbf{B}_{ii}^d}\) and the matrix of the eigenvectors is

\[\begin{split}\mathbf{X} = \mathbf{Q}^1\mathbf{Q}^2\dots\mathbf{Q}^k\dots \begin{bmatrix} \frac{1}{\sqrt{\mathbf{A}_{11}^d}} & & & \\ & \frac{1}{\sqrt{\mathbf{A}_{22}^d}} & & \\ & \dots & \dots & \\ & & & \frac{1}{\sqrt{\mathbf{A}_{nn}^d}} \end{bmatrix}\end{split}\]

Each matrix \(\mathbf{Q}^k\) is chosen so that a diagonal and non-zero term \((i,j)\) of \(\mathbf{A}^k\) or of \(\mathbf{B}^k\) is zero after the transformation.

13.1.15.2. Some choices

In this algorithm, we realize that the important points are the following:

  • How to choose the terms to cancel?

  • How to measure the diagonal character of matrices when \(k \rightarrow \infty\)?

  • How to measure convergence?

13.1.15.2.1. Non-diagonal terms to be canceled

Several methods can be used to choose terms to cancel:

  • At each step \(k\) choose the largest element in magnitude not on the diagonal of matrix \(\mathbf{A}_k\) or \(\mathbf{B}_k\) and cancel it by rotating. This choice ensures the convergence of the Jacobi method but is relatively expensive (due to the search for the maximum).

  • Successively cancel all the non-diagonal elements of these matrices following the natural order \(\mathbf{A}_{13}^k, \dots, \mathbf{A}_{1n}^k, \mathbf{A}_{23}^k, \dots\). When we arrive at \(\mathbf{A}_{n-1,n}^k\), we start the cycle again. This method converges slowly.

A variant of this method, while following the natural order of the terms, cancels only those which are greater than a given tolerance \(\varepsilon_k\). At the end of a cycle, one decreases the value of this criterion and starts again. It is this strategy which is used in Code_Aster.

13.1.15.2.2. Convergence test

To test the convergence and the diagonal character of the matrices, we check that all coupling factors defined by:

\[f_{\mathbf{A}_{ij}} = \frac{|\mathbf{A}_{ij}|} {\sqrt{|\mathbf{A}_{ii}\mathbf{A}_{jj}|}} \quad f_{\mathbf{B}_{ij}} = \frac{|\mathbf{B}_{ij}|} {\sqrt{|\mathbf{B}_{ii}\mathbf{B}_{jj}|}} \quad i \ne j\]

are less than a given precision (diagonal character of the matrices). We also control the convergence of eigenvalues ​​via the indicator

\[f_\lambda = \underset{i}{\text{max}} \frac{|\lambda_i^k-\lambda_i^{k-1}|}{\lambda_i^{k-1}}\]

so that it remains lower than a given precision \(\varepsilon_{\text{jaco}}\).

13.1.15.2.3. Algorithm implemented in Code_Aster

The algorithm implemented in Code_Aster is described below:

Algorithm 3.2. Jacobi method implemented in Code_Aster.

  • Initialize the matrix of eigen vectors to the identity matrix

  • Initialize the eigenvalues

  • Define the required dynamic convergence tolerance

  • Define the global tolerance \(\varepsilon_{\text{glob}}\)

  • For each cycle \(k = 1, n^{\text{max}}_{\text{jaco}}\) do

    • Define the dynamic tolerance: \(\varepsilon_k = (\varepsilon_{\text{dyn}})^k\)

    • \(l = 0\)

    • For each row \(i = 1, n\) do

      • For each column \(j = 1, n\) do

        • Calculate the coupling factors,

          \[f_{\mathbf{A}_{ij}^{k,l}} = \frac{|\mathbf{A}_{ij}^{k,l}|} {\sqrt{|\mathbf{A}_{ii}^{k,l}\mathbf{A}_{jj}^{k,l}|}} \quad f_{\mathbf{B}_{ij}^{k,l}} = \frac{|\mathbf{B}_{ij}^{k,l}|} {\sqrt{|\mathbf{B}_{ii}^{k,l}\mathbf{B}_{jj}^{k,l}|}} \quad i \ne j\]

        • If \(f_{\mathbf{A}_{ij}^{k,l}}\) or \(f_{\mathbf{B}_{ij}^{k,l}} \ge \varepsilon_k\) then

          • Calculate the coefficients of the Givens rotation

          • Transform matrices \(\mathbf{A}_{k,l}\) and \(\mathbf{B}_{k,l}\)

          • Transform eigenvectors

          • \(l = l+1\)

        • End If

      • End For

    • End For

    • Calculate the eigenvalues

      \[\lambda_i^k = \frac{\mathbf{A}_{ii}^{k,l}}{\mathbf{B}_{ii}^{k,l}}\]

    • Calculate

      \[f_\lambda = \underset{i}{\text{max}} \frac{|\lambda_i^k-\lambda_i^{k-1}|}{\lambda_i^{k-1}}\]

    • Calculate global coupling factors

      \[f_{\mathbf{A}} = \underset{i, j, i \ne j}{\text{max}} \frac{|\mathbf{A}_{ij}^{k,l}|} {\sqrt{|\mathbf{A}_{ii}^{k,l}\mathbf{A}_{jj}^{k,l}|}} \quad f_{\mathbf{B}} = \underset{i, j, i\ne j}{\text{max}} \frac{|\mathbf{B}_{ij}^{k,l}|} {\sqrt{|\mathbf{B}_{ii}^{k,l}\mathbf{B}_{jj}^{k,l}|}}\]

    • If \(f_{\mathbf{A}} \le \varepsilon_{\text{glob}}\) and \(f_{\mathbf{B}} \le \varepsilon_{\text{glob}}\) and \(f_{\lambda} \le \varepsilon_{\text{glob}}\) then

      • Correct the eigen modes (divide by \(\sqrt{\mathbf{B}_{ii}^{k}}\))

      • Exit

    • End if

  • End loop

Note

  • \(n^{\text{max}}_{\text{jaco}}\) is the maximum number of iterations allowed.

  • The matrices \(\mathbf{A}\) and \(\mathbf{B}\) being symmetrical, only half of the terms is stored in the form of a vector.

13.1.15.2.4. Givens rotation matrix

We seek to cancel the terms in positions \(i\) and \(j\) \(( i \ne j )\) in matrices \(\mathbf{A}^{k,l}\) and \(\mathbf{B}^{k,l}\) by multiplying them by a rotation matrix which has the form:

\[\mathbf{G}_{kk} = 1, ~,~~k = 1, n;~~ \mathbf{G}_{ij} = a;~~ \mathbf{G}_{ji} = b\]

the other terms being zero. We want

\[\mathbf{A}_{ij}^{k, l+1} = \mathbf{B}_{ij}^{k, l+1} = 0\]

which leads to:

\[\begin{split}a \mathbf{A}_{ii}^{k,l} + (1+ab) \mathbf{A}_{ij}^{k,l} + b \mathbf{A}_{jj}^{k,l} & = 0\\ a \mathbf{B}_{ii}^{k,l} + (1+ab) \mathbf{B}_{ij}^{k,l} + b \mathbf{B}_{jj}^{k,l} & = 0\end{split}\]

because

\[\mathbf{A}^{k, l+1} = \mathbf{G}^T \mathbf{A}^{k,l} \mathbf{G} \quad \text{and} \quad \mathbf{B}^{k, l+1} = \mathbf{G}^T \mathbf{B}^{k,l} \mathbf{G}\]

If the two equations are proportional, we have:

\[a = 0 \quad \text{and} \quad b = -\frac{\mathbf{A}_{ij}^{k,l}}{\mathbf{A}_{jj}^{k,l}}\]

Otherwise:

\[a = \frac{c_2}{d} \quad \text{and} \quad b = -\frac{c_1}{d}\]

with:

\[\begin{split}c_1 &= \mathbf{A}_{ii}^{k,l} \mathbf{B}_{ij}^{k,l} - \mathbf{B}_{ii}^{k,l} \mathbf{A}_{jj}^{k,l} \\ c_2 &= \mathbf{A}_{jj}^{k,l} \mathbf{B}_{ij}^{k,l} - \mathbf{B}_{jj}^{k,l} \mathbf{A}_{ij}^{k,l} \\ c_3 &= \mathbf{A}_{ii}^{k,l} \mathbf{B}_{jj}^{k,l} - \mathbf{B}_{ii}^{k,l} \mathbf{A}_{jj}^{k,l} \\ d &= \frac{c_3}{2}+ \text{sign}(c_3) \sqrt{\left(\frac{c_3}{2}\right)^2+ c_1 c_2}\end{split}\]

Note

If \(d = 0\), we are in the proportional case.