ARPACK utilities (structsolve.arpackutils)#

Utilities for ARPACK, against the Intel MKL bug of dsteqr

The symmetric drivers of ARPACK, used by scipy.sparse.linalg.eigsh(), compute the Ritz vectors with the LAPACK routine dsteqr applied to the ncv x ncv tridiagonal matrix of the Lanczos process. The dsteqr of Intel MKL 2024.2.0 to 2025.0.0 returns wrong eigenvectors for matrices larger than 32 x 32, such that scipy.sparse.linalg.eigsh() returns wrong eigenpairs, without any error, or raises ArpackError -8, whenever ncv > 32, e.g. with SciPy’s default ncv = 2 k + 1 for k >= 16. MKL 2025.0.1 or newer is not affected, neither are the non-symmetric drivers of scipy.sparse.linalg.eigs(), which do not call dsteqr.

capped_eigsh() keeps ncv <= 32 when arpack_max_ncv() says so, computing more than 15 eigenpairs in several calls with deflation, whose eigenvectors are combined with a Rayleigh-Ritz projection such that they are orthonormal also across the calls.

The cap is decided by arpack_max_ncv() and can be overridden with the module attribute ARPACK_MAX_NCV or the environment variable STRUCTSOLVE_ARPACK_MAX_NCV, e.g. to lift it with a fixed MKL that cannot be identified as such:

import structsolve.arpackutils
structsolve.arpackutils.ARPACK_MAX_NCV = 0  # no cap
structsolve.arpackutils.ARPACK_MAX_NCV = None#

Maximum ncv of capped_eigsh(). None, the default, reads the environment variable STRUCTSOLVE_ARPACK_MAX_NCV and, when it is not defined, uses 'auto': the cap DSTEQR_MAX_SIZE is applied when the LAPACK of SciPy is detected as faulty, see arpack_max_ncv(). 0 (or 'none') disables the cap, an integer of at least 3 sets it.

structsolve.arpackutils.DSTEQR_MAX_SIZE = 32#

Largest tridiagonal matrix for which the faulty dsteqr is correct

structsolve.arpackutils.MKL_FAULTY_VERSIONS = ((2024, 2, 0), (2025, 0, 0))#

First and last Intel MKL versions with the faulty dsteqr

structsolve.arpackutils.arpack_max_ncv()#

Maximum ncv of capped_eigsh(), None when not capped

The value of ARPACK_MAX_NCV is used, or, when it is None, the environment variable STRUCTSOLVE_ARPACK_MAX_NCV. Both accept 'auto', 0 or 'none' (no cap), or an integer of at least 3. In the 'auto' mode, the default, ncv is capped at DSTEQR_MAX_SIZE when:

The detection runs once per process.

structsolve.arpackutils.arpack_ncv(k, n)#

ncv of scipy.sparse.linalg.eigsh() for k eigenpairs of a problem of size n

SciPy’s default min(n, max(2 k + 1, 20)), capped at arpack_max_ncv(). With a cap, k should not exceed (max_ncv - 1)//2, see capped_eigsh().

structsolve.arpackutils.capped_eigsh(A, k, M, sigma, which, mode='normal', tol=0)#

scipy.sparse.linalg.eigsh() with ncv from arpack_ncv()

Solves \([A]\{x\} = w [M]\{x\}\) with a positive definite M about the shift sigma, in the shift-invert (mode='normal') or Cayley (mode='cayley') mode, selecting the eigenvalues by increasing \(w\), e.g. the Cayley mode with which='SM' and sigma larger than all \(|w|\), or the shift-invert mode with which='LM' and sigma lower than all \(w\).

When ncv is capped at max_ncv, see arpack_max_ncv(), and k is larger than (max_ncv - 1)//2, i.e. 15, the eigenpairs are computed in passes of at most 15 eigenpairs, each with ncv=max_ncv. With k close to ncv ARPACK converges very slowly, not at all, or misses eigenvalues, e.g. k=31 and ncv=32. The first pass is a call in the requested mode. The following passes use the shift-invert mode about \(\sigma_d\), the middle of the largest gap between the last five eigenvalues found, with the eigenvectors \([V]\) of the previous passes deflated from the operator:

\[\{y\} = ([I] - [V][V]^T[M])([A] - \sigma_d [M])^{-1}\{x\}\]

such that the eigenvalues of \([V]\), already found, become null and are not found again. Each pass therefore returns the eigenvalues nearest \(\sigma_d\) among those not found yet, i.e. the next ones by increasing \(w\). A shift close to them, instead of sigma, separates them better: in the study of the panels cylinders it was 30 to 45% faster and gave eigenvalues 5 to 10 times more accurate. Each pass factorizes \([A] - \sigma_d [M]\) with SuperLU, and every call starts from start_vector(). The eigenvectors of all the passes are finally combined with a Rayleigh-Ritz projection, see _rayleigh_ritz(), such that they are \([M]\)-orthonormal also across the passes.

Returns:
w, xndarray

The eigenvalues and the \([M]\)-orthonormal eigenvectors, by increasing \(w\) when computed in passes.

structsolve.arpackutils.dsteqr_is_faulty()#

Tell whether the LAPACK dsteqr of SciPy returns wrong eigenvectors

scipy.linalg.lapack.dstev computes the eigenvectors of a symmetric tridiagonal matrix with dsteqr. It is called for three random matrices of sizes 33, 40 and 64, whose eigenpairs must have relative residuals and an orthogonality error below 1e-10. With Intel MKL 2025.0.0 every size above 32 fails, with rare exceptions for sizes close to 32, hence the three sizes.

Returns:
faultybool or None

None when dstev is not available.

structsolve.arpackutils.lapack_is_mkl()#

Tell whether the LAPACK of SciPy is Intel MKL

scipy.show_config(mode='dicts') is read first, NumPy’s build configuration only when SciPy’s does not name the LAPACK library. Older versions of SciPy without mode='dicts' are checked through scipy.__config__.

Returns:
is_mklbool or None

None when the build configuration does not tell.

structsolve.arpackutils.mkl_version()#

Version of the Intel MKL loaded in the process, as a tuple (major, minor, update)

The build configuration of SciPy only gives the major version, e.g. '2025', therefore the version is read at runtime from the optional packages mkl-service (mkl.get_version()) or threadpoolctl.

Returns:
versiontuple or None

None when neither package is installed or no MKL is loaded.

structsolve.arpackutils.release_memory()#

Collect the reference cycles of SciPy’s ARPACK wrappers in the browser

The objects of scipy.sparse.linalg.eigsh() and scipy.sparse.linalg.eigs() keep the SuperLU factorizations of the shift-invert modes in reference cycles, freed only by the cyclic garbage collector. Natively this is harmless, but the WebAssembly heap of Pyodide (sys.platform == 'emscripten') is exhausted before the collector runs, e.g. after about ten eigs calls on a problem with 1,000 dofs, raising MemoryError. gc.collect() is then called after each ARPACK call; it does nothing on other platforms.

structsolve.arpackutils.start_vector(n)#

Fixed starting vector of ARPACK, uniform in [-1, 1] as ARPACK’s random one, such that the results are reproducible