# ___________________________________________________________________________ # # Pyomo: Python Optimization Modeling Objects # Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC # Under the terms of Contract DE-NA0003525 with National Technology and # Engineering Solutions of Sandia, LLC, the U.S. Government retains certain # rights in this software. # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ from scipy.sparse import isspmatrix_coo, coo_matrix import numpy as np from pyomo.common.dependencies import attempt_import mumps, mumps_available = attempt_import( 'mumps', error_message="Error importing mumps. PyNumero's " "mumps_interface requires pymumps; install it with, e.g., " "'conda install -c conda-forge pymumps'") from pyomo.contrib.pynumero.sparse import BlockVector class MumpsCentralizedAssembledLinearSolver(object): """ A thin wrapper around pymumps which uses the centralized assembled matrix format. In other words ICNTL(5) = 0 and ICNTL(18) = 0. Solve matrix * x = rhs for x. See the Mumps documentation for descriptions of the parameters. The section numbers listed below refer to the Mumps documentation for version 5.2.1. Parameters ---------- sym: int, optional See section 5.2.1 of the Mumps documentation par: int, optional See section 5.1.3 comm: mpi4py comm, optional See section 5.1.3 cntl_options: dict, optional See section 6.2 icntl_options: dict, optional See section 6.1 """ def __init__(self, sym=0, par=1, comm=None, cntl_options=None, icntl_options=None): self._nnz = None self._dim = None self._mumps = mumps.DMumpsContext(sym=sym, par=par, comm=comm) self._mumps.set_silent() self._icntl_options = dict() self._cntl_options = dict() if cntl_options is None: cntl_options = dict() if icntl_options is None: icntl_options = dict() for k, v in cntl_options.items(): self.set_cntl(k, v) for k, v in icntl_options.items(): self.set_icntl(k, v) def _init(self): """ The purpose of this method is to address issue #12 from pymumps """ self._mumps.run(job=-1) self._mumps.set_silent() for k, v in self._cntl_options.items(): self.set_cntl(k, v) for k, v in self._icntl_options.items(): self.set_icntl(k, v) def do_symbolic_factorization(self, matrix): """ Perform Mumps analysis. Parameters ---------- matrix: scipy.sparse.spmatrix or pyomo.contrib.pynumero.sparse.BlockMatrix This matrix must have the same nonzero structure as the matrix passed into do_numeric_factorization. The matrix will be converted to coo format if it is not already in coo format. If sym is 1 or 2, the matrix must be lower or upper triangular. """ self._init() if type(matrix) == np.ndarray: matrix = coo_matrix(matrix) if not isspmatrix_coo(matrix): matrix = matrix.tocoo() nrows, ncols = matrix.shape if nrows != ncols: raise ValueError('matrix is not square') self._dim = nrows self._nnz = matrix.nnz self._mumps.set_shape(nrows) self._mumps.set_centralized_assembled_rows_cols(matrix.row + 1, matrix.col + 1) self._mumps.run(job=1) def do_numeric_factorization(self, matrix): """ Perform Mumps factorization. Note that do_symbolic_factorization should be called before do_numeric_factorization. Parameters ---------- matrix: scipy.sparse.spmatrix or pyomo.contrib.pynumero.sparse.BlockMatrix This matrix must have the same nonzero structure as the matrix passed into do_symbolic_factorization. The matrix will be converted to coo format if it is not already in coo format. If sym is 1 or 2, the matrix must be lower or upper triangular. """ if self._nnz is None: raise RuntimeError('Call do_symbolic_factorization first.') if type(matrix) == np.ndarray: matrix = coo_matrix(matrix) if not isspmatrix_coo(matrix): matrix = matrix.tocoo() nrows, ncols = matrix.shape if nrows != ncols: raise ValueError('matrix is not square') if self._dim != nrows: raise ValueError('The shape of the matrix changed between symbolic and numeric factorization') if self._nnz != matrix.nnz: raise ValueError('The number of nonzeros changed between symbolic and numeric factorization') self._mumps.set_centralized_assembled_values(matrix.data) self._mumps.run(job=2) def do_back_solve(self, rhs): """ Perform back solve with Mumps. Note that both do_symbolic_factorization and do_numeric_factorization should be called before do_back_solve. Parameters ---------- rhs: numpy.ndarray or pyomo.contrib.pynumero.sparse.BlockVector The right hand side in matrix * x = rhs. Returns ------- result: numpy.ndarray or pyomo.contrib.pynumero.sparse.BlockVector The x in matrix * x = rhs. If rhs is a BlockVector, then, result will be a BlockVector with the same block structure as rhs. """ if isinstance(rhs, BlockVector): _rhs = rhs.flatten() result = _rhs else: result = rhs.copy() self._mumps.set_rhs(result) self._mumps.run(job=3) if isinstance(rhs, BlockVector): _result = rhs.copy_structure() _result.copyfrom(result) result = _result return result def __del__(self): self._mumps.destroy() def set_icntl(self, key, value): self._icntl_options[key] = value self._mumps.set_icntl(key, value) def set_cntl(self, key, value): self._cntl_options[key] = value self._mumps.id.cntl[key - 1] = value def solve(self, matrix, rhs): self.do_symbolic_factorization(matrix) self.do_numeric_factorization(matrix) return self.do_back_solve(rhs) def get_icntl(self, key): return self._mumps.id.icntl[key - 1] def get_cntl(self, key): return self._mumps.id.cntl[key - 1] def get_info(self, key): return self._mumps.id.info[key - 1] def get_infog(self, key): return self._mumps.id.infog[key - 1] def get_rinfo(self, key): return self._mumps.id.rinfo[key - 1] def get_rinfog(self, key): return self._mumps.id.rinfog[key - 1]