Utilities for ARPACK, against the Intel MKL bug of dsteqr (pyfe3d.solver.arpackutils)#

With Intel MKL 2024.2.0 to 2025.0.0, scipy.sparse.linalg.eigsh() returns wrong eigenpairs, without any error, or raises ArpackError -8, whenever its number of Lanczos vectors is ncv > 32, e.g. with SciPy’s default ncv = 2 k + 1 for k >= 16 eigenpairs. MKL 2025.0.1 or newer is not affected, neither is scipy.sparse.linalg.eigs(). The detection and the cap depend only on NumPy and SciPy.

capped_eigsh() keeps ncv <= 32 when arpack_max_ncv() says so, computing more than 15 eigenpairs in several calls with a deflation of the eigenvectors already found, which are combined with a Rayleigh-Ritz projection such that they are orthonormal also across the calls. It is used by linear_buckling() and natural_frequency() for symmetric matrices.

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

import pyfe3d.solver.arpackutils
pyfe3d.solver.arpackutils.ARPACK_MAX_NCV = 0  # no cap

Notes#

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, hence the wrong eigenpairs whenever ncv > 32. The non-symmetric drivers of scipy.sparse.linalg.eigs() do not call dsteqr.

pyfe3d.solver.arpackutils.ARPACK_MAX_NCV = None#

Maximum ncv of capped_eigsh(). None, the default, reads the environment variable PYFE3D_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.

pyfe3d.solver.arpackutils.DSTEQR_MAX_SIZE = 32#

Largest tridiagonal matrix for which the faulty dsteqr is correct

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

First and last Intel MKL versions with the faulty dsteqr

pyfe3d.solver.arpackutils._arpack_takes_rng()#

Tell whether scipy.sparse.linalg.eigs() and scipy.sparse.linalg.eigsh() accept the argument rng, i.e. SciPy 1.17 or newer

pyfe3d.solver.arpackutils._auto_max_ncv()#

Cap of ncv in the 'auto' mode of arpack_max_ncv(), DSTEQR_MAX_SIZE or None, decided once per process

pyfe3d.solver.arpackutils._deflation_shift(w, far_shift)#

Shift of the passes after the first one of capped_eigsh()

The middle of the largest gap between the last five eigenvalues found, sorted by increasing value. far_shift when they are all equal, i.e. when their largest gap is not larger than 1e-8 times their largest modulus.

Notes

The shift-invert mode finds first the eigenvalues nearest the shift. The middle of the largest gap is close to the next eigenvalues and away from the deflated ones.

pyfe3d.solver.arpackutils._lapack_name_from_config(module)#

Name of the LAPACK library in the build configuration of a module

Parameters:
modulemodule

scipy or numpy.

Returns:
namestr or None

The name in lower case, None when unknown. Builds with numpy.distutils give 'mkl' when their configuration has an MKL entry.

pyfe3d.solver.arpackutils._parse_mkl_version(text)#

Version (major, minor, update) in a text such as '2025.0.1', None when there is none; a missing update is 0

pyfe3d.solver.arpackutils._rayleigh_ritz(A, M, V)#

Rayleigh-Ritz projection of A x = w M x onto the columns of V

Returns the Ritz values and the M-orthonormal Ritz vectors of the eigenvalue problem of the symmetric A and the positive definite M in the subspace spanned by the columns of V, which must be linearly independent. Used by capped_eigsh() to combine the eigenvectors of its passes. The projected problem is solved with the LAPACK driver dsyevr (driver='evr') of scipy.linalg.eigh().

Returns:
w, xndarray

The Ritz values by increasing value and the Ritz vectors.

Notes

The eigenvectors of different passes of capped_eigsh() are only \(\boldsymbol{M}\)-orthogonal to the accuracy of ARPACK, e.g. \(10^{-6}\) between close eigenvalues of different passes. \(\boldsymbol{V}\) is \(\boldsymbol{M}\)-orthonormalized with the Cholesky factor \(\boldsymbol{L}\) of \(\boldsymbol{V}^\top \boldsymbol{M} \boldsymbol{V}\), and the projected problem \(\boldsymbol{W}^\top \boldsymbol{A} \boldsymbol{W} \boldsymbol{y} = w \boldsymbol{y}\), with \(\boldsymbol{W} = \boldsymbol{V} \boldsymbol{L}^{-\top}\), is solved with the LAPACK driver dsyevr (driver='evr'), which does not use the dsteqr of Intel MKL broken above 32 x 32 matrices, see dsteqr_is_faulty(). The returned eigenvectors \(\boldsymbol{W} \boldsymbol{y}\) are \(\boldsymbol{M}\)-orthonormal to the rounding errors, and the Ritz values are at least as accurate as the eigenvalues of the passes. \(\boldsymbol{V}^\top \boldsymbol{M} \boldsymbol{V}\) and \(\boldsymbol{W}^\top \boldsymbol{A} \boldsymbol{W}\) are symmetrized, \(\frac{1}{2}(\boldsymbol{B} + \boldsymbol{B}^\top)\), before the Cholesky factorization and the eigenvalue solution.

The Ritz values are the optimal approximations of the eigenvalues that can be obtained from the subspace spanned by the columns of \(\boldsymbol{V}\) ([1], Chapter 11, “Approximations from a Subspace”, Section 11.4, “Optimality”, Eq. (11.5) and Theorem 11.4.1, p. 236): the \(j\)-th Ritz value is the minimum over the \(j\)-dimensional subspaces of the span of the maximum of the Rayleigh quotient, so it can only improve when the span grows, and the span contains the eigenvectors of every pass. The theorem is stated for the standard problem; it holds for the pencil in the \(\boldsymbol{M}\) inner product, through the Cholesky factor of \(\boldsymbol{V}^\top \boldsymbol{M} \boldsymbol{V}\).

References

[1]

Parlett, B. N., The Symmetric Eigenvalue Problem, SIAM, Philadelphia, 1998. https://doi.org/10.1137/1.9781611971163

pyfe3d.solver.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 PYFE3D_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.

Returns:
max_ncvint or None

The cap of ncv, at least 3, or None when ncv is not capped.

Raises:
ValueError

When the value given is an integer lower than 3 and not 0.

pyfe3d.solver.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().

Parameters:
kint

Number of eigenpairs.

nint

Size of the problem, the number of rows of the matrices.

Returns:
ncvint

Number of Lanczos vectors.

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

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

Solves \(\boldsymbol{A} \boldsymbol{x} = w \boldsymbol{M} \boldsymbol{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. The first pass is a call in the requested mode. The following passes use the shift-invert mode about the shift sigma_d given by _deflation_shift(), with the eigenvectors of the previous passes deflated from the operator, and each of them returns the next eigenvalues by increasing \(w\). Each pass factorizes A - sigma_d*M with SuperLU, and every call starts from start_kwargs(). 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.

Parameters:
Asparse matrix

Square matrix of the left-hand side, of size n x n.

kint

Number of eigenpairs.

Msparse matrix

Positive definite matrix of the right-hand side, of size n x n.

sigmafloat

Shift of the shift-invert or Cayley mode, in the units of \(w\).

whichstr

Eigenvalues to find, passed to scipy.sparse.linalg.eigsh(), e.g. 'LM' or 'SM'.

modestr, optional

'normal' (shift-invert, default) or 'cayley', passed to scipy.sparse.linalg.eigsh().

tolfloat, optional

Relative accuracy of the eigenvalues, passed to scipy.sparse.linalg.eigsh(); 0 (default) means the machine precision.

Returns:
w, xndarray

The eigenvalues and the \(\boldsymbol{M}\)-orthonormal eigenvectors, by increasing \(w\) when computed in passes. w is a 1-D array of size k and x an array of shape (n, k).

See also

arpack_max_ncv, arpack_ncv, start_kwargs, start_vector

Notes

In the shift-invert mode ARPACK computes the eigenvalues \(\nu = 1/(w - \sigma)\) of largest modulus of the operator \((\boldsymbol{A} - \sigma \boldsymbol{M})^{-1} \boldsymbol{M}\), i.e. the \(w\) nearest \(\sigma\) ([1], Chapter 3, Section 3.2, Eq. (3.2.2)). The passes keep k <= (max_ncv - 1)//2, i.e. ncv >= 2 k + 1. With k close to ncv ARPACK converges very slowly, not at all, or misses eigenvalues, e.g. k=31 and ncv=32.

The passes after the first one use the shift-invert mode about \(\sigma_d\), the middle of the largest gap between the last five eigenvalues found, with the eigenvectors \(\boldsymbol{V}\) of the previous passes deflated from the operator:

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

such that the eigenvalues of \(\boldsymbol{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. ARPACK applies the operator to \(\boldsymbol{M} \boldsymbol{x}\); for an eigenvector \(\boldsymbol{v}_i\) of \(\boldsymbol{V}\), \((\boldsymbol{A} - \sigma_d \boldsymbol{M})^{-1} \boldsymbol{M} \boldsymbol{v}_i = \boldsymbol{v}_i/(w_i - \sigma_d)\) and, \(\boldsymbol{V}\) being \(\boldsymbol{M}\)-orthonormal, \((\boldsymbol{I} - \boldsymbol{V} \boldsymbol{V}^\top \boldsymbol{M}) \boldsymbol{v}_i = \boldsymbol{0}\).

References

[1]

Lehoucq, R. B., Sorensen, D. C. and Yang, C., ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM, Philadelphia, 1998. https://doi.org/10.1137/1.9780898719628

pyfe3d.solver.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.

pyfe3d.solver.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.

pyfe3d.solver.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.

pyfe3d.solver.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.

pyfe3d.solver.arpackutils.start_kwargs(n)#

Arguments of scipy.sparse.linalg.eigs() and scipy.sparse.linalg.eigsh() such that the results are reproducible

v0 is start_vector() and, with SciPy 1.17 or newer, rng is a numpy.random.Generator with the seed 42, from which ARPACK draws a new starting vector whenever the Krylov subspace is invariant, e.g. when the operator is nilpotent.

Parameters:
nint

Size of the problem.

Returns:
kwargsdict

{'v0': start_vector(n)}, with the key 'rng' when supported.

Notes

Since SciPy 1.17, ARPACK asks the caller for a random vector when the Arnoldi or Lanczos process breaks down, i.e. finds an invariant subspace, and SciPy draws it from rng. With rng=None it comes from the entropy of the operating system: with the same v0 the results then differ from one call to the next. Older versions of SciPy draw it inside ARPACK, with a seed kept from one call to the next in the same process, which rng cannot fix.

pyfe3d.solver.arpackutils.start_vector(n)#

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

The vector is drawn from numpy.random.RandomState with the seed 42, so that the same n always gives the same vector. The calls of ARPACK use start_kwargs(), which also fixes the vectors drawn when ARPACK restarts.

Parameters:
nint

Size of the problem.

Returns:
v0ndarray

1-D array of size n.