Source code for sfepy.solvers.eigen

from __future__ import absolute_import

import numpy as nm
import scipy.sparse as sps

from sfepy.base.base import get_default, try_imports, Struct
from sfepy.base.timing import Timer
from sfepy.solvers.solvers import Solver, EigenvalueSolver
import six
from six.moves import range

[docs] def eig(mtx_a, mtx_b=None, n_eigs=None, eigenvectors=True, return_time=None, solver_kind='eig.scipy', **ckwargs): """ Utility function that constructs an eigenvalue solver given by `solver_kind`, calls it and returns its output. """ kwargs = {'name' : 'aux', 'kind' : solver_kind} kwargs.update(ckwargs) conf = Struct(**kwargs) solver = Solver.any_from_conf(conf) status = {} out = solver(mtx_a, mtx_b, n_eigs, eigenvectors, status) if return_time is not None: return_time[0] = status['time'] return out
[docs] def standard_call(call): """ Decorator handling argument preparation and timing for eigensolvers. """ def _standard_call(self, mtx_a, mtx_b=None, n_eigs=None, eigenvectors=None, status=None, conf=None, **kwargs): timer = Timer(start=True) conf = get_default(conf, self.conf) mtx_a = get_default(mtx_a, self.mtx_a) mtx_b = get_default(mtx_b, self.mtx_b) n_eigs = get_default(n_eigs, self.n_eigs) eigenvectors = get_default(eigenvectors, self.eigenvectors) status = get_default(status, self.status) if n_eigs == 0: result = self._ret_zero(mtx_a, eigenvectors=eigenvectors) else: result = call(self, mtx_a, mtx_b, n_eigs, eigenvectors, status, conf, **kwargs) elapsed = timer.stop() if status is not None: status['time'] = elapsed return result return _standard_call
[docs] class ScipyEigenvalueSolver(EigenvalueSolver): """ SciPy-based solver for both dense and sparse problems. The problem is consirered sparse if `n_eigs` argument is not None. """ name = 'eig.scipy' _parameters = [ ('method', "{'eig', 'eigh', 'eigs', 'eigsh'}", 'eigs', False, """The method for solving general or symmetric eigenvalue problems: for dense problems :func:`eig()` or :func:`eigh()` can be used, for sparse problems :func:`eigs()` or :func:`eigsh()` should be used."""), ('which', "'LM' | 'SM' | 'LR' | 'SR' | 'LI' | 'SI'", 'SM', False, """Which eigenvectors and eigenvalues to find, see :func:`scipy.sparse.linalg.eigs()` or :func:`scipy.sparse.linalg.eigsh()`. For dense problmes, only 'LM' and 'SM' can be used"""), ('linear_solver', "({'ls.cholesky', 'ls.mumps', ...}, ls_conf)", None, False, """The method used to construct an inverse linear operator. If None, the eigenvalue solver will solve the linear system internally."""), ('*', '*', None, False, 'Additional parameters supported by the method.'), ] def __init__(self, conf, **kwargs): EigenvalueSolver.__init__(self, conf, **kwargs) import scipy.linalg as sla self.sla = sla aux = try_imports(['import scipy.sparse.linalg as ssla'], 'cannot import scipy sparse eigenvalue solvers!') self.ssla = aux['ssla'] @standard_call def __call__(self, mtx_a, mtx_b=None, n_eigs=None, eigenvectors=None, status=None, conf=None): kwargs = self.build_solver_kwargs(conf) if conf.method in ('eig', 'eigh'): mtx_a, mtx_b = self._to_array(mtx_a, mtx_b) if conf.method == 'eig': out = self.sla.eig(mtx_a, mtx_b, right=eigenvectors, **kwargs) else: out = self.sla.eigh(mtx_a, mtx_b, eigvals_only=not eigenvectors, **kwargs) else: eig = self.ssla.eigs if conf.method == 'eigs' else self.ssla.eigsh if conf.linear_solver is not None: import sfepy.solvers.ls as ls from sfepy.solvers import solver_table import scipy.sparse as sps from sfepy.base.base import Struct ls_solvers = {} for k, v in solver_table.items(): if not k.startswith('ls'): continue aux = [par[0] == 'use_presolve' for par in v._parameters] if nm.any(aux): ls_solvers[k] = v fake_mtx_a = sps.csc_matrix(mtx_a.shape) solver_name, solver_conf = conf.linear_solver ls_conf = Struct(use_presolve=True, **solver_conf) ls = ls_solvers[solver_name](ls_conf) ls.mtx = mtx_a ls.presolve(mtx_a) matvec = ls.__call__ inv_op_a = self.ssla.LinearOperator(mtx_a.shape, matvec=matvec) out = eig(fake_mtx_a, M=mtx_b, k=n_eigs, which=conf.which, OPinv=inv_op_a, return_eigenvectors=eigenvectors, **kwargs) else: out = eig(mtx_a, M=mtx_b, k=n_eigs, which=conf.which, return_eigenvectors=eigenvectors, **kwargs) if eigenvectors: eigs = out[0] else: eigs = out if nm.iscomplexobj(eigs): ii = nm.argsort(nm.linalg.norm(eigs[:, None], axis=1)) else: ii = nm.argsort(eigs) if n_eigs is not None and (conf.method in ('eig', 'eigh')): if conf.which == 'SM': ii = ii[:n_eigs] elif conf.which == 'LM': ii = ii[:-n_eigs-1:-1] else: raise ValueError("only 'LM' or 'SM' can be used with dense" " problems! (%s)" % conf.which) if eigenvectors: mtx_ev = out[1][:, ii] out = (eigs[ii], mtx_ev) else: out = eigs[ii] return out
[docs] class ScipySGEigenvalueSolver(EigenvalueSolver): """ SciPy-based solver for dense symmetric problems. """ name = 'eig.sgscipy' def __init__(self, conf, **kwargs): EigenvalueSolver.__init__(self, conf, **kwargs) try: import scipy.linalg.lapack as llapack except ImportError: import scipy.lib.lapack as llapack self.llapack = llapack @standard_call def __call__(self, mtx_a, mtx_b=None, n_eigs=None, eigenvectors=None, status=None, conf=None): """ Notes ----- Eigenvectors argument is ignored, as they are computed always. """ ll = self.llapack mtx_a, mtx_b = self._to_array(mtx_a, mtx_b) if nm.iscomplexobj(mtx_a): if mtx_b is None: fun = ll.get_lapack_funcs(['heev'], arrays=(mtx_a,))[0] else: fun = ll.get_lapack_funcs(['hegv'], arrays=(mtx_a,))[0] else: if mtx_b is None: fun = ll.get_lapack_funcs(['syev'], arrays=(mtx_a,))[0] else: fun = ll.get_lapack_funcs(['sygv'], arrays=(mtx_a,))[0] if mtx_b is None: out = fun(mtx_a) else: out = fun(mtx_a, mtx_b) # Fix output order of scipy.linalg.lapack functions. if out[0].ndim == 2: out = (out[1], out[0]) + out[2:] if not eigenvectors: if n_eigs is None: out = out[0] else: out = out[0][:n_eigs] else: if n_eigs is None: out = out[:-1] else: out = (out[0][:n_eigs], out[1][:, :n_eigs]) return out
[docs] class LOBPCGEigenvalueSolver(EigenvalueSolver): """ SciPy-based LOBPCG solver for sparse symmetric problems. """ name = 'eig.scipy_lobpcg' _parameters = [ ('i_max', 'int', 20, False, 'The maximum number of iterations.'), ('eps_a', 'float', None, False, 'The absolute tolerance for the convergence.'), ('largest', 'bool', True, False, 'If True, solve for the largest eigenvalues, otherwise the smallest.'), ('precond', '{dense matrix, sparse matrix, LinearOperator}', None, False, 'The preconditioner.'), ] def __init__(self, conf, **kwargs): EigenvalueSolver.__init__(self, conf, **kwargs) from scipy.sparse.linalg import lobpcg self.lobpcg = lobpcg @standard_call def __call__(self, mtx_a, mtx_b=None, n_eigs=None, eigenvectors=None, status=None, conf=None): if n_eigs is None: n_eigs = mtx_a.shape[0] else: n_eigs = min(n_eigs, mtx_a.shape[0]) rng = nm.random.default_rng(12345) x = rng.normal(size=(mtx_a.shape[0], n_eigs)) out = self.lobpcg(mtx_a, x, mtx_b, M=conf.precond, tol=conf.eps_a, maxiter=conf.i_max, largest=conf.largest, verbosityLevel=conf.verbose) if not eigenvectors: out = out[0] return out
[docs] def init_slepc_args(): try: import sys, slepc4py except ImportError: return argv = [arg for arg in sys.argv if arg not in ['-h', '--help']] slepc4py.init(argv)
[docs] class SLEPcEigenvalueSolver(EigenvalueSolver): """ General SLEPc eigenvalue problem solver. """ name = 'eig.slepc' _parameters = [ ('method', 'str', 'krylovschur', False, 'The actual solver to use.'), ('problem', 'str', 'gnhep', False, """The problem type: Hermitian (hep), non-Hermitian (nhep), generalized Hermitian (ghep), generalized non-Hermitian (gnhep), generalized non-Hermitian with positive semi-definite B (pgnhep), and generalized Hermitian-indefinite (ghiep)."""), ('i_max', 'int', 20, False, 'The maximum number of iterations.'), ('eps', 'float', None, False, 'The convergence tolerance.'), ('conv_test', '{"abs", "rel", "norm", "user"}, ', 'abs', False, 'The type of convergence test.'), ('which', """{'largest_magnitude', 'smallest_magnitude', 'largest_real', 'smallest_real', 'largest_imaginary', 'smallest_imaginary', 'target_magnitude', 'target_real', 'target_imaginary', 'all', 'which_user'}""", 'largest_magnitude', False, 'Which eigenvectors and eigenvalues to find.'), ('*', '*', None, False, 'Additional parameters supported by the method.'), ] def __init__(self, conf, comm=None, context=None, **kwargs): if comm is None: init_slepc_args() from petsc4py import PETSc as petsc from slepc4py import SLEPc as slepc EigenvalueSolver.__init__(self, conf, petsc=petsc, slepc=slepc, comm=comm, context=context, **kwargs)
[docs] def create_eps(self, options=None, comm=None): optDB = self.petsc.Options() if options is not None: for key, val in six.iteritems(options): optDB[key] = val es = self.slepc.EPS() es.create(comm) return es
[docs] def create_petsc_matrix(self, mtx, comm=None): if mtx is None or isinstance(mtx, self.petsc.Mat): pmtx = mtx else: mtx = sps.csr_matrix(mtx) pmtx = self.petsc.Mat() pmtx.createAIJ(mtx.shape, csr=(mtx.indptr, mtx.indices, mtx.data), comm=comm) return pmtx
@standard_call def __call__(self, mtx_a, mtx_b=None, n_eigs=None, eigenvectors=None, status=None, conf=None, comm=None, context=None): solver_kwargs = self.build_solver_kwargs(conf) pmtx_a = self.create_petsc_matrix(mtx_a, comm=comm) pmtx_b = self.create_petsc_matrix(mtx_b, comm=comm) es = self.create_eps(options=solver_kwargs, comm=comm) es.setType(conf.method) es.setProblemType(getattr(es.ProblemType, conf.problem.upper())) es.setDimensions(nev=n_eigs) es.setTolerances(tol=conf.eps, max_it=conf.i_max) es.setOperators(pmtx_a, pmtx_b) es.setConvergenceTest(getattr(es.Conv, conf.conv_test.upper())) es.setWhichEigenpairs(getattr(es.Which, conf.which.upper())) es.setFromOptions() es.solve() n_converged = es.getConverged() if status is not None: status['n_iter'] = es.getIterationNumber() status['n_converged'] = n_converged if not eigenvectors: out = nm.array([es.getEigenvalue(ii) for ii in range(n_converged)]) else: vr, vi = pmtx_a.createVecs() eigs = [] vrs, vis = [], [] is_real = True for ii in range(n_converged): val = es.getEigenpair(ii, vr, vi) eigs.append(val if val.imag != 0 else val.real) vrs.append(vr.getArray().copy()) vis.append(vi.getArray().copy()) if is_real and nm.sum(nm.abs(vis[-1])) > 0.0: is_real = False eigs = nm.array(eigs) vecs = nm.array(vrs).T if not is_real: vecs += 1j * nm.array(vis) out = (eigs, vecs) return out
[docs] class MatlabEigenvalueSolver(EigenvalueSolver): """ Matlab eigenvalue problem solver. """ name = 'eig.matlab' _parameters = [ ('method', """{'eig', 'eigs', None}""", 'eigs', False, """The solution method. Note that eig() function cannot be used for all inputs. If `n_eigs` is not None, eigs() is used regardless of this parameter."""), ('balance', """{'balance', 'nobalance'}""", 'balance', False, 'The balance option for eig().'), ('algorithm', """{'chol', 'qz'}""", 'chol', False, 'The algorithm option for eig().'), ('which', """{'lm', 'sm', 'la', 'sa', 'be' 'lr', 'sr', 'li', 'si', sigma}""", 'lm', False, 'Which eigenvectors and eigenvalues to find with eigs().'), ('*', '*', None, False, 'Additional parameters supported by eigs().'), ] def __init__(self, conf, comm=None, context=None, **kwargs): import matlab.engine as me EigenvalueSolver.__init__(self, conf, me=me, context=context, **kwargs)
[docs] def solver_call(self, mtx_filename, eigs_filename, which=None): import os import scipy.io as sio eng = self.me.start_matlab() eng.cd(os.path.dirname(__file__)) eng.matlab_eig(mtx_filename, eigs_filename) eng.quit() return sio.loadmat(eigs_filename)
@standard_call def __call__(self, mtx_a, mtx_b=None, n_eigs=None, eigenvectors=None, status=None, conf=None, comm=None, context=None): import os import shutil import tempfile import scipy.io as sio solver_kwargs = self.build_solver_kwargs(conf) dirname = tempfile.mkdtemp() mtx_filename = os.path.join(dirname, 'matrices.mat') eigs_filename = os.path.join(dirname, 'eigs.mat') sio.savemat(mtx_filename, { 'A' : mtx_a, 'B' : mtx_b if mtx_b is not None else 'None', 'n_eigs' : n_eigs if n_eigs is not None else 'None', 'eigenvectors' : (eigenvectors if eigenvectors is not None else False), 'method' : conf.method, 'balance' : conf.balance, 'algorithm' : conf.algorithm, 'which' : conf.which, 'verbose' : conf.verbose, 'eigs_options' : solver_kwargs, }) evp = self.solver_call(mtx_filename, eigs_filename, conf.which) shutil.rmtree(dirname) out = evp['vals'][:, 0] if eigenvectors: out = (out, evp['vecs']) return out
[docs] class OctaveEigenvalueSolver(MatlabEigenvalueSolver): """ Octave eigenvalue problem solver. """ name = 'eig.octave' def __init__(self, conf, comm=None, context=None, **kwargs): from oct2py import octave as oe EigenvalueSolver.__init__(self, conf, oe=oe, context=context, **kwargs)
[docs] def solver_call(self, mtx_filename, eigs_filename, which): import os import scipy.io as sio self.oe.addpath(os.path.dirname(__file__)) self.oe.matlab_eig(mtx_filename, eigs_filename) self.oe.exit() evp = sio.loadmat(eigs_filename) evals = evp['vals'] which = which.lower() sort_funs = { 'lm': (nm.abs, -1), 'sm': (nm.abs, 1), 'la': (nm.real, -1), 'sa': (nm.real, 1), 'lr': (nm.real, -1), 'sr': (nm.real, 1), 'li': (nm.imag, -1), 'si': (nm.imag, 1), 'be': (nm.real, 1), } sfun, order = sort_funs[which] idxs = nm.argsort(sfun(evals), axis=0) idxs = idxs.ravel()[::order] if which == 'be': idxs = nm.hstack([idxs[::-1], idxs]) k2 = len(idxs) // 2 idxs = idxs[k2:(k2 + len(idxs))] out = {'vals': evals[idxs, :]} if 'vecs' in evp: out['vecs'] = evp['vecs'][:, idxs] return out
[docs] class PrimmeEigenvalueSolver(EigenvalueSolver): """ PRIMME eigenvalue problem solver. https://github.com/primme/primme Installation: pip install primme """ name = 'eig.primme' _parameters = [ ('which', "{'LM', 'SM', 'LA', 'SA', 'CLT', 'CGT'}", 'LM', False, 'Which eigenvectors and eigenvalues to find.'), ('sigma', 'float', None, False, 'Find eigenvalues near sigma.'), ('maxiter', 'int', None, False, 'Maximum number of iterations.'), ('tol', 'float', 0, False, 'Tolerance for eigenpairs (stopping criterion).'), ('*', '*', None, False, 'Additional parameters supported by eigsh().'), ] def __init__(self, conf, comm=None, context=None, **kwargs): import primme EigenvalueSolver.__init__(self, conf, primme=primme, context=context, **kwargs) @standard_call def __call__(self, mtx_a, mtx_b=None, n_eigs=None, eigenvectors=None, status=None, conf=None, comm=None, context=None): if n_eigs is None: n_eigs = mtx_a.shape[0] solver_kwargs = self.build_solver_kwargs(conf) out = self.primme.eigsh(mtx_a, n_eigs, M=mtx_b, which=conf.which.upper(), tol=conf.tol, maxiter=conf.maxiter, sigma=conf.sigma, return_eigenvectors=eigenvectors, **solver_kwargs) return out