Eigenvalue Solver

The jdsym Module

The jdsym module provides an implementation of the JDSYM algorithm, that is conveniently callable from Python. JDSYM is an eigenvalue solver to compute eigenpairs of a generalised matrix eigenvalue problem of the form

(1)\[\mathbf{A} \mathbf{x} = \lambda \mathbf{M} \mathbf{x}\]

or a standard eigenvalue problem of the form

(2)\[\mathbf{A} \mathbf{x} = \lambda \mathbf{x}\]

where \(\mathbf{A}\) is symmetric and \(\mathbf{M}\) is symmetric positive definite.

The module exports a single function:

jdsym.jdsym(A, M, K, kmax, tau, jdtol, itmax, linsolver, **kwargs)

Implements Jacobi-Davidson iterative method to identify a given number of eigenvalues near a target value.

Parameters :

A:the matrix \(\mathbf{A}\) in (1) or (2). A must provide the shape attribute and the matvec and matvec_transp methods. M:the matrix \(\mathbf{M}\) in (1). \(\mathbf{M}\) must provide the shape attribute and the matvec and matvec_transp methods. If the standard eigenvalue problem (2) is to be solved, M should be set to None. K:a preconditioner object that supplies the shape attribute and the precon method. If no preconditioner is to used, then the None value can be passed for this parameter. kmax:the number of eigenpairs to be computed. tau:the target value \(\tau\). Eigenvalues in the vicinity of \(\tau\) will be computed. jdtol:the convergence tolerance for eigenpairs \((\lambda,\mathbf{x})\). The converged eigenpairs have a residual \(\|\mathbf{A} \mathbf{x} - \lambda \mathbf{M} \mathbf{x}\|_2\) less than jdtol. itmax:an integer that specifies the maximum number of Jacobi-Davidson iterations to perform. linsolver:a function that implements an iterative method for solving linear systems of equations. The function linsolver is required to conform to the standards mentioned in Iterative Solvers.

Keywords :

jmax:the maximum dimension of the search subspace (default: 25). jmin:the dimension of the search subspace after a restart (default: 10). blksize:the block size used in the JDSYM algorithm (default: 1). blkwise:is an integer that affects the convergence criterion if blksize is larger than 1 (default: 0). V0:a NumPy array of rank one or two. It specifies the initial search subspace (default: a randomly generated initial search subspace). optype:is an integer specifying the operator type used in the correction equation. If optype=1, the non-symmetric version is used. If optype=2, the symmetric version is used (default: 2). linitmax:the maximum number of steps taken in the inner iteration (iterative linear solver) (default: 200). eps_tr:the tracking parameter (default: 1.0e-3). toldecay:is a float value that influences the dynamic adaptation of the stopping criterion of the inner iteration (default: 1.5). clvl:verbosity level. The higher the clvl parameter, the more output is sent to the standard output. clvl=0 produces no output (default: 0). strategy:is an integer specifying shifting and sorting strategy of JDSYM. strategy=0 enables the default JDSYM algorithm. strategy=1 enables JDSYM to avoid convergence to eigenvalues smaller than \(\tau\) (default: 0). projector:is used to keep the search subspace and the eigenvectors in a certain subspace. The parameter projector can be any Python object that has a shape attribute and a project method. The project method takes a vector (a rank-1 NumPy array) as its sole argument and projects that vector in-place. This parameter can be used to implement the DIRPROJ and SAUG methods (default: no projection).

Returns:

kconv:the number of converged eigenpairs. lambda:a rank-1 NumPy array containing the converged eigenvalues. Q:a rank-2 NumPy array containing the converged eigenvectors. The i-th eigenvector is accessed by Q[:,i]. it:an integer indicating the number of Jacobi-Davidson steps (outer iteration steps) performed.

Example: Maxwell Problem

Todo

Update the timings below.

Warning

The timings below are Roman’s old benchmarks. We should run them again.

The following code illustrates the use of the jdsym module. Two matrices \(\mathbf{A}\) and \(\mathbf{M}\) are read from files. A Jacobi preconditioner from \(\mathbf{A} - \tau\mathbf{M}\) is built. Then the JDSYM eigensolver is called, calculating 5 eigenvalues near 25.0 and the associated eigenvalues to an accuracy of \(10^{-10}\). We set strategy=1 to avoid convergence to the high-dimensional null space of (\(\mathbf{A}\), \(\mathbf{M}\)):

from pysparse.sparse import spmatrix
from pysparse.itsolvers import krylov
from pysparse.eig import jdsym
from pysparse.precon import precon

A = spmatrix.ll_mat_from_mtx('edge6x3x5_A.mtx')
M = spmatrix.ll_mat_from_mtx('edge6x3x5_B.mtx')
tau = 25.0

Atau = A.copy()
Atau.shift(-tau, M)
K = precon.jacobi(Atau)

A = A.to_sss(); M = M.to_sss()
k_conv, lmbd, Q, it  = jdsym.jdsym(A, M, K, 5, tau,
                                   1e-10, 150, krylov.qmrs,
                                   jmin=5, jmax=10, clvl=1, strategy=1)

This code takes 33.71 seconds to compute the five eigenpairs. A native C version, using the same computational kernels, takes 35.64 for the same task. We expected the Python version to be slower due to the overhead generated when calling the matrix-vector multiplication and the preconditioner, but surprisingly the Python code was even a bit faster.