[Scipy-svn] r3861 - in trunk/scipy: . linalg linalg/tests sandbox/arpack/tests splinalg splinalg/isolve splinalg/isolve/tests

scipy-svn@scip... scipy-svn@scip...
Fri Jan 25 14:01:12 CST 2008


Author: wnbell
Date: 2008-01-25 14:00:52 -0600 (Fri, 25 Jan 2008)
New Revision: 3861

Added:
   trunk/scipy/splinalg/
   trunk/scipy/splinalg/eigen/
   trunk/scipy/splinalg/isolve/
   trunk/scipy/splinalg/isolve/iterative.py
   trunk/scipy/splinalg/isolve/iterative/
   trunk/scipy/splinalg/isolve/tests/
   trunk/scipy/splinalg/isolve/tests/test_iterative.py
Removed:
   trunk/scipy/linalg/iterative.py
   trunk/scipy/linalg/iterative/
   trunk/scipy/linalg/tests/test_iterative.py
Modified:
   trunk/scipy/linalg/__init__.py
   trunk/scipy/linalg/setup.py
   trunk/scipy/sandbox/arpack/tests/test_speigs.py
   trunk/scipy/setup.py
Log:
moved iterative solvers to splinalg.isolve


Modified: trunk/scipy/linalg/__init__.py
===================================================================
--- trunk/scipy/linalg/__init__.py	2008-01-24 22:11:36 UTC (rev 3860)
+++ trunk/scipy/linalg/__init__.py	2008-01-25 20:00:52 UTC (rev 3861)
@@ -1,5 +1,5 @@
 #
-# linalg - Linear algebra routines
+# linalg - Dense Linear Algebra routines
 #
 
 from info import __doc__
@@ -9,8 +9,11 @@
 from decomp import *
 from matfuncs import *
 from blas import *
-from iterative import *
 
+#from iterative import *
+# TODO remove this
+from scipy.splinalg.isolve import *
+
 __all__ = filter(lambda s:not s.startswith('_'),dir())
 
 from numpy.dual import register_func

Deleted: trunk/scipy/linalg/iterative.py
===================================================================
--- trunk/scipy/linalg/iterative.py	2008-01-24 22:11:36 UTC (rev 3860)
+++ trunk/scipy/linalg/iterative.py	2008-01-25 20:00:52 UTC (rev 3861)
@@ -1,807 +0,0 @@
-## Automatically adapted for scipy Oct 18, 2005 by
-
-
-# Iterative methods using reverse-communication raw material
-#   These methods solve
-#   Ax = b  for x
-
-#   where A must have A.matvec(x,*args) defined
-#    or be a numeric array
-
-
-__all__ = ['bicg','bicgstab','cg','cgs','gmres','qmr']
-from scipy.linalg import _iterative
-import numpy as sb
-import copy
-
-try:
-    False, True
-except NameError:
-    False, True = 0, 1
-
-_type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
-
-_coerce_rules = {('f','f'):'f', ('f','d'):'d', ('f','F'):'F',
-                 ('f','D'):'D', ('d','f'):'d', ('d','d'):'d',
-                 ('d','F'):'D', ('d','D'):'D', ('F','f'):'F',
-                 ('F','d'):'D', ('F','F'):'F', ('F','D'):'D',
-                 ('D','f'):'D', ('D','d'):'D', ('D','F'):'D',
-                 ('D','D'):'D'}
-
-class get_matvec:
-    methname = 'matvec'
-    def __init__(self, obj, *args):
-        self.obj = obj
-        self.args = args
-        if isinstance(obj, sb.matrix):
-            self.callfunc = self.type1m
-            return
-        if isinstance(obj, sb.ndarray):
-            self.callfunc = self.type1
-            return
-        meth = getattr(obj,self.methname,None)
-        if not callable(meth):
-            raise ValueError, "Object must be an array "\
-                  "or have a callable %s attribute." % (self.methname,)
-
-        self.obj = meth
-        self.callfunc = self.type2
-
-    def __call__(self, x):
-        return self.callfunc(x)
-
-    def type1(self, x):
-        return sb.dot(self.obj, x)
-
-    def type1m(self, x):
-        return sb.dot(self.obj.A, x)
-
-    def type2(self, x):
-        return self.obj(x,*self.args)
-
-class get_rmatvec(get_matvec):
-    methname = 'rmatvec'
-    def type1(self, x):
-        return sb.dot(x, self.obj)
-    def type1m(self, x):
-        return sb.dot(x, self.obj.A)
-
-class get_psolve:
-    methname = 'psolve'
-    def __init__(self, obj, *args):
-        self.obj = obj
-        self.args = args
-        meth = getattr(obj,self.methname,None)
-        if meth is None:  # no preconditiong available
-            self.callfunc = self.type1
-            return
-
-        if not callable(meth):
-            raise ValueError, "Preconditioning method %s "\
-                  "must be callable." % (self.methname,)
-
-        self.obj = meth
-        self.callfunc = self.type2
-
-    def __call__(self, x):
-        return self.callfunc(x)
-
-    def type1(self, x):
-        return x
-
-    def type2(self, x):
-        return self.obj(x,*self.args)
-
-class get_rpsolve(get_psolve):
-    methname = 'rpsolve'
-
-class get_psolveq(get_psolve):
-
-    def __call__(self, x, which):
-        return self.callfunc(x, which)
-
-    def type1(self, x, which):
-        return x
-
-    def type2(self, x, which):
-        return self.obj(x,which,*self.args)
-
-class get_rpsolveq(get_psolveq):
-    methname = 'rpsolve'
-
-def bicg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
-    """Use BIConjugate Gradient iteration to solve A x = b
-
-    Inputs:
-
-    A --   An array or an object with matvec(x) and rmatvec(x) methods
-           to represent A * x and A^H * x respectively.  May also have
-           psolve(b) and rpsolve(b) methods for representing solutions
-           to the preconditioning equations M * x = b and
-           M^H * x = b respectively.
-    b --   An n-length vector
-
-    Outputs:
-
-    x  --  The converged solution
-    info -- output result
-            0  : successful exit
-            >0 : convergence to tolerance not achieved, number of iterations
-            <0 : illegal input or breakdown
-
-    Optional Inputs:
-
-    x0  -- (0) default starting guess.
-    tol -- (1e-5) relative tolerance to achieve
-    maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be
-               determined from A.dtype.char and b.  If A does not have a
-               typecode method then it will compute A.matvec(x0) to get a
-               typecode.   To save the extra computation when A does not
-               have a typecode attribute use xtype=0 for the same type as
-               b or use xtype='f','d','F',or 'D'
-    callback -- an optional user-supplied function to call after each
-                iteration.  It is called as callback(xk), where xk is the
-                current parameter vector.
-    """
-    b = sb.asarray(b)+0.0
-    n = len(b)
-    if maxiter is None:
-        maxiter = n*10
-
-    if x0 is None:
-        x = sb.zeros(n)
-    else:
-        x = copy.copy(x0)
-
-    if xtype is None:
-        try:
-            atyp = A.dtype.char
-        except AttributeError:
-            atyp = None
-        if atyp is None:
-            atyp = A.matvec(x).dtype.char
-        typ = _coerce_rules[b.dtype.char,atyp]
-    elif xtype == 0:
-        typ = b.dtype.char
-    else:
-        typ = xtype
-        if typ not in 'fdFD':
-            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
-
-    x = sb.asarray(x,typ)
-    b = sb.asarray(b,typ)
-
-    matvec, psolve, rmatvec, rpsolve = (None,)*4
-    ltr = _type_conv[typ]
-    revcom = _iterative.__dict__[ltr+'bicgrevcom']
-    stoptest = _iterative.__dict__[ltr+'stoptest2']
-
-    resid = tol
-    ndx1 = 1
-    ndx2 = -1
-    work = sb.zeros(6*n,typ)
-    ijob = 1
-    info = 0
-    ftflag = True
-    bnrm2 = -1.0
-    iter_ = maxiter
-    while True:
-        olditer = iter_
-        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
-           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
-        if callback is not None and iter_ > olditer:
-            callback(x)
-        slice1 = slice(ndx1-1, ndx1-1+n)
-        slice2 = slice(ndx2-1, ndx2-1+n)
-        if (ijob == -1):
-            break
-        elif (ijob == 1):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(work[slice1])
-        elif (ijob == 2):
-            if rmatvec is None:
-                rmatvec = get_rmatvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*rmatvec(work[slice1])
-        elif (ijob == 3):
-            if psolve is None:
-                psolve = get_psolve(A)
-            work[slice1] = psolve(work[slice2])
-        elif (ijob == 4):
-            if rpsolve is None:
-                rpsolve = get_rpsolve(A)
-            work[slice1] = rpsolve(work[slice2])
-        elif (ijob == 5):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(x)
-        elif (ijob == 6):
-            if ftflag:
-                info = -1
-                ftflag = False
-            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
-        ijob = 2
-
-    return x, info
-
-
-def bicgstab(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
-    """Use BIConjugate Gradient STABilized iteration to solve A x = b
-
-    Inputs:
-
-    A --   An array or an object with a matvec(x) method
-           to represent A * x.  May also have a psolve(b) methods for
-           representing solution to the preconditioning equation
-           M * x = b.
-    b --   An n-length vector
-
-    Outputs:
-
-    x  --  The converged solution
-    info -- output result
-            0  : successful exit
-            >0 : convergence to tolerance not achieved, number of iterations
-            <0 : illegal input or breakdown
-
-    Optional Inputs:
-
-    x0  -- (0) default starting guess.
-    tol -- (1e-5) relative tolerance to achieve
-    maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be
-               determined from A.dtype.char and b.  If A does not have a
-               typecode method then it will compute A.matvec(x0) to get a
-               typecode.   To save the extra computation when A does not
-               have a typecode attribute use xtype=0 for the same type as
-               b or use xtype='f','d','F',or 'D'
-    callback -- an optional user-supplied function to call after each
-                iteration.  It is called as callback(xk), where xk is the
-                current parameter vector.
-    """
-    b = sb.asarray(b)+0.0
-    n = len(b)
-    if maxiter is None:
-        maxiter = n*10
-
-    if x0 is None:
-        x = sb.zeros(n)
-    else:
-        x = copy.copy(x0)
-
-
-    if xtype is None:
-        try:
-            atyp = A.dtype.char
-        except AttributeError:
-            atyp = None
-        if atyp is None:
-            atyp = A.matvec(x).dtype.char
-        typ = _coerce_rules[b.dtype.char,atyp]
-    elif xtype == 0:
-        typ = b.dtype.char
-    else:
-        typ = xtype
-        if typ not in 'fdFD':
-            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
-
-    x = sb.asarray(x,typ)
-    b = sb.asarray(b,typ)
-
-    matvec, psolve = (None,)*2
-    ltr = _type_conv[typ]
-    revcom = _iterative.__dict__[ltr+'bicgstabrevcom']
-    stoptest = _iterative.__dict__[ltr+'stoptest2']
-
-    resid = tol
-    ndx1 = 1
-    ndx2 = -1
-    work = sb.zeros(7*n,typ)
-    ijob = 1
-    info = 0
-    ftflag = True
-    bnrm2 = -1.0
-    iter_ = maxiter
-    while True:
-        olditer = iter_
-        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
-           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
-        if callback is not None and iter_ > olditer:
-            callback(x)
-        slice1 = slice(ndx1-1, ndx1-1+n)
-        slice2 = slice(ndx2-1, ndx2-1+n)
-        if (ijob == -1):
-            break
-        elif (ijob == 1):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(work[slice1])
-        elif (ijob == 2):
-            if psolve is None:
-                psolve = get_psolve(A)
-            work[slice1] = psolve(work[slice2])
-        elif (ijob == 3):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(x)
-        elif (ijob == 4):
-            if ftflag:
-                info = -1
-                ftflag = False
-            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
-        ijob = 2
-
-    return x, info
-
-
-def cg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
-    """Use Conjugate Gradient iteration to solve A x = b (A^H = A)
-
-    Inputs:
-
-    A --   An array or an object with a matvec(x) method
-           to represent A * x.  May also have a psolve(b) methods for
-           representing solution to the preconditioning equation
-           M * x = b.
-    b --   An n-length vector
-
-
-    Outputs:
-
-    x  --  The converged solution
-    info -- output result
-            0  : successful exit
-            >0 : convergence to tolerance not achieved, number of iterations
-            <0 : illegal input or breakdown
-
-    Optional Inputs:
-
-    x0  -- (0) default starting guess.
-    tol -- (1e-5) relative tolerance to achieve
-    maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be
-               determined from A.dtype.char and b.  If A does not have a
-               typecode method then it will compute A.matvec(x0) to get a
-               typecode.   To save the extra computation when A does not
-               have a typecode attribute use xtype=0 for the same type as
-               b or use xtype='f','d','F',or 'D'
-    callback -- an optional user-supplied function to call after each
-                iteration.  It is called as callback(xk), where xk is the
-                current parameter vector.
-    """
-    b = sb.asarray(b)+0.0
-    n = len(b)
-    if maxiter is None:
-        maxiter = n*10
-
-    if x0 is None:
-        x = sb.zeros(n)
-    else:
-        x = copy.copy(x0)
-
-
-    if xtype is None:
-        try:
-            atyp = A.dtype.char
-        except AttributeError:
-            atyp = None
-        if atyp is None:
-            atyp = A.matvec(x).dtype.char
-        typ = _coerce_rules[b.dtype.char,atyp]
-    elif xtype == 0:
-        typ = b.dtype.char
-    else:
-        typ = xtype
-        if typ not in 'fdFD':
-            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
-
-    x = sb.asarray(x,typ)
-    b = sb.asarray(b,typ)
-
-    matvec, psolve = (None,)*2
-    ltr = _type_conv[typ]
-    revcom = _iterative.__dict__[ltr+'cgrevcom']
-    stoptest = _iterative.__dict__[ltr+'stoptest2']
-
-    resid = tol
-    ndx1 = 1
-    ndx2 = -1
-    work = sb.zeros(4*n,typ)
-    ijob = 1
-    info = 0
-    ftflag = True
-    bnrm2 = -1.0
-    iter_ = maxiter
-    while True:
-        olditer = iter_
-        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
-           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
-        if callback is not None and iter_ > olditer:
-            callback(x)
-        slice1 = slice(ndx1-1, ndx1-1+n)
-        slice2 = slice(ndx2-1, ndx2-1+n)
-        if (ijob == -1):
-            break
-        elif (ijob == 1):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(work[slice1])
-        elif (ijob == 2):
-            if psolve is None:
-                psolve = get_psolve(A)
-            work[slice1] = psolve(work[slice2])
-        elif (ijob == 3):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(x)
-        elif (ijob == 4):
-            if ftflag:
-                info = -1
-                ftflag = False
-            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
-        ijob = 2
-
-    return x, info
-
-
-def cgs(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
-    """Use Conjugate Gradient Squared iteration to solve A x = b
-
-    Inputs:
-
-    A --   An array or an object with a matvec(x) method
-           to represent A * x.  May also have a psolve(b) methods for
-           representing solution to the preconditioning equation
-           M * x = b.
-    b --   An n-length vector
-
-
-    Outputs:
-
-    x  --  The converged solution
-    info -- output result
-            0  : successful exit
-            >0 : convergence to tolerance not achieved, number of iterations
-            <0 : illegal input or breakdown
-
-    Optional Inputs:
-
-    x0  -- (0) default starting guess.
-    tol -- (1e-5) relative tolerance to achieve
-    maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be
-               determined from A.dtype.char and b.  If A does not have a
-               typecode method then it will compute A.matvec(x0) to get a
-               typecode.   To save the extra computation when A does not
-               have a typecode attribute use xtype=0 for the same type as
-               b or use xtype='f','d','F',or 'D'
-    callback -- an optional user-supplied function to call after each
-                iteration.  It is called as callback(xk), where xk is the
-                current parameter vector.
-    """
-    b = sb.asarray(b) + 0.0
-    n = len(b)
-    if maxiter is None:
-        maxiter = n*10
-
-    if x0 is None:
-        x = sb.zeros(n)
-    else:
-        x = copy.copy(x0)
-
-    if xtype is None:
-        try:
-            atyp = A.dtype.char
-        except AttributeError:
-            atyp = None
-        if atyp is None:
-            atyp = A.matvec(x).dtype.char
-        typ = _coerce_rules[b.dtype.char,atyp]
-    elif xtype == 0:
-        typ = b.dtype.char
-    else:
-        typ = xtype
-        if typ not in 'fdFD':
-            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
-
-    x = sb.asarray(x,typ)
-    b = sb.asarray(b,typ)
-
-    matvec, psolve = (None,)*2
-    ltr = _type_conv[typ]
-    revcom = _iterative.__dict__[ltr+'cgsrevcom']
-    stoptest = _iterative.__dict__[ltr+'stoptest2']
-
-    resid = tol
-    ndx1 = 1
-    ndx2 = -1
-    work = sb.zeros(7*n,typ)
-    ijob = 1
-    info = 0
-    ftflag = True
-    bnrm2 = -1.0
-    iter_ = maxiter
-    while True:
-        olditer = iter_
-        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
-           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
-        if callback is not None and iter_ > olditer:
-            callback(x)
-        slice1 = slice(ndx1-1, ndx1-1+n)
-        slice2 = slice(ndx2-1, ndx2-1+n)
-        if (ijob == -1):
-            break
-        elif (ijob == 1):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(work[slice1])
-        elif (ijob == 2):
-            if psolve is None:
-                psolve = get_psolve(A)
-            work[slice1] = psolve(work[slice2])
-        elif (ijob == 3):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(x)
-        elif (ijob == 4):
-            if ftflag:
-                info = -1
-                ftflag = False
-            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
-        ijob = 2
-
-    return x, info
-
-
-def gmres(A, b, x0=None, tol=1e-5, restrt=None, maxiter=None, xtype=None, callback=None):
-    """Use Generalized Minimal RESidual iteration to solve A x = b
-
-    Inputs:
-
-    A --   An array or an object with a matvec(x) method
-           to represent A * x.  May also have a psolve(b) methods for
-           representing solution to the preconditioning equation
-           M * x = b.
-    b --   An n-length vector
-
-
-    Outputs:
-
-    x  --  The converged solution
-    info -- output result
-            0  : successful exit
-            >0 : convergence to tolerance not achieved, number of iterations
-            <0 : illegal input or breakdown
-
-    Optional Inputs:
-
-    x0  -- (0) default starting guess.
-    tol -- (1e-5) relative tolerance to achieve
-    restrt -- (n) When to restart (change this to get faster performance -- but
-                   may not converge).
-    maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be
-               determined from A.dtype.char and b.  If A does not have a
-               typecode method then it will compute A.matvec(x0) to get a
-               typecode.   To save the extra computation when A does not
-               have a typecode attribute use xtype=0 for the same type as
-               b or use xtype='f','d','F',or 'D'
-    callback -- an optional user-supplied function to call after each
-                iteration.  It is called as callback(xk), where xk is the
-                current parameter vector.
-    """
-    b = sb.asarray(b)+0.0
-    n = len(b)
-    if maxiter is None:
-        maxiter = n*10
-
-    if x0 is None:
-        x = sb.zeros(n)
-    else:
-        x = copy.copy(x0)
-
-    if xtype is None:
-        try:
-            atyp = A.dtype.char
-        except AttributeError:
-            atyp = A.matvec(x).dtype.char
-        typ = _coerce_rules[b.dtype.char,atyp]
-    elif xtype == 0:
-        typ = b.dtype.char
-    else:
-        typ = xtype
-        if typ not in 'fdFD':
-            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
-
-    x = sb.asarray(x,typ)
-    b = sb.asarray(b,typ)
-
-    matvec, psolve = (None,)*2
-    ltr = _type_conv[typ]
-    revcom = _iterative.__dict__[ltr+'gmresrevcom']
-    stoptest = _iterative.__dict__[ltr+'stoptest2']
-
-    if restrt is None:
-        restrt = n
-    resid = tol
-    ndx1 = 1
-    ndx2 = -1
-    work = sb.zeros((6+restrt)*n,typ)
-    work2 = sb.zeros((restrt+1)*(2*restrt+2),typ)
-    ijob = 1
-    info = 0
-    ftflag = True
-    bnrm2 = -1.0
-    iter_ = maxiter
-    while True:
-        olditer = iter_
-        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
-           revcom(b, x, restrt, work, work2, iter_, resid, info, ndx1, ndx2, ijob)
-        if callback is not None and iter_ > olditer:
-            callback(x)
-        slice1 = slice(ndx1-1, ndx1-1+n)
-        slice2 = slice(ndx2-1, ndx2-1+n)
-        if (ijob == -1):
-            break
-        elif (ijob == 1):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(x)
-        elif (ijob == 2):
-            if psolve is None:
-                psolve = get_psolve(A)
-            work[slice1] = psolve(work[slice2])
-        elif (ijob == 3):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(work[slice1])
-        elif (ijob == 4):
-            if ftflag:
-                info = -1
-                ftflag = False
-            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
-        ijob = 2
-
-    return x, info
-
-
-def qmr(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
-    """Use Quasi-Minimal Residual iteration to solve A x = b
-
-    Inputs:
-
-    A --   An array or an object with matvec(x) and rmatvec(x) methods
-           to represent A * x and A^H * x respectively.  May also have
-           psolve(b,<which>) and rpsolve(b,<which>) methods for
-           representing solutions to the preconditioning equations
-           M * x = b and M^H * x = b respectively.   The <which> argument
-           may be given to specify 'left' or 'right' preconditioning.
-    b --   An n-length vector
-
-    Outputs:
-
-    x  --  The converged solution
-    info -- output result
-            0  : successful exit
-            >0 : convergence to tolerance not achieved, number of iterations
-            <0 : illegal input or breakdown
-
-    Optional Inputs:
-
-    x0  -- (0) default starting guess.
-    tol -- (1e-5) relative tolerance to achieve
-    maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be
-               determined from A.dtype.char and b.  If A does not have a
-               typecode method then it will compute A.matvec(x0) to get a
-               typecode.   To save the extra computation when A does not
-               have a typecode attribute use xtype=0 for the same type as
-               b or use xtype='f','d','F',or 'D'
-    callback -- an optional user-supplied function to call after each
-                iteration.  It is called as callback(xk), where xk is the
-                current parameter vector.
-    """
-    b = sb.asarray(b)+0.0
-    n = len(b)
-    if maxiter is None:
-        maxiter = n*10
-
-    if x0 is None:
-        x = sb.zeros(n)
-    else:
-        x = copy.copy(x0)
-
-    if xtype is None:
-        try:
-            atyp = A.dtype.char
-        except AttributeError:
-            atyp = None
-        if atyp is None:
-            atyp = A.matvec(x).dtype.char
-        typ = _coerce_rules[b.dtype.char,atyp]
-    elif xtype == 0:
-        typ = b.dtype.char
-    else:
-        typ = xtype
-        if typ not in 'fdFD':
-            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
-
-    x = sb.asarray(x,typ)
-    b = sb.asarray(b,typ)
-
-
-    matvec, psolve, rmatvec, rpsolve = (None,)*4
-    ltr = _type_conv[typ]
-    revcom = _iterative.__dict__[ltr+'qmrrevcom']
-    stoptest = _iterative.__dict__[ltr+'stoptest2']
-
-    resid = tol
-    ndx1 = 1
-    ndx2 = -1
-    work = sb.zeros(11*n,typ)
-    ijob = 1
-    info = 0
-    ftflag = True
-    bnrm2 = -1.0
-    iter_ = maxiter
-    while True:
-        olditer = iter_
-        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
-           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
-        if callback is not None and iter_ > olditer:
-            callback(x)
-        slice1 = slice(ndx1-1, ndx1-1+n)
-        slice2 = slice(ndx2-1, ndx2-1+n)
-        if (ijob == -1):
-            break
-        elif (ijob == 1):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(work[slice1])
-        elif (ijob == 2):
-            if rmatvec is None:
-                rmatvec = get_rmatvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*rmatvec(work[slice1])
-        elif (ijob == 3):
-            if psolve is None:
-                psolve = get_psolveq(A)
-            work[slice1] = psolve(work[slice2],'left')
-        elif (ijob == 4):
-            if psolve is None:
-                psolve = get_psolveq(A)
-            work[slice1] = psolve(work[slice2],'right')
-        elif (ijob == 5):
-            if rpsolve is None:
-                rpsolve = get_rpsolveq(A)
-            work[slice1] = rpsolve(work[slice2],'left')
-        elif (ijob == 6):
-            if rpsolve is None:
-                rpsolve = get_rpsolveq(A)
-            work[slice1] = rpsolve(work[slice2],'right')
-        elif (ijob == 7):
-            if matvec is None:
-                matvec = get_matvec(A)
-            work[slice2] *= sclr2
-            work[slice2] += sclr1*matvec(x)
-        elif (ijob == 8):
-            if ftflag:
-                info = -1
-                ftflag = False
-            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
-        ijob = 2
-
-    return x, info

Modified: trunk/scipy/linalg/setup.py
===================================================================
--- trunk/scipy/linalg/setup.py	2008-01-24 22:11:36 UTC (rev 3860)
+++ trunk/scipy/linalg/setup.py	2008-01-25 20:00:52 UTC (rev 3861)
@@ -183,23 +183,23 @@
                          extra_info = lapack_opt
                          )
 
-    # iterative methods
-    methods = ['BiCGREVCOM.f.src',
-               'BiCGSTABREVCOM.f.src',
-               'CGREVCOM.f.src',
-               'CGSREVCOM.f.src',
-#               'ChebyREVCOM.f.src',
-               'GMRESREVCOM.f.src',
-#               'JacobiREVCOM.f.src',
-               'QMRREVCOM.f.src',
-#               'SORREVCOM.f.src'
-               ]
-    Util = ['STOPTEST2.f.src','getbreak.f.src']
-    sources = Util + methods + ['_iterative.pyf.src']
-    config.add_extension('_iterative',
-                         sources = [join('iterative',x) for x in sources],
-                         extra_info = lapack_opt
-                         )
+#    # iterative methods
+#    methods = ['BiCGREVCOM.f.src',
+#               'BiCGSTABREVCOM.f.src',
+#               'CGREVCOM.f.src',
+#               'CGSREVCOM.f.src',
+##               'ChebyREVCOM.f.src',
+#               'GMRESREVCOM.f.src',
+##               'JacobiREVCOM.f.src',
+#               'QMRREVCOM.f.src',
+##               'SORREVCOM.f.src'
+#               ]
+#    Util = ['STOPTEST2.f.src','getbreak.f.src']
+#    sources = Util + methods + ['_iterative.pyf.src']
+#    config.add_extension('_iterative',
+#                         sources = [join('iterative',x) for x in sources],
+#                         extra_info = lapack_opt
+#                         )
 
     config.add_data_dir('tests')
 

Deleted: trunk/scipy/linalg/tests/test_iterative.py
===================================================================
--- trunk/scipy/linalg/tests/test_iterative.py	2008-01-24 22:11:36 UTC (rev 3860)
+++ trunk/scipy/linalg/tests/test_iterative.py	2008-01-25 20:00:52 UTC (rev 3861)
@@ -1,83 +0,0 @@
-#!/usr/bin/env python
-#
-# Created by: Ed Schofield, Jan 2007
-#
-""" Test functions for the linalg.iterative module
-"""
-__usage__ = """
-Build linalg:
-  python setup_linalg.py build
-Run tests if scipy is installed:
-  python -c 'import scipy;scipy.linalg.test(<level>)'
-Run tests if linalg is not installed:
-  python tests/test_iterative.py [<level>]
-"""
-
-import sys
-
-from numpy import zeros, dot, diag, ones
-from scipy.testing import *
-from numpy.random import rand
-#from numpy import arange, add, array, dot, zeros, identity, conjugate, transpose
-
-from scipy.linalg import iterative, norm, cg, cgs, bicg, bicgstab, gmres, qmr
-
-
-def callback(x):
-    global A, b
-    res = b-dot(A,x)
-    #print "||A.x - b|| = " + str(norm(dot(A,x)-b))
-
-class TestIterativeSolvers(TestCase):
-    def __init__(self, *args, **kwds):
-        TestCase.__init__(self, *args, **kwds)
-        self.setUp()
-    def setUp (self):
-        global A, b
-        n = 10
-        self.tol = 1e-5
-        self.x0 = zeros(n, float)
-        A = rand(n, n)+diag(4*ones(n))
-        self.A = 0.5 * (A+A.T)
-        A = self.A
-        self.b = rand(n)
-        b = self.b
-
-    def test_cg(self):
-        bx0 = self.x0.copy()
-        x, info = cg(self.A, self.b, self.x0, callback=callback)
-        assert_array_equal(bx0, self.x0)
-        assert norm(dot(self.A, x) - self.b) < 5*self.tol
-
-    def test_bicg(self):
-        bx0 = self.x0.copy()
-        x, info = bicg(self.A, self.b, self.x0, callback=callback)
-        assert_array_equal(bx0, self.x0)
-        assert norm(dot(self.A, x) - self.b) < 5*self.tol
-
-    def test_cgs(self):
-        bx0 = self.x0.copy()
-        x, info = cgs(self.A, self.b, self.x0, callback=callback)
-        assert_array_equal(bx0, self.x0)
-        assert norm(dot(self.A, x) - self.b) < 5*self.tol
-
-    def test_bicgstab(self):
-        bx0 = self.x0.copy()
-        x, info = bicgstab(self.A, self.b, self.x0, callback=callback)
-        assert_array_equal(bx0, self.x0)
-        assert norm(dot(self.A, x) - self.b) < 5*self.tol
-
-    def test_gmres(self):
-        bx0 = self.x0.copy()
-        x, info = gmres(self.A, self.b, self.x0, callback=callback)
-        assert_array_equal(bx0, self.x0)
-        assert norm(dot(self.A, x) - self.b) < 5*self.tol
-
-    def test_qmr(self):
-        bx0 = self.x0.copy()
-        x, info = qmr(self.A, self.b, self.x0, callback=callback)
-        assert_array_equal(bx0, self.x0)
-        assert norm(dot(self.A, x) - self.b) < 5*self.tol
-
-if __name__ == "__main__":
-    nose.run(argv=['', __file__])

Modified: trunk/scipy/sandbox/arpack/tests/test_speigs.py
===================================================================
--- trunk/scipy/sandbox/arpack/tests/test_speigs.py	2008-01-24 22:11:36 UTC (rev 3860)
+++ trunk/scipy/sandbox/arpack/tests/test_speigs.py	2008-01-25 20:00:52 UTC (rev 3861)
@@ -22,7 +22,7 @@
         vals = vals[uv_sortind]
         vecs = vecs[:,uv_sortind]
 
-        from scipy.linalg.iterative import get_matvec
+        from scipy.splinalg.isolve.iterative import get_matvec
         matvec = get_matvec(A)
         #= lambda x: N.asarray(A*x)[0]
         nev=4

Modified: trunk/scipy/setup.py
===================================================================
--- trunk/scipy/setup.py	2008-01-24 22:11:36 UTC (rev 3860)
+++ trunk/scipy/setup.py	2008-01-25 20:00:52 UTC (rev 3861)
@@ -18,6 +18,7 @@
     config.add_subpackage('signal')
     config.add_subpackage('sparse')
     config.add_subpackage('special')
+    config.add_subpackage('splinalg')
     config.add_subpackage('stats')
     config.add_subpackage('ndimage')
     config.add_subpackage('stsci')

Copied: trunk/scipy/splinalg/isolve/iterative (from rev 3860, trunk/scipy/linalg/iterative)

Copied: trunk/scipy/splinalg/isolve/iterative.py (from rev 3860, trunk/scipy/linalg/iterative.py)
===================================================================
--- trunk/scipy/linalg/iterative.py	2008-01-24 22:11:36 UTC (rev 3860)
+++ trunk/scipy/splinalg/isolve/iterative.py	2008-01-25 20:00:52 UTC (rev 3861)
@@ -0,0 +1,808 @@
+## Automatically adapted for scipy Oct 18, 2005 by
+
+
+# Iterative methods using reverse-communication raw material
+#   These methods solve
+#   Ax = b  for x
+
+#   where A must have A.matvec(x,*args) defined
+#    or be a numeric array
+
+
+__all__ = ['bicg','bicgstab','cg','cgs','gmres','qmr']
+
+import _iterative
+import numpy as sb
+import copy
+
+try:
+    False, True
+except NameError:
+    False, True = 0, 1
+
+_type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
+
+_coerce_rules = {('f','f'):'f', ('f','d'):'d', ('f','F'):'F',
+                 ('f','D'):'D', ('d','f'):'d', ('d','d'):'d',
+                 ('d','F'):'D', ('d','D'):'D', ('F','f'):'F',
+                 ('F','d'):'D', ('F','F'):'F', ('F','D'):'D',
+                 ('D','f'):'D', ('D','d'):'D', ('D','F'):'D',
+                 ('D','D'):'D'}
+
+class get_matvec:
+    methname = 'matvec'
+    def __init__(self, obj, *args):
+        self.obj = obj
+        self.args = args
+        if isinstance(obj, sb.matrix):
+            self.callfunc = self.type1m
+            return
+        if isinstance(obj, sb.ndarray):
+            self.callfunc = self.type1
+            return
+        meth = getattr(obj,self.methname,None)
+        if not callable(meth):
+            raise ValueError, "Object must be an array "\
+                  "or have a callable %s attribute." % (self.methname,)
+
+        self.obj = meth
+        self.callfunc = self.type2
+
+    def __call__(self, x):
+        return self.callfunc(x)
+
+    def type1(self, x):
+        return sb.dot(self.obj, x)
+
+    def type1m(self, x):
+        return sb.dot(self.obj.A, x)
+
+    def type2(self, x):
+        return self.obj(x,*self.args)
+
+class get_rmatvec(get_matvec):
+    methname = 'rmatvec'
+    def type1(self, x):
+        return sb.dot(x, self.obj)
+    def type1m(self, x):
+        return sb.dot(x, self.obj.A)
+
+class get_psolve:
+    methname = 'psolve'
+    def __init__(self, obj, *args):
+        self.obj = obj
+        self.args = args
+        meth = getattr(obj,self.methname,None)
+        if meth is None:  # no preconditiong available
+            self.callfunc = self.type1
+            return
+
+        if not callable(meth):
+            raise ValueError, "Preconditioning method %s "\
+                  "must be callable." % (self.methname,)
+
+        self.obj = meth
+        self.callfunc = self.type2
+
+    def __call__(self, x):
+        return self.callfunc(x)
+
+    def type1(self, x):
+        return x
+
+    def type2(self, x):
+        return self.obj(x,*self.args)
+
+class get_rpsolve(get_psolve):
+    methname = 'rpsolve'
+
+class get_psolveq(get_psolve):
+
+    def __call__(self, x, which):
+        return self.callfunc(x, which)
+
+    def type1(self, x, which):
+        return x
+
+    def type2(self, x, which):
+        return self.obj(x,which,*self.args)
+
+class get_rpsolveq(get_psolveq):
+    methname = 'rpsolve'
+
+def bicg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
+    """Use BIConjugate Gradient iteration to solve A x = b
+
+    Inputs:
+
+    A --   An array or an object with matvec(x) and rmatvec(x) methods
+           to represent A * x and A^H * x respectively.  May also have
+           psolve(b) and rpsolve(b) methods for representing solutions
+           to the preconditioning equations M * x = b and
+           M^H * x = b respectively.
+    b --   An n-length vector
+
+    Outputs:
+
+    x  --  The converged solution
+    info -- output result
+            0  : successful exit
+            >0 : convergence to tolerance not achieved, number of iterations
+            <0 : illegal input or breakdown
+
+    Optional Inputs:
+
+    x0  -- (0) default starting guess.
+    tol -- (1e-5) relative tolerance to achieve
+    maxiter -- (10*n) maximum number of iterations
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
+    """
+    b = sb.asarray(b)+0.0
+    n = len(b)
+    if maxiter is None:
+        maxiter = n*10
+
+    if x0 is None:
+        x = sb.zeros(n)
+    else:
+        x = copy.copy(x0)
+
+    if xtype is None:
+        try:
+            atyp = A.dtype.char
+        except AttributeError:
+            atyp = None
+        if atyp is None:
+            atyp = A.matvec(x).dtype.char
+        typ = _coerce_rules[b.dtype.char,atyp]
+    elif xtype == 0:
+        typ = b.dtype.char
+    else:
+        typ = xtype
+        if typ not in 'fdFD':
+            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
+
+    x = sb.asarray(x,typ)
+    b = sb.asarray(b,typ)
+
+    matvec, psolve, rmatvec, rpsolve = (None,)*4
+    ltr = _type_conv[typ]
+    revcom = _iterative.__dict__[ltr+'bicgrevcom']
+    stoptest = _iterative.__dict__[ltr+'stoptest2']
+
+    resid = tol
+    ndx1 = 1
+    ndx2 = -1
+    work = sb.zeros(6*n,typ)
+    ijob = 1
+    info = 0
+    ftflag = True
+    bnrm2 = -1.0
+    iter_ = maxiter
+    while True:
+        olditer = iter_
+        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
+           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
+        slice1 = slice(ndx1-1, ndx1-1+n)
+        slice2 = slice(ndx2-1, ndx2-1+n)
+        if (ijob == -1):
+            break
+        elif (ijob == 1):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(work[slice1])
+        elif (ijob == 2):
+            if rmatvec is None:
+                rmatvec = get_rmatvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*rmatvec(work[slice1])
+        elif (ijob == 3):
+            if psolve is None:
+                psolve = get_psolve(A)
+            work[slice1] = psolve(work[slice2])
+        elif (ijob == 4):
+            if rpsolve is None:
+                rpsolve = get_rpsolve(A)
+            work[slice1] = rpsolve(work[slice2])
+        elif (ijob == 5):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(x)
+        elif (ijob == 6):
+            if ftflag:
+                info = -1
+                ftflag = False
+            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
+        ijob = 2
+
+    return x, info
+
+
+def bicgstab(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
+    """Use BIConjugate Gradient STABilized iteration to solve A x = b
+
+    Inputs:
+
+    A --   An array or an object with a matvec(x) method
+           to represent A * x.  May also have a psolve(b) methods for
+           representing solution to the preconditioning equation
+           M * x = b.
+    b --   An n-length vector
+
+    Outputs:
+
+    x  --  The converged solution
+    info -- output result
+            0  : successful exit
+            >0 : convergence to tolerance not achieved, number of iterations
+            <0 : illegal input or breakdown
+
+    Optional Inputs:
+
+    x0  -- (0) default starting guess.
+    tol -- (1e-5) relative tolerance to achieve
+    maxiter -- (10*n) maximum number of iterations
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
+    """
+    b = sb.asarray(b)+0.0
+    n = len(b)
+    if maxiter is None:
+        maxiter = n*10
+
+    if x0 is None:
+        x = sb.zeros(n)
+    else:
+        x = copy.copy(x0)
+
+
+    if xtype is None:
+        try:
+            atyp = A.dtype.char
+        except AttributeError:
+            atyp = None
+        if atyp is None:
+            atyp = A.matvec(x).dtype.char
+        typ = _coerce_rules[b.dtype.char,atyp]
+    elif xtype == 0:
+        typ = b.dtype.char
+    else:
+        typ = xtype
+        if typ not in 'fdFD':
+            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
+
+    x = sb.asarray(x,typ)
+    b = sb.asarray(b,typ)
+
+    matvec, psolve = (None,)*2
+    ltr = _type_conv[typ]
+    revcom = _iterative.__dict__[ltr+'bicgstabrevcom']
+    stoptest = _iterative.__dict__[ltr+'stoptest2']
+
+    resid = tol
+    ndx1 = 1
+    ndx2 = -1
+    work = sb.zeros(7*n,typ)
+    ijob = 1
+    info = 0
+    ftflag = True
+    bnrm2 = -1.0
+    iter_ = maxiter
+    while True:
+        olditer = iter_
+        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
+           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
+        slice1 = slice(ndx1-1, ndx1-1+n)
+        slice2 = slice(ndx2-1, ndx2-1+n)
+        if (ijob == -1):
+            break
+        elif (ijob == 1):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(work[slice1])
+        elif (ijob == 2):
+            if psolve is None:
+                psolve = get_psolve(A)
+            work[slice1] = psolve(work[slice2])
+        elif (ijob == 3):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(x)
+        elif (ijob == 4):
+            if ftflag:
+                info = -1
+                ftflag = False
+            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
+        ijob = 2
+
+    return x, info
+
+
+def cg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
+    """Use Conjugate Gradient iteration to solve A x = b (A^H = A)
+
+    Inputs:
+
+    A --   An array or an object with a matvec(x) method
+           to represent A * x.  May also have a psolve(b) methods for
+           representing solution to the preconditioning equation
+           M * x = b.
+    b --   An n-length vector
+
+
+    Outputs:
+
+    x  --  The converged solution
+    info -- output result
+            0  : successful exit
+            >0 : convergence to tolerance not achieved, number of iterations
+            <0 : illegal input or breakdown
+
+    Optional Inputs:
+
+    x0  -- (0) default starting guess.
+    tol -- (1e-5) relative tolerance to achieve
+    maxiter -- (10*n) maximum number of iterations
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
+    """
+    b = sb.asarray(b)+0.0
+    n = len(b)
+    if maxiter is None:
+        maxiter = n*10
+
+    if x0 is None:
+        x = sb.zeros(n)
+    else:
+        x = copy.copy(x0)
+
+
+    if xtype is None:
+        try:
+            atyp = A.dtype.char
+        except AttributeError:
+            atyp = None
+        if atyp is None:
+            atyp = A.matvec(x).dtype.char
+        typ = _coerce_rules[b.dtype.char,atyp]
+    elif xtype == 0:
+        typ = b.dtype.char
+    else:
+        typ = xtype
+        if typ not in 'fdFD':
+            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
+
+    x = sb.asarray(x,typ)
+    b = sb.asarray(b,typ)
+
+    matvec, psolve = (None,)*2
+    ltr = _type_conv[typ]
+    revcom = _iterative.__dict__[ltr+'cgrevcom']
+    stoptest = _iterative.__dict__[ltr+'stoptest2']
+
+    resid = tol
+    ndx1 = 1
+    ndx2 = -1
+    work = sb.zeros(4*n,typ)
+    ijob = 1
+    info = 0
+    ftflag = True
+    bnrm2 = -1.0
+    iter_ = maxiter
+    while True:
+        olditer = iter_
+        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
+           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
+        slice1 = slice(ndx1-1, ndx1-1+n)
+        slice2 = slice(ndx2-1, ndx2-1+n)
+        if (ijob == -1):
+            break
+        elif (ijob == 1):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(work[slice1])
+        elif (ijob == 2):
+            if psolve is None:
+                psolve = get_psolve(A)
+            work[slice1] = psolve(work[slice2])
+        elif (ijob == 3):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(x)
+        elif (ijob == 4):
+            if ftflag:
+                info = -1
+                ftflag = False
+            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
+        ijob = 2
+
+    return x, info
+
+
+def cgs(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
+    """Use Conjugate Gradient Squared iteration to solve A x = b
+
+    Inputs:
+
+    A --   An array or an object with a matvec(x) method
+           to represent A * x.  May also have a psolve(b) methods for
+           representing solution to the preconditioning equation
+           M * x = b.
+    b --   An n-length vector
+
+
+    Outputs:
+
+    x  --  The converged solution
+    info -- output result
+            0  : successful exit
+            >0 : convergence to tolerance not achieved, number of iterations
+            <0 : illegal input or breakdown
+
+    Optional Inputs:
+
+    x0  -- (0) default starting guess.
+    tol -- (1e-5) relative tolerance to achieve
+    maxiter -- (10*n) maximum number of iterations
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
+    """
+    b = sb.asarray(b) + 0.0
+    n = len(b)
+    if maxiter is None:
+        maxiter = n*10
+
+    if x0 is None:
+        x = sb.zeros(n)
+    else:
+        x = copy.copy(x0)
+
+    if xtype is None:
+        try:
+            atyp = A.dtype.char
+        except AttributeError:
+            atyp = None
+        if atyp is None:
+            atyp = A.matvec(x).dtype.char
+        typ = _coerce_rules[b.dtype.char,atyp]
+    elif xtype == 0:
+        typ = b.dtype.char
+    else:
+        typ = xtype
+        if typ not in 'fdFD':
+            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
+
+    x = sb.asarray(x,typ)
+    b = sb.asarray(b,typ)
+
+    matvec, psolve = (None,)*2
+    ltr = _type_conv[typ]
+    revcom = _iterative.__dict__[ltr+'cgsrevcom']
+    stoptest = _iterative.__dict__[ltr+'stoptest2']
+
+    resid = tol
+    ndx1 = 1
+    ndx2 = -1
+    work = sb.zeros(7*n,typ)
+    ijob = 1
+    info = 0
+    ftflag = True
+    bnrm2 = -1.0
+    iter_ = maxiter
+    while True:
+        olditer = iter_
+        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
+           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
+        slice1 = slice(ndx1-1, ndx1-1+n)
+        slice2 = slice(ndx2-1, ndx2-1+n)
+        if (ijob == -1):
+            break
+        elif (ijob == 1):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(work[slice1])
+        elif (ijob == 2):
+            if psolve is None:
+                psolve = get_psolve(A)
+            work[slice1] = psolve(work[slice2])
+        elif (ijob == 3):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(x)
+        elif (ijob == 4):
+            if ftflag:
+                info = -1
+                ftflag = False
+            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
+        ijob = 2
+
+    return x, info
+
+
+def gmres(A, b, x0=None, tol=1e-5, restrt=None, maxiter=None, xtype=None, callback=None):
+    """Use Generalized Minimal RESidual iteration to solve A x = b
+
+    Inputs:
+
+    A --   An array or an object with a matvec(x) method
+           to represent A * x.  May also have a psolve(b) methods for
+           representing solution to the preconditioning equation
+           M * x = b.
+    b --   An n-length vector
+
+
+    Outputs:
+
+    x  --  The converged solution
+    info -- output result
+            0  : successful exit
+            >0 : convergence to tolerance not achieved, number of iterations
+            <0 : illegal input or breakdown
+
+    Optional Inputs:
+
+    x0  -- (0) default starting guess.
+    tol -- (1e-5) relative tolerance to achieve
+    restrt -- (n) When to restart (change this to get faster performance -- but
+                   may not converge).
+    maxiter -- (10*n) maximum number of iterations
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
+    """
+    b = sb.asarray(b)+0.0
+    n = len(b)
+    if maxiter is None:
+        maxiter = n*10
+
+    if x0 is None:
+        x = sb.zeros(n)
+    else:
+        x = copy.copy(x0)
+
+    if xtype is None:
+        try:
+            atyp = A.dtype.char
+        except AttributeError:
+            atyp = A.matvec(x).dtype.char
+        typ = _coerce_rules[b.dtype.char,atyp]
+    elif xtype == 0:
+        typ = b.dtype.char
+    else:
+        typ = xtype
+        if typ not in 'fdFD':
+            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
+
+    x = sb.asarray(x,typ)
+    b = sb.asarray(b,typ)
+
+    matvec, psolve = (None,)*2
+    ltr = _type_conv[typ]
+    revcom = _iterative.__dict__[ltr+'gmresrevcom']
+    stoptest = _iterative.__dict__[ltr+'stoptest2']
+
+    if restrt is None:
+        restrt = n
+    resid = tol
+    ndx1 = 1
+    ndx2 = -1
+    work = sb.zeros((6+restrt)*n,typ)
+    work2 = sb.zeros((restrt+1)*(2*restrt+2),typ)
+    ijob = 1
+    info = 0
+    ftflag = True
+    bnrm2 = -1.0
+    iter_ = maxiter
+    while True:
+        olditer = iter_
+        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
+           revcom(b, x, restrt, work, work2, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
+        slice1 = slice(ndx1-1, ndx1-1+n)
+        slice2 = slice(ndx2-1, ndx2-1+n)
+        if (ijob == -1):
+            break
+        elif (ijob == 1):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(x)
+        elif (ijob == 2):
+            if psolve is None:
+                psolve = get_psolve(A)
+            work[slice1] = psolve(work[slice2])
+        elif (ijob == 3):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(work[slice1])
+        elif (ijob == 4):
+            if ftflag:
+                info = -1
+                ftflag = False
+            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
+        ijob = 2
+
+    return x, info
+
+
+def qmr(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
+    """Use Quasi-Minimal Residual iteration to solve A x = b
+
+    Inputs:
+
+    A --   An array or an object with matvec(x) and rmatvec(x) methods
+           to represent A * x and A^H * x respectively.  May also have
+           psolve(b,<which>) and rpsolve(b,<which>) methods for
+           representing solutions to the preconditioning equations
+           M * x = b and M^H * x = b respectively.   The <which> argument
+           may be given to specify 'left' or 'right' preconditioning.
+    b --   An n-length vector
+
+    Outputs:
+
+    x  --  The converged solution
+    info -- output result
+            0  : successful exit
+            >0 : convergence to tolerance not achieved, number of iterations
+            <0 : illegal input or breakdown
+
+    Optional Inputs:
+
+    x0  -- (0) default starting guess.
+    tol -- (1e-5) relative tolerance to achieve
+    maxiter -- (10*n) maximum number of iterations
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
+    """
+    b = sb.asarray(b)+0.0
+    n = len(b)
+    if maxiter is None:
+        maxiter = n*10
+
+    if x0 is None:
+        x = sb.zeros(n)
+    else:
+        x = copy.copy(x0)
+
+    if xtype is None:
+        try:
+            atyp = A.dtype.char
+        except AttributeError:
+            atyp = None
+        if atyp is None:
+            atyp = A.matvec(x).dtype.char
+        typ = _coerce_rules[b.dtype.char,atyp]
+    elif xtype == 0:
+        typ = b.dtype.char
+    else:
+        typ = xtype
+        if typ not in 'fdFD':
+            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
+
+    x = sb.asarray(x,typ)
+    b = sb.asarray(b,typ)
+
+
+    matvec, psolve, rmatvec, rpsolve = (None,)*4
+    ltr = _type_conv[typ]
+    revcom = _iterative.__dict__[ltr+'qmrrevcom']
+    stoptest = _iterative.__dict__[ltr+'stoptest2']
+
+    resid = tol
+    ndx1 = 1
+    ndx2 = -1
+    work = sb.zeros(11*n,typ)
+    ijob = 1
+    info = 0
+    ftflag = True
+    bnrm2 = -1.0
+    iter_ = maxiter
+    while True:
+        olditer = iter_
+        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
+           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
+        slice1 = slice(ndx1-1, ndx1-1+n)
+        slice2 = slice(ndx2-1, ndx2-1+n)
+        if (ijob == -1):
+            break
+        elif (ijob == 1):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(work[slice1])
+        elif (ijob == 2):
+            if rmatvec is None:
+                rmatvec = get_rmatvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*rmatvec(work[slice1])
+        elif (ijob == 3):
+            if psolve is None:
+                psolve = get_psolveq(A)
+            work[slice1] = psolve(work[slice2],'left')
+        elif (ijob == 4):
+            if psolve is None:
+                psolve = get_psolveq(A)
+            work[slice1] = psolve(work[slice2],'right')
+        elif (ijob == 5):
+            if rpsolve is None:
+                rpsolve = get_rpsolveq(A)
+            work[slice1] = rpsolve(work[slice2],'left')
+        elif (ijob == 6):
+            if rpsolve is None:
+                rpsolve = get_rpsolveq(A)
+            work[slice1] = rpsolve(work[slice2],'right')
+        elif (ijob == 7):
+            if matvec is None:
+                matvec = get_matvec(A)
+            work[slice2] *= sclr2
+            work[slice2] += sclr1*matvec(x)
+        elif (ijob == 8):
+            if ftflag:
+                info = -1
+                ftflag = False
+            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
+        ijob = 2
+
+    return x, info

Copied: trunk/scipy/splinalg/isolve/tests/test_iterative.py (from rev 3860, trunk/scipy/linalg/tests/test_iterative.py)
===================================================================
--- trunk/scipy/linalg/tests/test_iterative.py	2008-01-24 22:11:36 UTC (rev 3860)
+++ trunk/scipy/splinalg/isolve/tests/test_iterative.py	2008-01-25 20:00:52 UTC (rev 3861)
@@ -0,0 +1,73 @@
+#!/usr/bin/env python
+""" Test functions for the splinalg.isolve.iterative module
+"""
+
+import sys
+
+from numpy import zeros, dot, diag, ones
+from scipy.testing import *
+from numpy.random import rand
+#from numpy import arange, add, array, dot, zeros, identity, conjugate, transpose
+
+from scipy.linalg import norm
+from scipy.splinalg.isolve import cg, cgs, bicg, bicgstab, gmres, qmr
+
+
+def callback(x):
+    global A, b
+    res = b-dot(A,x)
+    #print "||A.x - b|| = " + str(norm(dot(A,x)-b))
+
+class TestIterativeSolvers(TestCase):
+    def __init__(self, *args, **kwds):
+        TestCase.__init__(self, *args, **kwds)
+        self.setUp()
+    def setUp (self):
+        global A, b
+        n = 10
+        self.tol = 1e-5
+        self.x0 = zeros(n, float)
+        A = rand(n, n)+diag(4*ones(n))
+        self.A = 0.5 * (A+A.T)
+        A = self.A
+        self.b = rand(n)
+        b = self.b
+
+    def test_cg(self):
+        bx0 = self.x0.copy()
+        x, info = cg(self.A, self.b, self.x0, callback=callback)
+        assert_array_equal(bx0, self.x0)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def test_bicg(self):
+        bx0 = self.x0.copy()
+        x, info = bicg(self.A, self.b, self.x0, callback=callback)
+        assert_array_equal(bx0, self.x0)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def test_cgs(self):
+        bx0 = self.x0.copy()
+        x, info = cgs(self.A, self.b, self.x0, callback=callback)
+        assert_array_equal(bx0, self.x0)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def test_bicgstab(self):
+        bx0 = self.x0.copy()
+        x, info = bicgstab(self.A, self.b, self.x0, callback=callback)
+        assert_array_equal(bx0, self.x0)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def test_gmres(self):
+        bx0 = self.x0.copy()
+        x, info = gmres(self.A, self.b, self.x0, callback=callback)
+        assert_array_equal(bx0, self.x0)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def test_qmr(self):
+        bx0 = self.x0.copy()
+        x, info = qmr(self.A, self.b, self.x0, callback=callback)
+        assert_array_equal(bx0, self.x0)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+if __name__ == "__main__":
+    nose.run(argv=['', __file__])



More information about the Scipy-svn mailing list