Skip to content

Commit

Permalink
Utils: PEP8 changes.
Browse files Browse the repository at this point in the history
  • Loading branch information
dpo committed Oct 5, 2015
1 parent 4ca77bc commit 572898d
Showing 1 changed file with 18 additions and 13 deletions.
31 changes: 18 additions & 13 deletions pykrylov/tools/utils.py
Original file line number Diff line number Diff line change
@@ -1,22 +1,26 @@
# Various utilities.
"""Utilities related to linear operators."""

import numpy as np
from math import copysign, sqrt


def machine_epsilon():
"""Return the machine epsilon in double precision."""
return np.finfo(np.double).eps


def roots_quadratic(q2, q1, q0, tol=1.0e-8, nitref=1):
"""
"""Find the real roots of a quadratic.
Find the real roots of the quadratic q(x) = q2 * x^2 + q1 * x + q0.
The numbers a0, a1 and a0 must be real.
This function is written after the GALAHAD function of the same name.
See http://galahad.rl.ac.uk.
"""
a2 = float(q2) ; a1 = float(q1) ; a0 = float(q0)
a2 = float(q2)
a1 = float(q1)
a0 = float(q0)

# Case of a linear function.
if a2 == 0.0:
Expand All @@ -26,25 +30,25 @@ def roots_quadratic(q2, q1, q0, tol=1.0e-8, nitref=1):
else:
return []
else:
roots = [-a0/a1]
roots = [-a0 / a1]
else:
# Case of a quadratic.
rhs = tol * a1 * a1
if abs(a0*a2) > rhs:
if abs(a0 * a2) > rhs:
rho = a1 * a1 - 4.0 * a2 * a0
if rho < 0.0:
return []
# There are two real roots.
d = -0.5 * (a1 + copysign(sqrt(rho), a1))
roots = [d/a2, a0/d]
roots = [d / a2, a0 / d]
else:
# Ill-conditioned quadratic.
roots = [-a1/a2, 0.0]
roots = [-a1 / a2, 0.0]

# Perform a few Newton iterations to improve accuracy.
new_roots = []
for root in roots:
for it in range(nitref):
for _ in range(nitref):
val = (a2 * root + a1) * root + a0
der = 2.0 * a2 * root + a1
if der == 0.0:
Expand All @@ -57,7 +61,8 @@ def roots_quadratic(q2, q1, q0, tol=1.0e-8, nitref=1):


def check_symmetric(op, repeats=10):
"""
"""Cheap symmetry check.
Cheap check that a linear operator is symmetric. Supply `op`, a callable
linear operator. A set of `repeats` random vectors will be generated.
This function returns `True` or `False`.
Expand All @@ -67,7 +72,7 @@ def check_symmetric(op, repeats=10):
return False
eps = machine_epsilon()
np.random.seed(1)
for k in xrange(repeats):
for _ in xrange(repeats):
x = np.random.random(n)
w = op * x
r = op * w
Expand All @@ -80,6 +85,6 @@ def check_symmetric(op, repeats=10):
return True


if __name__ == '__main__':
roots = roots_quadratic(2.0e+20, .1, -4)
print 'Received: ', roots
# if __name__ == '__main__':
# roots = roots_quadratic(2.0e+20, .1, -4)
# print 'Received: ', roots

0 comments on commit 572898d

Please sign in to comment.