X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=dunshire%2Fmatrices.py;h=f35827d16be600cec3ee5d1c2d0fbcea27a9484d;hb=0274de467062ab29d2a41d2a91ec0b28fcd95c8d;hp=94f841d39f888392624077afd12bcdcb7e0c3c19;hpb=cb34e104961ce1da242017a0e5ebf60c3fe18a1d;p=dunshire.git

diff --git a/dunshire/matrices.py b/dunshire/matrices.py
index 94f841d..f35827d 100644
--- a/dunshire/matrices.py
+++ b/dunshire/matrices.py
@@ -18,8 +18,11 @@ def append_col(left, right):
     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
     -------
@@ -57,8 +60,11 @@ def append_row(top, bottom):
     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
     -------
@@ -265,8 +271,11 @@ def inner_product(vec1, vec2):
     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
     -------
@@ -330,17 +339,53 @@ def norm(matrix_or_vector):
     --------
 
         >>> 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)
-        2.0
+        2.0...
 
     """
     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``.
@@ -422,21 +467,25 @@ def condition_number(mat):
     Examples
     --------
 
-    >>> condition_number(identity(3, typecode='d'))
+    >>> condition_number(identity(3))
     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
 
-    >>> 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')
-    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]