X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2Fdunshire%2Fmatrices.py;h=ef3c4c0ba6b8926c67332b242a19c070032319ee;hb=2bf0987eea59afe20d1391babcd9b8562f7ec5b0;hp=3670d8f1e07db9cb0de23bce7b9a3fdbea9b0dfb;hpb=e708b535dd508993c2c58dbf552d9acc47d23e49;p=dunshire.git

diff --git a/src/dunshire/matrices.py b/src/dunshire/matrices.py
index 3670d8f..ef3c4c0 100644
--- a/src/dunshire/matrices.py
+++ b/src/dunshire/matrices.py
@@ -5,7 +5,7 @@ class:`cvxopt.base.matrix` class).
 
 from math import sqrt
 from cvxopt import matrix
-from cvxopt.lapack import syev
+from cvxopt.lapack import gees, syevr
 
 import options
 
@@ -95,6 +95,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
     ----------
 
@@ -136,10 +140,75 @@ def eigenvalues(symmat):
 
     domain_dim = symmat.size[0]
     eigs = matrix(0, (domain_dim, 1), tc='d')
-    syev(symmat, eigs)
+    syevr(symmat, eigs)
     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.
@@ -156,6 +225,12 @@ def identity(domain_dim):
     matrix
         A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
 
+    Raises
+    ------
+
+    ValueError
+        If you ask for the identity on zero or fewer dimensions.
+
     Examples
     --------
 
@@ -191,6 +266,12 @@ def inner_product(vec1, vec2):
     float
         The inner product of ``vec1`` and ``vec2``.
 
+    Raises
+    ------
+
+    TypeError
+        If the lengths of ``vec1`` and ``vec2`` differ.
+
     Examples
     --------