"""
Interface to the UMFPACK library.
--
Author: Robert Cimrman
"""
#from base import Struct, pause
import numpy as np
import scipy.sparse as sp
import re
try: # Silence import error.
import _umfpack as _um
except:
_um = None
assumeSortedIndices = False
##
# 10.01.2006, c
def configure( **kwargs ):
"""
Valid keyword arguments with defaults (other ignored):
assumeSortedIndices = False
Umfpack requires a CSR/CSC matrix to have sorted column/row indices. If
sure that the matrix fulfills this, pass assumeSortedIndices =
True to gain some speed.
"""
if 'assumeSortedIndices' in kwargs:
globals()['assumeSortedIndices'] = kwargs['assumeSortedIndices']
##
# 30.11.2005, c
def updateDictWithVars( adict, module, pattern, group = None ):
match = re.compile( pattern ).match
for name in [ii for ii in vars( module ).keys()
if match( ii )]:
if group is not None:
outName = match( name ).group( group )
else:
outName = name
adict[outName] = module.__dict__[name]
return adict
##
# How to list these automagically?
umfControls = [
'UMFPACK_PRL',
'UMFPACK_DENSE_ROW',
'UMFPACK_DENSE_COL',
'UMFPACK_BLOCK_SIZE',
'UMFPACK_STRATEGY',
'UMFPACK_2BY2_TOLERANCE',
'UMFPACK_FIXQ',
'UMFPACK_AMD_DENSE',
'UMFPACK_AGGRESSIVE',
'UMFPACK_PIVOT_TOLERANCE',
'UMFPACK_ALLOC_INIT',
'UMFPACK_SYM_PIVOT_TOLERANCE',
'UMFPACK_SCALE',
'UMFPACK_FRONT_ALLOC_INIT',
'UMFPACK_DROPTOL',
'UMFPACK_IRSTEP',
'UMFPACK_COMPILED_WITH_BLAS',
'UMFPACK_COMPILED_FOR_MATLAB',
'UMFPACK_COMPILED_WITH_GETRUSAGE',
'UMFPACK_COMPILED_IN_DEBUG_MODE',
'UMFPACK_STRATEGY_AUTO',
'UMFPACK_STRATEGY_UNSYMMETRIC',
'UMFPACK_STRATEGY_2BY2',
'UMFPACK_STRATEGY_SYMMETRIC',
'UMFPACK_SCALE_NONE',
'UMFPACK_SCALE_SUM',
'UMFPACK_SCALE_MAX',
]
umfInfo = [
'UMFPACK_STATUS',
'UMFPACK_NROW',
'UMFPACK_NCOL',
'UMFPACK_NZ',
'UMFPACK_SIZE_OF_UNIT',
'UMFPACK_SIZE_OF_INT',
'UMFPACK_SIZE_OF_LONG',
'UMFPACK_SIZE_OF_POINTER',
'UMFPACK_SIZE_OF_ENTRY',
'UMFPACK_NDENSE_ROW',
'UMFPACK_NEMPTY_ROW',
'UMFPACK_NDENSE_COL',
'UMFPACK_NEMPTY_COL',
'UMFPACK_SYMBOLIC_DEFRAG',
'UMFPACK_SYMBOLIC_PEAK_MEMORY',
'UMFPACK_SYMBOLIC_SIZE',
'UMFPACK_SYMBOLIC_TIME',
'UMFPACK_SYMBOLIC_WALLTIME',
'UMFPACK_STRATEGY_USED',
'UMFPACK_ORDERING_USED',
'UMFPACK_QFIXED',
'UMFPACK_DIAG_PREFERRED',
'UMFPACK_PATTERN_SYMMETRY',
'UMFPACK_NZ_A_PLUS_AT',
'UMFPACK_NZDIAG',
'UMFPACK_SYMMETRIC_LUNZ',
'UMFPACK_SYMMETRIC_FLOPS',
'UMFPACK_SYMMETRIC_NDENSE',
'UMFPACK_SYMMETRIC_DMAX',
'UMFPACK_2BY2_NWEAK',
'UMFPACK_2BY2_UNMATCHED',
'UMFPACK_2BY2_PATTERN_SYMMETRY',
'UMFPACK_2BY2_NZ_PA_PLUS_PAT',
'UMFPACK_2BY2_NZDIAG',
'UMFPACK_COL_SINGLETONS',
'UMFPACK_ROW_SINGLETONS',
'UMFPACK_N2',
'UMFPACK_S_SYMMETRIC',
'UMFPACK_NUMERIC_SIZE_ESTIMATE',
'UMFPACK_PEAK_MEMORY_ESTIMATE',
'UMFPACK_FLOPS_ESTIMATE',
'UMFPACK_LNZ_ESTIMATE',
'UMFPACK_UNZ_ESTIMATE',
'UMFPACK_VARIABLE_INIT_ESTIMATE',
'UMFPACK_VARIABLE_PEAK_ESTIMATE',
'UMFPACK_VARIABLE_FINAL_ESTIMATE',
'UMFPACK_MAX_FRONT_SIZE_ESTIMATE',
'UMFPACK_MAX_FRONT_NROWS_ESTIMATE',
'UMFPACK_MAX_FRONT_NCOLS_ESTIMATE',
'UMFPACK_NUMERIC_SIZE',
'UMFPACK_PEAK_MEMORY',
'UMFPACK_FLOPS',
'UMFPACK_LNZ',
'UMFPACK_UNZ',
'UMFPACK_VARIABLE_INIT',
'UMFPACK_VARIABLE_PEAK',
'UMFPACK_VARIABLE_FINAL',
'UMFPACK_MAX_FRONT_SIZE',
'UMFPACK_MAX_FRONT_NROWS',
'UMFPACK_MAX_FRONT_NCOLS',
'UMFPACK_NUMERIC_DEFRAG',
'UMFPACK_NUMERIC_REALLOC',
'UMFPACK_NUMERIC_COSTLY_REALLOC',
'UMFPACK_COMPRESSED_PATTERN',
'UMFPACK_LU_ENTRIES',
'UMFPACK_NUMERIC_TIME',
'UMFPACK_UDIAG_NZ',
'UMFPACK_RCOND',
'UMFPACK_WAS_SCALED',
'UMFPACK_RSMIN',
'UMFPACK_RSMAX',
'UMFPACK_UMIN',
'UMFPACK_UMAX',
'UMFPACK_ALLOC_INIT_USED',
'UMFPACK_FORCED_UPDATES',
'UMFPACK_NUMERIC_WALLTIME',
'UMFPACK_NOFF_DIAG',
'UMFPACK_ALL_LNZ',
'UMFPACK_ALL_UNZ',
'UMFPACK_NZDROPPED',
'UMFPACK_IR_TAKEN',
'UMFPACK_IR_ATTEMPTED',
'UMFPACK_OMEGA1',
'UMFPACK_OMEGA2',
'UMFPACK_SOLVE_FLOPS',
'UMFPACK_SOLVE_TIME',
'UMFPACK_SOLVE_WALLTIME',
'UMFPACK_ORDERING_COLAMD',
'UMFPACK_ORDERING_AMD',
'UMFPACK_ORDERING_GIVEN',
]
if _um:
##
# Export UMFPACK constants from _um.
umfDefines = updateDictWithVars( {}, _um, 'UMFPACK_.*' )
locals().update( umfDefines )
umfStatus = {
UMFPACK_OK : 'UMFPACK_OK',
UMFPACK_WARNING_singular_matrix : 'UMFPACK_WARNING_singular_matrix',
UMFPACK_WARNING_determinant_underflow : 'UMFPACK_WARNING_determinant_underflow',
UMFPACK_WARNING_determinant_overflow : 'UMFPACK_WARNING_determinant_overflow',
UMFPACK_ERROR_out_of_memory : 'UMFPACK_ERROR_out_of_memory',
UMFPACK_ERROR_invalid_Numeric_object : 'UMFPACK_ERROR_invalid_Numeric_object',
UMFPACK_ERROR_invalid_Symbolic_object : 'UMFPACK_ERROR_invalid_Symbolic_object',
UMFPACK_ERROR_argument_missing : 'UMFPACK_ERROR_argument_missing',
UMFPACK_ERROR_n_nonpositive : 'UMFPACK_ERROR_n_nonpositive',
UMFPACK_ERROR_invalid_matrix : 'UMFPACK_ERROR_invalid_matrix',
UMFPACK_ERROR_different_pattern : 'UMFPACK_ERROR_different_pattern',
UMFPACK_ERROR_invalid_system : 'UMFPACK_ERROR_invalid_system',
UMFPACK_ERROR_invalid_permutation : 'UMFPACK_ERROR_invalid_permutation',
UMFPACK_ERROR_internal_error : 'UMFPACK_ERROR_internal_error',
UMFPACK_ERROR_file_IO : 'UMFPACK_ERROR_file_IO',
}
umfSys = [
UMFPACK_A,
UMFPACK_At,
UMFPACK_Aat,
UMFPACK_Pt_L,
UMFPACK_L,
UMFPACK_Lt_P,
UMFPACK_Lat_P,
UMFPACK_Lt,
UMFPACK_U_Qt,
UMFPACK_U,
UMFPACK_Q_Ut,
UMFPACK_Q_Uat,
UMFPACK_Ut,
UMFPACK_Uat,
]
# Real, complex.
umfSys_transposeMap = [
{UMFPACK_A : UMFPACK_At,
UMFPACK_At : UMFPACK_A,
UMFPACK_Aat : UMFPACK_A},
{UMFPACK_A : UMFPACK_Aat,
UMFPACK_Aat : UMFPACK_A}
]
umfFamilyTypes = {'di' : int, 'dl' : long, 'zi' : int, 'zl' : long}
umfRealTypes = ('di', 'dl')
umfComplexTypes = ('zi', 'zl')
##
# 02.01.2005
class Struct( object ):
# 03.10.2005, c
# 26.10.2005
def __init__( self, **kwargs ):
if kwargs:
self.__dict__.update( kwargs )
# 08.03.2005
def __str__( self ):
ss = "%s\n" % self.__class__
for key, val in self.__dict__.iteritems():
if (issubclass( self.__dict__[key].__class__, Struct )):
ss += " %s:\n %s\n" % (key, self.__dict__[key].__class__)
else:
aux = "\n" + str( val )
aux = aux.replace( "\n", "\n " );
ss += " %s:\n%s\n" % (key, aux[1:])
return( ss.rstrip() )
##
# 30.11.2005, c
class UmfpackContext( Struct ):
##
# 30.11.2005, c
# 01.12.2005
# 21.12.2005
# 01.03.2006
def __init__( self, family = 'di', **kwargs ):
"""
Arguments:
family .. family of UMFPACK functions ('di', 'dl', 'zi', 'zl')
Keyword arguments:
maxCond .. if extimated condition number is greater than maxCond,
a warning is printed (default: 1e12)"""
if _um is None:
raise ImportError('Scipy was built without UMFPACK support. '
'You need to install the UMFPACK library and '
'header files before building scipy.')
self.maxCond = 1e12
Struct.__init__( self, **kwargs )
if family not in umfFamilyTypes.keys():
raise TypeError, 'wrong family: %s' % family
self.family = family
self.control = np.zeros( (UMFPACK_CONTROL, ), dtype = np.double )
self.info = np.zeros( (UMFPACK_INFO, ), dtype = np.double )
self._symbolic = None
self._numeric = None
self.mtx = None
self.isReal = self.family in umfRealTypes
##
# Functions corresponding to <family> are stored in self.funs.
pattern = 'umfpack_' + family + '_(.*)'
fn = updateDictWithVars( {}, _um, pattern, group = 1 )
self.funs = Struct( **fn )
self.funs.defaults( self.control )
self.control[UMFPACK_PRL] = 3
##
# 30.11.2005, c
def strControl( self ):
maxLen = max( [len( name ) for name in umfControls] )
format = '%%-%ds : %%d' % maxLen
aux = [format % (name, self.control[umfDefines[name]])
for name in umfControls if name in umfDefines]
return '\n'.join( aux )
##
# 01.12.2005, c
def strInfo( self ):
maxLen = max( [len( name ) for name in umfInfo] )
format = '%%-%ds : %%d' % maxLen
aux = [format % (name, self.info[umfDefines[name]])
for name in umfInfo if name in umfDefines]
return '\n'.join( aux )
##
# 30.11.2005, c
# 01.12.2005
# 14.12.2005
# 01.03.2006
def _getIndx( self, mtx ):
if sp.isspmatrix_csc( mtx ):
indx = mtx.indices
self.isCSR = 0
elif sp.isspmatrix_csr( mtx ):
indx = mtx.indices
self.isCSR = 1
else:
raise TypeError, 'must be a CSC/CSR matrix (is %s)' % mtx.__class__
##
# Should check types of indices to correspond to familyTypes.
if self.family[1] == 'i':
if (indx.dtype != np.dtype('i')) \
or mtx.indptr.dtype != np.dtype('i'):
raise ValueError, 'matrix must have int indices'
else:
if (indx.dtype != np.dtype('l')) \
or mtx.indptr.dtype != np.dtype('l'):
raise ValueError, 'matrix must have long indices'
if self.isReal:
if mtx.data.dtype != np.dtype('<f8'):
raise ValueError, 'matrix must have float64 values'
else:
if mtx.data.dtype != np.dtype('<c16'):
raise ValueError, 'matrix must have complex128 values'
return indx
##
# 30.11.2005, c
# last revision: 10.01.2007
def symbolic( self, mtx ):
"""Symbolic object (symbolic LU decomposition) computation for a given
sparsity pattern."""
self.free_symbolic()
indx = self._getIndx( mtx )
if not assumeSortedIndices:
# row/column indices cannot be assumed to be sorted
mtx.sort_indices()
if self.isReal:
status, self._symbolic\
= self.funs.symbolic( mtx.shape[0], mtx.shape[1],
mtx.indptr, indx, mtx.data,
self.control, self.info )
else:
real, imag = mtx.data.real.copy(), mtx.data.imag.copy()
status, self._symbolic\
= self.funs.symbolic( mtx.shape[0], mtx.shape[1],
mtx.indptr, indx,
real, imag,
self.control, self.info )
## print status, self._symbolic
if status != UMFPACK_OK:
raise RuntimeError, '%s failed with %s' % (self.funs.symbolic,
umfStatus[status])
self.mtx = mtx
##
# 30.11.2005, c
# 01.12.2005
# 02.12.2005
# 01.03.2006
def numeric( self, mtx ):
"""Numeric object (LU decomposition) computation using the
symbolic decomposition. The symbolic decomposition is (re)computed
if necessary."""
self.free_numeric()
if self._symbolic is None:
self.symbolic( mtx )
indx = self._getIndx( mtx )
failCount = 0
while 1:
if self.isReal:
status, self._numeric\
= self.funs.numeric( mtx.indptr, indx, mtx.data,
self._symbolic,
self.control, self.info )
else:
real, imag = mtx.data.real.copy(), mtx.data.imag.copy()
status, self._numeric\
= self.funs.numeric( mtx.indptr, indx,
real, imag,
self._symbolic,
self.control, self.info )
## print status, self._numeric
if status != UMFPACK_OK:
if status == UMFPACK_WARNING_singular_matrix:
print 'warning: singular matrix'
break
elif status in (UMFPACK_ERROR_different_pattern,
UMFPACK_ERROR_invalid_Symbolic_object):
# Try again.
print 'warning: recomputing symbolic'
self.symbolic( mtx )
failCount += 1
else:
failCount += 100
else:
break
if failCount >= 2:
raise RuntimeError, '%s failed with %s' % (self.funs.numeric,
umfStatus[status])
##
# 14.12.2005, c
def report_symbolic( self ):
"""Print information about the symbolic object. Output depends on
self.control[UMFPACK_PRL]."""
self.funs.report_symbolic( self._symbolic, self.control )
##
# 14.12.2005, c
def report_numeric( self ):
"""Print information about the numeric object. Output depends on
self.control[UMFPACK_PRL]."""
self.funs.report_numeric( self._numeric, self.control )
##
# 14.12.2005, c
def report_control( self ):
"""Print control values."""
self.funs.report_control( self.control )
##
# 14.12.2005, c
def report_info( self ):
"""Print all status information. Output depends on
self.control[UMFPACK_PRL]."""
self.funs.report_info( self.control, self.info )
##
# 30.11.2005, c
# 01.12.2005
def free_symbolic( self ):
if self._symbolic is not None:
self.funs.free_symbolic( self._symbolic )
self._symbolic = None
self.mtx = None
##
# 30.11.2005, c
# 01.12.2005
def free_numeric( self ):
if self._numeric is not None:
self.funs.free_numeric( self._numeric )
self._numeric = None
self.free_symbolic()
##
# 30.11.2005, c
def free( self ):
self.free_symbolic()
self.free_numeric()
##
# 30.11.2005, c
# 01.12.2005
# 02.12.2005
# 21.12.2005
# 01.03.2006
def solve( self, sys, mtx, rhs, autoTranspose = False ):
"""
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
if autoTranspose and self.isCSR:
##
# UMFPACK uses CSC internally...
if self.family in umfRealTypes: ii = 0
else: ii = 1
if sys in umfSys_transposeMap[ii]:
sys = umfSys_transposeMap[ii][sys]
else:
raise RuntimeError, 'autoTranspose ambiguous, switch it off'
if self._numeric is not None:
if self.mtx is not mtx:
raise ValueError, 'must be called with same matrix as numeric()'
else:
raise RuntimeError, 'numeric() not called'
indx = self._getIndx( mtx )
if self.isReal:
rhs = rhs.astype( np.float64 )
sol = np.zeros( (mtx.shape[1],), dtype = np.float64 )
status = self.funs.solve( sys, mtx.indptr, indx, mtx.data, sol, rhs,
self._numeric, self.control, self.info )
else:
rhs = rhs.astype( np.complex128 )
sol = np.zeros( (mtx.shape[1],), dtype = np.complex128 )
mreal, mimag = mtx.data.real.copy(), mtx.data.imag.copy()
sreal, simag = sol.real.copy(), sol.imag.copy()
rreal, rimag = rhs.real.copy(), rhs.imag.copy()
status = self.funs.solve( sys, mtx.indptr, indx,
mreal, mimag, sreal, simag, rreal, rimag,
self._numeric, self.control, self.info )
sol.real, sol.imag = sreal, simag
#self.funs.report_info( self.control, self.info )
#pause()
if status != UMFPACK_OK:
if status == UMFPACK_WARNING_singular_matrix:
## Change inf, nan to zeros.
print 'zeroing nan and inf entries...'
sol[~np.isfinite( sol )] = 0.0
else:
raise RuntimeError, '%s failed with %s' % (self.funs.solve,
umfStatus[status])
econd = 1.0 / self.info[UMFPACK_RCOND]
if econd > self.maxCond:
print 'warning: (almost) singular matrix! '\
+ '(estimated cond. number: %.2e)' % econd
return sol
##
# 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.
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
if self._numeric is None:
self.numeric( mtx )
else:
if self.mtx is not mtx:
self.numeric( mtx )
sol = self.solve( sys, mtx, rhs, autoTranspose )
self.free_numeric()
return sol
##
# 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.
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 = np.zeros( (n_row+1,), dtype = i_type )
Lj = np.zeros( (lnz,), dtype = i_type )
Lx = np.zeros( (lnz,), dtype = np.double )
Up = np.zeros( (n_col+1,), dtype = i_type )
Ui = np.zeros( (unz,), dtype = i_type )
Ux = np.zeros( (unz,), dtype = np.double )
P = np.zeros( (n_row,), dtype = i_type )
Q = np.zeros( (n_col,), dtype = i_type )
Dx = np.zeros( (min(n_row,n_col),), dtype = np.double )
Rs = np.zeros( (n_row,), dtype = np.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,bool(do_recip))
else:
#allocate additional storage for imaginary parts
Lz = np.zeros( (lnz,), dtype = np.double )
Uz = np.zeros( (unz,), dtype = np.double )
Dz = np.zeros( (min(n_row,n_col),), dtype = np.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 = np.zeros( (lnz,), dtype = np.complex128 )
Uxz = np.zeros( (unz,), dtype = np.complex128 )
Dxz = np.zeros( (min(n_row,n_col),), dtype = np.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,bool(do_recip))
|