From 31dc7ac9c96ef265f4d22fd6e66f91283432a40a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Wed, 17 Jan 2018 19:29:57 +0100 Subject: [PATCH 01/11] utils: fix int/float warning in NormalizedRootsPolynomial --- krypy/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/krypy/utils.py b/krypy/utils.py index 5a73cde..b7d56be 100644 --- a/krypy/utils.py +++ b/krypy/utils.py @@ -2114,7 +2114,7 @@ def __call__(self, points): for j in range(vals.shape[1]): sort_tmp = numpy.argsort(numpy.abs(vals[:, j])) sort = numpy.zeros((n,), dtype=numpy.int) - mid = numpy.ceil(float(n)/2) + mid = int(numpy.ceil(float(n)/2)) sort[::2] = sort_tmp[:mid] sort[1::2] = sort_tmp[mid:][::-1] vals[:, j] = vals[sort, j] From d9bfb5871a8d2aea4c102100696be578c698a46f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Wed, 17 Jan 2018 19:41:11 +0100 Subject: [PATCH 02/11] utils: fix "object() takes no parameters" error in BoundMinres --- krypy/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/krypy/utils.py b/krypy/utils.py index b7d56be..ab7517d 100644 --- a/krypy/utils.py +++ b/krypy/utils.py @@ -1989,7 +1989,7 @@ def __new__(cls, evals): pos = True if pos: return BoundCG(evals) - return super(BoundMinres, cls).__new__(cls, evals) + return super(BoundMinres, cls).__new__(cls) def __init__(self, evals): '''Initialize with array/list of eigenvalues or Intervals object.''' From 8ce56ed27d2be36ac992447bba53e846ba8ba3ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Wed, 17 Jan 2018 19:44:27 +0100 Subject: [PATCH 03/11] utils: fix dead code in norm() fixes #61 --- krypy/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/krypy/utils.py b/krypy/utils.py index ab7517d..4539f2f 100644 --- a/krypy/utils.py +++ b/krypy/utils.py @@ -188,12 +188,12 @@ def norm(x, y=None, ip_B=None): Compute :math:`\sqrt{\langle x,y\rangle}` where the inner product is defined via ``ip_B``. ''' - if y is None: - y = x # Euclidean inner product? if y is None and (ip_B is None or isinstance(ip_B, IdentityLinearOperator)): return numpy.linalg.norm(x, 2) + if y is None: + y = x ip = inner(x, y, ip_B=ip_B) nrm_diag = numpy.linalg.norm(numpy.diag(ip), 2) nrm_diag_imag = numpy.linalg.norm(numpy.imag(numpy.diag(ip)), 2) From 9da95d1860679d6a8053426ae8beba6c0bb0e435 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Wed, 17 Jan 2018 19:55:09 +0100 Subject: [PATCH 04/11] travis: install numpy and scipy via addons --- .travis.yml | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/.travis.yml b/.travis.yml index 49e548f..423b413 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,8 +4,11 @@ python: - "3.2" virtualenv: system_site_packages: true -before_install: - - sudo add-apt-repository -y ppa:pylab/stable - - sudo apt-get update -qq - - sudo apt-get install -qq python-numpy python-scipy python3-numpy python3-scipy +addons: + apt: + packages: + - python-numpy + - python-scipy + - python3-numpy + - python3-scipy script: nosetests krypy/tests From 420b86c1c7b7fc0770db944e52897c7f1650ede4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Wed, 17 Jan 2018 20:07:49 +0100 Subject: [PATCH 05/11] travis: use python 3.6 --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 423b413..2d79a6d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,7 +1,7 @@ language: python python: - "2.7" - - "3.2" + - "3.6" virtualenv: system_site_packages: true addons: From e1d2390ce40daa825b4038f92f57fd4ccf0c64f8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Wed, 17 Jan 2018 20:18:31 +0100 Subject: [PATCH 06/11] travis: use non-sudo setup --- .travis.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index 2d79a6d..dfb8e1f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,9 +1,8 @@ +sudo: false language: python python: - "2.7" - "3.6" -virtualenv: - system_site_packages: true addons: apt: packages: From 8caa646dd1873bbeb4ecbf43ade557b9b897e2c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Wed, 17 Jan 2018 20:30:41 +0100 Subject: [PATCH 07/11] travis: install dependencies via pip and cache them --- .travis.yml | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/.travis.yml b/.travis.yml index dfb8e1f..17684b8 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,13 +1,9 @@ sudo: false language: python +cache: pip python: - "2.7" - "3.6" -addons: - apt: - packages: - - python-numpy - - python-scipy - - python3-numpy - - python3-scipy +install: + - travis_retry pip install -r requirements.txt script: nosetests krypy/tests From d8fd40a8c5135b4bc3dfb62631b67b0f53a619c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Wed, 17 Jan 2018 21:03:06 +0100 Subject: [PATCH 08/11] utils: remove hessenberg() --- krypy/utils.py | 93 -------------------------------------------------- 1 file changed, 93 deletions(-) diff --git a/krypy/utils.py b/krypy/utils.py index 4539f2f..c5ede84 100644 --- a/krypy/utils.py +++ b/krypy/utils.py @@ -365,99 +365,6 @@ def matrix(self): return numpy.eye(n, n) - self.beta * numpy.dot(self.v, self.v.T.conj()) -def hessenberg(a, calc_q=False, overwrite_a=False, check_finite=True): - """ - Compute Hessenberg form of a matrix. - - The Hessenberg decomposition is:: - - A = Q H Q^H - - where `Q` is unitary/orthogonal and `H` has only zero elements below - the first sub-diagonal. - - Parameters - ---------- - a : (M, M) array_like - Matrix to bring into Hessenberg form. - calc_q : bool, optional - Whether to compute the transformation matrix. Default is False. - overwrite_a : bool, optional - Whether to overwrite `a`; may improve performance. - Default is False. - check_finite : boolean, optional - Whether to check that the input matrix contains only finite numbers. - Disabling may give a performance gain, but may result in problems - (crashes, non-termination) if the inputs do contain infinities or NaNs. - - Returns - ------- - H : (M, M) ndarray - Hessenberg form of `a`. - Q : (M, M) ndarray - Unitary/orthogonal similarity transformation matrix ``A = Q H Q^H``. - Only returned if ``calc_q=True``. - - """ - from numpy import asarray_chkfinite, asarray - from scipy.linalg import calc_lwork - from scipy.linalg.misc import _datacopied - from scipy.linalg.lapack import get_lapack_funcs - from scipy.linalg.blas import get_blas_funcs - - if check_finite: - a1 = asarray_chkfinite(a) - else: - a1 = asarray(a) - if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]): - raise ValueError('expected square matrix') - overwrite_a = overwrite_a or (_datacopied(a1, a)) - - # if 2x2 or smaller: already in Hessenberg form - if a1.shape[0] <= 2: - if calc_q: - return a1, numpy.eye(a1.shape[0]) - return a1 - - gehrd, gebal = get_lapack_funcs(('gehrd', 'gebal'), (a1,)) - ba, lo, hi, pivscale, info = gebal(a1, permute=0, overwrite_a=overwrite_a) - if info < 0: - raise ValueError('illegal value in %d-th argument of internal gebal ' - '(hessenberg)' % -info) - n = len(a1) - lwork = calc_lwork.gehrd(gehrd.typecode, n, lo, hi) - hq, tau, info = gehrd(ba, lo=lo, hi=hi, lwork=lwork, overwrite_a=1) - if info < 0: - raise ValueError('illegal value in %d-th argument of internal gehrd ' - '(hessenberg)' % -info) - - if not calc_q: - for i in range(lo, hi): - hq[i+2:hi+1, i] = 0.0 - return hq - - # XXX: Use ORGHR routines to compute q. - typecode = hq.dtype - ger, gemm = get_blas_funcs(('ger', 'gemm'), dtype=typecode) - q = None - for i in range(lo, hi): - if tau[i] == 0.0: - continue - v = numpy.zeros(n, dtype=typecode) - v[i+1] = 1.0 - v[i+2:hi+1] = hq[i+2:hi+1, i] - hq[i+2:hi+1, i] = 0.0 - h = ger(-tau[i], v, v, - a=numpy.diag(numpy.ones(n, dtype=typecode)), overwrite_a=1) - if q is None: - q = h - else: - q = gemm(1.0, q, h) - if q is None: - q = numpy.diag(numpy.ones(n, dtype=typecode)) - return hq, q - - class Givens: def __init__(self, x): """Compute Givens rotation for provided vector x. From 060c6acccc25d92a47527a3cadba7d327ba8342c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Wed, 17 Jan 2018 21:03:16 +0100 Subject: [PATCH 09/11] deflation: use scipy's hessenberg() --- krypy/deflation.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/krypy/deflation.py b/krypy/deflation.py index a039cc1..ad8e462 100644 --- a/krypy/deflation.py +++ b/krypy/deflation.py @@ -398,10 +398,8 @@ def get(self, Wt, full=False): # Arnoldify WtoQ = Q.apply(Wto.T.conj()).T.conj() - # TODO: replace with scipy's Hessenberg when bug is fixed - #from scipy.linalg import hessenberg - #Hh, T = hessenberg( - Hh, T = utils.hessenberg( + from scipy.linalg import hessenberg + Hh, T = hessenberg( Q.apply(Wto.T.conj().dot(self.J).dot(Pt*(self.L.dot(WtoQ)))), calc_q=True) From 060ee887d515974c01c706b8ff03e6c1fb23e040 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Thu, 18 Jan 2018 11:12:05 +0100 Subject: [PATCH 10/11] setup: update dependency versions --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 0c01e4a..f00a0c0 100644 --- a/setup.py +++ b/setup.py @@ -29,7 +29,7 @@ def read(fname): author='André Gaul', author_email='gaul@web-yard.de', url='https://github.com/andrenarchy/krypy', - install_requires=['numpy (>=1.7)', 'scipy (>=0.13)'], + install_requires=['numpy (>=1.11)', 'scipy (>=0.17)'], classifiers=[ 'Development Status :: 4 - Beta', 'Intended Audience :: Science/Research', From 9bc13396a35f38c306d7f0017155172031366c8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Gaul?= Date: Thu, 18 Jan 2018 12:09:41 +0100 Subject: [PATCH 11/11] docs: update bibliography and intro --- docs/bibliography.rst | 4 ++++ docs/index.rst | 19 ++++++++++++++++--- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/docs/bibliography.rst b/docs/bibliography.rst index 5a1e6a0..1d72299 100644 --- a/docs/bibliography.rst +++ b/docs/bibliography.rst @@ -1,6 +1,10 @@ Bibliography ============ +.. [Gau14] A. Gaul. Recycling Krylov subspace methods for sequences of + linear systems. PhD thesis. TU Berlin, 2014. + https://dx.doi.org/10.14279/depositonce-4147 + .. [GauGLN13] A. Gaul, M. H. Gutknecht, J. Liesen, and R. Nabben. A framework for deflated and augmented Krylov subspace methods. SIAM J. Matrix Anal. Appl., 34 (2013), pp. 495-518. diff --git a/docs/index.rst b/docs/index.rst index 0c9d36b..365526a 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -12,7 +12,20 @@ What is KryPy and where is the code? KryPy is a Krylov subspace methods package for Python. If you're looking for the source code or bug reports, take a look at `KryPy's github page `_. These pages provide the -documentation of KryPy's API. +documentation of KryPy's API. The project was initiated by André Gaul while +researching Krylov subspace methods. The theoretical background as well as +applications of this software package can be found in the PhD thesis [Gau14]_. + +KryPy allows Python users to easily use Krylov subspace methods, e.g., for solving +linear systems or eigenvalue problems. With its built-in deflation and +recycling capabilities it is suitable for advanced applications of Krylov +subspace methods (see :doc:`krypy.deflation` and :doc:`krypy.recycling`). +It is also ideal for experimenting with Krylov subspaces since you have access to +all data that is generated (e.g., Arnoldi/Lanczos relations), you can use +different orthogonalization algorithms (Lanczos short recurrences, modified +Gram-Schmidt, double modified Gram-Schmidt, Householder), compare subspaces via +angles, and much more. And if you need more: KryPy is free software, it's easy +to extend, and pull requests are more than welcome! Contents -------- @@ -64,10 +77,10 @@ strategy. Just for illustration, the same linear system is solved twice in the following code:: from krypy.recycling import RecyclingMinres - + # get recycling solver with approximate Krylov subspace strategy rminres = RecyclingMinres(vector_factory='RitzApproxKrylov') - + # solve twice rsolver1 = rminres.solve(linear_system) rsolver2 = rminres.solve(linear_system)