Merge pull request #62 from andrenarchy/fix/61-dead-norm-code

Fix dead code in norm()

.travis.yml

@@ -1,11 +1,9 @@
+sudo: false
 language: python
+cache: pip
 python:
     - "2.7"
-    - "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
+    - "3.6"
+install:
+    - travis_retry pip install -r requirements.txt
 script: nosetests krypy/tests

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.

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
 <https://github.com/andrenarchy/krypy>`_. 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)

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)
 

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)
@@ -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.
@@ -1989,7 +1896,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.'''
@@ -2114,7 +2021,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]

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',