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).

  • 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 \([K]^{-1}[K_G]\), retried with a 10 times larger shift. The dense fallback condenses out the degrees of freedom where \([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([K] + s [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.

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 \([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.