Modal eigenvalue problem and shift-invert Lanczos
=================================================

Free-vibration modal analysis solves the
**generalised symmetric-definite** eigenvalue problem

.. math::
    :label: gevp

    \mathbf{K}\, \boldsymbol{\phi} \;=\; \omega^{2}\, \mathbf{M}\, \boldsymbol{\phi},

with :math:`\mathbf{K}` the assembled stiffness matrix
(:doc:`assembly`) and :math:`\mathbf{M}` the mass matrix.  Both
are sparse, symmetric, and positive (semi-)definite —
:math:`\mathbf{K}` becomes positive-definite once the rigid-
body modes are removed by Dirichlet BCs (:doc:`bc_elimination`).

This chapter develops the algorithm
:func:`scipy.sparse.linalg.eigsh` runs underneath
femorph-solver's :meth:`Model.modal_solve`: an Arnoldi /
Lanczos iteration on the **shift-inverted** operator that
converts the smallest-eigenvalue search into the largest-
magnitude search ARPACK was designed for.

Why not just solve :math:`\mathbf{K}^{-1}\, \mathbf{M}`?
--------------------------------------------------------

The straightforward power iteration on
:math:`\mathbf{K}^{-1}\, \mathbf{M}` extracts the **largest**
eigenvalue of that operator — equivalently the **smallest**
:math:`\omega^{2}` of :eq:`gevp`.  But power iteration
converges at a rate set by the spectral gap
:math:`|\lambda_2 / \lambda_1|`, which for the lowest few
modes is typically near 1.  Lanczos accelerates this by
projecting onto a Krylov subspace and computing all extreme
eigenvalues simultaneously; ARPACK is the canonical
implementation (Lehoucq, Sorensen & Yang 1998).

Shift-invert Lanczos
--------------------

Define the **shifted-inverted operator**

.. math::
    :label: shift-invert

    \mathbf{S}_{\sigma}
    \;=\;
    (\mathbf{K} - \sigma\, \mathbf{M})^{-1}\, \mathbf{M}.

The eigenvalues of :math:`\mathbf{S}_{\sigma}` are
:math:`\mu = 1 / (\omega^{2} - \sigma)`.  Eigenvalues
:math:`\omega^{2}` near :math:`\sigma` map to the
**largest-magnitude** eigenvalues of :math:`\mathbf{S}_{\sigma}`
— exactly what Lanczos converges to fastest.  The default
shift in femorph-solver is :math:`\sigma = 0`, so
:math:`\mathbf{S}_{0} = \mathbf{K}^{-1}\, \mathbf{M}` and the
lowest :math:`n_\mathrm{modes}` are returned by ascending
:math:`\omega^{2}`.

Implementation:

1. **Factor** :math:`\mathbf{K} - \sigma\, \mathbf{M}` once
   through the registered linear backend (Pardiso / CHOLMOD /
   MUMPS / SuperLU auto-chain — see
   :func:`femorph_solver.solvers.linear.list_linear_solvers`).
2. Each Lanczos iteration applies
   :math:`\mathbf{S}_{\sigma}\, \mathbf{v}` by
   solving :math:`(\mathbf{K} - \sigma \mathbf{M})\, \mathbf{w}
   = \mathbf{M}\, \mathbf{v}` against the cached factor.
3. ARPACK's reverse-communication interface returns
   :math:`(\mu, \mathbf{q})` pairs; femorph-solver back-
   transforms via :math:`\omega^{2} = \sigma + 1/\mu` and
   returns the original eigenvectors :math:`\boldsymbol{\phi}
   = \mathbf{q}` (Lanczos eigenvectors of
   :math:`\mathbf{S}_{\sigma}` and of the original GEVP
   :eq:`gevp` are the same up to back-transformation of the
   eigenvalue).

The whole flow factors :math:`\mathbf{K} - \sigma \mathbf{M}`
**once** and re-uses the factor for every Lanczos
matrix-vector product — a critical optimisation because the
factor cost dominates a sparse direct solve.  See Ericsson &
Ruhe (1980) for the original treatment.

Cluster handling near :math:`\sigma`
------------------------------------

Eigenvalues clustered near the shift produce ill-conditioned
:math:`\mathbf{S}_{\sigma}` factors and slow Lanczos
convergence.  The standard remedy is **shift-and-invert with a
slightly off-cluster shift**: pick :math:`\sigma` near but not
*at* the cluster.  ARPACK exposes this via the ``sigma``
parameter to ``eigsh``; femorph-solver currently uses
:math:`\sigma = 0` (suitable for static-stress-free models
where rigid-body modes have already been removed by BCs) but
the API allows an explicit override.  When a frequency-shifted
sweep is needed (e.g., to extract modes near 1 kHz on a model
whose first 50 modes are below 100 Hz), set
``Model.modal_solve(n_modes=N, sigma=2 * pi * 1000)``.

Real-symmetric vs complex-Hermitian
-----------------------------------

For an isotropic linear-elastic structure with conforming
mesh, both :math:`\mathbf{K}` and :math:`\mathbf{M}` are real
symmetric.  ``scipy.sparse.linalg.eigsh`` (the real-symmetric
Lanczos) is the natural choice and femorph-solver's modal
solver dispatches there.

The cyclic-symmetry path (:doc:`assembly`) reduces the global
GEVP to a per-harmonic complex-Hermitian sub-problem
:math:`\mathbf{K}_{k}\, \boldsymbol{\phi}_{k} = \omega^{2}\,
\mathbf{M}_{k}\, \boldsymbol{\phi}_{k}`.  Rather than calling
the complex ``eigs`` (which is slower and lacks
positive-definite support), femorph-solver applies the
**Grimes–Lewis–Simon (1994) real 2n-augmentation**:

.. math::

    \tilde{\mathbf{K}}_{k}
    =
    \begin{bmatrix}
        \mathrm{Re}(\mathbf{K}_{k}) & -\mathrm{Im}(\mathbf{K}_{k}) \\
        \mathrm{Im}(\mathbf{K}_{k}) & \mathrm{Re}(\mathbf{K}_{k})
    \end{bmatrix},
    \qquad
    \tilde{\mathbf{M}}_{k}
    =
    \begin{bmatrix}
        \mathrm{Re}(\mathbf{M}_{k}) & -\mathrm{Im}(\mathbf{M}_{k}) \\
        \mathrm{Im}(\mathbf{M}_{k}) & \mathrm{Re}(\mathbf{M}_{k})
    \end{bmatrix},

which are real-symmetric and **share the same eigenvalue
spectrum** (each appearing with multiplicity 2).  The same
``eigsh`` shift-invert path runs unchanged on the doubled
matrices; the eigenvectors recombine into a complex pair on
the way out.

Mode normalisation
------------------

Two conventions for normalising :math:`\boldsymbol{\phi}` see
common use:

* **Mass-orthonormal** —
  :math:`\boldsymbol{\phi}_{i}^{\!\top}\, \mathbf{M}\,
  \boldsymbol{\phi}_{j} = \delta_{ij}`.  Used by Bathe §8 and
  by femorph-solver's modal API.
* **Unit-amplitude** — :math:`\max_n |\phi_{i,n}| = 1`.  Used
  by some legacy MAPDL pipelines.

mass-orthonormal is the natural choice for modal-superposition
post-processing because the modal mass matrix becomes the
identity and the harmonic / transient response decouples per
mode.

What ``Model.modal_solve`` returns
----------------------------------

:meth:`Model.modal_solve` populates a
:class:`ModalResult` instance with:

* ``frequency`` — :math:`(\omega / 2\pi)` in Hz, ascending.
* ``omega_sq`` — :math:`\omega^{2}` in (rad/s)², matching the
  GEVP eigenvalue directly.
* ``mode_shapes`` — :math:`\boldsymbol{\phi}` arrays
  mass-orthonormalised to ``Model.dof_map()`` ordering.

Convergence diagnostics (Lanczos iteration count, ARPACK
``info`` flag) are surfaced as attributes when ``progress=True``.

References
----------

* Parlett, B. N. (1998) *The Symmetric Eigenvalue Problem*,
  SIAM, §13 (Lanczos), §14 (shift-invert).
* Lehoucq, R. B., Sorensen, D. C. and Yang, C. (1998)
  *ARPACK Users' Guide: Solution of Large-Scale Eigenvalue
  Problems with Implicitly Restarted Arnoldi Methods*, SIAM
  Software Environments and Tools 6.
* Ericsson, T. and Ruhe, A. (1980) "The spectral
  transformation Lanczos method for the numerical solution
  of large sparse generalized symmetric eigenvalue problems,"
  *Mathematics of Computation* 35, 1251–1268 (the original
  shift-invert formulation).
* Saad, Y. (2011) *Numerical Methods for Large Eigenvalue
  Problems*, 2nd ed., SIAM, §6 (Lanczos for symmetric
  problems), §13 (shift-and-invert strategies).
* Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed.,
  §11 (eigensolution methods), §12 (mode-superposition
  dynamics).
* Wilkinson, J. H. (1965) *The Algebraic Eigenvalue Problem*,
  Oxford UP — foundational treatment of symmetric-definite
  eigenproblem stability.
* Grimes, R. G., Lewis, J. G. and Simon, H. D. (1994) "A
  shifted block Lanczos algorithm for solving sparse symmetric
  generalized eigenproblems," *SIAM J. Matrix Analysis and
  Applications* 15, 228–272 (real 2n augmentation for the
  complex-Hermitian cyclic-symmetry sub-problem).