Add eigenvalues_re() to compute real parts of eigenvalues of any matrix.

diff --git a/src/dunshire/matrices.py b/src/dunshire/matrices.py
index bf8ae03ee416e443245ed45b76fc1ce614d8990a..384b11fdddcb836f2bef455df8210d0a92b07d47 100644 (file)
--- a/src/dunshire/matrices.py
+++ b/src/dunshire/matrices.py
@@ -3,9 +3,10 @@ Utility functions for working with CVXOPT matrices (instances of the
 class:`cvxopt.base.matrix` class).
 """
 
+from copy import copy
 from math import sqrt
 from cvxopt import matrix
-from cvxopt.lapack import syevr
+from cvxopt.lapack import gees, syevr
 
 import options
 
@@ -95,6 +96,10 @@ def eigenvalues(symmat):
     """
     Return the eigenvalues of the given symmetric real matrix.
 
+    On the surface, this appears redundant to the :func:`eigenvalues_re`
+    function. However, if we know in advance that our input is
+    symmetric, a better algorithm can be used.
+
     Parameters
     ----------
 
@@ -140,6 +145,71 @@ def eigenvalues(symmat):
     return list(eigs)
 
 
+def eigenvalues_re(anymat):
+    """
+    Return the real parts of the eigenvalues of the given square matrix.
+
+    Parameters
+    ----------
+
+    anymat : matrix
+        The square matrix whose eigenvalues you want.
+
+    Returns
+    -------
+
+    list of float
+       A list of the real parts (in no particular order) of the
+       eigenvalues of ``anymat``.
+
+    Raises
+    ------
+
+    TypeError
+        If the input matrix is not square.
+
+    Examples
+    --------
+
+    This is symmetric and has two real eigenvalues:
+
+        >>> A = matrix([[2,1],[1,2]], tc='d')
+        >>> sorted(eigenvalues_re(A))
+        [1.0, 3.0]
+
+    But this rotation matrix has eigenvalues `i` and `-i`, both of whose
+    real parts are zero:
+
+        >>> A = matrix([[0,-1],[1,0]])
+        >>> eigenvalues_re(A)
+        [0.0, 0.0]
+
+    If the input matrix is not square, it doesn't have eigenvalues::
+
+        >>> A = matrix([[1,2],[3,4],[5,6]])
+        >>> eigenvalues_re(A)
+        Traceback (most recent call last):
+        ...
+        TypeError: input matrix must be square
+
+    """
+    if not anymat.size[0] == anymat.size[1]:
+        raise TypeError('input matrix must be square')
+
+    domain_dim = anymat.size[0]
+    eigs = matrix(0, (domain_dim, 1), tc='z')
+
+    # Create a copy of ``anymat`` here for two reasons:
+    #
+    #   1. ``gees`` clobbers its input.
+    #   2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'.
+    #
+    dummy = matrix(anymat, anymat.size, tc='d')
+
+    gees(dummy, eigs)
+    return [eig.real for eig in eigs]
+
+
 def identity(domain_dim):
     """
     Create an identity matrix of the given dimensions.