4. Subspace method (METHODE =” TRI_DIAG “/” SORENSEN “)#

4.1. From QEP to SEP#

4.1.1. Passage QEP/GEP: linearization#

As explained in §2.4, the transition from QEP to GEP takes place via the following linearization

(L2) » \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{-}\mathrm{K}& {0}_{n}\\ {0}_{n}& {\mathrm{M}}_{R}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{-}\mathrm{C}& \mathrm{M}\\ \mathrm{M}& {0}_{n}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (4.1-1)

The GEP obtained is structurally symmetric real and inherits the properties of the matrices \(\mathrm{K}\), \(\mathrm{M}\) and \(\mathrm{C}\). As soon as one of them is non-symmetric or complex (only possible with \(\mathrm{K}\)), GEP becomes non-symmetric or complex. The symmetric or not aspect of GEP is in fact not blocking because each method will work to find the « work operator/scalar product » pair that allows it to work. Besides, only Lanczos needs symmetric torque. Can the IRAM method work very well when non-symmetric and often with more robustness [21] _ .

4.1.2. Passage GEP/SEP: spectral transformation#

Then a spectral transformation (cf. [Boi09] §3.7/§5) makes it possible to complete the QEP pre-treatment phase by transforming the previous GEP \(({L}_{2})’’\) into an SEP. When manipulating complex modes, two scenarios are possible:

Work completely in**complex arithmetic*,

Remain**in real arithmetic* most often by isolating the real and imaginary contributions of the classical spectral transformation.

The first strategy classically uses the so-called « shift-and-invert » spectral transformation and metamorphose \(({L}_{2})’’\) into

\(({S}_{1})\) \(\begin{array}{c}\underset{{\mathrm{A}}_{\sigma }}{\underset{\underbrace{}}{{(\mathrm{A}\mathrm{-}\sigma \mathrm{B})}^{\mathrm{-}1}\mathrm{B}}}\mathrm{w}\mathrm{=}\mu \mathrm{w}\\ \text{avec}\mu \mathrm{:}\mathrm{=}\frac{1}{\lambda \mathrm{-}\sigma },\mathrm{w}\mathrm{:}\mathrm{=}\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\end{array}\) (4.1-2)

It is implemented by activating the complex approach of IRAM (“SORENSEN “+ APPROCHE =” COMPLEXE”). It is the only alternative when manipulating a complex \(\mathrm{K}\) matrix [22] _ . It can be used with \(\mathrm{K},\mathrm{M}\) and \(\mathrm{C}\) symmetric or not.

The second strategy is based on the « double shift sum » spectral transformation introduced in 1987 by B.N.Parlett & Y.Saad [PS87], which leads to the manipulation of two types of SEP

\(({S}_{2})\) \(\begin{array}{c}\text{Re}\underset{{\mathrm{A}}_{\sigma }^{\text{+}}}{\underset{\underbrace{}}{({(\mathrm{A}\mathrm{-}\sigma \mathrm{B})}^{\mathrm{-}1}\mathrm{B})}}\mathrm{w}\mathrm{=}{\mu }^{\text{+}}\mathrm{w}\\ \text{avec}{\mu }^{\text{+}}\mathrm{:}\mathrm{=}\frac{1}{2}(\frac{1}{\lambda \mathrm{-}\sigma }+\frac{1}{\lambda \mathrm{-}\stackrel{ˉ}{\sigma }}),\mathrm{w}\mathrm{:}\mathrm{=}\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\end{array}\) (4.1-3)

\(({S}_{3})\) \(\begin{array}{c}\text{Im}\underset{{\mathrm{A}}_{\sigma }^{\text{-}}}{\underset{\underbrace{}}{({(\mathrm{A}\mathrm{-}\sigma \mathrm{B})}^{\mathrm{-}1}\mathrm{B})}}\mathrm{w}\mathrm{=}{\mu }^{\text{-}}\mathrm{w}\\ \text{avec}{\mu }^{\text{-}}\mathrm{:}\mathrm{=}\frac{1}{\mathrm{2i}}(\frac{1}{\lambda \mathrm{-}\sigma }+\frac{1}{\lambda \mathrm{-}\stackrel{ˉ}{\sigma }}),\mathrm{w}\mathrm{:}\mathrm{=}\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\end{array}\) (4.1-4)

These two variants are implemented with Lanczos and IRAM via the values, respectively, “REEL” and “IMAG” of the keyword APPROCHE. With Lanczos, the scope of use is limited to classic QEP (real symmetric), while IRAM is less restricted and also allows non-symmetric ones.

Notes:

  • The origin of the term « double shift sum » of the spectral transformation above appears more clearly if we reformulate the operators as follows

\(\begin{array}{c}{\mathrm{A}}_{\sigma }^{\text{+}}\mathrm{:}\mathrm{=}\frac{1}{2}\left[{(\mathrm{A}\mathrm{-}\sigma \mathrm{B})}^{\mathrm{-}1}+{(\mathrm{A}\mathrm{-}\stackrel{ˉ}{\sigma }\mathrm{B})}^{\mathrm{-}1}\right]\mathrm{B}\\ {\mathrm{A}}_{\sigma }^{\text{-}}\mathrm{:}\mathrm{=}\frac{1}{\mathrm{2i}}\left[{(\mathrm{A}\mathrm{-}\sigma \mathrm{B})}^{\mathrm{-}1}+{(\mathrm{A}\mathrm{-}\stackrel{ˉ}{\sigma }\mathrm{B})}^{\mathrm{-}1}\right]\mathrm{B}\end{array}\) (4.1-5)

  • The « sum » terminology is in contrast to the « double shift product » originally proposed by J.C.Francis and used, moreover, secretly, by the QZ method and the « implicit » restarts of IRAM (cf. [Boi09] Appendix 1):

\({\mathrm{A}}_{\sigma }\mathrm{:}\mathrm{=}{\left[(\mathrm{A}\mathrm{-}\sigma \mathrm{B})+(\mathrm{A}\mathrm{-}\stackrel{ˉ}{\sigma }\mathrm{B})\right]}^{\mathrm{-}1}\mathrm{B}\) (4.1-6)

This has the major disadvantage of producing very full working matrices (therefore more expensive to handle and store) compared to the Parlett & Saad strategy.

Once the work operator has been determined, it remains to choose his associated dot product. For IRAM, we can be satisfied with a non-symmetric « operator/dot product » pair (cf. [Boi09] §6), but for Lanczos, this one must necessarily be symmetric (cf. [Boi09] §7). Several choices of matrix (nickname) scalar products allow this symmetry:

\(\begin{array}{c}{(\mathrm{w},\mathrm{z})}_{\mathrm{B}}\mathrm{:}\mathrm{=}{\mathrm{z}}^{T}\mathrm{Bw}\\ {(\mathrm{w},\mathrm{z})}_{\text{+}}\mathrm{:}\mathrm{=}{\mathrm{z}}^{T}\text{Re}{\left[{(\mathrm{A}\mathrm{-}\sigma \mathrm{B})}^{\mathrm{-}1}\right]}^{\mathrm{-}1}\mathrm{w}\\ {(\mathrm{w},\mathrm{z})}_{\text{-}}\mathrm{:}\mathrm{=}{\mathrm{z}}^{T}\text{Im}{\left[{(\mathrm{A}\mathrm{-}\sigma \mathrm{B})}^{\mathrm{-}1}\right]}^{\mathrm{-}1}\mathrm{w}\end{array}\) (4.1-7)

The first works with both formulations (S2) and (S3), while the second (resp. the third) only makes \({\mathrm{A}}_{\sigma }^{\text{+}}\) symmetric (resp. \({\mathrm{A}}_{\sigma }^{\text{-}}\)). In the effective implementation in Code_Aste r, it is these last two that are retained for Lanczos’s real and imaginary approaches.

Note:

  • These scalar products are the « natural » extensions of the (nickname) matrix dot product used by the Newmann-Pipano variant in GEP:

\((\mathrm{w},\mathrm{z})\mathrm{:}\mathrm{=}{\mathrm{z}}^{T}(\mathrm{A}\mathrm{-}\sigma \mathrm{B})\mathrm{w}\) (4.1-8)

4.2. Implementation in Code_Aster#

4.2.1. Choice of spectral shift#

In GEP there are four ways to choose this shift (cf. [Boi09] §5.4). In QEP, the “BANDE” option being proscribed due to the absence of an adapted Sturm test, this figure is reduced to three:

  • \(\sigma =0,\) we are looking for the smallest eigenvalues of the original problem. This corresponds to OPTION =” PLUS_PETITE “under the keyword factor CALC_FREQ.

  • \(\sigma ={\sigma }_{0}\text{avec}{\sigma }_{0}={(2\pi {f}_{0})}^{2},\) we are looking for frequencies close to the frequency FREQ = \({f}_{0}\) (OPTION =” CENTRE “).

  • \(\sigma ={\sigma }_{0}\text{avec}{\sigma }_{0}={(2\pi {f}_{0})}^{2}\text{et}{\sigma }_{1}=\eta {(2\pi {f}_{0})}^{2}/2\) we are looking for frequencies close to the frequency FREQ = \({f}_{0}\) (OPTION =” CENTRE “) and the reduced damping AMOR_REDUIT = \(\eta\) (OPTION =” CENTRE”).

The number of frequencies to be calculated is generally given by the user using NMAX_FREQ under the keyword factor CALC_FREQ.

4.2.2. Calculation of the operator and the dot product of work#

In the following table, we summarize all the possible pairs (work operator, dot product), according to the options chosen by the user.

APPROCHE/METHODE

“TRI_DIAG”

“SORENSEN” (by default)

“COMPLEXE”

Unavailable

\(({A}_{\sigma },{L}^{2})\)

“REEL” ( by default)

\(({A}_{\sigma }^{\text{+}},{()}_{\text{+}})\)

\(({A}_{\sigma }^{\text{+}},{L}^{2})\)

“IMAG”

\(({A}_{\sigma }^{\text{-}},{()}_{\text{-}})\)

\(({A}_{\sigma }^{\text{-}},{L}^{2})\)

Table 4.2-1. Set of possible couples (work operator, scalar product) according to the options chosen by the user.

The calculation of work operators \({\mathrm{A}}_{\sigma }^{\mathrm{\pm }}\) can be reduced to simple matrix-vector products and to an inversion of the dynamic matrix \(\mathrm{Q}(\lambda )\) (previously factored) via the following (clever) formulation

\({(\mathrm{A}\mathrm{-}\sigma \mathrm{B})}^{\mathrm{-}1}\mathrm{B}\mathrm{=}\left[\begin{array}{cc}{0}_{n}& {0}_{n}\\ {\mathrm{M}}_{R}^{\mathrm{-}1}\mathrm{M}& {0}_{n}\end{array}\right]\mathrm{-}\left[\begin{array}{cc}\mathrm{Q}{(\sigma )}^{\mathrm{-}1}& {0}_{n}\\ {0}_{n}& \sigma {\mathrm{M}}_{R}^{\mathrm{-}1}\mathrm{M}\mathrm{Q}{(\sigma )}^{\mathrm{-}1}\end{array}\right]\left[\begin{array}{cc}\mathrm{C}+\sigma \mathrm{M}& \mathrm{M}\\ \mathrm{C}+\sigma \mathrm{M}& \mathrm{M}\end{array}\right]\) (4.2-1)

On the other hand, this formulation has the « good taste » of allowing treatments in pure real arithmetic and of maintaining the sparse structure of the matrices. The following operating mode can thus be constructed:

  • Preparation in \(C\) [23] _

(only once, during the initialization phase of the algorithm, cf.

§2.5):

Former \(\mathrm{Q}(\sigma )\mathrm{:}\mathrm{=}({\sigma }^{2}\mathrm{M}+\sigma \mathrm{C}+\mathrm{K})\)

Factorize \(\mathrm{Q}(\sigma )\) in the form LDLT.

  • Calculation of \(\begin{array}{c}\\ {\mathrm{A}}_{{\sigma }^{\mathrm{\pm }}}\underset{\mathrm{w}}{\underset{\underbrace{}}{\left[\begin{array}{c}\mathrm{u}\\ \mathrm{v}\end{array}\right]}}\end{array}\) (several times per Lanczos or IRAM iteration):

Train in \(R{\mathrm{u}}_{1}\mathrm{:}\mathrm{=}\mathrm{Cu},{\mathrm{u}}_{2}\mathrm{:}\mathrm{=}\mathrm{Mu},{\mathrm{u}}_{3}\mathrm{:}\mathrm{=}\mathrm{Mv},\)

Calculate in \(C{\mathrm{u}}_{4}\mathrm{:}\mathrm{=}\mathrm{Q}{(\sigma )}^{\mathrm{-}1}+\sigma {\mathrm{u}}_{2}+{\mathrm{u}}_{3}\)

  • Depending on the choice of approach \(\begin{array}{c}{\mathrm{A}}_{\sigma }^{\text{+}}\mathrm{w}\mathrm{=}\left[\begin{array}{c}\mathrm{-}\text{Re}({\mathrm{u}}_{4})\\ {\mathrm{M}}_{R}^{\mathrm{-}1}\mathrm{M}(\mathrm{u}\mathrm{-}\text{Re}(\sigma {\mathrm{u}}_{4}))\end{array}\right]\\ {\mathrm{A}}_{\sigma }^{\text{-}}\mathrm{w}\mathrm{=}\left[\begin{array}{c}\mathrm{-}\text{Im}({\mathrm{u}}_{4})\\ \mathrm{-}{\mathrm{M}}_{R}^{\mathrm{-}1}\mathrm{M}(\mathrm{u}\mathrm{-}\text{Im}(\sigma {\mathrm{u}}_{4}))\end{array}\right]\end{array}\)

Algorithm 3. Procedure for calculating work operators \({\mathrm{A}}_{\sigma }^{\mathrm{\pm }}\) .

In a complex approach, the calculation is done directly via using the formula (4.1-2) in complex arithmetic.

With regard to matrix scalar products, their manipulations in Code_Aster are based on the following formulations:

\(\begin{array}{c}{(\mathrm{w},\mathrm{z})}_{\text{+}}\mathrm{:}\mathrm{=}{\mathrm{z}}^{T}\left[(\mathrm{A}\mathrm{-}\text{Re}(\sigma )\mathrm{B})+{\text{Im}}^{2}(\sigma )\mathrm{B}{(\mathrm{A}\mathrm{-}\text{Re}(\sigma )\mathrm{B})}^{\mathrm{-}1}\mathrm{B}\right]\mathrm{w}\\ {(\mathrm{w},\mathrm{z})}_{\text{-}}\mathrm{:}\mathrm{=}{\mathrm{z}}^{T}\left[\frac{1}{\text{Im}(\sigma )}(\mathrm{A}\mathrm{-}\text{Re}(\sigma )\mathrm{B}){\mathrm{B}}^{\mathrm{-}1}(\mathrm{A}\mathrm{-}\text{Re}(\sigma )\mathrm{B})+{\text{Im}}^{2}(\sigma )\mathrm{B}\right]\mathrm{w}\end{array}\) (4.2-2)

Note:

  • The calculation of matrix dot products can only be done in real arithmetic. For the first dot product we need an additional dynamic matrix \(\mathrm{Q}(\text{Re}(\sigma ))\) and its factorized (resulting in additional calculation and memory costs) .

4.3. Scope of use#

QEP with real symmetric matries for Lanczos, possibly nonsymmetric (and \(\mathrm{K}\) may be complex symmetric) for IRAM.

The user can only specify TYPE_RESU =” DYNAMIQUE “as the class to which their calculation belongs (no buckling). The vector FREQ is then optionally entered (and not CHAR_CRIT).

4.4. Display in the message file#

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

THE NOMBRE OF DDL

TOTAL EST: 130

FROM LAGRANGE EST: 16

THE NOMBRE OF DDL ACTIFS EST: 106

CALCUL MODAL: METHODE OF ITERATION SIMULTANEE

METHODE OF SORENSEN

NUMERO FREQUENCE (HZ) AMORTISSEMENT NORME D’ERREUR

1 1.24163E+02 -1.89229E-02 1.16550E-09

2 1.24164E+02 1.89229E-02 1.04882E-09

…

6 1.07321E+03 1.63495E-01 1.36092E-13

NORME from ERREUR MOYENNE: 0.36947E-09

Example 4. Trace of CALC_MODES with OPTIONparmi [“CENTRE”, “PLUS_PETITE”]*and METHODE =” SORENSEN “(QEP) in the message file. Excerpt from the sdll123b test case.

For each eigenvalue (represented in the form FREQUENCE = \(\frac{\text{Im}({\lambda }_{i})}{2\pi }\) and AMORTISSEMENT = and = \(\frac{-\text{Re}({\lambda }_{i})}{\mid {\lambda }_{i}\mid }\) cf §2.5), the error norm of the residue determined according to algorithm No. 1 (§2.5) is plotted.

4.5. Summary of the settings#

Let’s now summarize the setting of the CALC_MODES operator with OPTION among [“CENTRE”, “PLUS_PETITE”] for the Lanczos and Sorensen methods.

Operand

Keyword

Default Value

References

TYPE_RESU =” DYNAMIQUE “

“ DYNAMIQUE “

“ “

§2.1

OPTION =” PLUS_PETITE “

“ PLUS_PETITE “

“ “

§4.2

“CENTRE”

§4.2

CALC_FREQ

FREQ (si” CENTRE “)

§4.2

AMOR_REDUIT (si” CENTRE “)

§4.2

NMAX_FREQ

10

§4.2

SOLVEUR_MODAL

METHODE ** =” TRI_DIAG “

“ SORENSEN “

[Boi09] §6

“SORENSEN”

[Boi09] §7

APPROCHE = “REEL”

“REEL”

“”

§4.1

“IMAG”

§4.1

“ COMPLEXE “(if” SORENSEN “)

§4.1

DIM_SOUS_ESPACE COEF_DIM_ESPACE

Calculated

[Boi09] §6/7

…

Usual setting of “TRI_DIAG”/”SORENSEN”

[Boi09] §6/7

Table 4.5-1. Summary of the settings of CALC_MODES with OPTIONparmi [“CENTRE”, “PLUS_PETITE”] and* METHODE =” TRI_DIAG “or “ SORENSEN “ * ** (QEP) * .

Notes:

  • During the first few passes, it is strongly recommended to modify only the main parameters noted in bold. The others are more about the arcana of the algorithm and they have been initialized empirically to standard values.

  • In particular, to improve the quality of a mode, the fundamental parameter is the dimension of the projection subspace of SEP, DIM_SOUS_ESPACE .