[Scipy-svn] r2208 - in trunk/Lib/linsolve/umfpack: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Sep 21 05:27:30 EDT 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)
+
More information about the Scipy-svn
mailing list