[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