Appendix 3. Jacobi method
===========================




Principle
--------

The Jacobi method [LT86] makes it possible to calculate all the eigenvalues of a generalized problem whose matrices are definite positive and symmetric (the matrices obtained at each iteration by the Bathe &Wilson method verify these properties; cf. § :ref:`8 <Ref487282183>`). It consists in transforming the matrices :math:`A` and :math:`B` of the GEP :math:`\mathrm{A}\mathrm{u}\mathrm{=}\lambda \mathrm{B}\mathrm{u}` into diagonal matrices, using successively similar orthogonal transformations (Givens rotation matrices). The process can be summarised as follows:


.. image:: images/Object_437.svg
    :width: 396
    :height: 106


.. _RefImage_Object_437.svg:


.. image:: images/Object_438.svg
    :width: 396
    :height: 106


.. _RefImage_Object_438.svg:

**Algorithm 3.1. Jacobi process**

The eigenvalues are given by :math:`\lambda \mathrm{=}{\mathrm{A}}^{d}{({\mathrm{B}}^{d})}^{\text{-1}}` i.e. :math:`\lambda \mathrm{=}\frac{{\mathrm{A}}_{\mathit{ii}}^{d}}{{\mathrm{B}}_{\mathit{ii}}^{d}}` and the eigenvector matrix verifies:


.. image:: images/Object_441.svg
    :width: 396
    :height: 106


.. _RefImage_Object_441.svg:

Each matrix :math:`{\mathrm{Q}}_{k}` is chosen so that a diagonal and non-zero :math:`(i,j)` term of :math:`{\mathrm{A}}_{k}` or :math:`{\mathrm{B}}_{k}` is zero after the transformation.




A few choices
--------------

In this algorithm, we realize that the important points are the following:


* How do I choose which terms to cancel?


* How do you measure the diagonal character of matrices when :math:`k` approaches infinity?


* How do you measure convergence?




.. _Ref487282278:


Non-diagonal terms to cancel
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

To choose the terms to cancel, there are several methods:


* The first consists in each step :math:`k`, choosing the largest element in non-diagonal modulus from the :math:`{\mathrm{A}}_{k}` or :math:`{\mathrm{B}}_{k}` matrix and cancelling it by performing a rotation. This choice ensures the convergence of the Jacobi method but is relatively expensive (search for the maximum element).


* The second solution consists in cancelling all the non-diagonal elements of these matrices in succession by following the natural order :math:`{\mathrm{A}}_{13}^{k},\mathrm{....},{\mathrm{A}}_{\mathrm{1n}}^{k},{\mathrm{A}}_{23}^{k},\mathrm{....}`. When you get to :math:`{\mathrm{A}}_{n\mathrm{-}\mathrm{1,}n}^{k}`, you start the cycle again. This method is converging slowly.

A variant of this method, while following the natural order of the terms, consists in cancelling out only those that are greater than a given :math:`{\varepsilon }_{k}` precision. At the end of a cycle, the value of this criterion is reduced and the value of this criterion is repeated. It is this strategy that is used in*Code_Aster*.




Convergence test
~~~~~~~~~~~~~~~~~~~~~

To test the convergence and the diagonal nature of the matrices, we operate in this way. It is verified that all the coupling factors defined by:

:math:`{f}_{{\mathrm{A}}_{\mathit{ij}}}\mathrm{=}\frac{∣{\mathrm{A}}_{\mathit{ij}}∣}{\sqrt{∣{\mathrm{A}}_{\mathit{ii}}{\mathrm{A}}_{\mathit{jj}}∣}}` :math:`{f}_{{\mathrm{B}}_{\mathit{ij}}}\mathrm{=}\frac{∣{\mathrm{B}}_{\mathit{ij}}∣}{\sqrt{∣{\mathrm{B}}_{\mathit{ii}}{\mathrm{B}}_{\mathit{jj}}∣}}` :math:`i\ne j`

are less than a given precision (diagonal nature of the matrices). The convergence of eigenvalues is also checked via the indicator.

:math:`{f}_{\lambda }=\underset{i}{\mathrm{max}}\frac{∣{\lambda }_{i}^{k}-{\lambda }_{i}^{k-1}∣}{{\lambda }_{i}^{k-1}}`

so that it stays below a given precision :math:`{\varepsilon }_{\mathrm{jaco}}`.




Algorithm implemented in Code_Aster
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The algorithm implemented in*Code_Aster* comes down to:

Initialization of the vector matrix specific to the identity matrix.

Initialization of eigenvalues.

Define the required dynamic convergence accuracy.

Define global accuracy :math:`{\varepsilon }_{\mathrm{glob}}`.

For each :math:`k\mathrm{=}\mathrm{1,}{\text{n\_max}}_{\mathit{jaco}}` cycle

Define dynamic tolerance: :math:`{\varepsilon }_{k}={({\varepsilon }_{\mathrm{dyn}})}^{k}`,

:math:`l=0`,

For each line :math:`i=\mathrm{1,}n`

For each column :math:`j=\mathrm{1,}n`

Calculating coupling factors, :math:`{f}_{{\mathrm{A}}_{\mathit{ij}}^{k,l}}\mathrm{=}\frac{∣{\mathrm{A}}_{\mathit{ij}}^{k,l}∣}{\sqrt{∣{\mathrm{A}}_{\mathit{ii}}^{k,l}{\mathrm{A}}_{\mathit{jj}}^{k,l}∣}}` :math:`{f}_{{\mathrm{B}}_{\mathit{ij}}^{k,l}}\mathrm{=}\frac{∣{\mathrm{B}}_{\mathit{ij}}^{k,l}∣}{\sqrt{∣{\mathrm{B}}_{\mathit{ii}}^{k,l}{\mathrm{B}}_{\mathit{jj}}^{k,l}∣}}` :math:`i\ne j`

If :math:`{f}_{{A}_{\mathrm{ij}}^{k,l}}\text{ou}{f}_{{B}_{\mathrm{ij}}^{k,l}}\ge {\varepsilon }_{k}` then

Calculation of the Givens rotation coefficients,

Transformation of matrices :math:`{\mathrm{A}}^{k\mathrm{,}\mathrm{l}}` and :math:`{\mathrm{B}}^{k\mathrm{,}\mathrm{l}}`,

Transformation of eigenvectors,

:math:`l=l+1`

End yes.

End loop.

End loop.

Calculating eigenvalues :math:`{\lambda }_{i}^{k}\mathrm{=}\frac{{\mathrm{A}}_{\mathit{ii}}^{k,l}}{{\mathrm{B}}_{\mathit{ii}}^{k,l}}`.

Calculation of :math:`{f}_{\lambda }=\underset{i}{\mathrm{max}}\frac{∣{\lambda }_{i}^{k}-{\lambda }_{i}^{k-1}∣}{{\lambda }_{i}^{k-1}}`.

Calculating global coupling factors

:math:`{f}_{\mathrm{A}}\mathrm{=}\underset{\begin{array}{c}i,j\\ i\mathrm{\ne }j\end{array}}{\mathit{Max}}\frac{∣{\mathrm{A}}_{\mathit{ij}}^{k,l}∣}{\sqrt{∣{\mathrm{A}}_{\mathit{ii}}^{k,l}{\mathrm{A}}_{\mathit{jj}}^{k,l}∣}}` :math:`{f}_{\mathrm{B}}\mathrm{=}\underset{\begin{array}{c}i,j\\ i\mathrm{\ne }j\end{array}}{\mathit{Max}}\frac{∣{\mathrm{B}}_{\mathit{ij}}^{k,l}∣}{\sqrt{∣{\mathrm{B}}_{\mathit{ii}}^{k,l}{\mathrm{B}}_{\mathit{jj}}^{k,l}∣}}`

If :math:`{f}_{\text{A}}\le {\varepsilon }_{\mathrm{glob}}\text{et}{f}_{\text{B}}\le {\varepsilon }_{\mathrm{glob}}\text{et}{f}_{\lambda }\le {\varepsilon }_{\mathrm{glob}}` then

Correcting the natural modes (let's divide by :math:`\sqrt{{\mathrm{B}}_{\mathit{ii}}^{k}}`),

Exit;

End yes.

End loop.

**Algorithm 3.1. Jacobi method implemented in Code_Aster.**

Note :math:`{\text{n\_max}}_{\mathrm{jaco}}` as the maximum number of iterations allowed.

Note:


* Since matrices :math:`\mathrm{A}` and :math:`B` are symmetric, only half of the terms are stored as a vector.




.. _Ref487282319:


Givens rotation matrix
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

At each step, we try to cancel the terms found in position :math:`i` and :math:`j` :math:`(i\ne j)` of the matrices :math:`{A}^{k,l}` and :math:`{B}^{k,l}` by multiplying them by a rotation matrix which has the following form:

:math:`{\mathrm{G}}_{\text{ll}}\mathrm{=}1l\mathrm{=}\mathrm{1,}n;{\mathrm{G}}_{\mathit{ij}}\mathrm{=}a;{\mathrm{G}}_{\mathit{ji}}\mathrm{=}b`

the other terms being null. We want to have :math:`{\mathrm{A}}_{\mathit{ij}}^{k,l+1}\mathrm{=}{\mathrm{B}}_{\mathit{ij}}^{k,l+1}\mathrm{=}0` which leads to:

:math:`\mathrm{\{}\begin{array}{c}a{\mathrm{A}}_{\mathit{ii}}^{k,l}+(1+\mathit{ab}){\mathrm{A}}_{\mathit{ij}}^{k,l}+b{\mathrm{A}}_{\mathit{jj}}^{k,l}\mathrm{=}0\\ a{\mathrm{B}}_{\mathit{ii}}^{k,l}+(1+\mathit{ab}){\mathrm{B}}_{\mathit{ij}}^{k,l}+b{\mathrm{B}}_{\mathit{jj}}^{k,l}\mathrm{=}0\end{array}`

because :math:`{\mathrm{A}}^{k,l+1}\mathrm{=}{\mathrm{G}}^{T}{\mathrm{A}}^{k,l}\mathrm{G}` and :math:`{\mathrm{B}}^{k,l+1}\mathrm{=}{\mathrm{G}}^{T}{\mathrm{B}}^{k,l}\mathrm{G}`. If the two equations are proportional we have:

:math:`a=0` and :math:`b\mathrm{=}\frac{\mathrm{-}{\mathrm{A}}_{\mathit{ij}}^{k,l}}{{\mathrm{A}}_{\mathit{jj}}^{k,l}}`

otherwise:

:math:`a=\frac{{c}_{2}}{d}b=\frac{-{c}_{1}}{d}`

with:

:math:`{c}_{1}\mathrm{=}{\mathrm{A}}_{\mathit{ii}}^{k,l}{\mathrm{B}}_{\mathit{ij}}^{k,l}\mathrm{-}{\mathrm{B}}_{\mathit{ii}}^{k,l}{\mathrm{A}}_{\mathit{jj}}^{k,l}` :math:`{c}_{2}\mathrm{=}{\mathrm{A}}_{\mathit{jj}}^{k,l}{\mathrm{B}}_{\mathit{ij}}^{k,l}\mathrm{-}{\mathrm{B}}_{\mathit{jj}}^{k,l}{\mathrm{A}}_{\mathit{ij}}^{k,l}`

:math:`{c}_{3}\mathrm{=}{\mathrm{A}}_{\mathit{ii}}^{k,l}{\mathrm{B}}_{\mathit{jj}}^{k,l}\mathrm{-}{\mathrm{B}}_{\mathit{ii}}^{k,l}{\mathrm{A}}_{\mathit{jj}}^{k,l}` :math:`d=\frac{{c}_{3}}{2}+\mathrm{sign}({c}_{3})\sqrt{{(\frac{{c}_{3}}{2})}^{2}+{c}_{1}{c}_{2}}`

Note:


* If :math:`d=0` then we are in the proportional case.