Remove unused import.

[dunshire.git] / src / dunshire / matrices.py

diff --git a/src/dunshire/matrices.py b/src/dunshire/matrices.py
index 38d4afe78efaad7665b71d936c8ba9e8309de581..ef3c4c0ba6b8926c67332b242a19c070032319ee 100644 (file)
--- a/src/dunshire/matrices.py
+++ b/src/dunshire/matrices.py
@@ -5,19 +5,28 @@ 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
 
+
 def append_col(left, right):
     """
-    Append the matrix ``right`` to the right side of the matrix ``left``.
+    Append two matrices side-by-side.
 
     Parameters
     ----------
+
     left, right : matrix
         The two matrices to append to one another.
 
+    Returns
+    -------
+
+    matrix
+        A new matrix consisting of ``right`` appended to the right
+        of ``left``.
+
     Examples
     --------
 
@@ -39,15 +48,23 @@ def append_col(left, right):
     """
     return matrix([left.trans(), right.trans()]).trans()
 
+
 def append_row(top, bottom):
     """
-    Append the matrix ``bottom`` to the bottom of the matrix ``top``.
+    Append two matrices top-to-bottom.
 
     Parameters
     ----------
+
     top, bottom : matrix
         The two matrices to append to one another.
 
+    Returns
+    -------
+
+    matrix
+        A new matrix consisting of ``bottom`` appended below ``top``.
+
     Examples
     --------
 
@@ -78,15 +95,25 @@ 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
     ----------
 
     symmat : matrix
         The real symmetric matrix whose eigenvalues you want.
 
+    Returns
+    -------
+
+    list of float
+       A list of the eigenvalues (in no particular order) of ``symmat``.
 
     Raises
     ------
+
     TypeError
         If the input matrix is not symmetric.
 
@@ -113,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.
@@ -129,9 +221,16 @@ def identity(domain_dim):
 
     Returns
     -------
+
     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
     --------
 
@@ -163,9 +262,16 @@ def inner_product(vec1, vec2):
 
     Returns
     -------
+
     float
         The inner product of ``vec1`` and ``vec2``.
 
+    Raises
+    ------
+
+    TypeError
+        If the lengths of ``vec1`` and ``vec2`` differ.
+
     Examples
     --------