[Scipy-svn] r3841 - in trunk/scipy: linsolve sandbox/multigrid sparse/tests
scipy-svn@scip...
scipy-svn@scip...
Wed Jan 16 06:27:47 CST 2008
Author: wnbell
Date: 2008-01-16 06:27:45 -0600 (Wed, 16 Jan 2008)
New Revision: 3841
Modified:
trunk/scipy/linsolve/linsolve.py
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sparse/tests/test_base.py
Log:
cleaned up linsolve
Modified: trunk/scipy/linsolve/linsolve.py
===================================================================
--- trunk/scipy/linsolve/linsolve.py 2008-01-16 12:05:14 UTC (rev 3840)
+++ trunk/scipy/linsolve/linsolve.py 2008-01-16 12:27:45 UTC (rev 3841)
@@ -1,7 +1,11 @@
-from scipy.sparse import isspmatrix_csc, isspmatrix_csr, isspmatrix, spdiags
-import _superlu
+from warnings import warn
+
from numpy import asarray
+from scipy.sparse import isspmatrix_csc, isspmatrix_csr, isspmatrix, \
+ SparseEfficiencyWarning, csc_matrix
+import _superlu
+
import umfpack
if hasattr( umfpack, 'UMFPACK_OK' ):
isUmfpack = True
@@ -10,6 +14,7 @@
isUmfpack = False
useUmfpack = True
+
#convert numpy char to superLU char
superLU_transtabl = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
@@ -34,30 +39,10 @@
if isUmfpack:
umfpack.configure( **kwargs )
-def _toCS_superLU( A ):
- if hasattr(A, 'tocsc') and not isspmatrix_csr( A ):
- mat = A.tocsc()
- csc = 1
- elif hasattr(A, 'tocsr'):
- mat = A.tocsr()
- csc = 0
- else:
- raise ValueError, "matrix cannot be converted to CSC/CSR"
- return mat, csc
-def _toCS_umfpack( A ):
- if isspmatrix_csr( A ) or isspmatrix_csc( A ):
- mat = A
- else:
- if hasattr(A, 'tocsc'):
- mat = A.tocsc()
- elif hasattr(A, 'tocsr'):
- mat = A.tocsr()
- else:
- raise ValueError, "matrix cannot be converted to CSC/CSR"
- return mat
-
def spsolve(A, b, permc_spec=2):
+ """Solve the sparse linear system Ax=b
+ """
if isspmatrix( b ):
b = b.toarray()
@@ -67,47 +52,44 @@
else:
raise ValueError, "rhs must be a vector (has shape %s)" % (b.shape,)
- if not isspmatrix(A):
- raise TypeError,'expected sparse matrix'
+ if not (isspmatrix_csc(A) or isspmatrix_csr(A)):
+ A = csc_matrix(A)
+ warn('spsolve requires CSC or CSR matrix format', SparseEfficiencyWarning)
+ A.sort_indices()
A = A.asfptype() #upcast to a floating point format
- if not hasattr(A, 'tocsr') and not hasattr(A, 'tocsc'):
- raise ValueError, "sparse matrix must be able to return CSC format--"\
- "A.tocsc()--or CSR format--A.tocsr()"
- if not hasattr(A, 'shape'):
- raise ValueError, "sparse matrix must be able to return shape" \
- " (rows, cols) = A.shape"
M, N = A.shape
if (M != N):
- raise ValueError, "matrix must be square (has shape %s)" % (A.shape,)
+ raise ValueError, "matrix must be square (has shape %s)" % (M,N)
if M != b.size:
raise ValueError, "matrix - rhs size mismatch (%s - %s)"\
- % (A.shape, b.shape)
+ % (A.shape, b.size)
-
if isUmfpack and useUmfpack:
- mat = _toCS_umfpack( A )
-
- if mat.dtype.char not in 'dD':
+ if A.dtype.char not in 'dD':
raise ValueError, "convert matrix data to double, please, using"\
" .astype(), or set linsolve.useUmfpack = False"
family = {'d' : 'di', 'D' : 'zi'}
- umf = umfpack.UmfpackContext( family[mat.dtype.char] )
- return umf.linsolve( umfpack.UMFPACK_A, mat, b,
+ umf = umfpack.UmfpackContext( family[A.dtype.char] )
+ return umf.linsolve( umfpack.UMFPACK_A, A, b,
autoTranspose = True )
else:
- mat, csc = _toCS_superLU( A )
- ftype = superLU_transtabl[mat.dtype.char]
- index0 = mat.indices
- lastel, data, index1 = mat.nnz, mat.data, mat.indptr
+ if isspmatrix_csc(A):
+ flag = 1 # CSC format
+ else:
+ flag = 0 # CSR format
+
+ ftype = superLU_transtabl[A.dtype.char]
+
gssv = eval('_superlu.' + ftype + 'gssv')
- b = asarray(b, dtype=data.dtype)
- return gssv(N, lastel, data, index0, index1, b, csc, permc_spec)[0]
+ b = asarray(b, dtype=A.dtype)
+ return gssv(N, A.nnz, A.data, A.indices, A.indptr, b, flag, permc_spec)[0]
+
def splu(A, permc_spec=2, diag_pivot_thresh=1.0,
drop_tol=0.0, relax=1, panel_size=10):
"""
@@ -118,19 +100,22 @@
See scipy.linsolve._superlu.dgstrf for more info.
"""
+
+ if not isspmatrix_csc(A):
+ A = csc_matrix(A)
+ warn('splu requires CSC matrix format', SparseEfficiencyWarning)
+
+ A.sort_indices()
+ A = A.asfptype() #upcast to a floating point format
+
M, N = A.shape
if (M != N):
- raise ValueError, "can only factor square matrices"
+ raise ValueError, "can only factor square matrices" #is this true?
-## if isUmfpack:
-## print "UMFPACK is present - try umfpack.numeric and umfpack.solve instead!"
+ ftype = superLU_transtabl[A.dtype.char]
- csc = A.tocsc()
- csc = csc.asfptype() #upcast to a floating point format
- ftype = superLU_transtabl[csc.dtype.char]
-
gstrf = eval('_superlu.' + ftype + 'gstrf')
- return gstrf(N, csc.nnz, csc.data, csc.indices, csc.indptr, permc_spec,
+ return gstrf(N, A.nnz, A.data, A.indices, A.indptr, permc_spec,
diag_pivot_thresh, drop_tol, relax, panel_size)
def factorized( A ):
@@ -143,29 +128,34 @@
x2 = solve( rhs2 ) # Uses again the LU factors.
"""
if isUmfpack and useUmfpack:
- mat = _toCS_umfpack( A )
+ if not isspmatrix_csc(A):
+ A = csc_matrix(A)
+ warn('splu requires CSC matrix format', SparseEfficiencyWarning)
- if mat.dtype.char not in 'dD':
+ A.sort_indices()
+ A = A.asfptype() #upcast to a floating point format
+
+ if A.dtype.char not in 'dD':
raise ValueError, "convert matrix data to double, please, using"\
" .astype(), or set linsolve.useUmfpack = False"
family = {'d' : 'di', 'D' : 'zi'}
- umf = umfpack.UmfpackContext( family[mat.dtype.char] )
+ umf = umfpack.UmfpackContext( family[A.dtype.char] )
# Make LU decomposition.
- umf.numeric( mat )
+ umf.numeric( A )
def solve( b ):
- return umf.solve( umfpack.UMFPACK_A, mat, b, autoTranspose = True )
+ return umf.solve( umfpack.UMFPACK_A, A, b, autoTranspose = True )
return solve
else:
return splu( A ).solve
def _testme():
- from scipy.sparse import csc_matrix
+ from scipy.sparse import csc_matrix, spdiags
from numpy import array
- from scipy.linsolve import spdiags, spsolve, use_solver
+ from scipy.linsolve import spsolve, use_solver
print "Inverting a sparse linear system:"
print "The sparse matrix (constructed from diagonals):"
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2008-01-16 12:05:14 UTC (rev 3840)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2008-01-16 12:27:45 UTC (rev 3841)
@@ -206,7 +206,7 @@
if len(self.As) == 1:
#TODO make spsolve preserve dimensions
- x[:] = spsolve(A,b).reshape(x.shape)
+ x[:] = spsolve(A.tocsc(),b).reshape(x.shape)
return
self.presmoother(A,x,b)
@@ -219,7 +219,7 @@
if lvl == len(self.As) - 2:
#use direct solver on coarsest level
#TODO reuse factors for efficiency?
- coarse_x[:] = spsolve(self.As[-1],coarse_b).reshape(coarse_x.shape)
+ coarse_x[:] = spsolve(self.As[-1].tocsc(),coarse_b).reshape(coarse_x.shape)
#coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0].reshape(coarse_x.shape)
#A_inv = asarray(scipy.linalg.pinv2(self.As[-1].todense()))
#coarse_x[:] = scipy.dot(A_inv,coarse_b)
Modified: trunk/scipy/sparse/tests/test_base.py
===================================================================
--- trunk/scipy/sparse/tests/test_base.py 2008-01-16 12:05:14 UTC (rev 3840)
+++ trunk/scipy/sparse/tests/test_base.py 2008-01-16 12:27:45 UTC (rev 3841)
@@ -13,6 +13,8 @@
python tests/test_sparse.py
"""
+import warnings
+
import numpy
from numpy import arange, zeros, array, dot, ones, matrix, asmatrix, \
asarray, vstack, ndarray, kron, transpose
@@ -27,7 +29,10 @@
from scipy.linsolve import splu
+warnings.simplefilter('ignore',SparseEfficiencyWarning)
+
+
#TODO test spmatrix( [[1,2],[3,4]] ) format
#TODO check that invalid shape in constructor raises exception
#TODO check that spmatrix( ... , copy=X ) is respected
@@ -485,8 +490,6 @@
class _TestGetSet:
def test_setelement(self):
- import warnings
- warnings.simplefilter('ignore',SparseEfficiencyWarning)
a = self.spmatrix((3,4))
a[1,2] = 4.0
a[0,1] = 3
More information about the Scipy-svn
mailing list