""" 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 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 def __del__(self): self.free() ## # 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))