[Scipy-svn] r2839 - in trunk/Lib/linsolve: . umfpack/tests

[scipy-svn at scipy.org]

Author: rc
Date: 2007-03-08 11:06:49 -0600 (Thu, 08 Mar 2007)
New Revision: 2839

Modified:
   trunk/Lib/linsolve/linsolve.py
   trunk/Lib/linsolve/umfpack/tests/test_umfpack.py
Log:
added linsolve.factorized()

Modified: trunk/Lib/linsolve/linsolve.py
===================================================================
--- trunk/Lib/linsolve/linsolve.py	2007-03-08 09:44:13 UTC (rev 2838)
+++ trunk/Lib/linsolve/linsolve.py	2007-03-08 17:06:49 UTC (rev 2839)
@@ -118,6 +118,35 @@
     return gstrf(N, csc.nnz, csc.data, csc.indices, csc.indptr, permc_spec,
                  diag_pivot_thresh, drop_tol, relax, panel_size)
 
+def factorized( A ):
+    """
+    Return a fuction for solving a linear system, with A pre-factorized.
+
+    Example:
+      solve = factorized( A ) # Makes LU decomposition.
+      x1 = solve( rhs1 ) # Uses the LU factors.
+      x2 = solve( rhs2 ) # Uses again the LU factors.
+    """
+    if isUmfpack and useUmfpack:
+        mat = _toCS_umfpack( A )
+
+        if mat.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] )
+
+        # Make LU decomposition.
+        umf.numeric( mat )
+
+        def solve( b ):
+            return umf.solve( umfpack.UMFPACK_A, mat, b, autoTranspose = True )
+            
+        return solve
+    else:
+        return splu( A ).solve
+
 def _testme():
     from scipy.sparse import csc_matrix
     from numpy import array

Modified: trunk/Lib/linsolve/umfpack/tests/test_umfpack.py
===================================================================
--- trunk/Lib/linsolve/umfpack/tests/test_umfpack.py	2007-03-08 09:44:13 UTC (rev 2838)
+++ trunk/Lib/linsolve/umfpack/tests/test_umfpack.py	2007-03-08 17:06:49 UTC (rev 2839)
@@ -73,11 +73,34 @@
         #print "Error: ", a*x-b
         assert_array_almost_equal(a*x, self.b)
 
+    def check_factorized_umfpack(self):
+        """Prefactorize (with UMFPACK) matrix for solving with multiple rhs"""
+        linsolve.use_solver( useUmfpack = True )
+        a = self.a.astype('d')
+        solve = linsolve.factorized( a )
+
+        x1 = solve( self.b )
+        assert_array_almost_equal(a*x1, self.b)
+        x2 = solve( self.b2 )
+        assert_array_almost_equal(a*x2, self.b2)
+
+    def check_factorized_without_umfpack(self):
+        """Prefactorize matrix for solving with multiple rhs"""
+        linsolve.use_solver( useUmfpack = False )
+        a = self.a.astype('d')
+        solve = linsolve.factorized( a )
+
+        x1 = solve( self.b )
+        assert_array_almost_equal(a*x1, self.b)
+        x2 = solve( self.b2 )
+        assert_array_almost_equal(a*x2, self.b2)
+
     def setUp(self):
         self.a = spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]], [0, 1], 5, 5)
         #print "The sparse matrix (constructed from diagonals):"
         #print self.a
         self.b = array([1, 2, 3, 4, 5])
+        self.b2 = array([5, 4, 3, 2, 1])