# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
from pymor.core.config import config
if config.HAVE_PYAMG:
import numpy as np
import pyamg
from pymor.algorithms.genericsolvers import _parse_options
from pymor.core.defaults import defaults
from pymor.core.exceptions import InversionError
from pymor.operators.numpy import NumpyMatrixOperator
[docs] @defaults('tol', 'maxiter', 'verb', 'rs_strength', 'rs_CF',
'rs_postsmoother', 'rs_max_levels', 'rs_max_coarse', 'rs_coarse_solver',
'rs_cycle', 'rs_accel', 'rs_tol', 'rs_maxiter',
'sa_symmetry', 'sa_strength', 'sa_aggregate', 'sa_smooth',
'sa_presmoother', 'sa_postsmoother', 'sa_improve_candidates', 'sa_max_levels',
'sa_max_coarse', 'sa_diagonal_dominance', 'sa_coarse_solver', 'sa_cycle',
'sa_accel', 'sa_tol', 'sa_maxiter')
def solver_options(tol=1e-5,
maxiter=400,
verb=False,
rs_strength=('classical', {'theta': 0.25}),
rs_CF='RS',
rs_presmoother=('gauss_seidel', {'sweep': 'symmetric'}),
rs_postsmoother=('gauss_seidel', {'sweep': 'symmetric'}),
rs_max_levels=10,
rs_max_coarse=500,
rs_coarse_solver='pinv2',
rs_cycle='V',
rs_accel=None,
rs_tol=1e-5,
rs_maxiter=100,
sa_symmetry='hermitian',
sa_strength='symmetric',
sa_aggregate='standard',
sa_smooth=('jacobi', {'omega': 4.0/3.0}),
sa_presmoother=('block_gauss_seidel', {'sweep': 'symmetric'}),
sa_postsmoother=('block_gauss_seidel', {'sweep': 'symmetric'}),
sa_improve_candidates=(('block_gauss_seidel', {'sweep': 'symmetric', 'iterations': 4}), None),
sa_max_levels=10,
sa_max_coarse=500,
sa_diagonal_dominance=False,
sa_coarse_solver='pinv2',
sa_cycle='V',
sa_accel=None,
sa_tol=1e-5,
sa_maxiter=100):
"""Returns available solvers with default |solver_options| for the PyAMG backend.
Parameters
----------
tol
Tolerance for `PyAMG <http://pyamg.github.io/>`_ blackbox solver.
maxiter
Maximum iterations for `PyAMG <http://pyamg.github.io/>`_ blackbox solver.
verb
Verbosity flag for `PyAMG <http://pyamg.github.io/>`_ blackbox solver.
rs_strength
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_CF
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_presmoother
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_postsmoother
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_max_levels
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_max_coarse
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_coarse_solver
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_cycle
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_accel
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_tol
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
rs_maxiter
Parameter for `PyAMG <http://pyamg.github.io/>`_ Ruge-Stuben solver.
sa_symmetry
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_strength
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_aggregate
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_smooth
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_presmoother
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_postsmoother
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_improve_candidates
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_max_levels
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_max_coarse
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_diagonal_dominance
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_coarse_solver
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_cycle
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_accel
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_tol
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
sa_maxiter
Parameter for `PyAMG <http://pyamg.github.io/>`_ Smoothed-Aggregation solver.
Returns
-------
A dict of available solvers with default |solver_options|.
"""
return {'pyamg_solve': {'type': 'pyamg_solve',
'tol': tol,
'maxiter': maxiter},
'pyamg_rs': {'type': 'pyamg_rs',
'strength': rs_strength,
'CF': rs_CF,
'presmoother': rs_presmoother,
'postsmoother': rs_postsmoother,
'max_levels': rs_max_levels,
'max_coarse': rs_max_coarse,
'coarse_solver': rs_coarse_solver,
'cycle': rs_cycle,
'accel': rs_accel,
'tol': rs_tol,
'maxiter': rs_maxiter},
'pyamg_sa': {'type': 'pyamg_sa',
'symmetry': sa_symmetry,
'strength': sa_strength,
'aggregate': sa_aggregate,
'smooth': sa_smooth,
'presmoother': sa_presmoother,
'postsmoother': sa_postsmoother,
'improve_candidates': sa_improve_candidates,
'max_levels': sa_max_levels,
'max_coarse': sa_max_coarse,
'diagonal_dominance': sa_diagonal_dominance,
'coarse_solver': sa_coarse_solver,
'cycle': sa_cycle,
'accel': sa_accel,
'tol': sa_tol,
'maxiter': sa_maxiter}}
[docs] @defaults('check_finite', 'default_solver')
def apply_inverse(op, V, initial_guess=None, options=None, least_squares=False,
check_finite=True, default_solver='pyamg_solve'):
"""Solve linear equation system.
Applies the inverse of `op` to the vectors in `V` using PyAMG.
Parameters
----------
op
The linear, non-parametric |Operator| to invert.
V
|VectorArray| of right-hand sides for the equation system.
initial_guess
|VectorArray| with the same length as `V` containing initial guesses
for the solution. Some implementations of `apply_inverse` may
ignore this parameter. If `None` a solver-dependent default is used.
options
The |solver_options| to use (see :func:`solver_options`).
least_squares
Must be `False`.
check_finite
Test if solution only contains finite values.
default_solver
Default solver to use (pyamg_solve, pyamg_rs, pyamg_sa).
Returns
-------
|VectorArray| of the solution vectors.
"""
assert V in op.range
assert initial_guess is None or initial_guess in op.source and len(initial_guess) == len(V)
if least_squares:
raise NotImplementedError
if isinstance(op, NumpyMatrixOperator):
matrix = op.matrix
else:
from pymor.algorithms.to_matrix import to_matrix
matrix = to_matrix(op)
options = _parse_options(options, solver_options(), default_solver, None, least_squares)
V = V.to_numpy()
promoted_type = np.promote_types(matrix.dtype, V.dtype)
R = np.empty((len(V), matrix.shape[1]), dtype=promoted_type)
if options['type'] == 'pyamg_solve':
if len(V) > 0:
V_iter = iter(enumerate(V))
R[0], ml = pyamg.solve(matrix, next(V_iter)[1],
tol=options['tol'],
maxiter=options['maxiter'],
return_solver=True)
for i, VV in V_iter:
R[i] = pyamg.solve(matrix, VV,
tol=options['tol'],
maxiter=options['maxiter'],
existing_solver=ml)
elif options['type'] == 'pyamg_rs':
ml = pyamg.ruge_stuben_solver(matrix,
strength=options['strength'],
CF=options['CF'],
presmoother=options['presmoother'],
postsmoother=options['postsmoother'],
max_levels=options['max_levels'],
max_coarse=options['max_coarse'],
coarse_solver=options['coarse_solver'])
for i, VV in enumerate(V):
R[i] = ml.solve(VV,
tol=options['tol'],
maxiter=options['maxiter'],
cycle=options['cycle'],
accel=options['accel'])
elif options['type'] == 'pyamg_sa':
ml = pyamg.smoothed_aggregation_solver(matrix,
symmetry=options['symmetry'],
strength=options['strength'],
aggregate=options['aggregate'],
smooth=options['smooth'],
presmoother=options['presmoother'],
postsmoother=options['postsmoother'],
improve_candidates=options['improve_candidates'],
max_levels=options['max_levels'],
max_coarse=options['max_coarse'],
diagonal_dominance=options['diagonal_dominance'])
for i, VV in enumerate(V):
R[i] = ml.solve(VV,
tol=options['tol'],
maxiter=options['maxiter'],
cycle=options['cycle'],
accel=options['accel'])
else:
raise ValueError('Unknown solver type')
if check_finite:
if not np.isfinite(np.sum(R)):
raise InversionError('Result contains non-finite values')
return op.source.from_numpy(R)