""" Interface to the UMFPACK library. -- Author: Robert Cimrman """ #from base import Struct, pause import numpy as nm import scipy.sparse as sp import re, imp 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 kwargs.has_key( 'assumeSortedIndices' ): 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)""" self.maxCond = 1e12 Struct.__init__( self, **kwargs ) if family not in umfFamilyTypes.keys(): raise TypeError, 'wrong family: %s' % family self.family = family self.control = nm.zeros( (UMFPACK_CONTROL, ), dtype = nm.double ) self.info = nm.zeros( (UMFPACK_INFO, ), dtype = nm.double ) self._symbolic = None self._numeric = None self.mtx = None self.isReal = self.family in umfRealTypes ## # Functions corresponding to 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 umfDefines.has_key( name )] 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 umfDefines.has_key( name )] 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 != nm.dtype('i')) \ or mtx.indptr.dtype != nm.dtype('i'): raise ValueError, 'matrix must have int indices' else: if (indx.dtype != nm.dtype('l')) \ or mtx.indptr.dtype != nm.dtype('l'): raise ValueError, 'matrix must have long indices' if self.isReal: if mtx.data.dtype != nm.dtype('= 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 umfSys_transposeMap[ii].has_key( sys ): 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( nm.float64 ) sol = nm.zeros( (mtx.shape[1],), dtype = nm.float64 ) status = self.funs.solve( sys, mtx.indptr, indx, mtx.data, sol, rhs, self._numeric, self.control, self.info ) else: rhs = rhs.astype( nm.complex128 ) sol = nm.zeros( (mtx.shape[1],), dtype = nm.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[~nm.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 = 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,bool(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,bool(do_recip))