Eigenvalue solvers#

The eigenvalue analyses lb() and freq() rely on sparse solvers and keep dense solvers only as a last resort. Sparse solvers scale with the number of non-null entries of the matrices, while a dense solution scales with the cube of the number of degrees of freedom. This matters in particular in the browser (Pyodide), where dense LAPACK is single-threaded, and for follower loads, whose load stiffness couples almost all the degrees of freedom.

Solver policy#

  • sparse_solver=True, the default, always runs the sparse solvers first, for every problem size, symmetric or unsymmetric. Every result is verified before it is returned: the relative residual of each eigenpair and, for symmetric matrices, the inertia of a shifted matrix, which tells whether an eigenvalue was missed (Sylvester’s law of inertia). The residual \(\|\boldsymbol{K} \boldsymbol{u} + \lambda \boldsymbol{B} \boldsymbol{u}\|\) is relative to \(\|\boldsymbol{K} \boldsymbol{u}\| + |\lambda| \, \|\boldsymbol{B} \boldsymbol{u}\|\), with \(\boldsymbol{B} = \boldsymbol{K}_G\) or \(\boldsymbol{M}\). Only for a rigid-body mode of freq(), whose eigenvalue and \(\boldsymbol{K} \boldsymbol{u}\) are both null up to the round-off, i.e. \(\|\boldsymbol{K} \boldsymbol{u}\| + |\lambda| \, \|\boldsymbol{M} \boldsymbol{u}\|\) not larger than RIGID_BODY_RTOL \(= 100 \epsilon\), \(\epsilon\) the machine precision, times the matrix-norm scale \((\|\boldsymbol{K}\|_F + |\lambda| \, \|\boldsymbol{M}\|_F)\|\boldsymbol{u}\|\), is the residual relative to this scale, see linear_buckling._relative_residuals(). An eigenpair of freq() whose relative residual is larger than check_rtol passes when its normwise backward error, the residual relative to the matrix-norm scale, is not larger than BACKWARD_ERROR_RTOL \(= \epsilon\), i.e. when it is accurate up to the round-off, e.g. an eigenvalue close to zero of a model pre-stressed close to buckling, see linear_buckling._backward_errors(). The check of the other eigenpairs, e.g. the inaccurate ones of models with stiff penalty springs, whose \(\|\boldsymbol{K}\|_F\) is many orders of magnitude larger than \(\|\boldsymbol{K} \boldsymbol{u}\|\), is not loosened, and lb() keeps the relative residual only.

  • When all the sparse solvers failed (verification, ARPACK error or no convergence, singular factorization, no positive real load multiplier for unsymmetric matrices…), a dense solver is used as a last resort, if the dense problem has at most max_dense_size degrees of freedom (default 2000, 0 disables the dense fallback). A DenseFallbackWarning gives the reasons of the failures, so that a dense solution is never silent, and the dense result is verified as well.

  • Otherwise a RuntimeError lists every attempted solver and the reason of its failure. No silently wrong result is returned.

  • sparse_solver=False explicitly solves the full problem with a dense solver, scipy.linalg.eigh() or scipy.linalg.eig(), without trying the sparse solvers, returning all the eigenvalues. It is a reference for validations and for the callers that need all of them.

The dense fallback makes structsolve robust in environments where ARPACK is unreliable, see below, while normal runs never pay for a dense solution. To turn the fallback into an error, e.g. in a test suite:

import warnings
from structsolve import DenseFallbackWarning
warnings.simplefilter('error', DenseFallbackWarning)

Linear buckling, symmetric matrices#

scipy.sparse.linalg.eigsh() in Cayley mode, with a shift estimated by power iterations on \(\boldsymbol{K}^{-1} \boldsymbol{K}_G\), retried with a 10 times larger shift. It gives the lowest positive load multipliers; when fewer than num_eigvalues exist, the same solver with \(-\boldsymbol{K}_G\) gives the negative ones closest to zero, i.e. the critical load multipliers of the reversed load, which are verified in the same way. The dense fallback condenses out the degrees of freedom where \(\boldsymbol{K}_G\) has null rows and solves the condensed problem with scipy.linalg.eigh().

Linear buckling, unsymmetric matrices#

With follower loads many negative or complex load multipliers may have a smaller modulus than the critical one, e.g. under an internal pressure, and the multipliers of smallest modulus may contain no positive real one. The sparse solver is a search of shifts \(s\) along the positive real axis: at each shift, scipy.sparse.linalg.eigs() finds the load multipliers nearest \(s\), within a disc that contains all of them. The shifts move until the lowest positive real multiplier found is inside the union of the discs that covers the positive real axis from zero, such that none was missed. A parity check of the sign of \(\det(\boldsymbol{K} + s \boldsymbol{K}_G)\) confirms that an even number, normally zero, of real multipliers lies below the critical one. See structsolve.linear_buckling._eigs_shift_search().

Frequency analysis#

scipy.sparse.linalg.eigs(), or scipy.sparse.linalg.eigsh() with symmetric=True, in shift-invert mode with a negative shift. The dense fallback returns the lowest frequencies, or the eigenvalues closest to the shift for unsymmetric matrices.

Factorizations with PARDISO#

linear_solver='pardiso' of lb() and freq() (default 'splu', SuperLU) factorizes the matrices of the sparse solvers of symmetric matrices with PARDISO of Intel oneMKL, see structsolve.pardiso, through static._factorize(): the shifted matrix of the Cayley or shift-invert mode, passed to scipy.sparse.linalg.eigsh() as OPinv, also in the passes of arpackutils.capped_eigsh(), the stiffness matrix of the estimates of the shifts, and the shifted matrices of the inertia counts, whose number of negative pivots is then that of the symmetric indefinite factorization of PARDISO, known also where SuperLU would need a row interchange. ARPACK remains the eigenvalue solver. A factorization of PARDISO is used only when the relative residual of its solution for a random right-hand side is not larger than arpackutils.PARDISO_ACCURACY_RTOL and, for the counts, no pivot remains perturbed; SuperLU factorizes that matrix otherwise, see arpackutils.pardiso_solve() and linear_buckling._negative_pivots(). The solvers of unsymmetric matrices and the dense fallbacks are unchanged.

Intel MKL and ARPACK#

The LAPACK routine dsteqr of Intel MKL 2024.2.0 to 2025.0.0 returns wrong eigenvectors for matrices larger than 32 x 32. ARPACK’s symmetric drivers, used by scipy.sparse.linalg.eigsh(), apply it to a matrix of size ncv, the dimension of the Krylov subspace, and then return wrong eigenpairs or raise ArpackError -8, e.g. with some Anaconda builds of SciPy. structsolve.arpackutils detects such a LAPACK and keeps ncv <= 32, computing more than 15 eigenpairs in several passes with deflation. The eigenvectors of all the passes are combined with a Rayleigh-Ritz projection, solved with the LAPACK driver dsyevr, which does not use dsteqr, such that they are \(\boldsymbol{M}\)-orthonormal also across the passes, where ARPACK alone gives them orthogonal only to its accuracy. The detection can be overridden, e.g. to lift the cap with MKL 2025.0.1 or newer when its version cannot be read:

STRUCTSOLVE_ARPACK_MAX_NCV=0 python my_script.py

or, in Python:

import structsolve.arpackutils
structsolve.arpackutils.ARPACK_MAX_NCV = 0

The non-symmetric drivers of scipy.sparse.linalg.eigs() are not affected.

In the browser#

In Pyodide (sys.platform == 'emscripten') SciPy uses OpenBLAS, ncv is not capped, and structsolve.arpackutils.release_memory() collects the reference cycles of SciPy’s ARPACK wrappers after each call, which otherwise keep the sparse factorizations alive until the WebAssembly heap is exhausted. The dense solvers of unsymmetric matrices call scipy.linalg.eig() with complex matrices, since the real drivers crash the Python process for some matrices with Intel MKL 2025.0.0.