Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update to current ecosystem #86

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 4 additions & 6 deletions krypy/deflation.py
Original file line number Diff line number Diff line change
Expand Up @@ -213,9 +213,7 @@ def estimate_time(self, nsteps, ndefl, deflweight=1.0):
"M": ndefl,
"Ml": ndefl,
"Mr": ndefl,
"ip_B": (
ndefl * (ndefl + 1) / 2 + ndefl ** 2 + 2 * ndefl * solver_ops["Ml"]
),
"ip_B": (ndefl * (ndefl + 1) / 2 + ndefl**2 + 2 * ndefl * solver_ops["Ml"]),
"axpy": (
ndefl * (ndefl + 1) / 2
+ ndefl * ndefl
Expand Down Expand Up @@ -670,7 +668,7 @@ def _auto():
)[0:-1]

def compute_pseudo(delta_log):
delta = 10 ** delta_log
delta = 10**delta_log
if ls_small.self_adjoint:
# pseudospectrum are intervals
pseudo_intervals = utils.Intervals(
Expand Down Expand Up @@ -701,7 +699,7 @@ def compute_pseudo(delta_log):
if pseudolen > 0:
polymax = numpy.max(numpy.abs(p(pseudo_path.vertices())))
else:
polymax = numpy.Inf
polymax = numpy.inf

# compute THE bound
return (
Expand Down Expand Up @@ -801,7 +799,7 @@ def __init__(self, deflated_solver, mode="ritz"):
self.values = numpy.zeros(m + n, dtype=sigmas.dtype)
zero = numpy.abs(sigmas) < numpy.finfo(float).eps
self.values[~zero] = 1.0 / sigmas[~zero]
self.values[zero] = numpy.Inf
self.values[zero] = numpy.inf
else:
raise utils.ArgumentError(
f"Invalid value '{mode}' for 'mode'. "
Expand Down
9 changes: 4 additions & 5 deletions krypy/linsys.py
Original file line number Diff line number Diff line change
Expand Up @@ -510,8 +510,7 @@ def operations(nsteps):
)

def _solve(self):
"""Abstract method that solves the linear system.
"""
"""Abstract method that solves the linear system."""
raise NotImplementedError(
"_solve has to be overridden by " "the derived solver class."
)
Expand Down Expand Up @@ -597,7 +596,7 @@ def _solve(self):
yk = numpy.zeros((N, 1), dtype=self.dtype)

# square of the old residual norm
self.rhos = rhos = [self.MMlr0_norm ** 2]
self.rhos = rhos = [self.MMlr0_norm**2]

# will be updated by _compute_rkn if explicit_residual is True
self.Mlrk = self.Mlr0.copy()
Expand Down Expand Up @@ -662,7 +661,7 @@ def _solve(self):

# compute norm and rho_new
MMlrk_norm = utils.norm(self.Mlrk, self.MMlrk, ip_B=self.linear_system.ip_B)
rhos.append(MMlrk_norm ** 2)
rhos.append(MMlrk_norm**2)

# compute Lanczos vector + new subdiagonal element
if self.store_arnoldi:
Expand All @@ -680,7 +679,7 @@ def _solve(self):
# update rho_new if it was updated in _compute_norms
if rkn is not None:
# new rho
rhos[-1] = rkn ** 2
rhos[-1] = rkn**2

self.iter += 1

Expand Down
4 changes: 2 additions & 2 deletions krypy/recycling/generators.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ def generate(self, ritz, remaining_subset):
class RitzSmall(_RitzSubsetsGenerator):
"""Successively returns the Ritz value of smallest magnitude."""

def __init__(self, max_vectors=numpy.Inf):
def __init__(self, max_vectors=numpy.inf):
self.max_vectors = max_vectors

def generate(self, ritz, remaining_subset):
Expand All @@ -34,7 +34,7 @@ class RitzExtremal(_RitzSubsetsGenerator):
smallest and largest magnitude are returned.
"""

def __init__(self, max_vectors=numpy.Inf):
def __init__(self, max_vectors=numpy.inf):
self.max_vectors = max_vectors

def generate(self, ritz, remaining_subset):
Expand Down
28 changes: 16 additions & 12 deletions krypy/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,6 @@
import scipy.linalg.blas as blas
from scipy.sparse import isspmatrix

# from scipy.sparse.linalg import LinearOperator, aslinearoperator
from scipy.sparse.sputils import isintlike

__all__ = [
"ArgumentError",
"AssumptionError",
Expand Down Expand Up @@ -143,6 +140,13 @@ def shape_vecs(*args):
return flat_vecs, ret_args


def isint(x):
try:
return int(x) == x
except Exception:
return False


def ip_euclid(X, Y):
"""Euclidean inner product.

Expand Down Expand Up @@ -233,7 +237,7 @@ def norm(x, y=None, ip_B=None):
raise InnerProductError(
"inner product defined by ip_B not positive "
"definite? ||diag(ip).imag||/||diag(ip)||="
f"{nrm_diag_imag/nrm_diag}"
f"{nrm_diag_imag / nrm_diag}"
)
return numpy.sqrt(numpy.linalg.norm(ip, 2))

Expand Down Expand Up @@ -355,7 +359,7 @@ def __init__(self, x):
alpha = 1 if gamma == 0 else gamma / xnorm
else:
sigma = numpy.linalg.norm(v[1:], 2)
xnorm = numpy.sqrt(numpy.abs(gamma) ** 2 + sigma ** 2)
xnorm = numpy.sqrt(numpy.abs(gamma) ** 2 + sigma**2)

# is x the multiple of first unit vector?
if sigma == 0:
Expand All @@ -372,7 +376,7 @@ def __init__(self, x):
alpha = -gamma / numpy.abs(gamma)

self.xnorm = xnorm
self.v = v / numpy.sqrt(numpy.abs(v[0]) ** 2 + sigma ** 2)
self.v = v / numpy.sqrt(numpy.abs(v[0]) ** 2 + sigma**2)
self.alpha = alpha
self.beta = beta

Expand Down Expand Up @@ -767,7 +771,7 @@ def angles(F, G, ip_B=None, compute_vectors=False):
else:
Y, s, Z = scipy.linalg.svd(inner(QF, QG, ip_B=ip_B))
Vcos = numpy.dot(QG, Z.T.conj())
n_large = numpy.flatnonzero((s ** 2) < 0.5).shape[0]
n_large = numpy.flatnonzero((s**2) < 0.5).shape[0]
n_small = s.shape[0] - n_large
theta = numpy.hstack(
[
Expand Down Expand Up @@ -1369,7 +1373,7 @@ class LinearOperator(object):
"""

def __init__(self, shape, dtype, dot=None, dot_adj=None):
if len(shape) != 2 or not isintlike(shape[0]) or not isintlike(shape[1]):
if len(shape) != 2 or not isint(shape[0]) or not isint(shape[1]):
raise LinearOperatorError("shape must be (m,n) with m and n " "integer")
self.shape = shape
self.dtype = numpy.dtype(dtype) # defaults to float64
Expand Down Expand Up @@ -1525,7 +1529,7 @@ def __init__(self, A, p):
raise LinearOperatorError("LinearOperator expected as A")
if A.shape[0] != A.shape[1]:
raise LinearOperatorError("square LinearOperator expected as A")
if not isintlike(p):
if not isint(p):
raise LinearOperatorError("integer expected as p")
self.args = (A, p)
super(_PowerLinearOperator, self).__init__(
Expand Down Expand Up @@ -1699,7 +1703,7 @@ def gap(lamda, sigma, mode="individual"):
# is a sigma value in lamda interval?
if not numpy.all(sigma_lo + sigma_hi):
return None
delta = numpy.Infinity
delta = numpy.inf
if numpy.any(sigma_lo):
delta = lamda_min - numpy.max(sigma[sigma_lo])
if numpy.any(sigma_hi):
Expand Down Expand Up @@ -1909,7 +1913,7 @@ def __init__(self, evals, exclude_zeros=False):

def eval_step(self, step):
"""Evaluate bound for given step."""
return 2 * self.base ** step
return 2 * self.base**step

def get_step(self, tol):
"""Return step at which bound falls below tolerance."""
Expand Down Expand Up @@ -1999,7 +2003,7 @@ def eval_step(self, step):
return 2 * self.base ** numpy.floor(step / 2.0)

def get_step(self, tol):
"""Return step at which bound falls below tolerance. """
"""Return step at which bound falls below tolerance."""
return 2 * numpy.log(tol / 2.0) / numpy.log(self.base)


Expand Down
2 changes: 1 addition & 1 deletion test/test_recycling.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ def test_RitzFactorySimple(Solver, which):
# ],
# )
# @pytest.mark.parametrize("self_adjoint", [True])
# @pytest.mark.parametrize("max_vectors", [numpy.Inf, 2])
# @pytest.mark.parametrize("max_vectors", [numpy.inf, 2])
# def test_SmallRitzGenerator(A, self_adjoint, max_vectors):
# ritz = _get_ritz(A, self_adjoint)
# small = krypy.recycling.generators.SmallRitz(max_vectors=max_vectors)
Expand Down
6 changes: 3 additions & 3 deletions test/test_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -504,7 +504,7 @@ def assert_arnoldi(
arnoldi_res = MAV - numpy.dot(V, H)
arnoldi_resn = krypy.utils.norm(arnoldi_res, ip_B=ip_B)
# inequality (2.3) in [1]
arnoldi_tol = arnoldi_const * k * (N ** 1.5) * eps * An
arnoldi_tol = arnoldi_const * k * (N**1.5) * eps * An
assert arnoldi_resn <= arnoldi_tol

# check orthogonality by measuring \| I - <V,V> \|_2
Expand All @@ -514,7 +514,7 @@ def assert_arnoldi(
ortho_res = numpy.eye(V.shape[1]) - krypy.utils.inner(V, V, ip_B=ip_B)
ortho_resn = numpy.linalg.norm(ortho_res, 2)
if ortho == "house":
ortho_tol = ortho_const * (k ** 1.5) * N * eps # inequality (2.4) in [1]
ortho_tol = ortho_const * (k**1.5) * N * eps # inequality (2.4) in [1]
else:
vAV_singvals = scipy.linalg.svd(
numpy.column_stack([V[:, [0]], (MAV[:, :-1] if invariant else MAV)]),
Expand All @@ -525,7 +525,7 @@ def assert_arnoldi(
else:
# inequality (2.5) in [1]
ortho_tol = (
ortho_const * (k ** 2) * N * eps * vAV_singvals[0] / vAV_singvals[-1]
ortho_const * (k**2) * N * eps * vAV_singvals[0] / vAV_singvals[-1]
)
# mgs or lanczos is not able to detect an invariant subspace reliably
if (ortho != "mgs" or N != k) and ortho != "lanczos":
Expand Down
Loading