Skip to content

Commit

Permalink
[pyamg] move pyamg glue code to bindings.pyamg
Browse files Browse the repository at this point in the history
  • Loading branch information
@sdrave
sdrave committed Jan 5, 2017
1 parent f480435 commit 57034d3
Show file tree
Hide file tree
Showing 2 changed files with 300 additions and 189 deletions.
273 changes: 273 additions & 0 deletions src/pymor/bindings/pyamg.py
View file
@@ -0,0 +1,273 @@
# -*- coding: utf-8 -*-
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2016 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:

from collections import OrderedDict

import numpy as np
import pyamg

from pymor.core.defaults import defaults
from pymor.core.exceptions import InversionError
from pymor.operators.numpy import NumpyMatrixOperator

@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',
sid_ignore=('verb',))
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 possible |solver_options| (with default values) 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 tuple of all possible |solver_options|.
"""

return OrderedDict(
(('pyamg', {'type': 'pyamg',
'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}))
)


@defaults('check_finite')
def apply_inverse(op, V, options=None, check_finite=True):
"""Solve linear equation system.
Applies the inverse of `op` to the row vectors in `V`.
See :func:`dense_options` for documentation of all possible options for
sparse matrices.
See :func:`sparse_options` for documentation of all possible options for
sparse matrices.
This method is called by :meth:`NumpyMatrixOperator.apply_inverse`
and usually should not be used directly.
Parameters
----------
op
The left-hand side of the linear equation systems to
solve as a |NumpyMatrixOperator|
V
2-dimensional |NumPy array| containing as row vectors
the right-hand sides of the linear equation systems to
solve.
options
The solver options to use. (See :func:`_options`.)
Returns
-------
|NumPy array| of the solution vectors.
"""

assert isinstance(op, NumpyMatrixOperator)
assert V in op.range
matrix = op._matrix
# TODO Use to_matrix(op) after sparse format issue has been resolved.
V = V.data

default_options = solver_options(matrix)

if options is None:
options = next(iter(default_options.values()))
elif isinstance(options, str):
if options == 'least_squares':
for k, v in default_options.items():
if k.startswith('least_squares'):
options = v
break
assert not isinstance(options, str)
else:
options = default_options[options]
else:
assert 'type' in options and options['type'] in default_options \
and options.keys() <= default_options[options['type']].keys()
user_options = options
options = default_options[user_options['type']]
options.update(user_options)

promoted_type = np.promote_types(matrix.dtype, V.dtype)
R = np.empty((len(V), matrix.shape[1]), dtype=promoted_type)

if options['type'] == 'pyamg':
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.make_array(R)

0 comments on commit 57034d3

Please sign in to comment.