[Scipy-svn] r6319 - in trunk: doc/source scipy/linalg scipy/linalg/tests

[scipy-svn at scipy.org]

Author: warren.weckesser
Date: 2010-04-09 20:47:20 -0500 (Fri, 09 Apr 2010)
New Revision: 6319

Added:
   trunk/scipy/linalg/decomp_cholesky.py
   trunk/scipy/linalg/decomp_lu.py
   trunk/scipy/linalg/decomp_qr.py
   trunk/scipy/linalg/decomp_schur.py
   trunk/scipy/linalg/decomp_svd.py
Modified:
   trunk/doc/source/linalg.rst
   trunk/scipy/linalg/__init__.py
   trunk/scipy/linalg/basic.py
   trunk/scipy/linalg/decomp.py
   trunk/scipy/linalg/generic_flapack.pyf
   trunk/scipy/linalg/info.py
   trunk/scipy/linalg/matfuncs.py
   trunk/scipy/linalg/misc.py
   trunk/scipy/linalg/tests/test_basic.py
   trunk/scipy/linalg/tests/test_decomp.py
Log:
linalg:
REF: Split basic.py and decomp.py into multiple modules.
     Removed duplicated versions of the functions cho_solve and lu_solve.
     Return value of solveh_banded no longer includes the Choleskey factor
     (ticket #676).
ENH: Added cho_solve_banded.
STY: Clean up python style in many places.


Modified: trunk/doc/source/linalg.rst
===================================================================
--- trunk/doc/source/linalg.rst	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/doc/source/linalg.rst	2010-04-10 01:47:20 UTC (rev 6319)
@@ -20,8 +20,8 @@
    pinv
    pinv2
 
-Eigenvalues and Decompositions
-==============================
+Eigenvalue Problem
+==================
 
 .. autosummary::
    :toctree: generated/
@@ -32,6 +32,13 @@
    eigvalsh
    eig_banded
    eigvals_banded
+
+Decompositions
+==============
+
+.. autosummary::
+   :toctree: generated/
+
    lu
    lu_factor
    lu_solve
@@ -43,6 +50,7 @@
    cholesky_banded
    cho_factor
    cho_solve
+   cho_solve_banded
    qr
    schur
    rsf2csf

Modified: trunk/scipy/linalg/__init__.py
===================================================================
--- trunk/scipy/linalg/__init__.py	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/scipy/linalg/__init__.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -5,12 +5,19 @@
 from info import __doc__
 from linalg_version import linalg_version as __version__
 
+from misc import *
 from basic import *
 from decomp import *
+from decomp_lu import *
+from decomp_cholesky import *
+from decomp_qr import *
+from decomp_svd import *
+from decomp_schur import *
 from matfuncs import *
 from blas import *
+from special_matrices import *
 
-__all__ = filter(lambda s:not s.startswith('_'),dir())
+__all__ = filter(lambda s: not s.startswith('_'), dir())
 
 from numpy.dual import register_func
 for k in ['norm', 'inv', 'svd', 'solve', 'det', 'eig', 'eigh', 'eigvals',

Modified: trunk/scipy/linalg/basic.py
===================================================================
--- trunk/scipy/linalg/basic.py	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/scipy/linalg/basic.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -1,108 +1,22 @@
-## Automatically adapted for scipy Oct 18, 2005 by
-
-## Automatically adapted for scipy Oct 18, 2005 by
-
 #
 # Author: Pearu Peterson, March 2002
 #
 # w/ additions by Travis Oliphant, March 2002
 
-__all__ = ['solve', 'inv', 'det', 'lstsq', 'pinv', 'pinv2',
-           'cholesky_banded', 'solveh_banded', 'lu_solve', 'cho_solve',
-           'solve_banded',
-           #  From special_matrices:
-           'tri','tril', 'triu', 'toeplitz', 'circulant', 'hankel', 'kron',
-           'hadamard', 'block_diag', 'all_mat',
-           #  From misc:
-           'LinAlgError', 'norm',
-            ]
+__all__ = ['solve', 'solveh_banded', 'solve_banded',
+            'inv', 'det', 'lstsq', 'pinv', 'pinv2']
 
 from numpy import asarray, zeros, sum, conjugate, dot, transpose, \
         asarray_chkfinite,  single
 import numpy
 
-#from blas import get_blas_funcs
 from flinalg import get_flinalg_funcs
 from lapack import get_lapack_funcs
-from misc import LinAlgError, norm
-from special_matrices import tri, tril, triu, toeplitz, circulant, hankel, \
-        hadamard, kron, block_diag, all_mat
+from misc import LinAlgError
 from scipy.linalg import calc_lwork
-import decomp
+import decomp_svd
 
 
-def lu_solve((lu, piv), b, trans=0, overwrite_b=False):
-    """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:
-
-        =====  =========
-        trans  system
-        =====  =========
-        0      a x   = b
-        1      a^T x = b
-        2      a^H x = b
-        =====  =========
-
-    Returns
-    -------
-    x : array
-        Solution to the system
-
-    See also
-    --------
-    lu_factor : LU factorize a matrix
-
-    """
-    b1 = asarray_chkfinite(b)
-    overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
-    if lu.shape[0] != b1.shape[0]:
-        raise ValueError, "incompatible dimensions."
-    getrs, = get_lapack_funcs(('getrs',),(lu,b1))
-    x,info = getrs(lu,piv,b1,trans=trans,overwrite_b=overwrite_b)
-    if info==0:
-        return x
-    raise ValueError(
-          'illegal value in %d-th argument of internal gesv|posv' % (-info))
-
-def cho_solve((c, lower), b, overwrite_b=False):
-    """Solve an equation system, a x = b, given the Cholesky factorization of a
-
-    Parameters
-    ----------
-    (c, lower)
-        Cholesky factorization of a, as given by cho_factor
-    b : array
-        Right-hand side
-
-    Returns
-    -------
-    x : array
-        The solution to the system a x = b
-
-    See also
-    --------
-    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__'))
-    if c.shape[0] != b1.shape[0]:
-        raise ValueError, "incompatible dimensions."
-    potrs, = get_lapack_funcs(('potrs',),(c,b1))
-    x,info = potrs(c,b1,lower=lower,overwrite_b=overwrite_b)
-    if info==0:
-        return x
-    raise ValueError(
-          'illegal value in %d-th argument of internal gesv|posv' % (-info))
-
 # Linear equations
 def solve(a, b, sym_pos=False, lower=False, overwrite_a=False, overwrite_b=False,
           debug=False):
@@ -141,25 +55,23 @@
         print 'solve:overwrite_a=',overwrite_a
         print 'solve:overwrite_b=',overwrite_b
     if sym_pos:
-        posv, = get_lapack_funcs(('posv',),(a1,b1))
-        c,x,info = posv(a1,b1,
-                        lower = lower,
+        posv, = get_lapack_funcs(('posv',), (a1,b1))
+        c, x, info = posv(a1, b1, lower=lower,
                         overwrite_a=overwrite_a,
                         overwrite_b=overwrite_b)
     else:
-        gesv, = get_lapack_funcs(('gesv',),(a1,b1))
-        lu,piv,x,info = gesv(a1,b1,
-                             overwrite_a=overwrite_a,
-                             overwrite_b=overwrite_b)
+        gesv, = get_lapack_funcs(('gesv',), (a1,b1))
+        lu, piv, x, info = gesv(a1, b1, overwrite_a=overwrite_a,
+                                            overwrite_b=overwrite_b)
 
-    if info==0:
+    if info == 0:
         return x
-    if info>0:
-        raise LinAlgError, "singular matrix"
-    raise ValueError(
-          'illegal value in %d-th argument of internal gesv|posv' % (-info))
+    if info > 0:
+        raise LinAlgError("singular matrix")
+    raise ValueError('illegal value in %d-th argument of internal gesv|posv'
+                                                                    % -info)
 
-def solve_banded((l,u), ab, b, overwrite_ab=False, overwrite_b=False,
+def solve_banded((l, u), ab, b, overwrite_ab=False, overwrite_b=False,
           debug=False):
     """Solve the equation a x = b for x, assuming a is banded matrix.
 
@@ -193,32 +105,29 @@
         The solution to the system a x = b
 
     """
-    a1, b1 = map(asarray_chkfinite,(ab,b))
+    a1, b1 = map(asarray_chkfinite, (ab, b))
 
     # Validate shapes.
     if a1.shape[-1] != b1.shape[0]:
         raise ValueError("shapes of ab and b are not compatible.")
-    if l+u+1 != a1.shape[0]:
+    if l + u + 1 != a1.shape[0]:
         raise ValueError("invalid values for the number of lower and upper diagonals:"
                 " l+u+1 (%d) does not equal ab.shape[0] (%d)" % (l+u+1, ab.shape[0]))
 
     overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
 
-    gbsv, = get_lapack_funcs(('gbsv',),(a1,b1))
-    a2 = zeros((2*l+u+1,a1.shape[1]), dtype=gbsv.dtype)
+    gbsv, = get_lapack_funcs(('gbsv',), (a1, b1))
+    a2 = zeros((2*l+u+1, a1.shape[1]), dtype=gbsv.dtype)
     a2[l:,:] = a1
-    lu,piv,x,info = gbsv(l,u,a2,b1,
-                         overwrite_ab=1,
-                         overwrite_b=overwrite_b)
-    if info==0:
+    lu, piv, x, info = gbsv(l, u, a2, b1, overwrite_ab=True,
+                                                overwrite_b=overwrite_b)
+    if info == 0:
         return x
-    if info>0:
-        raise LinAlgError, "singular matrix"
-    raise ValueError(
-          'illegal value in %d-th argument of internal gbsv' % (-info))
+    if info > 0:
+        raise LinAlgError("singular matrix")
+    raise ValueError('illegal value in %d-th argument of internal gbsv' % -info)
 
-def solveh_banded(ab, b, overwrite_ab=False, overwrite_b=False,
-                  lower=False):
+def solveh_banded(ab, b, overwrite_ab=False, overwrite_b=False, lower=False):
     """Solve equation a x = b. a is Hermitian positive-definite banded matrix.
 
     The matrix a is stored in ab either in lower diagonal or upper
@@ -256,81 +165,27 @@
 
     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))
+    ab, b = map(asarray_chkfinite, (ab, b))
 
     # Validate shapes.
     if ab.shape[-1] != b.shape[0]:
         raise ValueError("shapes of ab and b are not compatible.")
 
-    pbsv, = get_lapack_funcs(('pbsv',),(ab,b))
-    c,x,info = pbsv(ab,b,
-                    lower=lower,
-                    overwrite_ab=overwrite_ab,
-                    overwrite_b=overwrite_b)
-    if info==0:
-        return c, x
-    if info>0:
-        raise LinAlgError, "%d-th leading minor not positive definite" % info
-    raise ValueError(
-          'illegal value in %d-th argument of internal pbsv' % (-info))
+    pbsv, = get_lapack_funcs(('pbsv',), (ab, b))
+    c, x, info = pbsv(ab, b, lower=lower, overwrite_ab=overwrite_ab,
+                                            overwrite_b=overwrite_b)
+    if info > 0:
+        raise LinAlgError("%d-th leading minor not positive definite" % info)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal pbsv'
+                                                                    % -info)
+    return x
 
-def cholesky_banded(ab, overwrite_ab=False, lower=False):
-    """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 (u + 1, M)
-        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 (u+1, M)
-        Cholesky factorization of a, in the same banded format as ab
-
-    """
-    ab = asarray_chkfinite(ab)
-
-    pbtrf, = get_lapack_funcs(('pbtrf',),(ab,))
-    c,info = pbtrf(ab,
-                   lower=lower,
-                   overwrite_ab=overwrite_ab)
-
-    if info==0:
-        return c
-    if info>0:
-        raise LinAlgError, "%d-th leading minor not positive definite" % info
-    raise ValueError(
-          'illegal value in %d-th argument of internal pbtrf' % (-info))
-
-
 # matrix inversion
 def inv(a, overwrite_a=False):
     """Compute the inverse of a matrix.
@@ -360,7 +215,7 @@
     """
     a1 = asarray_chkfinite(a)
     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
-        raise ValueError, 'expected square matrix'
+        raise ValueError('expected square matrix')
     overwrite_a = overwrite_a or (a1 is not a and not hasattr(a,'__array__'))
     #XXX: I found no advantage or disadvantage of using finv.
 ##     finv, = get_flinalg_funcs(('inv',),(a1,))
@@ -371,21 +226,21 @@
 ##         if info>0: raise LinAlgError, "singular matrix"
 ##         if info<0: raise ValueError,\
 ##            'illegal value in %d-th argument of internal inv.getrf|getri'%(-info)
-    getrf,getri = get_lapack_funcs(('getrf','getri'),(a1,))
+    getrf, getri = get_lapack_funcs(('getrf','getri'), (a1,))
     #XXX: C ATLAS versions of getrf/i have rowmajor=1, this could be
     #     exploited for further optimization. But it will be probably
     #     a mess. So, a good testing site is required before trying
     #     to do that.
-    if getrf.module_name[:7]=='clapack'!=getri.module_name[:7]:
+    if getrf.module_name[:7] == 'clapack' != getri.module_name[:7]:
         # ATLAS 3.2.1 has getrf but not getri.
-        lu,piv,info = getrf(transpose(a1),
-                            rowmajor=0,overwrite_a=overwrite_a)
+        lu, piv, info = getrf(transpose(a1), rowmajor=0,
+                                                overwrite_a=overwrite_a)
         lu = transpose(lu)
     else:
-        lu,piv,info = getrf(a1,overwrite_a=overwrite_a)
-    if info==0:
+        lu, piv, info = getrf(a1, overwrite_a=overwrite_a)
+    if info == 0:
         if getri.module_name[:7] == 'flapack':
-            lwork = calc_lwork.getri(getri.prefix,a1.shape[0])
+            lwork = calc_lwork.getri(getri.prefix, a1.shape[0])
             lwork = lwork[1]
             # XXX: the following line fixes curious SEGFAULT when
             # benchmarking 500x500 matrix inverse. This seems to
@@ -393,14 +248,15 @@
             # minimal (when using lwork[0] instead of lwork[1]) then
             # all tests pass. Further investigation is required if
             # more such SEGFAULTs occur.
-            lwork = int(1.01*lwork)
-            inv_a,info = getri(lu,piv,
-                               lwork=lwork,overwrite_lu=1)
+            lwork = int(1.01 * lwork)
+            inv_a, info = getri(lu, piv, lwork=lwork, overwrite_lu=1)
         else: # clapack
-            inv_a,info = getri(lu,piv,overwrite_lu=1)
-    if info>0: raise LinAlgError, "singular matrix"
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal getrf|getri' % (-info))
+            inv_a, info = getri(lu, piv, overwrite_lu=1)
+    if info > 0:
+        raise LinAlgError("singular matrix")
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal '
+                                                    'getrf|getri' % -info)
     return inv_a
 
 
@@ -424,12 +280,13 @@
     """
     a1 = asarray_chkfinite(a)
     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
-        raise ValueError, 'expected square matrix'
+        raise ValueError('expected square matrix')
     overwrite_a = overwrite_a or (a1 is not a and not hasattr(a,'__array__'))
-    fdet, = get_flinalg_funcs(('det',),(a1,))
-    a_det,info = fdet(a1,overwrite_a=overwrite_a)
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal det.getrf' % (-info))
+    fdet, = get_flinalg_funcs(('det',), (a1,))
+    a_det, info = fdet(a1, overwrite_a=overwrite_a)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal '
+                                                        'det.getrf' % -info)
     return a_det
 
 ### Linear Least Squares
@@ -468,41 +325,47 @@
     Raises LinAlgError if computation does not converge
 
     """
-    a1, b1 = map(asarray_chkfinite,(a,b))
+    a1, b1 = map(asarray_chkfinite, (a, b))
     if len(a1.shape) != 2:
         raise ValueError, 'expected matrix'
-    m,n = a1.shape
-    if len(b1.shape)==2: nrhs = b1.shape[1]
-    else: nrhs = 1
+    m, n = a1.shape
+    if len(b1.shape) == 2:
+        nrhs = b1.shape[1]
+    else:
+        nrhs = 1
     if m != b1.shape[0]:
-        raise ValueError, 'incompatible dimensions'
-    gelss, = get_lapack_funcs(('gelss',),(a1,b1))
-    if n>m:
+        raise ValueError('incompatible dimensions')
+    gelss, = get_lapack_funcs(('gelss',), (a1, b1))
+    if n > m:
         # need to extend b matrix as it will be filled with
         # a larger solution matrix
-        b2 = zeros((n,nrhs), dtype=gelss.dtype)
-        if len(b1.shape)==2: b2[:m,:] = b1
-        else: b2[:m,0] = b1
+        b2 = zeros((n, nrhs), dtype=gelss.dtype)
+        if len(b1.shape) == 2:
+            b2[:m,:] = b1
+        else:
+            b2[:m,0] = b1
         b1 = b2
     overwrite_a = overwrite_a or (a1 is not a and not hasattr(a,'__array__'))
     overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
     if gelss.module_name[:7] == 'flapack':
-        lwork = calc_lwork.gelss(gelss.prefix,m,n,nrhs)[1]
-        v,x,s,rank,info = gelss(a1,b1,cond = cond,
-                                lwork = lwork,
-                                overwrite_a = overwrite_a,
-                                overwrite_b = overwrite_b)
+        lwork = calc_lwork.gelss(gelss.prefix, m, n, nrhs)[1]
+        v, x, s, rank, info = gelss(a1, b1, cond=cond, lwork=lwork,
+                                                overwrite_a=overwrite_a,
+                                                overwrite_b=overwrite_b)
     else:
-        raise NotImplementedError,'calling gelss from %s' % (gelss.module_name)
-    if info>0: raise LinAlgError, "SVD did not converge in Linear Least Squares"
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal gelss' % (-info))
+        raise NotImplementedError('calling gelss from %s' % gelss.module_name)
+    if info > 0:
+        raise LinAlgError("SVD did not converge in Linear Least Squares")
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal gelss'
+                                                                    % -info)
     resids = asarray([], dtype=x.dtype)
-    if n<m:
+    if n < m:
         x1 = x[:n]
-        if rank==n: resids = sum(x[n:]**2,axis=0)
+        if rank == n:
+            resids = sum(x[n:]**2, axis=0)
         x = x1
-    return x,resids,rank,s
+    return x, resids, rank, s
 
 
 def pinv(a, cond=None, rcond=None):
@@ -544,11 +407,6 @@
     return lstsq(a, b, cond=cond)[0]
 
 
-eps = numpy.finfo(float).eps
-feps = numpy.finfo(single).eps
-
-_array_precision = {'f': 0, 'd': 1, 'F': 0, 'D': 1}
-
 def pinv2(a, cond=None, rcond=None):
     """Compute the (Moore-Penrose) pseudo-inverse of a matrix.
 
@@ -585,15 +443,18 @@
 
     """
     a = asarray_chkfinite(a)
-    u, s, vh = decomp.svd(a)
+    u, s, vh = decomp_svd.svd(a)
     t = u.dtype.char
     if rcond is not None:
         cond = rcond
     if cond in [None,-1]:
+        eps = numpy.finfo(float).eps
+        feps = numpy.finfo(single).eps
+        _array_precision = {'f': 0, 'd': 1, 'F': 0, 'D': 1}
         cond = {0: feps*1e3, 1: eps*1e6}[_array_precision[t]]
-    m,n = a.shape
+    m, n = a.shape
     cutoff = cond*numpy.maximum.reduce(s)
-    psigma = zeros((m,n),t)
+    psigma = zeros((m, n), t)
     for i in range(len(s)):
         if s[i] > cutoff:
             psigma[i,i] = 1.0/conjugate(s[i])

Modified: trunk/scipy/linalg/decomp.py
===================================================================
--- trunk/scipy/linalg/decomp.py	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/scipy/linalg/decomp.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -1,5 +1,3 @@
-## Automatically adapted for scipy Oct 18, 2005 by
-
 #
 # Author: Pearu Peterson, March 2002
 #
@@ -9,83 +7,70 @@
 # additions by Bart Vandereycken, June 2006
 # additions by Andrew D Straw, May 2007
 # additions by Tiziano Zito, November 2008
-
+#
+# April 2010: Functions for LU, QR, SVD, Schur and Cholesky decompositions were
+# moved to their own files.  Still in this file are functions for eigenstuff
+# and for the Hessenberg form.
+ 
 __all__ = ['eig','eigh','eig_banded','eigvals','eigvalsh', 'eigvals_banded',
-           'lu','svd','svdvals','diagsvd','cholesky','qr','qr_old','rq',
-           'schur','rsf2csf','lu_factor','cho_factor','cho_solve','orth',
            'hessenberg']
 
-from misc import LinAlgError
-import misc
-import special_matrices
-
-from warnings import warn
-from lapack import get_lapack_funcs, find_best_lapack_type
-from blas import get_blas_funcs
-from flinalg import get_flinalg_funcs
-from scipy.linalg import calc_lwork
 import numpy
 from numpy import array, asarray_chkfinite, asarray, diag, zeros, ones, \
-        single, isfinite, inexact, complexfloating, nonzero, iscomplexobj
+        isfinite, inexact, nonzero, iscomplexobj, cast
 
-cast = numpy.cast
-r_ = numpy.r_
+# Local imports
+from scipy.linalg import calc_lwork
+from misc import LinAlgError, _datanotshared
+from lapack import get_lapack_funcs
+from blas import get_blas_funcs
 
+
 _I = cast['F'](1j)
-def _make_complex_eigvecs(w,vin,cmplx_tcode):
-    v = numpy.array(vin,dtype=cmplx_tcode)
+
+def _make_complex_eigvecs(w, vin, cmplx_tcode):
+    v = numpy.array(vin, dtype=cmplx_tcode)
     #ind = numpy.flatnonzero(numpy.not_equal(w.imag,0.0))
-    ind = numpy.flatnonzero(numpy.logical_and(numpy.not_equal(w.imag,0.0),
+    ind = numpy.flatnonzero(numpy.logical_and(numpy.not_equal(w.imag, 0.0),
                             numpy.isfinite(w)))
-    vnew = numpy.zeros((v.shape[0],len(ind)>>1),cmplx_tcode)
-    vnew.real = numpy.take(vin,ind[::2],1)
-    vnew.imag = numpy.take(vin,ind[1::2],1)
+    vnew = numpy.zeros((v.shape[0], len(ind)>>1), cmplx_tcode)
+    vnew.real = numpy.take(vin, ind[::2],1)
+    vnew.imag = numpy.take(vin, ind[1::2],1)
     count = 0
     conj = numpy.conjugate
     for i in range(len(ind)/2):
-        v[:,ind[2*i]] = vnew[:,count]
-        v[:,ind[2*i+1]] = conj(vnew[:,count])
+        v[:, ind[2*i]] = vnew[:, count]
+        v[:, ind[2*i+1]] = conj(vnew[:, count])
         count += 1
     return v
 
-
-
-def _datanotshared(a1,a):
-    if a1 is a:
-        return False
-    else:
-        #try comparing data pointers
-        try:
-            return a1.__array_interface__['data'][0] != a.__array_interface__['data'][0]
-        except:
-            return True
-
-
-def _geneig(a1,b,left,right,overwrite_a,overwrite_b):
+def _geneig(a1, b, left, right, overwrite_a, overwrite_b):
     b1 = asarray(b)
-    overwrite_b = overwrite_b or _datanotshared(b1,b)
+    overwrite_b = overwrite_b or _datanotshared(b1, b)
     if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
-        raise ValueError, 'expected square matrix'
-    ggev, = get_lapack_funcs(('ggev',),(a1,b1))
-    cvl,cvr = left,right
+        raise ValueError('expected square matrix')
+    ggev, = get_lapack_funcs(('ggev',), (a1, b1))
+    cvl, cvr = left, right
     if ggev.module_name[:7] == 'clapack':
-        raise NotImplementedError,'calling ggev from %s' % (ggev.module_name)
-    res = ggev(a1,b1,lwork=-1)
+        raise NotImplementedError('calling ggev from %s' % ggev.module_name)
+    res = ggev(a1, b1, lwork=-1)
     lwork = res[-2][0]
     if ggev.prefix in 'cz':
-        alpha,beta,vl,vr,work,info = ggev(a1,b1,cvl,cvr,lwork,
-                                          overwrite_a,overwrite_b)
+        alpha, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr, lwork,
+                                                    overwrite_a, overwrite_b)
         w = alpha / beta
     else:
-        alphar,alphai,beta,vl,vr,work,info = ggev(a1,b1,cvl,cvr,lwork,
-                                                  overwrite_a,overwrite_b)
-        w = (alphar+_I*alphai)/beta
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal ggev' % (-info))
-    if info>0: raise LinAlgError(
-       "generalized eig algorithm did not converge (info=%d)" % (info))
+        alphar, alphai, beta, vl, vr, work, info = ggev(a1, b1, cvl, cvr, lwork,
+                                                        overwrite_a,overwrite_b)
+        w = (alphar + _I * alphai) / beta
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal ggev'
+                                                                    % -info)
+    if info > 0:
+        raise LinAlgError("generalized eig algorithm did not converge (info=%d)"
+                                                                    % info)
 
-    only_real = numpy.logical_and.reduce(numpy.equal(w.imag,0.0))
+    only_real = numpy.logical_and.reduce(numpy.equal(w.imag, 0.0))
     if not (ggev.prefix in 'cz' or only_real):
         t = w.dtype.char
         if left:
@@ -100,7 +85,7 @@
         return w, vl
     return w, vr
 
-def eig(a,b=None, left=False, right=True, overwrite_a=False, overwrite_b=False):
+def eig(a, b=None, left=False, right=True, overwrite_a=False, overwrite_b=False):
     """Solve an ordinary or generalized eigenvalue problem of a square matrix.
 
     Find eigenvalues w and right or left eigenvectors of a general matrix::
@@ -153,49 +138,50 @@
     a1 = asarray_chkfinite(a)
     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
         raise ValueError('expected square matrix')
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
     if b is not None:
         b = asarray_chkfinite(b)
         if b.shape != a1.shape:
             raise ValueError('a and b must have the same shape')
-        return _geneig(a1,b,left,right,overwrite_a,overwrite_b)
-    geev, = get_lapack_funcs(('geev',),(a1,))
-    compute_vl,compute_vr=left,right
+        return _geneig(a1, b, left, right, overwrite_a, overwrite_b)
+    geev, = get_lapack_funcs(('geev',), (a1,))
+    compute_vl, compute_vr = left, right
     if geev.module_name[:7] == 'flapack':
-        lwork = calc_lwork.geev(geev.prefix,a1.shape[0],
-                                compute_vl,compute_vr)[1]
+        lwork = calc_lwork.geev(geev.prefix, a1.shape[0],
+                                    compute_vl, compute_vr)[1]
         if geev.prefix in 'cz':
-            w,vl,vr,info = geev(a1,lwork = lwork,
-                                compute_vl=compute_vl,
-                                compute_vr=compute_vr,
-                                overwrite_a=overwrite_a)
+            w, vl, vr, info = geev(a1, lwork=lwork,
+                                        compute_vl=compute_vl,
+                                        compute_vr=compute_vr,
+                                        overwrite_a=overwrite_a)
         else:
-            wr,wi,vl,vr,info = geev(a1,lwork = lwork,
-                                    compute_vl=compute_vl,
-                                    compute_vr=compute_vr,
-                                    overwrite_a=overwrite_a)
+            wr, wi, vl, vr, info = geev(a1, lwork=lwork,
+                                        compute_vl=compute_vl,
+                                        compute_vr=compute_vr,
+                                        overwrite_a=overwrite_a)
             t = {'f':'F','d':'D'}[wr.dtype.char]
-            w = wr+_I*wi
+            w = wr + _I * wi
     else: # 'clapack'
         if geev.prefix in 'cz':
-            w,vl,vr,info = geev(a1,
-                                compute_vl=compute_vl,
-                                compute_vr=compute_vr,
-                                overwrite_a=overwrite_a)
-        else:
-            wr,wi,vl,vr,info = geev(a1,
+            w, vl, vr, info = geev(a1,
                                     compute_vl=compute_vl,
                                     compute_vr=compute_vr,
                                     overwrite_a=overwrite_a)
+        else:
+            wr, wi, vl, vr, info = geev(a1,
+                                        compute_vl=compute_vl,
+                                        compute_vr=compute_vr,
+                                        overwrite_a=overwrite_a)
             t = {'f':'F','d':'D'}[wr.dtype.char]
-            w = wr+_I*wi
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal geev' % (-info))
-    if info>0: raise LinAlgError(
-       "eig algorithm did not converge (only eigenvalues "\
-       "with order >=%d have converged)" % (info))
+            w = wr + _I * wi
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal geev'
+                                                                    % -info)
+    if info > 0:
+        raise LinAlgError("eig algorithm did not converge (only eigenvalues "
+                            "with order >= %d have converged)" % info)
 
-    only_real = numpy.logical_and.reduce(numpy.equal(w.imag,0.0))
+    only_real = numpy.logical_and.reduce(numpy.equal(w.imag, 0.0))
     if not (geev.prefix in 'cz' or only_real):
         t = w.dtype.char
         if left:
@@ -281,17 +267,17 @@
     """
     a1 = asarray_chkfinite(a)
     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
-        raise ValueError, 'expected square matrix'
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
+        raise ValueError('expected square matrix')
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
     if iscomplexobj(a1):
         cplx = True
     else:
         cplx = False
     if b is not None:
         b1 = asarray_chkfinite(b)
-        overwrite_b = overwrite_b or _datanotshared(b1,b)
+        overwrite_b = overwrite_b or _datanotshared(b1, b)
         if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
-            raise ValueError, 'expected square matrix'
+            raise ValueError('expected square matrix')
 
         if b1.shape != a1.shape:
             raise ValueError("wrong b dimensions %s, should "
@@ -401,7 +387,7 @@
 
 def eig_banded(a_band, lower=False, eigvals_only=False, overwrite_a_band=False,
                select='a', select_range=None, max_ev = 0):
-    """Solve real symmetric or complex hermetian band matrix eigenvalue problem.
+    """Solve real symmetric or complex hermitian band matrix eigenvalue problem.
 
     Find eigenvalues w and optionally right eigenvectors v of a::
 
@@ -472,34 +458,34 @@
     """
     if eigvals_only or overwrite_a_band:
         a1 = asarray_chkfinite(a_band)
-        overwrite_a_band = overwrite_a_band or (_datanotshared(a1,a_band))
+        overwrite_a_band = overwrite_a_band or (_datanotshared(a1, a_band))
     else:
         a1 = array(a_band)
         if issubclass(a1.dtype.type, inexact) and not isfinite(a1).all():
-            raise ValueError, "array must not contain infs or NaNs"
+            raise ValueError("array must not contain infs or NaNs")
         overwrite_a_band = 1
 
     if len(a1.shape) != 2:
-        raise ValueError, 'expected two-dimensional array'
+        raise ValueError('expected two-dimensional array')
     if select.lower() not in [0, 1, 2, 'a', 'v', 'i', 'all', 'value', 'index']:
-        raise ValueError, 'invalid argument for select'
+        raise ValueError('invalid argument for select')
     if select.lower() in [0, 'a', 'all']:
         if a1.dtype.char in 'GFD':
-            bevd, = get_lapack_funcs(('hbevd',),(a1,))
+            bevd, = get_lapack_funcs(('hbevd',), (a1,))
             # FIXME: implement this somewhen, for now go with builtin values
             # FIXME: calc optimal lwork by calling ?hbevd(lwork=-1)
             #        or by using calc_lwork.f ???
             # lwork = calc_lwork.hbevd(bevd.prefix, a1.shape[0], lower)
             internal_name = 'hbevd'
         else: # a1.dtype.char in 'fd':
-            bevd, = get_lapack_funcs(('sbevd',),(a1,))
+            bevd, = get_lapack_funcs(('sbevd',), (a1,))
             # FIXME: implement this somewhen, for now go with builtin values
             #         see above
             # lwork = calc_lwork.sbevd(bevd.prefix, a1.shape[0], lower)
             internal_name = 'sbevd'
-        w,v,info = bevd(a1, compute_v = not eigvals_only,
-                        lower = lower,
-                        overwrite_ab = overwrite_a_band)
+        w,v,info = bevd(a1, compute_v=not eigvals_only,
+                        lower=lower,
+                        overwrite_ab=overwrite_a_band)
     if select.lower() in [1, 2, 'i', 'v', 'index', 'value']:
         # calculate certain range only
         if select.lower() in [2, 'i', 'index']:
@@ -517,31 +503,33 @@
             max_ev = 1
         # calculate optimal abstol for dsbevx (see manpage)
         if a1.dtype.char in 'fF':  # single precision
-            lamch, = get_lapack_funcs(('lamch',),(array(0, dtype='f'),))
+            lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='f'),))
         else:
-            lamch, = get_lapack_funcs(('lamch',),(array(0, dtype='d'),))
+            lamch, = get_lapack_funcs(('lamch',), (array(0, dtype='d'),))
         abstol = 2 * lamch('s')
         if a1.dtype.char in 'GFD':
-            bevx, = get_lapack_funcs(('hbevx',),(a1,))
+            bevx, = get_lapack_funcs(('hbevx',), (a1,))
             internal_name = 'hbevx'
         else: # a1.dtype.char in 'gfd'
-            bevx, = get_lapack_funcs(('sbevx',),(a1,))
+            bevx, = get_lapack_funcs(('sbevx',), (a1,))
             internal_name = 'sbevx'
         # il+1, iu+1: translate python indexing (0 ... N-1) into Fortran
         # indexing (1 ... N)
         w, v, m, ifail, info = bevx(a1, vl, vu, il+1, iu+1,
-                                    compute_v = not eigvals_only,
-                                    mmax = max_ev,
-                                    range = select, lower = lower,
-                                    overwrite_ab = overwrite_a_band,
+                                    compute_v=not eigvals_only,
+                                    mmax=max_ev,
+                                    range=select, lower=lower,
+                                    overwrite_ab=overwrite_a_band,
                                     abstol=abstol)
         # crop off w and v
         w = w[:m]
         if not eigvals_only:
             v = v[:, :m]
-    if info<0: raise ValueError(
-        'illegal value in %d-th argument of internal %s'%(-info, internal_name))
-    if info>0: raise LinAlgError,"eig algorithm did not converge"
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal %s'
+                                                    % (-info, internal_name))
+    if info > 0:
+        raise LinAlgError("eig algorithm did not converge")
 
     if eigvals_only:
         return w
@@ -580,7 +568,7 @@
     eigh : eigenvalues and eigenvectors of symmetric/Hermitean arrays.
 
     """
-    return eig(a,b=b,left=0,right=0,overwrite_a=overwrite_a)
+    return eig(a, b=b, left=0, right=0, overwrite_a=overwrite_a)
 
 def eigvalsh(a, b=None, lower=True, overwrite_a=False,
              overwrite_b=False, turbo=True, eigvals=None, type=1):
@@ -710,837 +698,13 @@
     eig : eigenvalues and right eigenvectors for non-symmetric arrays
 
     """
-    return eig_banded(a_band,lower=lower,eigvals_only=1,
+    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=False):
-    """Compute pivoted LU decomposition of a matrix.
-
-    The decomposition is::
-
-        A = P L U
-
-    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)
-
-    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].
-
-    See also
-    --------
-    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]):
-        raise ValueError, 'expected square matrix'
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
-    getrf, = get_lapack_funcs(('getrf',),(a1,))
-    lu, piv, info = getrf(a,overwrite_a=overwrite_a)
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal getrf (lu_factor)' % (-info))
-    if info>0: warn("Diagonal number %d is exactly zero. Singular matrix." % info,
-                    RuntimeWarning)
-    return lu, piv
-
-def lu_solve(a_lu_pivots, b):
-    """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
-
-    See also
-    --------
-    lu_factor : LU factorize a matrix
-
-    """
-    a_lu, pivots = a_lu_pivots
-    a_lu = asarray_chkfinite(a_lu)
-    pivots = asarray_chkfinite(pivots)
-    b = asarray_chkfinite(b)
-    _assert_squareness(a_lu)
-
-    getrs, = get_lapack_funcs(('getrs',),(a_lu,))
-    b, info = getrs(a_lu,pivots,b)
-    if info < 0:
-        msg = "Argument %d to lapack's ?getrs() has an illegal value." % info
-        raise TypeError, msg
-    if info > 0:
-        msg = "Unknown error occured int ?getrs(): error code = %d" % info
-        raise TypeError, msg
-    return b
-
-
-def lu(a, permute_l=False, overwrite_a=False):
-    """Compute pivoted LU decompostion of a matrix.
-
-    The decomposition is::
-
-        A = P L U
-
-    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)
-
-    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:
-        raise ValueError, 'expected matrix'
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
-    flu, = get_flinalg_funcs(('lu',),(a1,))
-    p,l,u,info = flu(a1,permute_l=permute_l,overwrite_a = overwrite_a)
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal lu.getrf' % (-info))
-    if permute_l:
-        return l,u
-    return p,l,u
-
-def svd(a, full_matrices=True, compute_uv=True, overwrite_a=False):
-    """Singular Value Decomposition.
-
-    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
-
-    For compute_uv = False, only s is returned.
-
-    Raises LinAlgError if SVD computation does not converge
-
-    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
-
-    See also
-    --------
-    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:
-        import numpy.linalg
-        return numpy.linalg.svd(a, full_matrices=0, compute_uv=compute_uv)
-    a1 = asarray_chkfinite(a)
-    if len(a1.shape) != 2:
-        raise ValueError, 'expected matrix'
-    m,n = a1.shape
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
-    gesdd, = get_lapack_funcs(('gesdd',),(a1,))
-    if gesdd.module_name[:7] == 'flapack':
-        lwork = calc_lwork.gesdd(gesdd.prefix,m,n,compute_uv)[1]
-        u,s,v,info = gesdd(a1,compute_uv = compute_uv, lwork = lwork,
-                      overwrite_a = overwrite_a)
-    else: # 'clapack'
-        raise NotImplementedError,'calling gesdd from %s' % (gesdd.module_name)
-    if info>0: raise LinAlgError, "SVD did not converge"
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal gesdd' % (-info))
-    if compute_uv:
-        return u,s,v
-    else:
-        return s
-
-def svdvals(a, overwrite_a=False):
-    """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
-
-    See also
-    --------
-    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):
-    """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)
-    if MorN == M:
-        return r_['-1',part,zeros((M,N-M),typ)]
-    elif MorN == N:
-        return r_[part,zeros((M-N,N),typ)]
-    else:
-        raise ValueError, "Length of s must be M or N."
-
-def cholesky(a, lower=False, overwrite_a=False):
-    """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]:
-        raise ValueError, 'expected square matrix'
-    overwrite_a = overwrite_a or _datanotshared(a1,a)
-    potrf, = get_lapack_funcs(('potrf',),(a1,))
-    c,info = potrf(a1,lower=lower,overwrite_a=overwrite_a,clean=1)
-    if info>0: raise LinAlgError(
-       "matrix not positive definite (leading minor of order %d"\
-       "is not positive definite)" % (info-1))
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal potrf' % (-info))
-    return c
-
-def cho_factor(a, lower=False, overwrite_a=False):
-    """Compute the Cholesky decomposition of a matrix, to use in cho_solve
-
-    Returns a matrix containing the Cholesky decomposition,
-    ``A = L L*`` or ``A = U* U`` of a Hermitian positive-definite matrix `a`.
-    The return value can be directly used as the first parameter to cho_solve.
-
-    .. warning::
-        The returned matrix also contains random data in the entries not
-        used by the Cholesky decomposition. If you need to zero these
-        entries, use the function `cholesky` instead.
-
-    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)
-        Matrix whose upper or lower triangle contains the Cholesky factor
-        of `a`. Other parts of the matrix contain random data.
-    lower : boolean
-        Flag indicating whether the factor is in the lower or upper triangle
-
-    Raises
-    ------
-    LinAlgError
-        Raised if decomposition fails.
-
-    """
-    a1 = asarray_chkfinite(a)
-    if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
-        raise ValueError, 'expected square matrix'
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
-    potrf, = get_lapack_funcs(('potrf',),(a1,))
-    c,info = potrf(a1,lower=lower,overwrite_a=overwrite_a,clean=0)
-    if info>0: raise LinAlgError(
-       "matrix not positive definite (leading minor of order %d"\
-       "is not positive definite)"%(info-1))
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal potrf' % (-info))
-    return c, lower
-
-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)
-    _assert_squareness(c)
-    b = asarray_chkfinite(b)
-    potrs, = get_lapack_funcs(('potrs',),(c,))
-    b, info = potrs(c,b,lower)
-    if info < 0:
-        msg = "Argument %d to lapack's ?potrs() has an illegal value." % info
-        raise TypeError, msg
-    if info > 0:
-        msg = "Unknown error occured int ?potrs(): error code = %d" % info
-        raise TypeError, msg
-    return b
-
-def qr(a, overwrite_a=False, lwork=None, econ=None, mode='qr'):
-    """Compute QR decomposition of a matrix.
-
-    Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
-    and R upper triangular.
-
-    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).
-        Default is False.
-    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
-
-    (for any mode)
-    R : double or complex array, shape (M, N) or (K, N) for econ==True
-        Size K = min(M, N)
-
-    Raises LinAlgError if decomposition fails
-
-    Notes
-    -----
-    This is an interface to the LAPACK routines dgeqrf, zgeqrf,
-    dorgqr, and zungqr.
-
-    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))
-
-    """
-    if econ is None:
-        econ = False
-    else:
-        warn("qr econ argument will be removed after scipy 0.7. "
-             "The economy transform will then be available through "
-             "the mode='economic' argument.", DeprecationWarning)
-
-    a1 = asarray_chkfinite(a)
-    if len(a1.shape) != 2:
-        raise ValueError("expected 2D array")
-    M, N = a1.shape
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
-
-    geqrf, = get_lapack_funcs(('geqrf',),(a1,))
-    if lwork is None or lwork == -1:
-        # get optimal work array
-        qr,tau,work,info = geqrf(a1,lwork=-1,overwrite_a=1)
-        lwork = work[0]
-
-    qr,tau,work,info = geqrf(a1,lwork=lwork,overwrite_a=overwrite_a)
-    if info<0:
-        raise ValueError("illegal value in %d-th argument of internal geqrf"
-            % -info)
-
-    if not econ or M<N:
-        R = special_matrices.triu(qr)
-    else:
-        R = special_matrices.triu(qr[0:N,0:N])
-
-    if mode=='r':
-        return R
-
-    if find_best_lapack_type((a1,))[0]=='s' or find_best_lapack_type((a1,))[0]=='d':
-        gor_un_gqr, = get_lapack_funcs(('orgqr',),(qr,))
-    else:
-        gor_un_gqr, = get_lapack_funcs(('ungqr',),(qr,))
-
-
-    if M<N:
-        # get optimal work array
-        Q,work,info = gor_un_gqr(qr[:,0:M],tau,lwork=-1,overwrite_a=1)
-        lwork = work[0]
-        Q,work,info = gor_un_gqr(qr[:,0:M],tau,lwork=lwork,overwrite_a=1)
-    elif econ:
-        # get optimal work array
-        Q,work,info = gor_un_gqr(qr,tau,lwork=-1,overwrite_a=1)
-        lwork = work[0]
-        Q,work,info = gor_un_gqr(qr,tau,lwork=lwork,overwrite_a=1)
-    else:
-        t = qr.dtype.char
-        qqr = numpy.empty((M,M),dtype=t)
-        qqr[:,0:N]=qr
-        # get optimal work array
-        Q,work,info = gor_un_gqr(qqr,tau,lwork=-1,overwrite_a=1)
-        lwork = work[0]
-        Q,work,info = gor_un_gqr(qqr,tau,lwork=lwork,overwrite_a=1)
-
-    if info < 0:
-        raise ValueError("illegal value in %d-th argument of internal gorgqr"
-            % -info)
-
-    return Q, R
-
-
-
-def qr_old(a, overwrite_a=False, lwork=None):
-    """Compute QR decomposition of a matrix.
-
-    Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
-    and R upper triangular.
-
-    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)
-
-    Raises LinAlgError if decomposition fails
-
-    """
-    a1 = asarray_chkfinite(a)
-    if len(a1.shape) != 2:
-        raise ValueError, 'expected matrix'
-    M,N = a1.shape
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
-    geqrf, = get_lapack_funcs(('geqrf',),(a1,))
-    if lwork is None or lwork == -1:
-        # get optimal work array
-        qr,tau,work,info = geqrf(a1,lwork=-1,overwrite_a=1)
-        lwork = work[0]
-    qr,tau,work,info = geqrf(a1,lwork=lwork,overwrite_a=overwrite_a)
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal geqrf' % (-info))
-    gemm, = get_blas_funcs(('gemm',),(qr,))
-    t = qr.dtype.char
-    R = special_matrices.triu(qr)
-    Q = numpy.identity(M,dtype=t)
-    ident = numpy.identity(M,dtype=t)
-    zeros = numpy.zeros
-    for i in range(min(M,N)):
-        v = zeros((M,),t)
-        v[i] = 1
-        v[i+1:M] = qr[i+1:M,i]
-        H = gemm(-tau[i],v,v,1+0j,ident,trans_b=2)
-        Q = gemm(1,Q,H)
-    return Q, R
-
-
-
-def rq(a, overwrite_a=False, lwork=None):
-    """Compute RQ decomposition of a square real matrix.
-
-    Calculate the decomposition :lm:`A = R Q` where Q is unitary/orthogonal
-    and R upper triangular.
-
-    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)
-    if len(a1.shape) != 2:
-        raise ValueError, 'expected matrix'
-    M,N = a1.shape
-    if M != N:
-        raise ValueError, 'expected square matrix'
-    if issubclass(a1.dtype.type,complexfloating):
-        raise ValueError, 'expected real (non-complex) matrix'
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
-    gerqf, = get_lapack_funcs(('gerqf',),(a1,))
-    if lwork is None or lwork == -1:
-        # get optimal work array
-        rq,tau,work,info = gerqf(a1,lwork=-1,overwrite_a=1)
-        lwork = work[0]
-    rq,tau,work,info = gerqf(a1,lwork=lwork,overwrite_a=overwrite_a)
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal geqrf' % (-info))
-    gemm, = get_blas_funcs(('gemm',),(rq,))
-    t = rq.dtype.char
-    R = special_matrices.triu(rq)
-    Q = numpy.identity(M,dtype=t)
-    ident = numpy.identity(M,dtype=t)
-    zeros = numpy.zeros
-
-    k = min(M,N)
-    for i in range(k):
-        v = zeros((M,),t)
-        v[N-k+i] = 1
-        v[0:N-k+i] = rq[M-k+i,0:N-k+i]
-        H = gemm(-tau[i],v,v,1+0j,ident,trans_b=2)
-        Q = gemm(1,Q,H)
-    return R, Q
-
 _double_precision = ['i','l','d']
 
-def schur(a, output='real', lwork=None, overwrite_a=False):
-    """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)
-
-    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.
-
-    See also
-    --------
-    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'"
-    a1 = asarray_chkfinite(a)
-    if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
-        raise ValueError, 'expected square matrix'
-    typ = a1.dtype.char
-    if output in ['complex','c'] and typ not in ['F','D']:
-        if typ in _double_precision:
-            a1 = a1.astype('D')
-            typ = 'D'
-        else:
-            a1 = a1.astype('F')
-            typ = 'F'
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
-    gees, = get_lapack_funcs(('gees',),(a1,))
-    if lwork is None or lwork == -1:
-        # get optimal work array
-        result = gees(lambda x: None,a,lwork=-1)
-        lwork = result[-2][0]
-    result = gees(lambda x: None,a,lwork=result[-2][0],overwrite_a=overwrite_a)
-    info = result[-1]
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal gees' % (-info))
-    elif info>0: raise LinAlgError, "Schur form not found.  Possibly ill-conditioned."
-    return result[0], result[-3]
-
-eps = numpy.finfo(float).eps
-feps = numpy.finfo(single).eps
-
-_array_kind = {'b':0, 'h':0, 'B': 0, 'i':0, 'l': 0, 'f': 0, 'd': 0, 'F': 1, 'D': 1}
-_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
-_array_type = [['f', 'd'], ['F', 'D']]
-
-def _commonType(*arrays):
-    kind = 0
-    precision = 0
-    for a in arrays:
-        t = a.dtype.char
-        kind = max(kind, _array_kind[t])
-        precision = max(precision, _array_precision[t])
-    return _array_type[kind][precision]
-
-def _castCopy(type, *arrays):
-    cast_arrays = ()
-    for a in arrays:
-        if a.dtype.char == type:
-            cast_arrays = cast_arrays + (a.copy(),)
-        else:
-            cast_arrays = cast_arrays + (a.astype(type),)
-    if len(cast_arrays) == 1:
-        return cast_arrays[0]
-    else:
-        return cast_arrays
-
-def _assert_squareness(*arrays):
-    for a in arrays:
-        if max(a.shape) != min(a.shape):
-            raise LinAlgError, 'Array must be square'
-
-def rsf2csf(T, Z):
-    """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
-
-    See also
-    --------
-    schur : Schur decompose a matrix
-
-    """
-    Z,T = map(asarray_chkfinite, (Z,T))
-    if len(Z.shape) !=2 or Z.shape[0] != Z.shape[1]:
-        raise ValueError, "matrix must be square."
-    if len(T.shape) !=2 or T.shape[0] != T.shape[1]:
-        raise ValueError, "matrix must be square."
-    if T.shape[0] != Z.shape[0]:
-        raise ValueError, "matrices must be same dimension."
-    N = T.shape[0]
-    arr = numpy.array
-    t = _commonType(Z, T, arr([3.0],'F'))
-    Z, T = _castCopy(t, Z, T)
-    conj = numpy.conj
-    dot = numpy.dot
-    r_ = numpy.r_
-    transp = numpy.transpose
-    for m in range(N-1,0,-1):
-        if abs(T[m,m-1]) > eps*(abs(T[m-1,m-1]) + abs(T[m,m])):
-            k = slice(m-1,m+1)
-            mu = eigvals(T[k,k]) - T[m,m]
-            r = misc.norm([mu[0], T[m,m-1]])
-            c = mu[0] / r
-            s = T[m,m-1] / r
-            G = r_[arr([[conj(c),s]],dtype=t),arr([[-s,c]],dtype=t)]
-            Gc = conj(transp(G))
-            j = slice(m-1,N)
-            T[k,j] = dot(G,T[k,j])
-            i = slice(0,m+1)
-            T[i,k] = dot(T[i,k], Gc)
-            i = slice(0,N)
-            Z[i,k] = dot(Z[i,k], Gc)
-        T[m,m-1] = 0.0;
-    return T, Z
-
-
-# Orthonormal decomposition
-
-def orth(A):
-    """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
-
-    See also
-    --------
-    svd : Singular value decomposition of a matrix
-
-    """
-    u,s,vh = svd(A)
-    M,N = A.shape
-    tol = max(M,N)*numpy.amax(s)*eps
-    num = numpy.sum(s > tol,dtype=int)
-    Q = u[:,:num]
-    return Q
-
 def hessenberg(a, calc_q=False, overwrite_a=False):
     """Compute Hessenberg form of a matrix.
 
@@ -1572,39 +736,41 @@
     """
     a1 = asarray(a)
     if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
-        raise ValueError, 'expected square matrix'
-    overwrite_a = overwrite_a or (_datanotshared(a1,a))
-    gehrd,gebal = get_lapack_funcs(('gehrd','gebal'),(a1,))
-    ba,lo,hi,pivscale,info = gebal(a,permute=1,overwrite_a = overwrite_a)
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal gebal (hessenberg)' % (-info))
+        raise ValueError('expected square matrix')
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
+    gehrd,gebal = get_lapack_funcs(('gehrd','gebal'), (a1,))
+    ba, lo, hi, pivscale, info = gebal(a, permute=1, overwrite_a=overwrite_a)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal gebal '
+                                                    '(hessenberg)' % -info)
     n = len(a1)
-    lwork = calc_lwork.gehrd(gehrd.prefix,n,lo,hi)
-    hq,tau,info = gehrd(ba,lo=lo,hi=hi,lwork=lwork,overwrite_a=1)
-    if info<0: raise ValueError(
-       'illegal value in %d-th argument of internal gehrd (hessenberg)' % (-info))
+    lwork = calc_lwork.gehrd(gehrd.prefix, n, lo, hi)
+    hq, tau, info = gehrd(ba, lo=lo, hi=hi, lwork=lwork, overwrite_a=1)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal gehrd '
+                                        '(hessenberg)' % -info)
 
     if not calc_q:
-        for i in range(lo,hi):
-            hq[i+2:hi+1,i] = 0.0
+        for i in range(lo, hi):
+            hq[i+2:hi+1, i] = 0.0
         return hq
 
     # XXX: Use ORGHR routines to compute q.
-    ger,gemm = get_blas_funcs(('ger','gemm'),(hq,))
+    ger,gemm = get_blas_funcs(('ger','gemm'), (hq,))
     typecode = hq.dtype.char
     q = None
-    for i in range(lo,hi):
+    for i in range(lo, hi):
         if tau[i]==0.0:
             continue
-        v = zeros(n,dtype=typecode)
+        v = zeros(n, dtype=typecode)
         v[i+1] = 1.0
-        v[i+2:hi+1] = hq[i+2:hi+1,i]
-        hq[i+2:hi+1,i] = 0.0
-        h = ger(-tau[i],v,v,a=diag(ones(n,dtype=typecode)),overwrite_a=1)
+        v[i+2:hi+1] = hq[i+2:hi+1, i]
+        hq[i+2:hi+1, i] = 0.0
+        h = ger(-tau[i], v, v,a=diag(ones(n, dtype=typecode)), overwrite_a=1)
         if q is None:
             q = h
         else:
-            q = gemm(1.0,q,h)
+            q = gemm(1.0, q, h)
     if q is None:
-        q = diag(ones(n,dtype=typecode))
-    return hq,q
+        q = diag(ones(n, dtype=typecode))
+    return hq, q

Added: trunk/scipy/linalg/decomp_cholesky.py
===================================================================
--- trunk/scipy/linalg/decomp_cholesky.py	                        (rev 0)
+++ trunk/scipy/linalg/decomp_cholesky.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -0,0 +1,217 @@
+"""Cholesky decomposition functions."""
+
+from numpy import asarray_chkfinite
+
+# Local imports
+from misc import LinAlgError, _datanotshared
+from lapack import get_lapack_funcs
+
+__all__ = ['cholesky', 'cho_factor', 'cho_solve', 'cholesky_banded',
+            'cho_solve_banded']
+
+
+def _cholesky(a, lower=False, overwrite_a=False, clean=True):
+    """Common code for cholesky() and cho_factor()."""
+
+    a1 = asarray_chkfinite(a)
+    if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
+        raise ValueError('expected square matrix')
+
+    overwrite_a = overwrite_a or _datanotshared(a1, a)
+    potrf, = get_lapack_funcs(('potrf',), (a1,))
+    c, info = potrf(a1, lower=lower, overwrite_a=overwrite_a, clean=clean)
+    if info > 0:
+        raise LinAlgError("%d-th leading minor not positive definite" % info)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal potrf'
+                                                                    % -info)
+    return c, lower
+
+def cholesky(a, lower=False, overwrite_a=False):
+    """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
+    -------
+    c : 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]])
+
+    """
+    c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True)
+    return c
+
+
+def cho_factor(a, lower=False, overwrite_a=False):
+    """Compute the Cholesky decomposition of a matrix, to use in cho_solve
+
+    Returns a matrix containing the Cholesky decomposition,
+    ``A = L L*`` or ``A = U* U`` of a Hermitian positive-definite matrix `a`.
+    The return value can be directly used as the first parameter to cho_solve.
+
+    .. warning::
+        The returned matrix also contains random data in the entries not
+        used by the Cholesky decomposition. If you need to zero these
+        entries, use the function `cholesky` instead.
+
+    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)
+        Matrix whose upper or lower triangle contains the Cholesky factor
+        of `a`. Other parts of the matrix contain random data.
+    lower : boolean
+        Flag indicating whether the factor is in the lower or upper triangle
+
+    Raises
+    ------
+    LinAlgError
+        Raised if decomposition fails.
+
+    See also
+    --------
+    cho_solve : Solve a linear set equations using the Cholesky factorization 
+                of a matrix.
+
+    """
+    c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=False)
+    return c, lower
+
+
+def cho_solve((c, lower), b, overwrite_b=False):
+    """Solve the linear equations A x = b, given the Cholesky factorization of A.
+
+    Parameters
+    ----------
+    (c, lower) : tuple, (array, bool)
+        Cholesky factorization of a, as given by cho_factor
+    b : array
+        Right-hand side
+
+    Returns
+    -------
+    x : array
+        The solution to the system A x = b
+
+    See also
+    --------
+    cho_factor : Cholesky factorization of a matrix
+
+    """
+
+    b1 = asarray_chkfinite(b)
+    c = asarray_chkfinite(c)
+    if c.ndim != 2 or c.shape[0] != c.shape[1]:
+        raise ValueError("The factored matrix c is not square.")
+    if c.shape[1] != b1.shape[0]:
+        raise ValueError("incompatible dimensions.")
+
+    overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
+
+    potrs, = get_lapack_funcs(('potrs',), (c, b1))
+    x, info = potrs(c, b1, lower=lower, overwrite_b=overwrite_b)
+    if info != 0:
+        raise ValueError('illegal value in %d-th argument of internal potrs'
+                                                                    % -info)
+    return x
+
+def cholesky_banded(ab, overwrite_ab=False, lower=False):
+    """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 (u + 1, M)
+        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 (u+1, M)
+        Cholesky factorization of a, in the same banded format as ab
+
+    """
+    ab = asarray_chkfinite(ab)
+
+    pbtrf, = get_lapack_funcs(('pbtrf',), (ab,))
+    c, info = pbtrf(ab, lower=lower, overwrite_ab=overwrite_ab)
+    if info > 0:
+        raise LinAlgError("%d-th leading minor not positive definite" % info)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal pbtrf'
+                                                                    % -info)
+    return c
+
+# my new function
+def cho_solve_banded((ab, lower), b, overwrite_b=False):
+    """To be written..."""
+
+    ab = asarray_chkfinite(ab)
+    b = asarray_chkfinite(b)
+
+    # Validate shapes.
+    if ab.shape[-1] != b.shape[0]:
+        raise ValueError("shapes of ab and b are not compatible.")
+
+    pbtrs, = get_lapack_funcs(('pbtrs',), (ab, b))
+    x, info = pbtrs(ab, b, lower=lower, overwrite_b=overwrite_b)
+    if info > 0:
+        raise LinAlgError("%d-th leading minor not positive definite" % info)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal pbtrs'
+                                                                    % -info)
+    return x

Added: trunk/scipy/linalg/decomp_lu.py
===================================================================
--- trunk/scipy/linalg/decomp_lu.py	                        (rev 0)
+++ trunk/scipy/linalg/decomp_lu.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -0,0 +1,159 @@
+"""LU decomposition functions."""
+
+from warnings import warn
+
+from numpy import asarray, asarray_chkfinite
+
+# Local imports
+from misc import _datanotshared
+from lapack import get_lapack_funcs
+from flinalg import get_flinalg_funcs
+
+
+def lu_factor(a, overwrite_a=False):
+    """Compute pivoted LU decomposition of a matrix.
+
+    The decomposition is::
+
+        A = P L U
+
+    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)
+
+    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].
+
+    See also
+    --------
+    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]):
+        raise ValueError, 'expected square matrix'
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
+    getrf, = get_lapack_funcs(('getrf',), (a1,))
+    lu, piv, info = getrf(a, overwrite_a=overwrite_a)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of '
+                                'internal getrf (lu_factor)' % -info)
+    if info > 0:
+        warn("Diagonal number %d is exactly zero. Singular matrix." % info,
+                    RuntimeWarning)
+    return lu, piv
+
+
+def lu_solve((lu, piv), b, trans=0, overwrite_b=False):
+    """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:
+
+        =====  =========
+        trans  system
+        =====  =========
+        0      a x   = b
+        1      a^T x = b
+        2      a^H x = b
+        =====  =========
+
+    Returns
+    -------
+    x : array
+        Solution to the system
+
+    See also
+    --------
+    lu_factor : LU factorize a matrix
+
+    """
+    b1 = asarray_chkfinite(b)
+    overwrite_b = overwrite_b or (b1 is not b and not hasattr(b, '__array__'))
+    if lu.shape[0] != b1.shape[0]:
+        raise ValueError("incompatible dimensions.")
+
+    getrs, = get_lapack_funcs(('getrs',), (lu, b1))
+    x,info = getrs(lu, piv, b1, trans=trans, overwrite_b=overwrite_b)
+    if info == 0:
+        return x
+    raise ValueError('illegal value in %d-th argument of internal gesv|posv'
+                                                                    % -info)
+
+
+def lu(a, permute_l=False, overwrite_a=False):
+    """Compute pivoted LU decompostion of a matrix.
+
+    The decomposition is::
+
+        A = P L U
+
+    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)
+
+    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:
+        raise ValueError('expected matrix')
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
+    flu, = get_flinalg_funcs(('lu',), (a1,))
+    p, l, u, info = flu(a1, permute_l=permute_l, overwrite_a=overwrite_a)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of '
+                                            'internal lu.getrf' % -info)
+    if permute_l:
+        return l, u
+    return p, l, u

Added: trunk/scipy/linalg/decomp_qr.py
===================================================================
--- trunk/scipy/linalg/decomp_qr.py	                        (rev 0)
+++ trunk/scipy/linalg/decomp_qr.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -0,0 +1,248 @@
+"""QR decomposition functions."""
+
+from warnings import warn
+
+import numpy
+from numpy import asarray_chkfinite, complexfloating
+
+# Local imports
+import special_matrices
+from blas import get_blas_funcs
+from lapack import get_lapack_funcs, find_best_lapack_type
+from misc import _datanotshared
+
+
+def qr(a, overwrite_a=False, lwork=None, econ=None, mode='qr'):
+    """Compute QR decomposition of a matrix.
+
+    Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
+    and R upper triangular.
+
+    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).
+        Default is False.
+    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
+
+    (for any mode)
+    R : double or complex array, shape (M, N) or (K, N) for econ==True
+        Size K = min(M, N)
+
+    Raises LinAlgError if decomposition fails
+
+    Notes
+    -----
+    This is an interface to the LAPACK routines dgeqrf, zgeqrf,
+    dorgqr, and zungqr.
+
+    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))
+
+    """
+    if econ is None:
+        econ = False
+    else:
+        warn("qr econ argument will be removed after scipy 0.7. "
+             "The economy transform will then be available through "
+             "the mode='economic' argument.", DeprecationWarning)
+
+    a1 = asarray_chkfinite(a)
+    if len(a1.shape) != 2:
+        raise ValueError("expected 2D array")
+    M, N = a1.shape
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
+
+    geqrf, = get_lapack_funcs(('geqrf',), (a1,))
+    if lwork is None or lwork == -1:
+        # get optimal work array
+        qr, tau, work, info = geqrf(a1, lwork=-1, overwrite_a=1)
+        lwork = work[0]
+
+    qr, tau, work, info = geqrf(a1, lwork=lwork, overwrite_a=overwrite_a)
+    if info < 0:
+        raise ValueError("illegal value in %d-th argument of internal geqrf"
+                                                                    % -info)
+    if not econ or M < N:
+        R = special_matrices.triu(qr)
+    else:
+        R = special_matrices.triu(qr[0:N, 0:N])
+
+    if mode == 'r':
+        return R
+
+    if find_best_lapack_type((a1,))[0] == 's' or \
+                find_best_lapack_type((a1,))[0] == 'd':
+        gor_un_gqr, = get_lapack_funcs(('orgqr',), (qr,))
+    else:
+        gor_un_gqr, = get_lapack_funcs(('ungqr',), (qr,))
+
+    if M < N:
+        # get optimal work array
+        Q, work, info = gor_un_gqr(qr[:,0:M], tau, lwork=-1, overwrite_a=1)
+        lwork = work[0]
+        Q, work, info = gor_un_gqr(qr[:,0:M], tau, lwork=lwork, overwrite_a=1)
+    elif econ:
+        # get optimal work array
+        Q, work, info = gor_un_gqr(qr, tau, lwork=-1, overwrite_a=1)
+        lwork = work[0]
+        Q, work, info = gor_un_gqr(qr, tau, lwork=lwork, overwrite_a=1)
+    else:
+        t = qr.dtype.char
+        qqr = numpy.empty((M, M), dtype=t)
+        qqr[:,0:N] = qr
+        # get optimal work array
+        Q, work, info = gor_un_gqr(qqr, tau, lwork=-1, overwrite_a=1)
+        lwork = work[0]
+        Q, work, info = gor_un_gqr(qqr, tau, lwork=lwork, overwrite_a=1)
+
+    if info < 0:
+        raise ValueError("illegal value in %d-th argument of internal gorgqr"
+                                                                    % -info)
+    return Q, R
+
+
+
+def qr_old(a, overwrite_a=False, lwork=None):
+    """Compute QR decomposition of a matrix.
+
+    Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
+    and R upper triangular.
+
+    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)
+
+    Raises LinAlgError if decomposition fails
+
+    """
+    a1 = asarray_chkfinite(a)
+    if len(a1.shape) != 2:
+        raise ValueError, 'expected matrix'
+    M,N = a1.shape
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
+    geqrf, = get_lapack_funcs(('geqrf',), (a1,))
+    if lwork is None or lwork == -1:
+        # get optimal work array
+        qr, tau, work, info = geqrf(a1, lwork=-1, overwrite_a=1)
+        lwork = work[0]
+    qr, tau, work, info = geqrf(a1, lwork=lwork, overwrite_a=overwrite_a)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal geqrf'
+                                                                    % -info)
+    gemm, = get_blas_funcs(('gemm',), (qr,))
+    t = qr.dtype.char
+    R = special_matrices.triu(qr)
+    Q = numpy.identity(M, dtype=t)
+    ident = numpy.identity(M, dtype=t)
+    zeros = numpy.zeros
+    for i in range(min(M, N)):
+        v = zeros((M,), t)
+        v[i] = 1
+        v[i+1:M] = qr[i+1:M, i]
+        H = gemm(-tau[i], v, v, 1+0j, ident, trans_b=2)
+        Q = gemm(1, Q, H)
+    return Q, R
+
+
+def rq(a, overwrite_a=False, lwork=None):
+    """Compute RQ decomposition of a square real matrix.
+
+    Calculate the decomposition :lm:`A = R Q` where Q is unitary/orthogonal
+    and R upper triangular.
+
+    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)
+    if len(a1.shape) != 2:
+        raise ValueError('expected matrix')
+    M,N = a1.shape
+    if M != N:
+        raise ValueError('expected square matrix')
+    if issubclass(a1.dtype.type, complexfloating):
+        raise ValueError('expected real (non-complex) matrix')
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
+    gerqf, = get_lapack_funcs(('gerqf',), (a1,))
+    if lwork is None or lwork == -1:
+        # get optimal work array
+        rq, tau, work, info = gerqf(a1, lwork=-1, overwrite_a=1)
+        lwork = work[0]
+    rq, tau, work, info = gerqf(a1, lwork=lwork, overwrite_a=overwrite_a)
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal geqrf'
+                                                                    % -info)
+    gemm, = get_blas_funcs(('gemm',), (rq,))
+    t = rq.dtype.char
+    R = special_matrices.triu(rq)
+    Q = numpy.identity(M, dtype=t)
+    ident = numpy.identity(M, dtype=t)
+    zeros = numpy.zeros
+
+    k = min(M, N)
+    for i in range(k):
+        v = zeros((M,), t)
+        v[N-k+i] = 1
+        v[0:N-k+i] = rq[M-k+i, 0:N-k+i]
+        H = gemm(-tau[i], v, v, 1+0j, ident, trans_b=2)
+        Q = gemm(1, Q, H)
+    return R, Q

Added: trunk/scipy/linalg/decomp_schur.py
===================================================================
--- trunk/scipy/linalg/decomp_schur.py	                        (rev 0)
+++ trunk/scipy/linalg/decomp_schur.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -0,0 +1,167 @@
+"""Schur decomposition functions."""
+
+import numpy
+from numpy import asarray_chkfinite, single
+
+# Local imports.
+import misc
+from misc import LinAlgError, _datanotshared
+from lapack import get_lapack_funcs
+from decomp import eigvals
+
+
+__all__ = ['schur', 'rsf2csf']
+
+_double_precision = ['i','l','d']
+
+def schur(a, output='real', lwork=None, overwrite_a=False):
+    """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)
+
+    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.
+
+    See also
+    --------
+    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'"
+    a1 = asarray_chkfinite(a)
+    if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
+        raise ValueError, 'expected square matrix'
+    typ = a1.dtype.char
+    if output in ['complex','c'] and typ not in ['F','D']:
+        if typ in _double_precision:
+            a1 = a1.astype('D')
+            typ = 'D'
+        else:
+            a1 = a1.astype('F')
+            typ = 'F'
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
+    gees, = get_lapack_funcs(('gees',), (a1,))
+    if lwork is None or lwork == -1:
+        # get optimal work array
+        result = gees(lambda x: None, a, lwork=-1)
+        lwork = result[-2][0]
+    result = gees(lambda x: None, a, lwork=result[-2][0], overwrite_a=overwrite_a)
+    info = result[-1]
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal gees'
+                                                                    % -info)
+    elif info > 0:
+        raise LinAlgError("Schur form not found.  Possibly ill-conditioned.")
+    return result[0], result[-3]
+
+
+eps = numpy.finfo(float).eps
+feps = numpy.finfo(single).eps
+
+_array_kind = {'b':0, 'h':0, 'B': 0, 'i':0, 'l': 0, 'f': 0, 'd': 0, 'F': 1, 'D': 1}
+_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
+_array_type = [['f', 'd'], ['F', 'D']]
+
+def _commonType(*arrays):
+    kind = 0
+    precision = 0
+    for a in arrays:
+        t = a.dtype.char
+        kind = max(kind, _array_kind[t])
+        precision = max(precision, _array_precision[t])
+    return _array_type[kind][precision]
+
+def _castCopy(type, *arrays):
+    cast_arrays = ()
+    for a in arrays:
+        if a.dtype.char == type:
+            cast_arrays = cast_arrays + (a.copy(),)
+        else:
+            cast_arrays = cast_arrays + (a.astype(type),)
+    if len(cast_arrays) == 1:
+        return cast_arrays[0]
+    else:
+        return cast_arrays
+
+
+def rsf2csf(T, Z):
+    """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
+
+    See also
+    --------
+    schur : Schur decompose a matrix
+
+    """
+    Z, T = map(asarray_chkfinite, (Z, T))
+    if len(Z.shape) != 2 or Z.shape[0] != Z.shape[1]:
+        raise ValueError("matrix must be square.")
+    if len(T.shape) != 2 or T.shape[0] != T.shape[1]:
+        raise ValueError("matrix must be square.")
+    if T.shape[0] != Z.shape[0]:
+        raise ValueError("matrices must be same dimension.")
+    N = T.shape[0]
+    arr = numpy.array
+    t = _commonType(Z, T, arr([3.0],'F'))
+    Z, T = _castCopy(t, Z, T)
+    conj = numpy.conj
+    dot = numpy.dot
+    r_ = numpy.r_
+    transp = numpy.transpose
+    for m in range(N-1, 0, -1):
+        if abs(T[m,m-1]) > eps*(abs(T[m-1,m-1]) + abs(T[m,m])):
+            k = slice(m-1, m+1)
+            mu = eigvals(T[k,k]) - T[m,m]
+            r = misc.norm([mu[0], T[m,m-1]])
+            c = mu[0] / r
+            s = T[m,m-1] / r
+            G = r_[arr([[conj(c), s]], dtype=t), arr([[-s, c]], dtype=t)]
+            Gc = conj(transp(G))
+            j = slice(m-1, N)
+            T[k,j] = dot(G, T[k,j])
+            i = slice(0, m+1)
+            T[i,k] = dot(T[i,k], Gc)
+            i = slice(0, N)
+            Z[i,k] = dot(Z[i,k], Gc)
+        T[m,m-1] = 0.0;
+    return T, Z

Added: trunk/scipy/linalg/decomp_svd.py
===================================================================
--- trunk/scipy/linalg/decomp_svd.py	                        (rev 0)
+++ trunk/scipy/linalg/decomp_svd.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -0,0 +1,173 @@
+"""SVD decomposition functions."""
+
+import numpy
+from numpy import asarray_chkfinite, zeros, r_, diag
+from scipy.linalg import calc_lwork
+
+# Local imports.
+from misc import LinAlgError, _datanotshared
+from lapack import get_lapack_funcs
+
+
+def svd(a, full_matrices=True, compute_uv=True, overwrite_a=False):
+    """Singular Value Decomposition.
+
+    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
+
+    For compute_uv = False, only s is returned.
+
+    Raises LinAlgError if SVD computation does not converge
+
+    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
+
+    See also
+    --------
+    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:
+        import numpy.linalg
+        return numpy.linalg.svd(a, full_matrices=0, compute_uv=compute_uv)
+    a1 = asarray_chkfinite(a)
+    if len(a1.shape) != 2:
+        raise ValueError('expected matrix')
+    m,n = a1.shape
+    overwrite_a = overwrite_a or (_datanotshared(a1, a))
+    gesdd, = get_lapack_funcs(('gesdd',), (a1,))
+    if gesdd.module_name[:7] == 'flapack':
+        lwork = calc_lwork.gesdd(gesdd.prefix, m, n, compute_uv)[1]
+        u,s,v,info = gesdd(a1,compute_uv = compute_uv, lwork = lwork,
+                                                overwrite_a = overwrite_a)
+    else: # 'clapack'
+        raise NotImplementedError('calling gesdd from %s' % gesdd.module_name)
+    if info > 0:
+        raise LinAlgError("SVD did not converge")
+    if info < 0:
+        raise ValueError('illegal value in %d-th argument of internal gesdd'
+                                                                    % -info)
+    if compute_uv:
+        return u, s, v
+    else:
+        return s
+
+def svdvals(a, overwrite_a=False):
+    """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
+
+    See also
+    --------
+    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):
+    """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)
+    if MorN == M:
+        return r_['-1', part, zeros((M, N-M), typ)]
+    elif MorN == N:
+        return r_[part, zeros((M-N,N), typ)]
+    else:
+        raise ValueError("Length of s must be M or N.")
+
+
+# Orthonormal decomposition
+
+def orth(A):
+    """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
+
+    See also
+    --------
+    svd : Singular value decomposition of a matrix
+
+    """
+    u, s, vh = svd(A)
+    M, N = A.shape
+    eps = numpy.finfo(float).eps
+    tol = max(M,N) * numpy.amax(s) * eps
+    num = numpy.sum(s > tol, dtype=int)
+    Q = u[:,:num]
+    return Q

Modified: trunk/scipy/linalg/generic_flapack.pyf
===================================================================
--- trunk/scipy/linalg/generic_flapack.pyf	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/scipy/linalg/generic_flapack.pyf	2010-04-10 01:47:20 UTC (rev 6319)
@@ -56,6 +56,57 @@
      
    end subroutine <tchar=c,z>pbtrf
 
+
+   subroutine <tchar=s,d>pbtrs(lower, n, kd, nrhs, ab, ldab, b, ldb, info)
+
+     ! Solve a system of linear equations A*X = B with a symmetric
+     ! positive definite band matrix A using the Cholesky factorization.
+     ! AB is the triangular factur U or L from the Cholesky factorization
+     ! previously computed with *PBTRF.
+     ! A = U^T * U, AB = U if lower = 0
+     ! A = L * L^T, AB = L if lower = 1
+
+     callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,&nrhs,ab,&ldab,b,&ldb,&info); 
+     callprotoargument char*,int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*
+
+     integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
+     integer intent(hide),depend(ab) :: n=shape(ab,1)
+     integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
+     integer intent(hide),depend(b) :: ldb=shape(b,0)
+     integer intent(hide),depend(b) :: nrhs=shape(b,1)
+     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0   
+
+     <type_in> dimension(ldb, nrhs),intent(in,out,copy,out=x) :: b
+     <type_in> dimension(ldab,n),intent(in) :: ab
+     integer intent(out) :: info
+     
+   end subroutine <tchar=s,d>pbtrs
+
+   subroutine <tchar=c,z>pbtrs(lower, n, kd, nrhs, ab, ldab, b, ldb, info)
+
+     ! Solve a system of linear equations A*X = B with a symmetric
+     ! positive definite band matrix A using the Cholesky factorization.
+     ! AB is the triangular factur U or L from the Cholesky factorization
+     ! previously computed with *PBTRF.
+     ! A = U^T * U, AB = U if lower = 0
+     ! A = L * L^T, AB = L if lower = 1
+
+     callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,&nrhs,ab,&ldab,b,&ldb,&info); 
+     callprotoargument char*,int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*
+
+     integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
+     integer intent(hide),depend(ab) :: n=shape(ab,1)
+     integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
+     integer intent(hide),depend(b) :: ldb=shape(b,0)
+     integer intent(hide),depend(b) :: nrhs=shape(b,1)
+     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
+
+     <type_in> dimension(ldb, nrhs),intent(in,out,copy,out=x) :: b
+     <type_in> dimension(ldab,n),intent(in) :: ab
+     integer intent(out) :: info
+     
+   end subroutine <tchar=c,z>pbtrs
+
    subroutine <tchar=s,d>pbsv(lower,n,kd,nrhs,ab,ldab,b,ldb,info)
    
      !

Modified: trunk/scipy/linalg/info.py
===================================================================
--- trunk/scipy/linalg/info.py	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/scipy/linalg/info.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -14,12 +14,17 @@
    pinv       --- Pseudo-inverse (Moore-Penrose) using lstsq
    pinv2      --- Pseudo-inverse using svd
 
-Eigenvalues and Decompositions::
+Eigenvalue Problem::
 
    eig        --- Find the eigenvalues and vectors of a square matrix
    eigvals    --- Find the eigenvalues of a square matrix
+   eigh       --- Find the eigenvalues and eigenvectors of a complex Hermitian or real symmetric matrix.
+   eigvalsh   --- Find the eigenvalues of a complex Hermitian or real symmetric matrix.
    eig_banded --- Find the eigenvalues and vectors of a band matrix
    eigvals_banded --- Find the eigenvalues of a band matrix
+
+Decompositions::
+
    lu         --- LU decomposition of a matrix
    lu_factor  --- LU decomposition returning unordered matrix and pivots
    lu_solve   --- solve Ax=b using back substitution with output of lu_factor

Modified: trunk/scipy/linalg/matfuncs.py
===================================================================
--- trunk/scipy/linalg/matfuncs.py	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/scipy/linalg/matfuncs.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -1,5 +1,3 @@
-## Automatically adapted for scipy Oct 18, 2005 by
-
 #
 # Author: Travis Oliphant, March 2002
 #
@@ -13,9 +11,15 @@
      isfinite, sqrt, identity, single
 from numpy import matrix as mat
 import numpy as np
-from basic import solve, inv, norm, triu, all_mat
-from decomp import eig, schur, rsf2csf, orth, svd
 
+# Local imports
+from misc import norm
+from basic import solve, inv
+from special_matrices import triu, all_mat
+from decomp import eig
+from decomp_svd import orth, svd
+from decomp_schur import schur, rsf2csf
+
 eps = np.finfo(float).eps
 feps = np.finfo(single).eps
 

Modified: trunk/scipy/linalg/misc.py
===================================================================
--- trunk/scipy/linalg/misc.py	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/scipy/linalg/misc.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -1,10 +1,21 @@
 import numpy as np
 from numpy.linalg import LinAlgError
 
-### Norm
+__all__ = ['LinAlgError', 'norm']
 
+
 def norm(a, ord=None):
     # Differs from numpy only in non-finite handling
     return np.linalg.norm(np.asarray_chkfinite(a), ord=ord)
 norm.__doc__ = np.linalg.norm.__doc__
 
+
+def _datanotshared(a1,a):
+    if a1 is a:
+        return False
+    else:
+        #try comparing data pointers
+        try:
+            return a1.__array_interface__['data'][0] != a.__array_interface__['data'][0]
+        except:
+            return True
\ No newline at end of file

Modified: trunk/scipy/linalg/tests/test_basic.py
===================================================================
--- trunk/scipy/linalg/tests/test_basic.py	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/scipy/linalg/tests/test_basic.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -107,7 +107,7 @@
         # with the RHS as a 1D array.
         ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]])
         b = array([1.0, 4.0, 1.0])
-        c, x = solveh_banded(ab, b)
+        x = solveh_banded(ab, b)
         assert_array_almost_equal(x, [0.0, 1.0, 0.0])
 
     def test_02_upper(self):
@@ -121,7 +121,7 @@
         b = array([[1.0, 4.0],
                    [4.0, 2.0],
                    [1.0, 4.0]])
-        c, x = solveh_banded(ab, b)
+        x = solveh_banded(ab, b)
         expected = array([[0.0, 1.0],
                           [1.0, 0.0],
                           [0.0, 1.0]])
@@ -135,7 +135,7 @@
         # with the RHS as a 2D array with shape (3,1).
         ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]])
         b = array([1.0, 4.0, 1.0]).reshape(-1,1)
-        c, x = solveh_banded(ab, b)
+        x = solveh_banded(ab, b)
         assert_array_almost_equal(x, array([0.0, 1.0, 0.0]).reshape(-1,1))
 
     def test_01_lower(self):
@@ -147,7 +147,7 @@
         ab = array([[4.0, 4.0, 4.0],
                     [1.0, 1.0, -99]])
         b = array([1.0, 4.0, 1.0])
-        c, x = solveh_banded(ab, b, lower=True)
+        x = solveh_banded(ab, b, lower=True)
         assert_array_almost_equal(x, [0.0, 1.0, 0.0])
 
     def test_02_lower(self):
@@ -161,7 +161,7 @@
         b = array([[1.0, 4.0],
                    [4.0, 2.0],
                    [1.0, 4.0]])
-        c, x = solveh_banded(ab, b, lower=True)
+        x = solveh_banded(ab, b, lower=True)
         expected = array([[0.0, 1.0],
                           [1.0, 0.0],
                           [0.0, 1.0]])
@@ -175,7 +175,7 @@
         #
         ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]], dtype=float32)
         b = array([1.0, 4.0, 1.0], dtype=float32)
-        c, x = solveh_banded(ab, b)
+        x = solveh_banded(ab, b)
         assert_array_almost_equal(x, [0.0, 1.0, 0.0])
 
     def test_02_float32(self):
@@ -189,7 +189,7 @@
         b = array([[1.0, 4.0],
                    [4.0, 2.0],
                    [1.0, 4.0]], dtype=float32)
-        c, x = solveh_banded(ab, b)
+        x = solveh_banded(ab, b)
         expected = array([[0.0, 1.0],
                           [1.0, 0.0],
                           [0.0, 1.0]])
@@ -203,7 +203,7 @@
         #
         ab = array([[-99, -1.0j, -1.0j], [4.0, 4.0, 4.0]])
         b = array([-1.0j, 4.0-1j, 4+1j])
-        c, x = solveh_banded(ab, b)
+        x = solveh_banded(ab, b)
         assert_array_almost_equal(x, [0.0, 1.0, 1.0])
 
     def test_02_complex(self):
@@ -217,7 +217,7 @@
         b = array([[   -1j,    4.0j],
                    [4.0-1j, -1.0-1j],
                    [4.0+1j,     4.0]])
-        c, x = solveh_banded(ab, b)
+        x = solveh_banded(ab, b)
         expected = array([[0.0, 1.0j],
                           [1.0,  0.0],
                           [1.0,  1.0]])

Modified: trunk/scipy/linalg/tests/test_decomp.py
===================================================================
--- trunk/scipy/linalg/tests/test_decomp.py	2010-04-09 17:57:46 UTC (rev 6318)
+++ trunk/scipy/linalg/tests/test_decomp.py	2010-04-10 01:47:20 UTC (rev 6319)
@@ -17,8 +17,8 @@
 import numpy as np
 from numpy.testing import *
 
-from scipy.linalg import eig,eigvals,lu,svd,svdvals,cholesky,qr, \
-     schur,rsf2csf, lu_solve,lu_factor,solve,diagsvd,hessenberg,rq, \
+from scipy.linalg import eig, eigvals, lu, svd, svdvals, cholesky, qr, \
+     schur, rsf2csf, lu_solve, lu_factor, solve, diagsvd, hessenberg, rq, \
      eig_banded, eigvals_banded, eigh
 from scipy.linalg.flapack import dgbtrf, dgbtrs, zgbtrf, zgbtrs, \
      dsbev, dsbevd, dsbevx, zhbevd, zhbevx
@@ -837,8 +837,8 @@
             c = transpose(c)
             a = dot(c,transpose(conjugate(c)))
             assert_array_almost_equal(cholesky(a,lower=1),c)
+        
 
-
 class TestQR(TestCase):
 
     def test_simple(self):