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

scipy-svn at scipy.org scipy-svn at scipy.org
Thu Sep 21 04:27:30 CDT 2006

```Author: rc
Date: 2006-09-21 04:26:58 -0500 (Thu, 21 Sep 2006)
New Revision: 2208

Modified:
trunk/Lib/linsolve/umfpack/info.py
trunk/Lib/linsolve/umfpack/tests/test_umfpack.py
trunk/Lib/linsolve/umfpack/umfpack.i
trunk/Lib/linsolve/umfpack/umfpack.py
Log:
merged lu() method added by Nathan Bell

Modified: trunk/Lib/linsolve/umfpack/info.py
===================================================================
--- trunk/Lib/linsolve/umfpack/info.py	2006-09-21 05:19:03 UTC (rev 2207)
+++ trunk/Lib/linsolve/umfpack/info.py	2006-09-21 09:26:58 UTC (rev 2208)
@@ -42,6 +42,7 @@
include_dirs = <dir>/UFsparse/UMFPACK/Include, <dir>/UFsparse/UFconfig
umfpack_libs = umfpack

+
Examples:
=========

@@ -97,6 +98,36 @@
# Print all statistics.
umfpack.report_info()

+-or-
+
+# Get LU factors and permutation matrices of a matrix.
+L, U, P, Q, R, do_recip = umfpack.lu( mtx )
+
+Then:
+           L - Lower triangular m-by-min(m,n) CSR matrix
+           U - Upper triangular min(m,n)-by-n CSC matrix
+           P - Vector of row permuations
+           Q - Vector of column permuations
+           R - Vector of diagonal row scalings
+           do_recip - boolean
+
+       For a given matrix A, the decomposition satisfies:
+               LU = PRAQ        when do_recip is true
+               LU = P(R^-1)AQ   when do_recip is false
+
+Description of arguments of UmfpackContext solution methods:
+=============================================
+This holds for: umfpack(), umfpack.linsolve(), umfpack.solve()
+
+ sys - one of UMFPACK system description constants, like
+       UMFPACK_A, UMFPACK_At, see umfSys list and UMFPACK
+       docs
+ mtx - sparse matrix (CSR or CSC)
+ rhs - right hand side vector
+ autoTranspose - automatically changes 'sys' to the
+       transposed type, if 'mtx' is in CSR, since UMFPACK
+       assumes CSC internally
+
Setting control parameters:
===========================
Assuming this module imported as um:
@@ -113,6 +144,8 @@

--
Author: Robert Cimrman
+
+Other contributors: Nathan Bell (lu() method wrappers)
"""

postpone_import = 1

Modified: trunk/Lib/linsolve/umfpack/tests/test_umfpack.py
===================================================================
--- trunk/Lib/linsolve/umfpack/tests/test_umfpack.py	2006-09-21 05:19:03 UTC (rev 2207)
+++ trunk/Lib/linsolve/umfpack/tests/test_umfpack.py	2006-09-21 09:26:58 UTC (rev 2208)
@@ -10,8 +10,12 @@
import random
from numpy.testing import *
set_package_path()
-from scipy import linsolve
+from scipy import linsolve, rand, matrix, diag, eye
from scipy.sparse import csc_matrix, dok_matrix, spdiags
+
+import numpy as nm
+import scipy.linsolve.umfpack as um
+
restore_path()

class test_solvers(ScipyTestCase):
@@ -67,5 +71,65 @@
self.b = array([1, 2, 3, 4, 5])

+
+class test_factorization(ScipyTestCase):
+    """Tests factorizing a sparse linear system"""
+
+    def check_complex_lu(self):
+        """Getting factors of complex matrix"""
+        umfpack = um.UmfpackContext("zi")
+
+        for A in self.complex_matrices:
+            umfpack.numeric(A)
+
+            (L,U,P,Q,R,do_recip) = umfpack.lu(A)
+
+            L = L.todense()
+            U = U.todense()
+            A = A.todense()
+            if not do_recip: R = 1.0/R
+            R = matrix(diag(R))
+            P = eye(A.shape[0])[P,:]
+            Q = eye(A.shape[1])[:,Q]
+
+            assert_array_almost_equal(P*R*A*Q,L*U)
+
+    def check_real_lu(self):
+        """Getting factors of real matrix"""
+        umfpack = um.UmfpackContext("di")
+
+        for A in self.real_matrices:
+            umfpack.numeric(A)
+
+            (L,U,P,Q,R,do_recip) = umfpack.lu(A)
+
+            L = L.todense()
+            U = U.todense()
+            A = A.todense()
+            if not do_recip: R = 1.0/R
+            R = matrix(diag(R))
+            P = eye(A.shape[0])[P,:]
+            Q = eye(A.shape[1])[:,Q]
+
+            assert_array_almost_equal(P*R*A*Q,L*U)
+
+
+    def setUp(self):
+        random.seed(0) #make tests repeatable
+        self.real_matrices = []
+        self.real_matrices.append(spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]],
+                                          [0, 1], 5, 5))
+        self.real_matrices.append(spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]],
+                                          [0, 1], 4, 5))
+        self.real_matrices.append(spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]],
+                                          [0, 2], 5, 5))
+        self.real_matrices.append(csc_matrix(rand(3,3)))
+        self.real_matrices.append(csc_matrix(rand(5,4)))
+        self.real_matrices.append(csc_matrix(rand(4,5)))
+
+        self.complex_matrices = [x.astype(nm.complex128)
+                                 for x in self.real_matrices]
+
+
if __name__ == "__main__":
ScipyTest().run()

Modified: trunk/Lib/linsolve/umfpack/umfpack.i
===================================================================
--- trunk/Lib/linsolve/umfpack/umfpack.i	2006-09-21 05:19:03 UTC (rev 2207)
+++ trunk/Lib/linsolve/umfpack/umfpack.i	2006-09-21 09:26:58 UTC (rev 2208)
@@ -220,3 +220,49 @@

%include <umfpack_free_symbolic.h>
%include <umfpack_free_numeric.h>
+
+
+
+/*
+ * wnbell - attempt to get L,U,P,Q out
+ */
+%include "typemaps.i"
+%apply int  *OUTPUT {
+    int *lnz,
+    int *unz,
+    int *n_row,
+    int *n_col,
+    int *nz_udiag
+};
+%apply long *OUTPUT {
+    long *lnz,
+    long *unz,
+    long *n_row,
+    long *n_col,
+    long *nz_udiag
+};
+%include <umfpack_get_lunz.h>
+
+
+ARRAY_IN( double, double, DOUBLE )
+%apply double *array {
+    double Lx [ ],
+    double Lz [ ],
+    double Ux [ ],
+    double Uz [ ],
+    double Dx [ ],
+    double Dz [ ],
+    double Rs [ ]
+};
+
+ARRAY_IN( int, int, INT )
+%apply int *array {
+    int Lp [ ],
+    int Lj [ ],
+    int Up [ ],
+    int Ui [ ],
+    int P [ ],
+    int Q [ ]
+};
+%apply int  *OUTPUT { int *do_recip};
+%include <umfpack_get_numeric.h>

Modified: trunk/Lib/linsolve/umfpack/umfpack.py
===================================================================
--- trunk/Lib/linsolve/umfpack/umfpack.py	2006-09-21 05:19:03 UTC (rev 2207)
+++ trunk/Lib/linsolve/umfpack/umfpack.py	2006-09-21 09:26:58 UTC (rev 2208)
@@ -467,7 +467,19 @@
# 21.12.2005
# 01.03.2006
def solve( self, sys, mtx, rhs, autoTranspose = False ):
-        """Solution of system of linear equation using the Numeric object."""
+        """
+        Solution of system of linear equation using the Numeric object.
+
+        Arguments:
+                sys - one of UMFPACK system description constants, like
+                      UMFPACK_A, UMFPACK_At, see umfSys list and UMFPACK
+                      docs
+                mtx - sparse matrix (CSR or CSC)
+                rhs - right hand side vector
+                autoTranspose - automatically changes 'sys' to the
+                      transposed type, if 'mtx' is in CSR, since UMFPACK
+                      assumes CSC internally
+        """
if sys not in umfSys:
raise ValueError, 'sys must be in' % umfSys

@@ -526,9 +538,21 @@
# 30.11.2005, c
# 01.12.2005
def linsolve( self, sys, mtx, rhs, autoTranspose = False ):
-        """One-shot solution of system of linear equation. Reuses Numeric
-        object if possible."""
+        """
+        One-shot solution of system of linear equation. Reuses Numeric object
+        if possible.

+        Arguments:
+                sys - one of UMFPACK system description constants, like
+                      UMFPACK_A, UMFPACK_At, see umfSys list and UMFPACK
+                      docs
+                mtx - sparse matrix (CSR or CSC)
+                rhs - right hand side vector
+                autoTranspose - automatically changes 'sys' to the
+                      transposed type, if 'mtx' is in CSR, since UMFPACK
+                      assumes CSC internally
+        """
+
#        print self.family
if sys not in umfSys:
raise ValueError, 'sys must be in' % umfSys
@@ -548,10 +572,116 @@
# 30.11.2005, c
# 01.12.2005
def __call__( self, sys, mtx, rhs, autoTranspose = False ):
-        """Uses solve() or linsolve() depending on the presence of
-        the Numeric object."""
+        """
+        Uses solve() or linsolve() depending on the presence of the Numeric
+        object.

+        Arguments:
+                sys - one of UMFPACK system description constants, like
+                      UMFPACK_A, UMFPACK_At, see umfSys list and UMFPACK
+                      docs
+                mtx - sparse matrix (CSR or CSC)
+                rhs - right hand side vector
+                autoTranspose - automatically changes 'sys' to the
+                      transposed type, if 'mtx' is in CSR, since UMFPACK
+                      assumes CSC internally
+        """
+
if self._numeric is not None:
return self.solve( sys, mtx, rhs, autoTranspose )
else:
return self.linsolve( sys, mtx, rhs, autoTranspose )
+
+    ##
+    # 21.09.2006, added by Nathan Bell
+    def lu( self, mtx ):
+        """
+        Returns an LU decomposition of an m-by-n matrix in the form
+        (L, U, P, Q, R, do_recip):
+
+            L - Lower triangular m-by-min(m,n) CSR matrix
+            U - Upper triangular min(m,n)-by-n CSC matrix
+            P - Vector of row permuations
+            Q - Vector of column permuations
+            R - Vector of diagonal row scalings
+            do_recip - boolean
+
+        For a given matrix A, the decomposition satisfies:
+                LU = PRAQ        when do_recip is true
+                LU = P(R^-1)AQ   when do_recip is false
+        """
+
+        #this should probably be changed
+        mtx = mtx.tocsc()
+        self.numeric( mtx )
+
+        #first find out how much space to reserve
+        (status, lnz, unz, n_row, n_col, nz_udiag)\
+                 = self.funs.get_lunz( self._numeric )
+
+        if status != UMFPACK_OK:
+            raise RuntimeError, '%s failed with %s' % (self.funs.get_lunz,
+                                                       umfStatus[status])
+
+        #allocate storage for decomposition data
+        i_type = mtx.indptr.dtype
+
+        Lp = nm.zeros( (n_row+1,), dtype = i_type )
+        Lj = nm.zeros( (lnz,), dtype = i_type )
+        Lx = nm.zeros( (lnz,), dtype = nm.double )
+
+        Up = nm.zeros( (n_col+1,), dtype = i_type )
+        Ui = nm.zeros( (unz,), dtype = i_type )
+        Ux = nm.zeros( (unz,), dtype = nm.double )
+
+        P  = nm.zeros( (n_row,), dtype = i_type )
+        Q  = nm.zeros( (n_col,), dtype = i_type )
+
+        Dx = nm.zeros( (min(n_row,n_col),), dtype = nm.double )
+
+        Rs = nm.zeros( (n_row,), dtype = nm.double )
+
+        if self.isReal:
+            (status,do_recip) = self.funs.get_numeric( Lp,Lj,Lx,Up,Ui,Ux,
+                                                       P,Q,Dx,Rs,
+                                                       self._numeric )
+
+            if status != UMFPACK_OK:
+                raise RuntimeError, '%s failed with %s'\
+                      % (self.funs.get_numeric, umfStatus[status])
+
+            L = sp.csr_matrix((Lx,Lj,Lp),(n_row,min(n_row,n_col)))
+            U = sp.csc_matrix((Ux,Ui,Up),(min(n_row,n_col),n_col))
+            R = Rs
+
+            return (L,U,P,Q,R,do_recip)
+
+        else:
+            #allocate additional storage for imaginary parts
+            Lz = nm.zeros( (lnz,), dtype = nm.double )
+            Uz = nm.zeros( (unz,), dtype = nm.double )
+            Dz = nm.zeros( (min(n_row,n_col),), dtype = nm.double )
+
+            (status,do_recip) = self.funs.get_numeric(Lp,Lj,Lx,Lz,Up,Ui,Ux,Uz,
+                                                      P,Q,Dx,Dz,Rs,
+                                                      self._numeric)
+
+            if status != UMFPACK_OK:
+                raise RuntimeError, '%s failed with %s'\
+                      % (self.funs.get_numeric, umfStatus[status])
+
+
+            Lxz = nm.zeros( (lnz,), dtype = nm.complex128 )
+            Uxz = nm.zeros( (unz,), dtype = nm.complex128 )
+            Dxz = nm.zeros( (min(n_row,n_col),), dtype = nm.complex128 )
+
+            Lxz.real,Lxz.imag = Lx,Lz
+            Uxz.real,Uxz.imag = Ux,Uz
+            Dxz.real,Dxz.imag = Dx,Dz
+
+            L = sp.csr_matrix((Lxz,Lj,Lp),(n_row,min(n_row,n_col)))
+            U = sp.csc_matrix((Uxz,Ui,Up),(min(n_row,n_col),n_col))
+            R = Rs
+
+            return (L,U,P,Q,R,do_recip)
+

```