A bunch more doc fixes.

[dunshire.git] / dunshire / matrices.py

diff --git a/dunshire/matrices.py b/dunshire/matrices.py
index 13e8150f5c807914782dbcb49733f36be513912c..f35827d16be600cec3ee5d1c2d0fbcea27a9484d 100644 (file)
--- a/dunshire/matrices.py
+++ b/dunshire/matrices.py
@@ -3,6 +3,7 @@ Utility functions for working with CVXOPT matrices (instances of the
 class:`cvxopt.base.matrix` class).
 """
 
 class:`cvxopt.base.matrix` class).
 """
 

+from copy import copy

 from math import sqrt
 from cvxopt import matrix
 from cvxopt.lapack import gees, gesdd, syevr
 from math import sqrt
 from cvxopt import matrix
 from cvxopt.lapack import gees, gesdd, syevr


@@ -17,8 +18,11 @@ def append_col(left, right):
     Parameters
     ----------
 
     Parameters
     ----------
 

-    left, right : matrix
-        The two matrices to append to one another.
+    left : matrix
+        The "original" matrix, the one that will wind up on the left.
+
+    right : matrix
+        The matrix to be appended on the right of ``left``.

 
     Returns
     -------
 
     Returns
     -------


@@ -56,8 +60,11 @@ def append_row(top, bottom):
     Parameters
     ----------
 
     Parameters
     ----------
 

-    top, bottom : matrix
-        The two matrices to append to one another.
+    top : matrix
+        The "original" matrix, the one that will wind up on top.
+
+    bottom : matrix
+        The matrix to be appended below ``top``.

 
     Returns
     -------
 
     Returns
     -------


@@ -140,7 +147,10 @@ def eigenvalues(symmat):
 
     domain_dim = symmat.size[0]
     eigs = matrix(0, (domain_dim, 1), tc='d')
 
     domain_dim = symmat.size[0]
     eigs = matrix(0, (domain_dim, 1), tc='d')

-    syevr(symmat, eigs)
+
+    # Create a copy of ``symmat`` here because ``syevr`` clobbers it.
+    dummy = copy(symmat)
+    syevr(dummy, eigs)

     return list(eigs)
 
 
     return list(eigs)
 
 


@@ -261,8 +271,11 @@ def inner_product(vec1, vec2):
     Parameters
     ----------
 
     Parameters
     ----------
 

-    vec1, vec2 : matrix
-        The two vectors whose inner product you want.
+    vec1 : matrix
+        The first vector, whose inner product with ``vec2`` you want.
+
+    vec2 : matrix
+        The second vector, whose inner product with ``vec1`` you want.

 
     Returns
     -------
 
     Returns
     -------


@@ -326,17 +339,53 @@ def norm(matrix_or_vector):
     --------
 
         >>> v = matrix([1,1])
     --------
 
         >>> v = matrix([1,1])

-        >>> print('{:.5f}'.format(norm(v)))
-        1.41421
+        >>> norm(v)
+        1.414...

 
         >>> A = matrix([1,1,1,1], (2,2))
         >>> norm(A)
 
         >>> A = matrix([1,1,1,1], (2,2))
         >>> norm(A)

-        2.0
+        2.0...

 
     """
     return sqrt(inner_product(matrix_or_vector, matrix_or_vector))
 
 
 
     """
     return sqrt(inner_product(matrix_or_vector, matrix_or_vector))
 
 

+def specnorm(mat):
+    """
+    Return the spectral norm of a matrix.
+
+    The spectral norm of a matrix is its largest singular value, and it
+    corresponds to the operator norm induced by the vector Euclidean norm.
+
+    Parameters
+    ----------
+
+    mat : matrix
+        The matrix whose spectral norm you want.
+
+    Examples:
+
+        >>> specnorm(identity(3))
+        1.0
+
+        >>> specnorm(5*identity(4))
+        5.0
+
+    """
+    num_eigs = min(mat.size)
+    eigs = matrix(0, (num_eigs, 1), tc='d')
+    typecode = 'd'
+    if any([isinstance(entry, complex) for entry in mat]):
+        typecode = 'z'
+    dummy = matrix(mat, mat.size, tc=typecode)
+    gesdd(dummy, eigs)
+
+    if len(eigs) > 0:
+        return eigs[0]
+    else:
+        return 0
+
+

 def vec(mat):
     """
     Create a long vector in column-major order from ``mat``.
 def vec(mat):
     """
     Create a long vector in column-major order from ``mat``.


@@ -418,25 +467,25 @@ def condition_number(mat):
     Examples
     --------
 
     Examples
     --------
 

-    >>> condition_number(identity(1, typecode='d'))
-    1.0
-    >>> condition_number(identity(2, typecode='d'))
-    1.0
-    >>> condition_number(identity(3, typecode='d'))
+    >>> condition_number(identity(3))

     1.0
 
     1.0
 

-    >>> A = matrix([[2,1],[1,2]], tc='d')
+    >>> A = matrix([[2,1],[1,2]])

     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
     True
 
     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
     True
 

-    >>> A = matrix([[2,1j],[-1j,2]], tc='z')
+    >>> A = matrix([[2,1j],[-1j,2]])

     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
     True
 
     """
     num_eigs = min(mat.size)
     eigs = matrix(0, (num_eigs, 1), tc='d')
     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
     True
 
     """
     num_eigs = min(mat.size)
     eigs = matrix(0, (num_eigs, 1), tc='d')

-    gesdd(mat, eigs)
+    typecode = 'd'
+    if any([isinstance(entry, complex) for entry in mat]):
+        typecode = 'z'
+    dummy = matrix(mat, mat.size, tc=typecode)
+    gesdd(dummy, eigs)

 
     if len(eigs) > 0:
         return eigs[0]/eigs[-1]
 
     if len(eigs) > 0:
         return eigs[0]/eigs[-1]