# [Scipy-svn] r3952 - trunk/scipy/linalg

scipy-svn@scip... scipy-svn@scip...
Wed Feb 20 19:57:49 CST 2008

```Author: wnbell
Date: 2008-02-20 19:57:37 -0600 (Wed, 20 Feb 2008)
New Revision: 3952

Modified:
trunk/scipy/linalg/basic.py
trunk/scipy/linalg/blas.py
trunk/scipy/linalg/decomp.py
trunk/scipy/linalg/matfuncs.py
Log:
applied patch from user pv
Reformat and ReSTify scipy.linalg docstrings
resolves ticket #596

Modified: trunk/scipy/linalg/basic.py
===================================================================
--- trunk/scipy/linalg/basic.py	2008-02-20 05:39:03 UTC (rev 3951)
+++ trunk/scipy/linalg/basic.py	2008-02-21 01:57:37 UTC (rev 3952)
@@ -25,22 +25,34 @@

def lu_solve((lu, piv), b, trans=0, overwrite_b=0):
-    """ lu_solve((lu, piv), b, trans=0, overwrite_b=0) -> x
+    """Solve an equation system, a x = b, given the LU factorization of a
+
+    Parameters
+    ----------
+    (lu, piv)
+        Factorization of the coefficient matrix a, as given by lu_factor
+    b : array
+        Right-hand side
+    trans : {0, 1, 2}
+        Type of system to solve:

-    Solve a system of equations given a previously factored matrix
+        =====  =========
+        trans  system
+        =====  =========
+        0      a x   = b
+        1      a^T x = b
+        2      a^H x = b
+        =====  =========

-    Inputs:
+    Returns
+    -------
+    x : array
+        Solution to the system

-      (lu,piv) -- The factored matrix, a (the output of lu_factor)
-      b        -- a set of right-hand sides
-      trans    -- type of system to solve:
-                  0 : a   * x = b   (no transpose)
-                  1 : a^T * x = b   (transpose)
-                  2   a^H * x = b   (conjugate transpose)
-
-    Outputs:
-
-       x -- the solution to the system
+    --------
+    lu_factor : LU factorize a matrix
+
"""
b1 = asarray_chkfinite(b)
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
@@ -54,18 +66,24 @@
'illegal value in %-th argument of internal gesv|posv'%(-info)

def cho_solve((c, lower), b, overwrite_b=0):
-    """ cho_solve((c, lower), b, overwrite_b=0) -> x
+    """Solve an equation system, a x = b, given the Cholesky factorization of a

-    Solve a system of equations given a previously cholesky factored matrix
+    Parameters
+    ----------
+    (c, lower)
+        Cholesky factorization of a, as given by cho_factor
+    b : array
+        Right-hand side

-    Inputs:
+    Returns
+    -------
+    x : array
+        The solution to the system a x = b

-      (c,lower) -- The factored matrix, a (the output of cho_factor)
-      b        -- a set of right-hand sides
-
-    Outputs:
-
-       x -- the solution to the system a*x = b
+    --------
+    cho_factor : Cholesky factorization of a matrix
+
"""
b1 = asarray_chkfinite(b)
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
@@ -81,22 +99,29 @@
# Linear equations
def solve(a, b, sym_pos=0, lower=0, overwrite_a=0, overwrite_b=0,
debug = 0):
-    """ solve(a, b, sym_pos=0, lower=0, overwrite_a=0, overwrite_b=0) -> x
+    """Solve the equation a x = b for x

-    Solve a linear system of equations a * x = b for x.
+    Parameters
+    ----------
+    a : array, shape (M, M)
+    b : array, shape (M,) or (M, N)
+    sym_pos : boolean
+        Assume a is symmetric and positive definite
+    lower : boolean
+        Use only data contained in the lower triangle of a, if sym_pos is true.
+        Default is to use upper triangle.
+    overwrite_a : boolean
+        Allow overwriting data in a (may enhance performance)
+    overwrite_b : boolean
+        Allow overwriting data in b (may enhance performance)
+
+    Returns
+    -------
+    x : array, shape (M,) or (M, N) depending on b
+        Solution to the system a x = b

-    Inputs:
-
-      a -- An N x N matrix.
-      b -- An N x nrhs matrix or N vector.
-      sym_pos -- Assume a is symmetric and positive definite.
-      lower -- Assume a is lower triangular, otherwise upper one.
-               Only used if sym_pos is true.
-      overwrite_y - Discard data in y, where y is a or b.
-
-    Outputs:
-
-      x -- The solution to the system a * x = b
+    Raises LinAlgError if a is singular
+
"""
a1, b1 = map(asarray_chkfinite,(a,b))
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -129,29 +154,37 @@

def solve_banded((l,u), ab, b, overwrite_ab=0, overwrite_b=0,
debug = 0):
-    """ solve_banded((l,u), ab, b, overwrite_ab=0, overwrite_b=0) -> x
+    """Solve the equation a x = b for x, assuming a is banded matrix.

-    Solve a linear system of equations a * x = b for x where
-    a is a banded matrix stored in diagonal orded form
+    The matrix a is stored in ab using the matrix diagonal orded form::

-     *   *     a1u
+        ab[u + i - j, j] == a[i,j]

-     *  a12 a23 ...
-    a11 a22 a33 ...
-    a21 a32 a43 ...
-    .
-    al1 ..         *
+    Example of ab (shape of a is (6,6), u=1, l=2)::

-    Inputs:
+        *    a01  a12  a23  a34  a45
+        a00  a11  a22  a33  a44  a55
+        a10  a21  a32  a43  a54   *
+        a20  a31  a42  a53   *    *

-      (l,u) -- number of non-zero lower and upper diagonals, respectively.
-      a -- An N x (l+u+1) matrix.
-      b -- An N x nrhs matrix or N vector.
-      overwrite_y - Discard data in y, where y is ab or b.
+    Parameters
+    ----------
+    (l, u) : (integer, integer)
+        Number of non-zero lower and upper diagonals
+    ab : array, shape (l+u+1, M)
+        Banded matrix
+    b : array, shape (M,) or (M, K)
+        Right-hand side
+    overwrite_ab : boolean
+        Discard data in ab (may enhance performance)
+    overwrite_b : boolean
+        Discard data in b (may enhance performance)

-    Outputs:
-
-      x -- The solution to the system a * x = b
+    Returns
+    -------
+    x : array, shape (M,) or (M, K)
+        The solution to the system a x = b
+
"""
a1, b1 = map(asarray_chkfinite,(ab,b))
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
@@ -171,34 +204,48 @@

def solveh_banded(ab, b, overwrite_ab=0, overwrite_b=0,
lower=0):
-    """ solveh_banded(ab, b, overwrite_ab=0, overwrite_b=0) -> c, x
+    """Solve equation a x = b. a is Hermitian positive-definite banded matrix.

-    Solve a linear system of equations a * x = b for x where
-    a is a banded symmetric or Hermitian positive definite
-    matrix stored in lower diagonal ordered form (lower=1)
+    The matrix a is stored in ab either in lower diagonal or upper
+    diagonal ordered form:

-    a11 a22 a33 a44 a55 a66
-    a21 a32 a43 a54 a65 *
-    a31 a42 a53 a64 *   *
+        ab[u + i - j, j] == a[i,j]        (if upper form; i <= j)
+        ab[    i - j, j] == a[i,j]        (if lower form; i >= j)

-    or upper diagonal ordered form
+    Example of ab (shape of a is (6,6), u=2)::

-    *   *   a31 a42 a53 a64
-    *   a21 a32 a43 a54 a65
-    a11 a22 a33 a44 a55 a66
+        upper form:
+        *   *   a02 a13 a24 a35
+        *   a01 a12 a23 a34 a45
+        a00 a11 a22 a33 a44 a55
+
+        lower form:
+        a00 a11 a22 a33 a44 a55
+        a10 a21 a32 a43 a54 *
+        a20 a31 a42 a53 *   *

-    Inputs:
+    Cells marked with * are not used.

-      ab -- An N x l
-      b -- An N x nrhs matrix or N vector.
-      overwrite_y - Discard data in y, where y is ab or b.
-      lower - is ab in lower or upper form?
+    Parameters
+    ----------
+    ab : array, shape (M, u + 1)
+        Banded matrix
+    b : array, shape (M,) or (M, K)
+        Right-hand side
+    overwrite_ab : boolean
+        Discard data in ab (may enhance performance)
+    overwrite_b : boolean
+        Discard data in b (may enhance performance)
+    lower : boolean
+        Is the matrix in the lower form. (Default is upper form)

-    Outputs:
-
-      c:  the Cholesky factorization of ab
-      x:  the solution to ab * x = b
-
+    Returns
+    -------
+    c : array, shape (M, u+1)
+        Cholesky factorization of a, in the same banded format as ab
+    x : array, shape (M,) or (M, K)
+        The solution to the system a x = b
+
"""
ab, b = map(asarray_chkfinite,(ab,b))

@@ -215,32 +262,40 @@
'illegal value in %d-th argument of internal pbsv'%(-info)

def cholesky_banded(ab, overwrite_ab=0, lower=0):
-    """ cholesky_banded(ab, overwrite_ab=0, lower=0) -> c
-
-    Compute the Cholesky decomposition of a
-    banded symmetric or Hermitian positive definite
-    matrix stored in lower diagonal ordered form (lower=1)
-
-    a11 a22 a33 a44 a55 a66
-    a21 a32 a43 a54 a65 *
-    a31 a42 a53 a64 *   *
-
-    or upper diagonal ordered form
-
-    *   *   a31 a42 a53 a64
-    *   a21 a32 a43 a54 a65
-    a11 a22 a33 a44 a55 a66
-
-    Inputs:
-
-      ab -- An N x l
-      overwrite_ab - Discard data in ab
-      lower - is ab in lower or upper form?
-
-    Outputs:
-
-      c:  the Cholesky factorization of ab
-
+    """Cholesky decompose a banded Hermitian positive-definite matrix
+
+    The matrix a is stored in ab either in lower diagonal or upper
+    diagonal ordered form:
+
+        ab[u + i - j, j] == a[i,j]        (if upper form; i <= j)
+        ab[    i - j, j] == a[i,j]        (if lower form; i >= j)
+
+    Example of ab (shape of a is (6,6), u=2)::
+
+        upper form:
+        *   *   a02 a13 a24 a35
+        *   a01 a12 a23 a34 a45
+        a00 a11 a22 a33 a44 a55
+
+        lower form:
+        a00 a11 a22 a33 a44 a55
+        a10 a21 a32 a43 a54 *
+        a20 a31 a42 a53 *   *
+
+    Parameters
+    ----------
+    ab : array, shape (M, u + 1)
+        Banded matrix
+    overwrite_ab : boolean
+        Discard data in ab (may enhance performance)
+    lower : boolean
+        Is the matrix in the lower form. (Default is upper form)
+
+    Returns
+    -------
+    c : array, shape (M, u+1)
+        Cholesky factorization of a, in the same banded format as ab
+
"""
ab = asarray_chkfinite(ab)

@@ -259,9 +314,30 @@

# matrix inversion
def inv(a, overwrite_a=0):
-    """ inv(a, overwrite_a=0) -> a_inv
+    """Compute the inverse of a matrix.
+
+    Parameters
+    ----------
+    a : array-like, shape (M, M)
+        Matrix to be inverted
+
+    Returns
+    -------
+    ainv : array-like, shape (M, M)
+        Inverse of the matrix a

-    Return inverse of square matrix a.
+    Raises LinAlgError if a is singular
+
+    Examples
+    --------
+    >>> a = array([[1., 2.], [3., 4.]])
+    >>> inv(a)
+    array([[-2. ,  1. ],
+           [ 1.5, -0.5]])
+    >>> dot(a, inv(a))
+    array([[ 1.,  0.],
+           [ 0.,  1.]])
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -312,35 +388,39 @@
## matrix and Vector norm
import decomp
def norm(x, ord=None):
-    """ norm(x, ord=None) -> n
+    """Matrix or vector norm.

-    Matrix or vector norm.
+    Parameters
+    ----------
+    x : array, shape (M,) or (M, N)
+    ord : number, or {None, 1, -1, 2, -2, inf, -inf, 'fro'}
+        Order of the norm:

-    Inputs:
+        =====  ============================  ==========================
+        ord    norm for matrices             norm for vectors
+        =====  ============================  ==========================
+        None   Frobenius norm                2-norm
+        'fro'  Frobenius norm                -
+        inf    max(sum(abs(x), axis=1))      max(abs(x))
+        -inf   min(sum(abs(x), axis=1))      min(abs(x))
+        1      max(sum(abs(x), axis=0))      as below
+        -1     min(sum(abs(x), axis=0))      as below
+        2      2-norm (largest sing. value)  as below
+        -2     smallest singular value       as below
+        other  -                             sum(abs(x)**ord)**(1./ord)
+        =====  ============================  ==========================
+
+    Returns
+    -------
+    n : float
+        Norm of the matrix or vector

-      x -- a rank-1 (vector) or rank-2 (matrix) array
-      ord -- the order of the norm.
-
-       For arrays of any rank, if ord is None:
-         calculate the square norm (Euclidean norm for vectors, Frobenius norm for matrices)
-
-       For vectors ord can be any real number including Inf or -Inf.
-         ord = Inf, computes the maximum of the magnitudes
-         ord = -Inf, computes minimum of the magnitudes
-         ord is finite, computes sum(abs(x)**ord,axis=0)**(1.0/ord)
-
-       For matrices ord can only be one of the following values:
-         ord = 2 computes the largest singular value
-         ord = -2 computes the smallest singular value
-         ord = 1 computes the largest column sum of absolute values
-         ord = -1 computes the smallest column sum of absolute values
-         ord = Inf computes the largest row sum of absolute values
-         ord = -Inf computes the smallest row sum of absolute values
-         ord = 'fro' computes the frobenius norm sqrt(sum(diag(X.H * X),axis=0))
-
-       For values ord < 0, the result is, strictly speaking, not a
-       mathematical 'norm', but it may still be useful for numerical purposes.
+    Notes
+    -----
+    For values ord < 0, the result is, strictly speaking, not a
+    mathematical 'norm', but it may still be useful for numerical
+    purposes.
+
"""
x = asarray_chkfinite(x)
if ord is None: # check the default case first and handle it immediately
@@ -382,9 +462,20 @@
### Determinant

def det(a, overwrite_a=0):
-    """ det(a, overwrite_a=0) -> d
+    """Compute the determinant of a matrix

-    Return determinant of a square matrix.
+    Parameters
+    ----------
+    a : array, shape (M, M)
+
+    Returns
+    -------
+    det : float or complex
+        Determinant of a
+
+    Notes
+    -----
+    The determinant is computed via LU factorization, LAPACK routine z/dgetrf.
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -399,25 +490,38 @@
### Linear Least Squares

def lstsq(a, b, cond=None, overwrite_a=0, overwrite_b=0):
-    """ lstsq(a, b, cond=None, overwrite_a=0, overwrite_b=0) -> x,resids,rank,s
-
-    Return least-squares solution of a * x = b.
-
-    Inputs:
-
-      a -- An M x N matrix.
-      b -- An M x nrhs matrix or M vector.
-      cond -- Used to determine effective rank of a.
-
-    Outputs:
-
-      x -- The solution (N x nrhs matrix) to the minimization problem:
-                  2-norm(| b - a * x |) -> min
-      resids -- The residual sum-of-squares for the solution matrix x
-                (only if M>N and rank==N).
-      rank -- The effective rank of a.
-      s -- Singular values of a in decreasing order. The condition number
-           of a is abs(s[0]/s[-1]).
+    """Compute least-squares solution to equation :m:`a x = b`
+
+    Compute a vector x such that the 2-norm :m:`|b - a x|` is minimised.
+
+    Parameters
+    ----------
+    a : array, shape (M, N)
+    b : array, shape (M,) or (M, K)
+    cond : float
+        Cutoff for 'small' singular values; used to determine effective
+        rank of a. Singular values smaller than rcond*largest_singular_value
+        are considered zero.
+    overwrite_a : boolean
+        Discard data in a (may enhance performance)
+    overwrite_b : boolean
+        Discard data in b (may enhance performance)
+
+    Returns
+    -------
+    x : array, shape (N,) or (N, K) depending on shape of b
+        Least-squares solution
+    residues : array, shape () or (1,) or (K,)
+        Sums of residues, squared 2-norm for each column in :m:`b - a x`
+        If rank of matrix a is < N or > M this is an empty array.
+        If b was 1-d, this is an (1,) shape array, otherwise the shape is (K,)
+    rank : integer
+        Effective rank of matrix a
+    s : array, shape (min(M,N),)
+        Singular values of a. The condition number of a is abs(s[0]/s[-1]).
+
+    Raises LinAlgError if computation does not converge
+
"""
a1, b1 = map(asarray_chkfinite,(a,b))
if len(a1.shape) != 2:
@@ -457,9 +561,36 @@

def pinv(a, cond=None, rcond=None):
-    """ pinv(a, rcond=None) -> a_pinv
+    """Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+
+    Calculate a generalized inverse of a matrix using a least-squares
+    solver.
+
+    Parameters
+    ----------
+    a : array, shape (M, N)
+        Matrix to be pseudo-inverted
+    cond, rcond : float
+        Cutoff for 'small' singular values in the least-squares solver.
+        Singular values smaller than rcond*largest_singular_value are
+        considered zero.
+
+    Returns
+    -------
+    B : array, shape (N, M)
+
+    Raises LinAlgError if computation does not converge

-    Compute generalized inverse of A using least-squares solver.
+    Examples
+    --------
+    >>> from numpy import *
+    >>> a = random.randn(9, 6)
+    >>> B = linalg.pinv(a)
+    >>> allclose(a, dot(a, dot(B, a)))
+    True
+    >>> allclose(B, dot(B, dot(a, B)))
+    True
+
"""
a = asarray_chkfinite(a)
b = numpy.identity(a.shape[0], dtype=a.dtype)
@@ -473,9 +604,39 @@

_array_precision = {'f': 0, 'd': 1, 'F': 0, 'D': 1}
def pinv2(a, cond=None, rcond=None):
-    """ pinv2(a, rcond=None) -> a_pinv
+    """Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+
+    Calculate a generalized inverse of a matrix using its
+    singular-value decomposition and including all 'large' singular
+    values.
+
+    Parameters
+    ----------
+    a : array, shape (M, N)
+        Matrix to be pseudo-inverted
+    cond, rcond : float or None
+        Cutoff for 'small' singular values.
+        Singular values smaller than rcond*largest_singular_value are
+        considered zero.

-    Compute the generalized inverse of A using svd.
+        If None or -1, suitable machine precision is used.
+
+    Returns
+    -------
+    B : array, shape (N, M)
+
+    Raises LinAlgError if SVD computation does not converge
+
+    Examples
+    --------
+    >>> from numpy import *
+    >>> a = random.randn(9, 6)
+    >>> B = linalg.pinv2(a)
+    >>> allclose(a, dot(a, dot(B, a)))
+    True
+    >>> allclose(B, dot(B, dot(a, B)))
+    True
+
"""
a = asarray_chkfinite(a)
u, s, vh = decomp.svd(a)
@@ -498,8 +659,37 @@
#-----------------------------------------------------------------------------

def tri(N, M=None, k=0, dtype=None):
-    """ returns a N-by-M matrix where all the diagonals starting from
-        lower left corner up to the k-th are all ones.
+    """Construct (N, M) matrix filled with ones at and below the k-th diagonal.
+
+    The matrix has A[i,j] == 1 for i <= j + k
+
+    Parameters
+    ----------
+    N : integer
+    M : integer
+        Size of the matrix. If M is None, M == N is assumed.
+    k : integer
+        Number of subdiagonal below which matrix is filled with ones.
+        k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
+    dtype : dtype
+        Data type of the matrix.
+
+    Returns
+    -------
+    A : array, shape (N, M)
+
+    Examples
+    --------
+    >>> from scipy.linalg import tri
+    >>> tri(3, 5, 2, dtype=int)
+    array([[1, 1, 1, 0, 0],
+           [1, 1, 1, 1, 0],
+           [1, 1, 1, 1, 1]])
+    >>> tri(3, 5, -1, dtype=int)
+    array([[0, 0, 0, 0, 0],
+           [1, 0, 0, 0, 0],
+           [1, 1, 0, 0, 0]])
+
"""
if M is None: M = N
if type(M) == type('d'):
@@ -514,8 +704,29 @@
return m.astype(dtype)

def tril(m, k=0):
-    """ returns the elements on and below the k-th diagonal of m.  k=0 is the
-        main diagonal, k > 0 is above and k < 0 is below the main diagonal.
+    """Construct a copy of a matrix with elements above the k-th diagonal zeroed.
+
+    Parameters
+    ----------
+    m : array
+        Matrix whose elements to return
+    k : integer
+        Diagonal above which to zero elements.
+        k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
+
+    Returns
+    -------
+    A : array, shape m.shape, dtype m.dtype
+
+    Examples
+    --------
+    >>> from scipy.linalg import tril
+    >>> tril([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
+    array([[ 0,  0,  0],
+           [ 4,  0,  0],
+           [ 7,  8,  0],
+           [10, 11, 12]])
+
"""
svsp = getattr(m,'spacesaver',lambda:0)()
m = asarray(m)
@@ -524,8 +735,29 @@
return out

def triu(m, k=0):
-    """ returns the elements on and above the k-th diagonal of m.  k=0 is the
-        main diagonal, k > 0 is above and k < 0 is below the main diagonal.
+    """Construct a copy of a matrix with elements below the k-th diagonal zeroed.
+
+    Parameters
+    ----------
+    m : array
+        Matrix whose elements to return
+    k : integer
+        Diagonal below which to zero elements.
+        k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
+
+    Returns
+    -------
+    A : array, shape m.shape, dtype m.dtype
+
+    Examples
+    --------
+    >>> from scipy.linalg import tril
+    >>> triu([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
+    array([[ 1,  2,  3],
+           [ 4,  5,  6],
+           [ 0,  8,  9],
+           [ 0,  0, 12]])
+
"""
svsp = getattr(m,'spacesaver',lambda:0)()
m = asarray(m)
@@ -534,16 +766,36 @@
return out

def toeplitz(c,r=None):
-    """ Construct a toeplitz matrix (i.e. a matrix with constant diagonals).
-
-        Description:
-
-           toeplitz(c,r) is a non-symmetric Toeplitz matrix with c as its first
-           column and r as its first row.
-
-           toeplitz(c) is a symmetric (Hermitian) Toeplitz matrix (r=c).
-
+    """Construct a Toeplitz matrix.
+
+    The Toepliz matrix has constant diagonals, c as its first column,
+    and r as its first row (if not given, r == c is assumed).
+
+    Parameters
+    ----------
+    c : array
+        First column of the matrix
+    r : array
+        First row of the matrix. If None, r == c is assumed.
+
+    Returns
+    -------
+    A : array, shape (len(c), len(r))
+        Constructed Toeplitz matrix.
+        dtype is the same as (c[0] + r[0]).dtype
+
+    Examples
+    --------
+    >>> from scipy.linalg import toeplitz
+    >>> toeplitz([1,2,3], [1,4,5,6])
+    array([[1, 4, 5, 6],
+           [2, 1, 4, 5],
+           [3, 2, 1, 4]])
+
+    --------
+    hankel : Hankel matrix
+
"""
isscalar = numpy.isscalar
if isscalar(c) or isscalar(r):
@@ -566,17 +818,37 @@

def hankel(c,r=None):
-    """ Construct a hankel matrix (i.e. matrix with constant anti-diagonals).
-
-        Description:
-
-          hankel(c,r) is a Hankel matrix whose first column is c and whose
-          last row is r.
-
-          hankel(c) is a square Hankel matrix whose first column is C.
-          Elements below the first anti-diagonal are zero.
-
+    """Construct a Hankel matrix.
+
+    The Hankel matrix has constant anti-diagonals, c as its first column,
+    and r as its last row (if not given, r == 0 os assumed).
+
+    Parameters
+    ----------
+    c : array
+        First column of the matrix
+    r : array
+        Last row of the matrix. If None, r == 0 is assumed.
+
+    Returns
+    -------
+    A : array, shape (len(c), len(r))
+        Constructed Hankel matrix.
+        dtype is the same as (c[0] + r[0]).dtype
+
+    Examples
+    --------
+    >>> from scipy.linalg import hankel
+    >>> hankel([1,2,3,4], [4,7,7,8,9])
+    array([[1, 2, 3, 4, 7],
+           [2, 3, 4, 7, 7],
+           [3, 4, 7, 7, 8],
+           [4, 7, 7, 8, 9]])
+
+    --------
+    toeplitz : Toeplitz matrix
+
"""
isscalar = numpy.isscalar
if isscalar(c) or isscalar(r):
@@ -599,12 +871,32 @@
return map(Matrix,args)

def kron(a,b):
-    """kronecker product of a and b
+    """Kronecker product of a and b.

-    Kronecker product of two matrices is block matrix
-    [[ a[ 0 ,0]*b, a[ 0 ,1]*b, ... , a[ 0 ,n-1]*b  ],
-     [ ...                                   ...   ],
-     [ a[m-1,0]*b, a[m-1,1]*b, ... , a[m-1,n-1]*b  ]]
+    The result is the block matrix::
+
+        a[0,0]*b    a[0,1]*b  ... a[0,-1]*b
+        a[1,0]*b    a[1,1]*b  ... a[1,-1]*b
+        ...
+        a[-1,0]*b   a[-1,1]*b ... a[-1,-1]*b
+
+    Parameters
+    ----------
+    a : array, shape (M, N)
+    b : array, shape (P, Q)
+
+    Returns
+    -------
+    A : array, shape (M*P, N*Q)
+        Kronecker product of a and b
+
+    Examples
+    --------
+    >>> from scipy import kron, array
+    >>> kron(array([[1,2],[3,4]]), array([[1,1,1]]))
+    array([[1, 1, 1, 2, 2, 2],
+           [3, 3, 3, 4, 4, 4]])
+
"""
if not a.flags['CONTIGUOUS']:
a = reshape(a, a.shape)

Modified: trunk/scipy/linalg/blas.py
===================================================================
--- trunk/scipy/linalg/blas.py	2008-02-20 05:39:03 UTC (rev 3951)
+++ trunk/scipy/linalg/blas.py	2008-02-21 01:57:37 UTC (rev 3952)
@@ -22,12 +22,15 @@
_inv_type_conv = {'s':'f','d':'d','c':'F','z':'D'}

def has_column_major_storage(arr):
+    """Is array stored in column-major format"""
return arr.flags['FORTRAN']

def get_blas_funcs(names,arrays=(),debug=0):
"""Return available BLAS function objects with names.
arrays are used to determine the optimal prefix of
-    BLAS routines."""
+    BLAS routines.
+
+    """
ordering = []
for i in range(len(arrays)):
t = arrays[i].dtype.char

Modified: trunk/scipy/linalg/decomp.py
===================================================================
--- trunk/scipy/linalg/decomp.py	2008-02-20 05:39:03 UTC (rev 3951)
+++ trunk/scipy/linalg/decomp.py	2008-02-21 01:57:37 UTC (rev 3952)
@@ -98,32 +98,54 @@
return w, vr

def eig(a,b=None, left=False, right=True, overwrite_a=False, overwrite_b=False):
-    """ Solve ordinary and generalized eigenvalue problem
-    of a square matrix.
+    """Solve an ordinary or generalized eigenvalue problem of a square matrix.

-    Inputs:
+    Find eigenvalues w and right or left eigenvectors of a general matrix::
+
+        a   vr[:,i] = w[i]        b   vr[:,i]
+        a.H vl[:,i] = w[i].conj() b.H vl[:,i]
+
+    where .H is the Hermitean conjugation.
+
+    Parameters
+    ----------
+    a : array, shape (M, M)
+        A complex or real matrix whose eigenvalues and eigenvectors
+        will be computed.
+    b : array, shape (M, M)
+        Right-hand side matrix in a generalized eigenvalue problem.
+        If omitted, identity matrix is assumed.
+    left : boolean
+        Whether to calculate and return left eigenvectors
+    right : boolean
+        Whether to calculate and return right eigenvectors
+
+    overwrite_a : boolean
+        Whether to overwrite data in a (may improve performance)
+    overwrite_b : boolean
+        Whether to overwrite data in b (may improve performance)
+
+    Returns
+    -------
+    w : double or complex array, shape (M,)
+        The eigenvalues, each repeated according to its multiplicity.

-      a     -- An N x N matrix.
-      b     -- An N x N matrix [default is identity(N)].
-      left  -- Return left eigenvectors [disabled].
-      right -- Return right eigenvectors [enabled].
-      overwrite_a, overwrite_b -- save space by overwriting the a and/or
-                                  b matrices (both False by default)
+    (if left == True)
+    vl : double or complex array, shape (M, M)
+        The normalized left eigenvector corresponding to the eigenvalue w[i]
+        is the column v[:,i].
+
+    (if right == True)
+    vr : double or complex array, shape (M, M)
+        The normalized right eigenvector corresponding to the eigenvalue w[i]
+        is the column vr[:,i].
+
+    Raises LinAlgError if eigenvalue computation does not converge

-    Outputs:
-
-      w      -- eigenvalues [left==right==False].
-      w,vr   -- w and right eigenvectors [left==False,right=True].
-      w,vl   -- w and left eigenvectors [left==True,right==False].
-      w,vl,vr  -- [left==right==True].
-
-    Definitions:
-
-      a * vr[:,i] = w[i] * b * vr[:,i]
-
-      a^H * vl[:,i] = conjugate(w[i]) * b^H * vl[:,i]
-
-    where a^H denotes transpose(conjugate(a)).
+    --------
+    eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -182,28 +204,44 @@
return w, vr

def eigh(a, lower=True, eigvals_only=False, overwrite_a=False):
-    """ Solve real symmetric or complex hermitian eigenvalue problem.
+    """Solve the eigenvalue problem for a Hermitian or real symmetric matrix.

-    Inputs:
+    Find eigenvalues w and optionally right eigenvectors v of a::
+
+        a v[:,i] = w[i] v[:,i]
+        v.H v    = identity
+
+    Parameters
+    ----------
+    a : array, shape (M, M)
+        A complex Hermitian or symmetric real matrix whose eigenvalues
+        and eigenvectors will be computed.
+    lower : boolean
+        Whether the pertinent array data is taken from the lower or upper
+        triangle of a. (Default: lower)
+    eigvals_only : boolean
+        Whether to calculate only eigenvalues and no eigenvectors.
+        (Default: both are calculated)
+    overwrite_a : boolean
+        Whether data in a is overwritten (may improve performance).
+
+    Returns
+    -------
+    w : double array, shape (M,)
+        The eigenvalues, in ascending order, each repeated according to its
+        multiplicity.
+
+    (if eigvals_only == False)
+    v : double or complex double array, shape (M, M)
+        The normalized eigenvector corresponding to the eigenvalue w[i] is
+        the column v[:,i].
+
+    Raises LinAlgError if eigenvalue computation does not converge

-      a            -- A hermitian N x N matrix.
-      lower        -- values in a are read from lower triangle
-                      [True: UPLO='L' (default) / False: UPLO='U']
-      eigvals_only -- don't compute eigenvectors.
-      overwrite_a  -- content of a may be destroyed
-
-    Outputs:
-
-      For eigvals_only == False (the default),
-      w,v     -- w: eigenvalues, v: eigenvectors
-      For eigvals_only == True,
-      w       -- eigenvalues
-
-    Definitions:
-
-      a * v[:,i] = w[i] * vr[:,i]
-      v.H * v = identity
-
+    --------
+    eig : eigenvalues and right eigenvectors for non-symmetric arrays
+
"""
if eigvals_only or overwrite_a:
a1 = asarray_chkfinite(a)
@@ -258,43 +296,74 @@

def eig_banded(a_band, lower=0, eigvals_only=0, overwrite_a_band=0,
select='a', select_range=None, max_ev = 0):
-    """ Solve real symmetric or complex hermetian band matrix problem.
+    """Solve real symmetric or complex hermetian band matrix eigenvalue problem.

-    Inputs:
+    Find eigenvalues w and optionally right eigenvectors v of a::
+
+        a v[:,i] = w[i] v[:,i]
+        v.H v    = identity
+
+    The matrix a is stored in ab either in lower diagonal or upper
+    diagonal ordered form:

-      a_band            -- A hermitian N x M matrix in 'packed storage'.
-                           Packed storage looks like this: ('upper form')
-                           [ ... (more off-diagonals) ...,
-                            [ *   *   a13 a24 a35 a46 ... a(n-2)(n)],
-                            [ *   a12 a23 a34 a45 a56 ... a(n-1)(n)],
-                            [ a11 a22 a33 a44 a55 a66 ... a(n)(n)  ]]
-                           The cells denoted with * may contain anything.
-      lower             -- a is in lower packed storage
-                           (default: upper packed form)
-      eigvals_only      -- if True, don't compute eigenvectors.
-      overwrite_a_band  -- content of a may be destroyed
-      select       -- 'a', 'all', 0   : return all eigenvalues/vectors
-                      'v', 'value', 1 : eigenvalues in the interval (min, max]
-                                        will be found
-                      'i', 'index', 2 : eigenvalues with index [min...max]
-                                        will be found
-      select_range -- select == 'v': eigenvalue limits as tuple (min, max)
-                      select == 'i': index limits as tuple (min, max)
-                      select == 'a': meaningless
-      max_ev       -- select == 'v': set to max. number of eigenvalues that is
-                                     expected. In doubt, leave untouched.
-                      select == 'i', 'a': meaningless
+        ab[u + i - j, j] == a[i,j]        (if upper form; i <= j)
+        ab[    i - j, j] == a[i,j]        (if lower form; i >= j)

-    Outputs:
+    Example of ab (shape of a is (6,6), u=2)::

-      w,v     -- w: eigenvalues, v: eigenvectors [for eigvals_only == False]
-      w       -- eigenvalues [for eigvals_only == True].
+        upper form:
+        *   *   a02 a13 a24 a35
+        *   a01 a12 a23 a34 a45
+        a00 a11 a22 a33 a44 a55
+
+        lower form:
+        a00 a11 a22 a33 a44 a55
+        a10 a21 a32 a43 a54 *
+        a20 a31 a42 a53 *   *

-    Definitions:
+    Cells marked with * are not used.

-      a_full * v[:,i] = w[i] * v[:,i]  (with full matrix corresponding to a)
-      v.H * v = identity
+    Parameters
+    ----------
+    a_band : array, shape (M, u+1)
+        Banded matrix whose eigenvalues to calculate
+    lower : boolean
+        Is the matrix in the lower form. (Default is upper form)
+    eigvals_only : boolean
+        Compute only the eigenvalues and no eigenvectors.
+        (Default: calculate also eigenvectors)
+    overwrite_a_band:
+        Discard data in a_band (may enhance performance)
+    select: {'a', 'v', 'i'}
+        Which eigenvalues to calculate

+        ======  ========================================
+        select  calculated
+        ======  ========================================
+        'a'     All eigenvalues
+        'v'     Eigenvalues in the interval (min, max]
+        'i'     Eigenvalues with indices min <= i <= max
+        ======  ========================================
+    select_range : (min, max)
+        Range of selected eigenvalues
+    max_ev : integer
+        For select=='v', maximum number of eigenvalues expected.
+        For other values of select, has no meaning.
+
+        In doubt, leave this parameter untouched.
+
+    Returns
+    -------
+    w : array, shape (M,)
+        The eigenvalues, in ascending order, each repeated according to its
+        multiplicity.
+
+    v : double or complex double array, shape (M, M)
+        The normalized eigenvector corresponding to the eigenvalue w[i] is
+        the column v[:,i].
+
+    Raises LinAlgError if eigenvalue computation does not converge
+
"""
if eigvals_only or overwrite_a_band:
a1 = asarray_chkfinite(a_band)
@@ -374,32 +443,180 @@
return w, v

def eigvals(a,b=None,overwrite_a=0):
-    """Return eigenvalues of square matrix."""
+    """Compute eigenvalues from an ordinary or generalized eigenvalue problem.
+
+    Find eigenvalues of a general matrix::
+
+        a   vr[:,i] = w[i]        b   vr[:,i]
+
+    Parameters
+    ----------
+    a : array, shape (M, M)
+        A complex or real matrix whose eigenvalues and eigenvectors
+        will be computed.
+    b : array, shape (M, M)
+        Right-hand side matrix in a generalized eigenvalue problem.
+        If omitted, identity matrix is assumed.
+    overwrite_a : boolean
+        Whether to overwrite data in a (may improve performance)
+
+    Returns
+    -------
+    w : double or complex array, shape (M,)
+        The eigenvalues, each repeated according to its multiplicity,
+        but not in any specific order.
+
+    Raises LinAlgError if eigenvalue computation does not converge
+
+    --------
+    eigvalsh : eigenvalues of symmetric or Hemitiean arrays
+    eig : eigenvalues and right eigenvectors of general arrays
+    eigh : eigenvalues and eigenvectors of symmetric/Hermitean arrays.
+
+    """
return eig(a,b=b,left=0,right=0,overwrite_a=overwrite_a)

def eigvalsh(a,lower=1,overwrite_a=0):
-    """Return eigenvalues of hermitean or real symmetric matrix."""
+    """Solve the eigenvalue problem for a Hermitian or real symmetric matrix.
+
+    Find eigenvalues w of a::
+
+        a v[:,i] = w[i] v[:,i]
+        v.H v    = identity
+
+    Parameters
+    ----------
+    a : array, shape (M, M)
+        A complex Hermitian or symmetric real matrix whose eigenvalues
+        and eigenvectors will be computed.
+    lower : boolean
+        Whether the pertinent array data is taken from the lower or upper
+        triangle of a. (Default: lower)
+    overwrite_a : boolean
+        Whether data in a is overwritten (may improve performance).
+
+    Returns
+    -------
+    w : double array, shape (M,)
+        The eigenvalues, in ascending order, each repeated according to its
+        multiplicity.
+
+    Raises LinAlgError if eigenvalue computation does not converge
+
+    --------
+    eigvals : eigenvalues of general arrays
+    eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
+    eig : eigenvalues and right eigenvectors for non-symmetric arrays
+
+    """
return eigh(a,lower=lower,eigvals_only=1,overwrite_a=overwrite_a)

def eigvals_banded(a_band,lower=0,overwrite_a_band=0,
select='a', select_range=None):
-    """Return eigenvalues of hermitean or real symmetric matrix."""
+    """Solve real symmetric or complex hermetian band matrix eigenvalue problem.
+
+    Find eigenvalues w of a::
+
+        a v[:,i] = w[i] v[:,i]
+        v.H v    = identity
+
+    The matrix a is stored in ab either in lower diagonal or upper
+    diagonal ordered form:
+
+        ab[u + i - j, j] == a[i,j]        (if upper form; i <= j)
+        ab[    i - j, j] == a[i,j]        (if lower form; i >= j)
+
+    Example of ab (shape of a is (6,6), u=2)::
+
+        upper form:
+        *   *   a02 a13 a24 a35
+        *   a01 a12 a23 a34 a45
+        a00 a11 a22 a33 a44 a55
+
+        lower form:
+        a00 a11 a22 a33 a44 a55
+        a10 a21 a32 a43 a54 *
+        a20 a31 a42 a53 *   *
+
+    Cells marked with * are not used.
+
+    Parameters
+    ----------
+    a_band : array, shape (M, u+1)
+        Banded matrix whose eigenvalues to calculate
+    lower : boolean
+        Is the matrix in the lower form. (Default is upper form)
+    overwrite_a_band:
+        Discard data in a_band (may enhance performance)
+    select: {'a', 'v', 'i'}
+        Which eigenvalues to calculate
+
+        ======  ========================================
+        select  calculated
+        ======  ========================================
+        'a'     All eigenvalues
+        'v'     Eigenvalues in the interval (min, max]
+        'i'     Eigenvalues with indices min <= i <= max
+        ======  ========================================
+    select_range : (min, max)
+        Range of selected eigenvalues
+
+    Returns
+    -------
+    w : array, shape (M,)
+        The eigenvalues, in ascending order, each repeated according to its
+        multiplicity.
+
+    Raises LinAlgError if eigenvalue computation does not converge
+
+    --------
+    eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian band matrices
+    eigvals : eigenvalues of general arrays
+    eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
+    eig : eigenvalues and right eigenvectors for non-symmetric arrays
+
+    """
return eig_banded(a_band,lower=lower,eigvals_only=1,
overwrite_a_band=overwrite_a_band, select=select,
select_range=select_range)

def lu_factor(a, overwrite_a=0):
-    """Return raw LU decomposition of a matrix and pivots, for use in solving
-    a system of linear equations.
+    """Compute pivoted LU decomposition of a matrix.
+
+    The decomposition is::

-    Inputs:
+        A = P L U

-      a  --- an NxN matrix
+    where P is a permutation matrix, L lower triangular with unit
+    diagonal elements, and U upper triangular.
+
+    Parameters
+    ----------
+    a : array, shape (M, M)
+        Matrix to decompose
+    overwrite_a : boolean
+        Whether to overwrite data in A (may increase performance)

-    Outputs:
+    Returns
+    -------
+    lu : array, shape (N, N)
+        Matrix containing U in its upper triangle, and L in its lower triangle.
+        The unit diagonal elements of L are not stored.
+    piv : array, shape (N,)
+        Pivot indices representing the permutation matrix P:
+        row i of matrix was interchanged with row piv[i].

-      lu --- the lu factorization matrix
-      piv --- an array of pivots
+    --------
+    lu_solve : solve an equation system using the LU factorization of a matrix
+
+    Notes
+    -----
+    This is a wrapper to the *GETRF routines from LAPACK.
+
"""
a1 = asarray(a)
if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
@@ -414,8 +631,24 @@
return lu, piv

def lu_solve(a_lu_pivots,b):
-    """Solve a previously factored system.  First input is a tuple (lu, pivots)
-    which is the output to lu_factor.  Second input is the right hand side.
+    """Solve an equation system, a x = b, given the LU factorization of a
+
+    Parameters
+    ----------
+    (lu, piv)
+        Factorization of the coefficient matrix a, as given by lu_factor
+    b : array
+        Right-hand side
+
+    Returns
+    -------
+    x : array
+        Solution to the system
+
+    --------
+    lu_factor : LU factorize a matrix
+
"""
a_lu, pivots = a_lu_pivots
a_lu = asarray_chkfinite(a_lu)
@@ -435,28 +668,46 @@

def lu(a,permute_l=0,overwrite_a=0):
-    """Return LU decompostion of a matrix.
+    """Compute pivoted LU decompostion of a matrix.

-    Inputs:
+    The decomposition is::

-      a     -- An M x N matrix.
-      permute_l  -- Perform matrix multiplication p * l [disabled].
+        A = P L U

-    Outputs:
+    where P is a permutation matrix, L lower triangular with unit
+    diagonal elements, and U upper triangular.
+
+    Parameters
+    ----------
+    a : array, shape (M, N)
+        Array to decompose
+    permute_l : boolean
+        Perform the multiplication P*L  (Default: do not permute)
+    overwrite_a : boolean
+        Whether to overwrite data in a (may improve performance)

-      p,l,u  -- LU decomposition matrices of a [permute_l=0]
-      pl,u   -- LU decomposition matrices of a [permute_l=1]
-
-    Definitions:
-
-      a = p * l * u    [permute_l=0]
-      a = pl * u       [permute_l=1]
-
-      p   -  An M x M permutation matrix
-      l   -  An M x K lower triangular or trapezoidal matrix
-             with unit-diagonal
-      u   -  An K x N upper triangular or trapezoidal matrix
-             K = min(M,N)
+    Returns
+    -------
+    (If permute_l == False)
+    p : array, shape (M, M)
+        Permutation matrix
+    l : array, shape (M, K)
+        Lower triangular or trapezoidal matrix with unit diagonal.
+        K = min(M, N)
+    u : array, shape (K, N)
+        Upper triangular or trapezoidal matrix
+
+    (If permute_l == True)
+    pl : array, shape (M, K)
+        Permuted L matrix.
+        K = min(M, N)
+    u : array, shape (K, N)
+        Upper triangular or trapezoidal matrix
+
+    Notes
+    -----
+    This is a LU factorization routine written for Scipy.
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2:
@@ -471,29 +722,60 @@
return p,l,u

def svd(a,full_matrices=1,compute_uv=1,overwrite_a=0):
-    """Compute singular value decomposition (SVD) of matrix a.
+    """Singular Value Decomposition.

-    Description:
+    Factorizes the matrix a into two unitary matrices U and Vh and
+    an 1d-array s of singular values (real, non-negative) such that
+    a == U S Vh  if S is an suitably shaped matrix of zeros whose
+    main diagonal is s.
+
+    Parameters
+    ----------
+    a : array, shape (M, N)
+        Matrix to decompose
+    full_matrices : boolean
+        If true,  U, Vh are shaped  (M,M), (N,N)
+        If false, the shapes are    (M,K), (K,N) where K = min(M,N)
+    compute_uv : boolean
+        Whether to compute also U, Vh in addition to s (Default: true)
+    overwrite_a : boolean
+        Whether data in a is overwritten (may improve performance)
+
+    Returns
+    -------
+    U:  array, shape (M,M) or (M,K) depending on full_matrices
+    s:  array, shape (K,)
+        The singular values, sorted so that s[i] >= s[i+1]. K = min(M, N)
+    Vh: array, shape (N,N) or (K,N) depending on full_matrices

-      Singular value decomposition of a matrix a is
-        a = u * sigma * v^H,
-      where v^H denotes conjugate(transpose(v)), u,v are unitary
-      matrices, sigma is zero matrix with a main diagonal containing
-      real non-negative singular values of the matrix a.
+    For compute_uv = False, only s is returned.

-    Inputs:
+    Raises LinAlgError if SVD computation does not converge

-      a -- An M x N matrix.
-      compute_uv -- If zero, then only the vector of singular values
-                    is returned.
+    Examples
+    --------
+    >>> from scipy import random, linalg, allclose, dot
+    >>> a = random.randn(9, 6) + 1j*random.randn(9, 6)
+    >>> U, s, Vh = linalg.svd(a)
+    >>> U.shape, Vh.shape, s.shape
+    ((9, 9), (6, 6), (6,))
+
+    >>> U, s, Vh = linalg.svd(a, full_matrices=False)
+    >>> U.shape, Vh.shape, s.shape
+    ((9, 6), (6, 6), (6,))
+    >>> S = linalg.diagsvd(s, 6, 6)
+    >>> allclose(a, dot(U, dot(S, Vh)))
+    True
+
+    >>> s2 = linalg.svd(a, compute_uv=False)
+    >>> allclose(s, s2)
+    True

-    Outputs:
-
-      u -- An M x M unitary matrix [compute_uv=1].
-      s -- An min(M,N) vector of singular values in descending order,
-           sigma = diagsvd(s).
-      vh -- An N x N unitary matrix [compute_uv=1], vh = v^H.
-
+    --------
+    svdvals : return singular values of a matrix
+    diagsvd : return the Sigma matrix, given the vector s
+
"""
# A hack until full_matrices == 0 support is fixed here.
if full_matrices == 0:
@@ -520,11 +802,47 @@
return s

def svdvals(a,overwrite_a=0):
-    """Return singular values of a matrix."""
+    """Compute singular values of a matrix.
+
+    Parameters
+    ----------
+    a : array, shape (M, N)
+        Matrix to decompose
+    overwrite_a : boolean
+        Whether data in a is overwritten (may improve performance)
+
+    Returns
+    -------
+    s:  array, shape (K,)
+        The singular values, sorted so that s[i] >= s[i+1]. K = min(M, N)
+
+    Raises LinAlgError if SVD computation does not converge
+
+    --------
+    svd : return the full singular value decomposition of a matrix
+    diagsvd : return the Sigma matrix, given the vector s
+
+    """
return svd(a,compute_uv=0,overwrite_a=overwrite_a)

def diagsvd(s,M,N):
-    """Return sigma from singular values and original size M,N."""
+    """Construct the sigma matrix in SVD from singular values and size M,N.
+
+    Parameters
+    ----------
+    s : array, shape (M,) or (N,)
+        Singular values
+    M : integer
+    N : integer
+        Size of the matrix whose singular values are s
+
+    Returns
+    -------
+    S : array, shape (M, N)
+        The S-matrix in the singular value decomposition
+
+    """
part = diag(s)
typ = part.dtype.char
MorN = len(s)
@@ -536,14 +854,40 @@
raise ValueError, "Length of s must be M or N."

def cholesky(a,lower=0,overwrite_a=0):
-    """Compute Cholesky decomposition of matrix.
-
-    Description:
-
-      For a hermitian positive-definite matrix a return the
-      upper-triangular (or lower-triangular if lower==1) matrix,
-      u such that u^H * u = a (or l * l^H = a).
-
+    """Compute the Cholesky decomposition of a matrix.
+
+    Returns the Cholesky decomposition, :lm:`A = L L^*` or :lm:`A = U^* U`
+    of a Hermitian positive-definite matrix :lm:`A`.
+
+    Parameters
+    ----------
+    a : array, shape (M, M)
+        Matrix to be decomposed
+    lower : boolean
+        Whether to compute the upper or lower triangular Cholesky factorization
+        (Default: upper-triangular)
+    overwrite_a : boolean
+        Whether to overwrite data in a (may improve performance)
+
+    Returns
+    -------
+    B : array, shape (M, M)
+        Upper- or lower-triangular Cholesky factor of A
+
+    Raises LinAlgError if decomposition fails
+
+    Examples
+    --------
+    >>> from scipy import array, linalg, dot
+    >>> a = array([[1,-2j],[2j,5]])
+    >>> L = linalg.cholesky(a, lower=True)
+    >>> L
+    array([[ 1.+0.j,  0.+0.j],
+           [ 0.+2.j,  1.+0.j]])
+    >>> dot(L, L.T.conj())
+    array([[ 1.+0.j,  0.-2.j],
+           [ 0.+2.j,  5.+0.j]])
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -557,8 +901,32 @@
return c

def cho_factor(a, lower=0, overwrite_a=0):
-    """ Compute Cholesky decomposition of matrix and return an object
-    to be used for solving a linear system using cho_solve.
+    """Compute the Cholesky decomposition of a matrix, to use in cho_solve
+
+    Returns the Cholesky decomposition, :lm:`A = L L^*` or :lm:`A = U^* U`
+    of a Hermitian positive-definite matrix :lm:`A`.
+
+    The return value can be directly used as the first parameter to cho_solve.
+
+    Parameters
+    ----------
+    a : array, shape (M, M)
+        Matrix to be decomposed
+    lower : boolean
+        Whether to compute the upper or lower triangular Cholesky factorization
+        (Default: upper-triangular)
+    overwrite_a : boolean
+        Whether to overwrite data in a (may improve performance)
+
+    Returns
+    -------
+    c : array, shape (M, M)
+        Upper- or lower-triangular Cholesky factor of A
+    lower : array, shape (M, M)
+        Flag indicating whether the factor is lower or upper triangular
+
+    Raises LinAlgError if decomposition fails
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -573,8 +941,29 @@

def cho_solve(clow, b):
"""Solve a previously factored symmetric system of equations.
+
+    The equation system is
+
+        A x = b,  A = U^H U = L L^H
+
+    and A is real symmetric or complex Hermitian.
+
+    Parameters
+    ----------
+    clow : tuple (c, lower)
+        Cholesky factor and a flag indicating whether it is lower triangular.
+        The return value from cho_factor can be used.
+    b : array
+        Right-hand side of the equation system
+
First input is a tuple (LorU, lower) which is the output to cho_factor.
Second input is the right-hand side.
+
+    Returns
+    -------
+    x : array
+        Solution to the equation system
+
"""
c, lower = clow
c = asarray_chkfinite(c)
@@ -591,29 +980,60 @@
return b

def qr(a,overwrite_a=0,lwork=None,econ=False,mode='qr'):
-    """QR decomposition of an M x N matrix a.
+    """Compute QR decomposition of a matrix.

-    Description:
+    Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
+    and R upper triangular.

-      Find a unitary (orthogonal) matrix, q, and an upper-triangular
-      matrix r such that q * r = a
+    Parameters
+    ----------
+    a : array, shape (M, N)
+        Matrix to be decomposed
+    overwrite_a : boolean
+        Whether data in a is overwritten (may improve performance)
+    lwork : integer
+        Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
+        is computed.
+    econ : boolean
+        Whether to compute the economy-size QR decomposition, making shapes
+        of Q and R (M, K) and (K, N) instead of (M,M) and (M,N). K=min(M,N)
+    mode : {'qr', 'r'}
+        Determines what information is to be returned: either both Q and R
+        or only R.
+
+    Returns
+    -------
+    (if mode == 'qr')
+    Q : double or complex array, shape (M, M) or (M, K) for econ==True

-    Inputs:
+    (for any mode)
+    R : double or complex array, shape (M, N) or (K, N) for econ==True
+        Size K = min(M, N)

-      a -- the matrix
-      overwrite_a=0 -- if non-zero then discard the contents of a,
-                     i.e. a is used as a work array if possible.
+    Raises LinAlgError if decomposition fails

-      lwork=None -- >= shape(a)[1]. If None (or -1) compute optimal
-                    work array size.
-      econ=False -- computes the skinny or economy-size QR decomposition
-                    only useful when M>N
-      mode='qr' -- if 'qr' then return both q and r; if 'r' then just return r
+    Notes
+    -----
+    This is an interface to the LAPACK routines dgeqrf, zgeqrf,
+    dorgqr, and zungqr.

-    Outputs:
-      q,r  - if mode=='qr'
-      r    - if mode=='r'
+    Examples
+    --------
+    >>> from scipy import random, linalg, dot
+    >>> a = random.randn(9, 6)
+    >>> q, r = linalg.qr(a)
+    >>> allclose(a, dot(q, r))
+    True
+    >>> q.shape, r.shape
+    ((9, 9), (9, 6))

+    >>> r2 = linalg.qr(a, mode='r')
+    >>> allclose(r, r2)
+
+    >>> q3, r3 = linalg.qr(a, econ=True)
+    >>> q3.shape, r3.shape
+    ((9, 6), (6, 6))
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2:
@@ -674,26 +1094,29 @@

def qr_old(a,overwrite_a=0,lwork=None):
-    """QR decomposition of an M x N matrix a.
+    """Compute QR decomposition of a matrix.

-    Description:
+    Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
+    and R upper triangular.

-      Find a unitary (orthogonal) matrix, q, and an upper-triangular
-      matrix r such that q * r = a
+    Parameters
+    ----------
+    a : array, shape (M, N)
+        Matrix to be decomposed
+    overwrite_a : boolean
+        Whether data in a is overwritten (may improve performance)
+    lwork : integer
+        Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
+        is computed.
+
+    Returns
+    -------
+    Q : double or complex array, shape (M, M)
+    R : double or complex array, shape (M, N)
+        Size K = min(M, N)

-    Inputs:
-
-      a -- the matrix
-      overwrite_a=0 -- if non-zero then discard the contents of a,
-                     i.e. a is used as a work array if possible.
-
-      lwork=None -- >= shape(a)[1]. If None (or -1) compute optimal
-                    work array size.
-
-    Outputs:
-
-      q, r -- matrices such that q * r = a
-
+    Raises LinAlgError if decomposition fails
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2:
@@ -725,26 +1148,30 @@

def rq(a,overwrite_a=0,lwork=None):
-    """RQ decomposition of an M x N matrix a.
+    """Compute RQ decomposition of a square real matrix.

-    Description:
+    Calculate the decomposition :lm:`A = R Q` where Q is unitary/orthogonal
+    and R upper triangular.

-      Find an upper-triangular matrix r and a unitary (orthogonal)
-      matrix q such that r * q = a
-
-    Inputs:
-
-      a -- the matrix
-      overwrite_a=0 -- if non-zero then discard the contents of a,
-                     i.e. a is used as a work array if possible.
-
-      lwork=None -- >= shape(a)[1]. If None (or -1) compute optimal
-                    work array size.
-
-    Outputs:
-
-      r, q -- matrices such that r * q = a
-
+    Parameters
+    ----------
+    a : array, shape (M, M)
+        Square real matrix to be decomposed
+    overwrite_a : boolean
+        Whether data in a is overwritten (may improve performance)
+    lwork : integer
+        Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
+        is computed.
+    econ : boolean
+
+    Returns
+    -------
+    R : double array, shape (M, N) or (K, N) for econ==True
+        Size K = min(M, N)
+    Q : double or complex array, shape (M, M) or (M, K) for econ==True
+
+    Raises LinAlgError if decomposition fails
+
"""
# TODO: implement support for non-square and complex arrays
a1 = asarray_chkfinite(a)
@@ -783,13 +1210,40 @@
_double_precision = ['i','l','d']

def schur(a,output='real',lwork=None,overwrite_a=0):
-    """Compute Schur decomposition of matrix a.
+    """Compute Schur decomposition of a matrix.
+
+    The Schur decomposition is
+
+        A = Z T Z^H
+
+    where Z is unitary and T is either upper-triangular, or for real
+    Schur decomposition (output='real'), quasi-upper triangular.  In
+    the quasi-triangular form, 2x2 blocks describing complex-valued
+    eigenvalue pairs may extrude from the diagonal.
+
+    Parameters
+    ----------
+    a : array, shape (M, M)
+        Matrix to decompose
+    output : {'real', 'complex'}
+        Construct the real or complex Schur decomposition (for real matrices).
+    lwork : integer
+        Work array size. If None or -1, it is automatically computed.
+    overwrite_a : boolean
+        Whether to overwrite data in a (may improve performance)

-    Description:
+    Returns
+    -------
+    T : array, shape (M, M)
+        Schur form of A. It is real-valued for the real Schur decomposition.
+    Z : array, shape (M, M)
+        An unitary Schur transformation matrix for A.
+        It is real-valued for the real Schur decomposition.

-      Return T, Z such that a = Z * T * (Z**H) where Z is a
-      unitary matrix and T is either upper-triangular or quasi-upper
-      triangular for output='real'
+    --------
+    rsf2csf : Convert real Schur form to complex Schur form
+
"""
if not output in ['real','complex','r','c']:
raise ValueError, "argument must be 'real', or 'complex'"
@@ -850,16 +1304,29 @@
raise LinAlgError, 'Array must be square'

def rsf2csf(T, Z):
-    """Convert real schur form to complex schur form.
-
-    Description:
-
-      If A is a real-valued matrix, then the real schur form is
-      quasi-upper triangular.  2x2 blocks extrude from the main-diagonal
-      corresponding to any complex-valued eigenvalues.
-
-      This function converts this real schur form to a complex schur form
-      which is upper triangular.
+    """Convert real Schur form to complex Schur form.
+
+    Convert a quasi-diagonal real-valued Schur form to the upper triangular
+    complex-valued Schur form.
+
+    Parameters
+    ----------
+    T : array, shape (M, M)
+        Real Schur form of the original matrix
+    Z : array, shape (M, M)
+        Schur transformation matrix
+
+    Returns
+    -------
+    T : array, shape (M, M)
+        Complex Schur form of the original matrix
+    Z : array, shape (M, M)
+        Schur transformation matrix corresponding to the complex form
+
+    --------
+    schur : Schur decompose a matrix
+
"""
Z,T = map(asarray_chkfinite, (Z,T))
if len(Z.shape) !=2 or Z.shape[0] != Z.shape[1]:
@@ -898,7 +1365,23 @@
# Orthonormal decomposition

def orth(A):
-    """Return an orthonormal basis for the range of A using svd"""
+    """Construct an orthonormal basis for the range of A using SVD
+
+    Parameters
+    ----------
+    A : array, shape (M, N)
+
+    Returns
+    -------
+    Q : array, shape (M, K)
+        Orthonormal basis for the range of A.
+        K = effective rank of A, as determined by automatic cutoff
+
+    --------
+    svd : Singular value decomposition of a matrix
+
+    """
u,s,vh = svd(A)
M,N = A.shape
tol = max(M,N)*numpy.amax(s)*eps
@@ -907,21 +1390,33 @@
return Q

def hessenberg(a,calc_q=0,overwrite_a=0):
-    """ Compute Hessenberg form of a matrix.
+    """Compute Hessenberg form of a matrix.
+
+    The Hessenberg decomposition is
+
+        A = Q H Q^H
+
+    where Q is unitary/orthogonal and H has only zero elements below the first
+    subdiagonal.
+
+    Parameters
+    ----------
+    a : array, shape (M,M)
+        Matrix to bring into Hessenberg form
+    calc_q : boolean
+        Whether to compute the transformation matrix
+    overwrite_a : boolean
+        Whether to ovewrite data in a (may improve performance)

-    Inputs:
+    Returns
+    -------
+    H : array, shape (M,M)
+        Hessenberg form of A

-      a -- the matrix
-      calc_q -- if non-zero then calculate unitary similarity
-                transformation matrix q.
-      overwrite_a=0 -- if non-zero then discard the contents of a,
-                     i.e. a is used as a work array if possible.
+    (If calc_q == True)
+    Q : array, shape (M,M)
+        Unitary/orthogonal similarity transformation matrix s.t. A = Q H Q^H

-    Outputs:
-
-      h    -- Hessenberg form of a                [calc_q=0]
-      h, q -- matrices such that a = q * h * q^T  [calc_q=1]
-
"""
a1 = asarray(a)
if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):

Modified: trunk/scipy/linalg/matfuncs.py
===================================================================
--- trunk/scipy/linalg/matfuncs.py	2008-02-20 05:39:03 UTC (rev 3951)
+++ trunk/scipy/linalg/matfuncs.py	2008-02-21 01:57:37 UTC (rev 3952)
@@ -20,7 +20,20 @@
feps = sb.finfo(single).eps

def expm(A,q=7):
-    """Compute the matrix exponential using Pade approximation of order q.
+    """Compute the matrix exponential using Pade approximation.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+        Matrix to be exponentiated
+    q : integer
+        Order of the Pade approximation
+
+    Returns
+    -------
+    expA : array, shape(M,M)
+        Matrix exponential of A
+
"""
A = asarray(A)
ss = True
@@ -61,6 +74,17 @@

def expm2(A):
"""Compute the matrix exponential using eigenvalue decomposition.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+        Matrix to be exponentiated
+
+    Returns
+    -------
+    expA : array, shape(M,M)
+        Matrix exponential of A
+
"""
A = asarray(A)
t = A.dtype.char
@@ -72,7 +96,20 @@
return dot(dot(vr,diag(exp(s))),vri).astype(t)

def expm3(A,q=20):
-    """Compute the matrix exponential using a Taylor series.of order q.
+    """Compute the matrix exponential using Taylor series.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+        Matrix to be exponentiated
+    q : integer
+        Order of the Taylor series
+
+    Returns
+    -------
+    expA : array, shape(M,M)
+        Matrix exponential of A
+
"""
A = asarray(A)
t = A.dtype.char
@@ -91,6 +128,16 @@
_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
def toreal(arr,tol=None):
"""Return as real array if imaginary part is small.
+
+    Parameters
+    ----------
+    arr : array
+    tol : float
+        Absolute tolerance
+
+    Returns
+    -------
+    arr : double or complex array
"""
if tol is None:
tol = {0:feps*1e3, 1:eps*1e6}[_array_precision[arr.dtype.char]]
@@ -100,7 +147,19 @@
return arr

def cosm(A):
-    """matrix cosine.
+    """Compute the matrix cosine.
+
+    This routine uses expm to compute the matrix exponentials.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+
+    Returns
+    -------
+    cosA : array, shape(M,M)
+        Matrix cosine of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D','G']:
@@ -110,7 +169,19 @@

def sinm(A):
-    """matrix sine.
+    """Compute the matrix sine.
+
+    This routine uses expm to compute the matrix exponentials.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+
+    Returns
+    -------
+    sinA : array, shape(M,M)
+        Matrix cosine of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D','G']:
@@ -119,7 +190,19 @@
return -0.5j*(expm(1j*A) - expm(-1j*A))

def tanm(A):
-    """matrix tangent.
+    """Compute the matrix tangent.
+
+    This routine uses expm to compute the matrix exponentials.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+
+    Returns
+    -------
+    tanA : array, shape(M,M)
+        Matrix tangent of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D','G']:
@@ -128,7 +211,19 @@
return solve(cosm(A), sinm(A))

def coshm(A):
-    """matrix hyperbolic cosine.
+    """Compute the hyperbolic matrix cosine.
+
+    This routine uses expm to compute the matrix exponentials.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+
+    Returns
+    -------
+    coshA : array, shape(M,M)
+        Hyperbolic matrix cosine of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D','G']:
@@ -137,7 +232,19 @@
return 0.5*(expm(A) + expm(-A))

def sinhm(A):
-    """matrix hyperbolic sine.
+    """Compute the hyperbolic matrix sine.
+
+    This routine uses expm to compute the matrix exponentials.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+
+    Returns
+    -------
+    sinhA : array, shape(M,M)
+        Hyperbolic matrix sine of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D']:
@@ -146,7 +253,19 @@
return 0.5*(expm(A) - expm(-A))

def tanhm(A):
-    """matrix hyperbolic tangent.
+    """Compute the hyperbolic matrix tangent.
+
+    This routine uses expm to compute the matrix exponentials.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+
+    Returns
+    -------
+    tanhA : array, shape(M,M)
+        Hyperbolic matrix tangent of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D']:
@@ -155,11 +274,32 @@
return solve(coshm(A), sinhm(A))

def funm(A,func,disp=1):
-    """matrix function for arbitrary callable object func.
+    """Evaluate a matrix function specified by a callable.
+
+    Returns the value of matrix-valued function f at A. The function f
+    is an extension of the scalar-valued function func to matrices.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+        Matrix at which to evaluate the function
+    func : callable
+        Callable object that evaluates a scalar function f.
+        Must be vectorized (eg. using vectorize).
+    disp : boolean
+        Print warning if error in the result is estimated large
+        instead of returning estimated error. (Default: True)
+
+    Returns
+    -------
+    fA : array, shape(M,M)
+        Value of the matrix function specified by func evaluated at A
+
+    (if disp == False)
+    errest : float
+        1-norm of the estimated error, ||err||_1 / ||A||_1
+
"""
-    # func should take a vector of arguments (see vectorize if
-    #  it needs wrapping.
-
# Perform Shur decomposition (lapack ?gees)
A = asarray(A)
if len(A.shape)!=2:
@@ -209,7 +349,28 @@
return F, err

def logm(A,disp=1):
-    """Matrix logarithm, inverse of expm."""
+    """Compute matrix logarithm.
+
+    The matrix logarithm is the inverse of expm: expm(logm(A)) == A
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+        Matrix whose logarithm to evaluate
+    disp : boolean
+        Print warning if error in the result is estimated large
+        instead of returning estimated error. (Default: True)
+
+    Returns
+    -------
+    logA : array, shape(M,M)
+        Matrix logarithm of A
+
+    (if disp == False)
+    errest : float
+        1-norm of the estimated error, ||err||_1 / ||A||_1
+
+    """
# Compute using general funm but then use better error estimator and
#   make one step in improving estimate using a rotation matrix.
A = mat(asarray(A))
@@ -239,7 +400,37 @@
return F, errest

def signm(a,disp=1):
-    """matrix sign"""
+    """Matrix sign function.
+
+    Extension of the scalar sign(x) to matrices.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+        Matrix at which to evaluate the sign function
+    disp : boolean
+        Print warning if error in the result is estimated large
+        instead of returning estimated error. (Default: True)
+
+    Returns
+    -------
+    sgnA : array, shape(M,M)
+        Value of the sign function at A
+
+    (if disp == False)
+    errest : float
+        1-norm of the estimated error, ||err||_1 / ||A||_1
+
+    Examples
+    --------
+    >>> from scipy.linalg import signm, eigvals
+    >>> a = [[1,2,3], [1,2,1], [1,1,1]]
+    >>> eigvals(a)
+    array([ 4.12488542+0.j, -0.76155718+0.j,  0.63667176+0.j])
+    >>> eigvals(signm(a))
+    array([-1.+0.j,  1.+0.j,  1.+0.j])
+
+    """
def rounded_sign(x):
rx = real(x)
if rx.dtype.char=='f':
@@ -286,12 +477,29 @@
return S0, errest

def sqrtm(A,disp=1):
-    """Matrix square root
+    """Matrix square root.
+
+    Parameters
+    ----------
+    A : array, shape(M,M)
+        Matrix whose square root to evaluate
+    disp : boolean
+        Print warning if error in the result is estimated large
+        instead of returning estimated error. (Default: True)
+
+    Returns
+    -------
+    sgnA : array, shape(M,M)
+        Value of the sign function at A

-    If disp is non-zero display warning if singular matrix.
-    If disp is zero then return residual ||A-X*X||_F / ||A||_F
+    (if disp == False)
+    errest : float
+        Frobenius norm of the estimated error, ||err||_F / ||A||_F

+    Notes
+    -----
Uses algorithm by Nicholas J. Higham
+
"""
A = asarray(A)
if len(A.shape)!=2:

```