From 7befdbc540ec9f5d4479c632cf4f66a34dafca3f Mon Sep 17 00:00:00 2001 From: Jonathan Guyer Date: Mon, 4 Dec 2023 17:41:14 -0500 Subject: [PATCH] Fix #963: PETSc 3.20.0 broke the world (#982) Fixes #963 (Actually, PETSC 3.19.0 broke the world.) This PR: - Assembles before `multTranspose` to prevent newly added exception - Renames `bandwidth` to `nonZerosPerRow` and removes `sizeHint` The two were confusingly redundant: - PySparse takes `sizeHint`, the number of non-zeros in the matrix. - PyTrilinos takes `NumEntriesPerRow`. - petsc4py didn't used to be clear what it took, but is now documented as number of non-zeros per row (of the local portion of the matrix, but we'll ignore that part). - scipy doesn't preallocate. - Linear algebra [defines bandwidth](https://en.wikipedia.org/wiki/Band_matrix#Bandwidth) as "number $k$ such that $a_{i,j}=0$ if $|i-j| > k$", which is roughly half the number of non-zeros per row (and only applies to a band-diagonal matrix). Better to be explicit about what we really mean. Now all take same parameter and PySparse adjusts as needed. `sizeHint` was introduced in @a15d696 (in 2006!) to "allow smaller preallocations", but it was never used that way. Now, `nonZerosPerRow` can take an array_like to specify row-by-row preallocations, which are directly supported by PyTrilinos and petsc4py, and can be simulated for PySparse. Added `exactNonZeros`, which may have performance benefits for PyTrilinos and petsc4py. Currently unused. - Uses `Term`'s knowledge of own stencil to preallocate more effectively. Still doesn't do a good job with vector equations, but that's a deeper change (the resolution of which might help #920). - Fixes(?) conda/mamba installs - Fixes(?) race condition --- .azure/pipelines.yml | 2 +- .azure/templates/install.yml | 17 +- .en-custom.txt | 3 + fipy/boundaryConditions/boundaryCondition.py | 4 +- fipy/boundaryConditions/fixedValue.py | 6 +- fipy/matrices/offsetSparseMatrix.py | 37 +++- fipy/matrices/petscMatrix.py | 182 +++++++++++------ fipy/matrices/pysparseMatrix.py | 148 +++++++++----- fipy/matrices/scipyMatrix.py | 81 ++++---- fipy/matrices/trilinosMatrix.py | 193 +++++++++++++------ fipy/meshes/gmshMesh.py | 4 + fipy/terms/abstractDiffusionTerm.py | 9 +- fipy/terms/cellTerm.py | 2 +- fipy/terms/faceTerm.py | 3 +- 14 files changed, 463 insertions(+), 228 deletions(-) diff --git a/.azure/pipelines.yml b/.azure/pipelines.yml index 0f8567d0ec..230c184b95 100644 --- a/.azure/pipelines.yml +++ b/.azure/pipelines.yml @@ -258,7 +258,7 @@ stages: condition: startsWith(variables.image, 'windows') - bash: | - mamba create --quiet --name wheelEnvironment --channel conda-forge python=3.10 mamba + conda create --quiet --name wheelEnvironment --channel conda-forge python=3.10 source activate wheelEnvironment mkdir tmp cd tmp diff --git a/.azure/templates/install.yml b/.azure/templates/install.yml index 2fda776e45..7e08479446 100644 --- a/.azure/templates/install.yml +++ b/.azure/templates/install.yml @@ -31,24 +31,11 @@ steps: - bash: | conda config --set always_yes yes --set changeps1 no conda config --remove channels defaults - if [[ ${{ parameters.python_version }} == "2.7" ]]; then - # mamba doesn't work in Py2.7 - conda create --quiet --name myEnvironment --channel conda-forge python=${{ parameters.python_version }} - else - # mamba needs to be installed in the Py3k environment for some reason - conda create --yes --quiet --name myEnvironment --channel conda-forge python=${{ parameters.python_version }} mamba - fi + conda config --set solver libmamba displayName: Create Anaconda environment - bash: | - if [[ ${{ parameters.python_version }} == "2.7" ]]; then - # mamba doesn't work in Py2.7 - conda install --quiet --name myEnvironment --channel conda-forge python=${{ parameters.python_version }} ${{ parameters.conda_packages }} - else - # do mamba installs from environment on Py3k - source activate myEnvironment - mamba install --quiet --channel conda-forge python=${{ parameters.python_version }} ${{ parameters.conda_packages }} - fi + conda create --yes --quiet --name myEnvironment --channel conda-forge python=${{ parameters.python_version }} ${{ parameters.conda_packages }} displayName: Install Anaconda packages - bash: | diff --git a/.en-custom.txt b/.en-custom.txt index 86fd9c879a..2779ba5111 100644 --- a/.en-custom.txt +++ b/.en-custom.txt @@ -271,6 +271,9 @@ pre-assembled pre-built pre-installed Pre-Installed +preallocated +preallocation +preallocating preconditioner preconditioners Pusztai diff --git a/fipy/boundaryConditions/boundaryCondition.py b/fipy/boundaryConditions/boundaryCondition.py index 0b781c0044..07691a9cc8 100644 --- a/fipy/boundaryConditions/boundaryCondition.py +++ b/fipy/boundaryConditions/boundaryCondition.py @@ -57,8 +57,8 @@ def _buildMatrix(self, SparseMatrix, Ncells, MaxFaces, coeff): Ncells : int Size of matrices MaxFaces : int - Maximum number of faces per cell (determines bandwidth of - :math:`\mathsf{L}`) + Maximum number of faces per cell (determines number of + non-zeros per row of :math:`\mathsf{L}`) coeff : list Contribution due to this face """ diff --git a/fipy/boundaryConditions/fixedValue.py b/fipy/boundaryConditions/fixedValue.py index a3eff96d69..f950d2cde0 100644 --- a/fipy/boundaryConditions/fixedValue.py +++ b/fipy/boundaryConditions/fixedValue.py @@ -51,15 +51,15 @@ def _buildMatrix(self, SparseMatrix, Ncells, MaxFaces, coeff): Ncells : int Size of matrices MaxFaces : int - Maximum number of faces per cell (determines bandwidth of - :math:`\mathsf{L}`) + Maximum number of faces per cell (determines number of + non-zeros per row of :math:`\mathsf{L}`) coeff : list Contribution to adjacent cell diagonal and :math:`\mathsf{b}` vector by this exterior face """ faces = self.faces.value - LL = SparseMatrix(mesh=self.faces.mesh, sizeHint=len(self.faces), bandwidth=1) + LL = SparseMatrix(mesh=self.faces.mesh, nonZerosPerRow=1) LL.addAt(coeff['cell 1 diag'][faces], self.adjacentCellIDs, self.adjacentCellIDs) ## The following has been commented out because diff --git a/fipy/matrices/offsetSparseMatrix.py b/fipy/matrices/offsetSparseMatrix.py index 33cab8c0fd..ce62dd42f4 100644 --- a/fipy/matrices/offsetSparseMatrix.py +++ b/fipy/matrices/offsetSparseMatrix.py @@ -14,16 +14,41 @@ class OffsetSparseMatrixClass(SparseMatrix): equationIndex = 0 varIndex = 0 - def __init__(self, mesh, bandwidth=0, sizeHint=None, - numberOfVariables=numberOfVariables, numberOfEquations=numberOfEquations): - SparseMatrix.__init__(self, mesh=mesh, bandwidth=bandwidth, sizeHint=sizeHint, - numberOfVariables=numberOfVariables, numberOfEquations=numberOfEquations) + def __init__(self, mesh, nonZerosPerRow=1, exactNonZeros=False, + numberOfVariables=numberOfVariables, + numberOfEquations=numberOfEquations): + if hasattr(nonZerosPerRow, "__iter__"): + # nonZerosPerRow is an iterable for each row. + # need to pad rows for other equations with zeros. + # can't compare to collections.abc.Iterable because PySparse. + before = self.equationIndex + after = numberOfEquations - self.equationIndex - 1 + N = len(nonZerosPerRow) + nonZerosPerRow = numerix.concatenate([[0] * N * before, + nonZerosPerRow, + [0] * N * after]) + else: + nonZerosPerRow //= numberOfEquations + SparseMatrix.__init__(self, + mesh=mesh, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, + numberOfVariables=numberOfVariables, + numberOfEquations=numberOfEquations) def put(self, vector, id1, id2): - SparseMatrix.put(self, vector, id1 + self.mesh.numberOfCells * self.equationIndex, id2 + self.mesh.numberOfCells * self.varIndex) + N = self.mesh.numberOfCells + SparseMatrix.put(self, + vector=vector, + id1=id1 + N * self.equationIndex, + id2=id2 + N * self.varIndex) def addAt(self, vector, id1, id2): - SparseMatrix.addAt(self, vector, id1 + self.mesh.numberOfCells * self.equationIndex, id2 + self.mesh.numberOfCells * self.varIndex) + N = self.mesh.numberOfCells + SparseMatrix.addAt(self, + vector=vector, + id1=id1 + N * self.equationIndex, + id2=id2 + N * self.varIndex) def addAtDiagonal(self, vector): if type(vector) in [type(1), type(1.)]: diff --git a/fipy/matrices/petscMatrix.py b/fipy/matrices/petscMatrix.py index bea752e65a..4dab46215f 100644 --- a/fipy/matrices/petscMatrix.py +++ b/fipy/matrices/petscMatrix.py @@ -5,6 +5,7 @@ __all__ = [] +from collections.abc import Iterable from petsc4py import PETSc from fipy.tools import numerix @@ -53,7 +54,7 @@ def __add__(self, other): """ Add two sparse matrices - >>> L = _PETScMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> print(L + _PETScIdentityMatrix(size=3)) 1.000000 10.000000 3.000000 @@ -103,9 +104,9 @@ def __mul__(self, other): """ Multiply a sparse matrix by another sparse matrix - >>> L1 = _PETScMatrixFromShape(rows=3, cols=3, bandwidth=2) + >>> L1 = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=2) >>> L1.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) - >>> L2 = _PETScIdentityMatrix(size=3, bandwidth=3) + >>> L2 = _PETScIdentityMatrix(size=3, nonZerosPerRow=3) >>> L2.put([4.38], [2], [1]) >>> L2.put([4.38,12357.2,1.1], [2,1,0], [1,0,2]) @@ -165,6 +166,8 @@ def __rmul__(self, other): x = PETSc.Vec().createMPI(N, comm=self.matrix.comm) y = x.duplicate() x[:] = other + self.matrix.assemble() + x.assemble() self.matrix.multTranspose(x, y) arr = numerix.asarray(y) x.destroy() @@ -185,7 +188,7 @@ def put(self, vector, id1, id2): """ Put elements of `vector` at positions of the matrix corresponding to (`id1`, `id2`) - >>> L = _PETScMatrixFromShape(rows=3, cols=3, bandwidth=2) + >>> L = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=2) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> print(L) --- 10.000000 3.000000 @@ -242,7 +245,7 @@ def putDiagonal(self, vector): """ Put elements of `vector` along diagonal of matrix - >>> L = _PETScMatrixFromShape(rows=3, cols=3, bandwidth=1) + >>> L = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=1) >>> L.putDiagonal([3.,10.,numerix.pi]) >>> print(L) 3.000000 --- --- @@ -279,7 +282,7 @@ def addAt(self, vector, id1, id2): """ Add elements of `vector` to the positions in the matrix corresponding to (`id1`,`id2`) - >>> L = _PETScMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> L.addAt([1.73,2.2,8.4,3.9,1.23], [1,2,0,0,1], [2,2,0,0,2]) >>> print(L) @@ -318,7 +321,7 @@ def CSR(self): Examples -------- - >>> L = _PETScMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> L.addAt([1.73,2.2,8.4,3.9,1.23], [1,2,0,0,1], [2,2,0,0,2]) >>> ptrs, cols, data = L.CSR @@ -348,7 +351,7 @@ def LIL(self): Examples -------- - >>> L = _PETScMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> L.addAt([1.73,2.2,8.4,3.9,1.23], [1,2,0,0,1], [2,2,0,0,2]) >>> rows, data = L.LIL @@ -480,7 +483,9 @@ def T(self): class _PETScMatrixFromShape(_PETScMatrix): - def __init__(self, rows, cols, bandwidth=0, sizeHint=None, matrix=None, comm=PETSc.COMM_SELF): + def __init__(self, rows, cols, + nonZerosPerRow=0, exactNonZeros=False, + matrix=None, comm=PETSc.COMM_SELF): """Instantiates and wraps a PETSc `Mat` matrix Parameters @@ -489,17 +494,20 @@ def __init__(self, rows, cols, bandwidth=0, sizeHint=None, matrix=None, comm=PET The number of local matrix rows. cols : int The number of local matrix columns - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~petsc4py.PETSc.Mat Pre-assembled PETSc matrix to use for storage. comm : ~PETSc.Comm The MPI communicator to use. """ - if (bandwidth == 0) and (sizeHint is not None): - bandwidth = sizeHint // max(rows, cols) if matrix is None: matrix = PETSc.Mat() matrix.create(comm) @@ -508,13 +516,19 @@ def __init__(self, rows, cols, bandwidth=0, sizeHint=None, matrix=None, comm=PET matrix.setSizes([[rows, None], [cols, None]]) matrix.setType('aij') # sparse matrix.setUp() -# matrix.setPreallocationNNZ(bandwidth) # FIXME: ??? None, bandwidth -# matrix.setOption(matrix.Option.NEW_NONZERO_ALLOCATION_ERR, False) + if isinstance(nonZerosPerRow, Iterable) or (nonZerosPerRow > 0): + if isinstance(nonZerosPerRow, Iterable): + nonZerosPerRow = numerix.asarray(nonZerosPerRow, dtype=PETSc.IntType) + matrix.setPreallocationNNZ(nonZerosPerRow) + if not exactNonZeros: + matrix.setOption(matrix.Option.NEW_NONZERO_ALLOCATION_ERR, False) super(_PETScMatrixFromShape, self).__init__(matrix=matrix) class _PETScBaseMeshMatrix(_PETScMatrixFromShape): - def __init__(self, mesh, rows, cols, m2m, bandwidth=0, sizeHint=None, matrix=None): + def __init__(self, mesh, rows, cols, m2m, + nonZerosPerRow=0, exactNonZeros=False, + matrix=None): """Creates a `_PETScMatrixFromShape` associated with a `Mesh`. Parameters @@ -527,10 +541,15 @@ def __init__(self, mesh, rows, cols, m2m, bandwidth=0, sizeHint=None, matrix=Non The number of local matrix columns. m2m : ~fipy.matrices.sparseMatrix._Mesh2Matrix Object to convert between mesh coordinates and matrix coordinates. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros. + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~petsc4py.PETSc.Mat Pre-assembled PETSc matrix to use for storage. """ @@ -539,14 +558,14 @@ def __init__(self, mesh, rows, cols, m2m, bandwidth=0, sizeHint=None, matrix=Non super(_PETScBaseMeshMatrix, self).__init__(rows=rows, cols=cols, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, comm=mesh.communicator.petsc4py_comm) def copy(self): tmp = super(_PETScBaseMeshMatrix, self).copy() - copy = self.__class__(mesh=self.mesh) # FIXME: ??? , bandwidth=self.bandwidth) + copy = self.__class__(mesh=self.mesh) # FIXME: ??? , nonZerosPerRow=self.nonZerosPerRow) copy.matrix = tmp.matrix return copy @@ -726,8 +745,8 @@ def addAt(self, vector, id1, id2, overlapping=False): super(_PETScBaseMeshMatrix, self).addAt(vector=vector, id1=id1, id2=id2) class _PETScRowMeshMatrix(_PETScBaseMeshMatrix): - def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, - sizeHint=None, matrix=None, m2m=None): + def __init__(self, mesh, cols, numberOfEquations=1, + nonZerosPerRow=0, exactNonZeros=False, matrix=None, m2m=None): """Creates a `_PETScMatrixFromShape` with rows associated with equations. Parameters @@ -739,10 +758,15 @@ def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, numberOfEquations : int The local rows of the matrix are determined by `numberOfEquations * mesh._localNonOverlappingCellIDs`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros. + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~petsc4py.PETSc.Mat Pre-assembled PETSc matrix to use for storage. m2m : ~fipy.matrices.sparseMatrix._RowMesh2Matrix @@ -756,12 +780,13 @@ def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, rows=numberOfEquations * len(mesh._localNonOverlappingCellIDs), cols=cols, m2m=m2m, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix) class _PETScColMeshMatrix(_PETScBaseMeshMatrix): - def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, sizeHint=None, matrix=None): + def __init__(self, mesh, rows, numberOfVariables=1, + nonZerosPerRow=0, exactNonZeros=False, matrix=None): """Creates a `_PETScMatrixFromShape` with columns associated with solution variables. Parameters @@ -773,10 +798,15 @@ def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, sizeHint=None, numberOfVariables : int The local columns of the matrix are determined by `numberOfVariables * mesh.globalNumberOfCells`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros. + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~petsc4py.PETSc.Mat Pre-assembled PETSc matrix to use for storage. """ @@ -787,13 +817,13 @@ def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, sizeHint=None, rows=rows, cols=numberOfVariables * len(mesh._localNonOverlappingCellIDs), m2m=m2m, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix) class _PETScMeshMatrix(_PETScRowMeshMatrix): def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, - bandwidth=0, sizeHint=None, matrix=None): + nonZerosPerRow=0, exactNonZeros=False, matrix=None): """Creates a `_PETScBaseMeshMatrix` associated with equations and variables. Parameters @@ -806,10 +836,15 @@ def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, numberOfEquations : int The local rows of the matrix are determined by `numberOfEquations * len(mesh._localNonOverlappingCellIDs)`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~petsc4py.PETSc.Mat Pre-assembled PETSc matrix to use for storage. """ @@ -820,17 +855,17 @@ def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, super(_PETScMeshMatrix, self).__init__(mesh=mesh, cols=numberOfVariables * len(mesh._localNonOverlappingCellIDs), numberOfEquations=numberOfEquations, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, m2m=m2m) def __mul__(self, other): """Multiply a sparse matrix by another sparse matrix - >>> L1 = _PETScMatrixFromShape(rows=3, cols=3, bandwidth=2) + >>> L1 = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=2) >>> L1.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) - >>> L2 = _PETScIdentityMatrix(size=3, bandwidth=3) + >>> L2 = _PETScIdentityMatrix(size=3, nonZerosPerRow=3) >>> L2.put([4.38], [2], [1]) >>> L2.put([4.38,12357.2,1.1], [2,1,0], [1,0,2]) @@ -905,7 +940,7 @@ def flush(self): def _test(self): """Tests - >>> m = _PETScMatrixFromShape(rows=3, cols=3, bandwidth=1) + >>> m = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=1) >>> m.addAt((1., 0., 2.), (0, 2, 1), (1, 2, 0)) >>> m.matrix.assemble() @@ -918,13 +953,52 @@ def _test(self): True >>> print(numerix.allclose(val, [1., 2., 0.])) True + + Storing more than preallocated is an error when `exactNonZeros` is set + + >>> m = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=1, exactNonZeros=True) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) # doctest: +IGNORE_EXCEPTION_DETAIL + Traceback (most recent call last): + ... + petsc4py.PETSc.Error: error code 63 + + This is also true if multiple values are accumulated into the + same matrix entry. + + >>> m = _PETScMatrixFromShape(rows=3, cols=3, nonZerosPerRow=1, exactNonZeros=True) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,0], [2,1,1,1]) # doctest: +IGNORE_EXCEPTION_DETAIL + Traceback (most recent call last): + ... + petsc4py.PETSc.Error: error code 63 + + Preallocation can be specified row-by-row + + >>> m = _PETScMatrixFromShape(rows=3, cols=3, + ... nonZerosPerRow=[2, 1, 1]) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) + + Preallocating on the wrong rows is not an error + + >>> m = _PETScMatrixFromShape(rows=3, cols=3, + ... nonZerosPerRow=[1, 2, 1]) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) + + but it is when `exactNonZeros` is specified. + + >>> m = _PETScMatrixFromShape(rows=3, cols=3, + ... nonZerosPerRow=[1, 2, 1], + ... exactNonZeros=True) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) # doctest: +IGNORE_EXCEPTION_DETAIL + Traceback (most recent call last): + ... + petsc4py.PETSc.Error: error code 63 """ pass class _PETScIdentityMatrix(_PETScMatrixFromShape): """Represents a sparse identity matrix for pysparse. """ - def __init__(self, size, bandwidth=1, comm=PETSc.COMM_SELF): + def __init__(self, size, nonZerosPerRow=1, comm=PETSc.COMM_SELF): """Create a sparse matrix with `1` in the diagonal >>> print(_PETScIdentityMatrix(size=3)) @@ -932,12 +1006,12 @@ def __init__(self, size, bandwidth=1, comm=PETSc.COMM_SELF): --- 1.000000 --- --- --- 1.000000 """ - _PETScMatrixFromShape.__init__(self, rows=size, cols=size, bandwidth=bandwidth, comm=comm) + _PETScMatrixFromShape.__init__(self, rows=size, cols=size, nonZerosPerRow=nonZerosPerRow, comm=comm) ids = numerix.arange(size) self.put(numerix.ones(size, 'd'), ids, ids) class _PETScIdentityMeshMatrix(_PETScIdentityMatrix): - def __init__(self, mesh, bandwidth=1): + def __init__(self, mesh, nonZerosPerRow=1): """Create a sparse matrix associated with a `Mesh` with `1` in the diagonal >>> from fipy import Grid1D @@ -948,7 +1022,7 @@ def __init__(self, mesh, bandwidth=1): --- 1.000000 --- --- --- 1.000000 """ - _PETScIdentityMatrix.__init__(self, size=mesh.numberOfCells, bandwidth=bandwidth, + _PETScIdentityMatrix.__init__(self, size=mesh.numberOfCells, nonZerosPerRow=nonZerosPerRow, comm=mesh.communicator.petsc4py_comm) def _test(): diff --git a/fipy/matrices/pysparseMatrix.py b/fipy/matrices/pysparseMatrix.py index 5833039310..b46cbada3a 100644 --- a/fipy/matrices/pysparseMatrix.py +++ b/fipy/matrices/pysparseMatrix.py @@ -289,7 +289,7 @@ def CSR(self): Examples -------- - >>> L = _PysparseMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _PysparseMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> L.addAt([1.73,2.2,8.4,3.9,1.23], [1,2,0,0,1], [2,2,0,0,2]) >>> ptrs, cols, data = L.CSR @@ -324,7 +324,7 @@ def LIL(self): Examples -------- - >>> L = _PysparseMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _PysparseMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> L.addAt([1.73,2.2,8.4,3.9,1.23], [1,2,0,0,1], [2,2,0,0,2]) >>> rows, data = L.LIL @@ -394,7 +394,9 @@ def T(self): class _PysparseMatrixFromShape(_PysparseMatrix): - def __init__(self, rows, cols, bandwidth=0, sizeHint=None, matrix=None, storeZeros=True): + def __init__(self, rows, cols, + nonZerosPerRow=0, exactNonZeros=False, + matrix=None, storeZeros=True): """Instantiates and wraps a Pysparse `ll_mat` matrix Parameters @@ -403,18 +405,24 @@ def __init__(self, rows, cols, bandwidth=0, sizeHint=None, matrix=None, storeZer The number of matrix rows cols : int The number of matrix columns - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + *ignored* matrix : ~pysparse.spmatrix.ll_mat Pre-assembled Pysparse matrix to use for storage storeZeros : bool Instructs pysparse to store zero values if possible. """ - sizeHint = sizeHint or max(rows, cols) * bandwidth if matrix is None: tmpMatrix = spmatrix.ll_mat(1, 1, 1) + try: + assert len(nonZerosPerRow) == rows + sizeHint = numerix.sum(numerix.asarray(nonZerosPerRow)) + except TypeError: + sizeHint = nonZerosPerRow * rows if hasattr(tmpMatrix, 'storeZeros'): matrix = spmatrix.ll_mat(rows, cols, sizeHint, storeZeros) else: @@ -423,7 +431,9 @@ def __init__(self, rows, cols, bandwidth=0, sizeHint=None, matrix=None, storeZer super(_PysparseMatrixFromShape, self).__init__(matrix=matrix) class _PysparseBaseMeshMatrix(_PysparseMatrixFromShape): - def __init__(self, mesh, rows, cols, bandwidth=0, sizeHint=None, matrix=None, storeZeros=True): + def __init__(self, mesh, rows, cols, + nonZerosPerRow=0, exactNonZeros=False, + matrix=None, storeZeros=True): """Creates a `_PysparseMatrixFromShape` associated with a `Mesh`. Parameters @@ -434,10 +444,15 @@ def __init__(self, mesh, rows, cols, bandwidth=0, sizeHint=None, matrix=None, st The number of local matrix rows. cols : int The number of local matrix columns. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros. + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~pysparse.spmatrix.ll_mat Pre-assembled SciPy matrix to use for storage. storeZeros : bool @@ -447,8 +462,8 @@ def __init__(self, mesh, rows, cols, bandwidth=0, sizeHint=None, matrix=None, st super(_PysparseBaseMeshMatrix, self).__init__(rows=rows, cols=cols, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, storeZeros=storeZeros) @@ -497,8 +512,9 @@ def addAt(self, vector, id1, id2, overlapping=False): super(_PysparseBaseMeshMatrix, self).addAt(vector=vector, id1=id1, id2=id2) class _PysparseRowMeshMatrix(_PysparseBaseMeshMatrix): - def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, - sizeHint=None, matrix=None, storeZeros=True): + def __init__(self, mesh, cols, numberOfEquations=1, + nonZerosPerRow=0, exactNonZeros=False, + matrix=None, storeZeros=True): """Creates a `_PysparseBaseMeshMatrix` with rows associated with equations. Parameters @@ -510,10 +526,15 @@ def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, numberOfEquations : int The rows of the matrix are determined by `numberOfEquations * mesh.numberOfCells`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~pysparse.spmatrix.ll_mat Pre-assembled Pysparse matrix to use for storage. storeZeros : bool @@ -524,14 +545,15 @@ def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, super(_PysparseRowMeshMatrix, self).__init__(mesh=mesh, rows=numberOfEquations * mesh.numberOfCells, cols=cols, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, storeZeros=storeZeros) class _PysparseColMeshMatrix(_PysparseBaseMeshMatrix): - def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, - sizeHint=None, matrix=None, storeZeros=True): + def __init__(self, mesh, rows, numberOfVariables=1, + nonZerosPerRow=0, exactNonZeros=False, + matrix=None, storeZeros=True): """Creates a `_PysparseBaseMeshMatrix` with columns associated with solution variables. Parameters @@ -543,10 +565,15 @@ def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, numberOfVariables : int The columns of the matrix are determined by `numberOfVariables * mesh.globalNumberOfCells`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~pysparse.spmatrix.ll_mat Pre-assembled Pysparse matrix to use for storage. storeZeros : bool @@ -557,14 +584,15 @@ def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, super(_PysparseColMeshMatrix, self).__init__(mesh=mesh, rows=rows, cols=numberOfVariables * mesh.numberOfCells, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, storeZeros=storeZeros) class _PysparseMeshMatrix(_PysparseRowMeshMatrix): def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, - bandwidth=0, sizeHint=None, matrix=None, storeZeros=True): + nonZerosPerRow=0, exactNonZeros=False, + matrix=None, storeZeros=True): """Creates a `_PysparseBaseMeshMatrix` with associated with equations and variables. Parameters @@ -577,10 +605,15 @@ def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, numberOfEquations : int The rows of the matrix are determined by `numberOfEquations * mesh.numberOfCells`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~pysparse.spmatrix.ll_mat Pre-assembled Pysparse matrix to use for storage storeZeros : bool @@ -591,8 +624,8 @@ def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, super(_PysparseMeshMatrix, self).__init__(mesh=mesh, cols=numberOfVariables * mesh.numberOfCells, numberOfEquations=numberOfEquations, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, storeZeros=storeZeros) @@ -617,12 +650,12 @@ def asTrilinosMeshMatrix(self): if not hasattr(self, 'trilinosMatrix'): if A.shape[0] == 0: - bandwidth = 0 + nonZerosPerRow = 0 else: - bandwidth = int(numerix.ceil(float(len(values)) / float(A.shape[0]))) - bandwidth = 1 + nonZerosPerRow = int(numerix.ceil(float(len(values)) / float(A.shape[0]))) + nonZerosPerRow = 1 from fipy.matrices.trilinosMatrix import _TrilinosMeshMatrixKeepStencil - tmp = _TrilinosMeshMatrixKeepStencil(mesh=self.mesh, bandwidth=bandwidth, + tmp = _TrilinosMeshMatrixKeepStencil(mesh=self.mesh, nonZerosPerRow=nonZerosPerRow, numberOfVariables=self.numberOfVariables, numberOfEquations=self.numberOfEquations) self.trilinosMatrix = tmp @@ -674,6 +707,35 @@ def _test(self): >>> print(numerix.allequal(list(m.matrix.values()), numerix.array([1.0, 2.0]))) True + Storing more than preallocated is not an error when `exactNonZeros` is set + + >>> m = _PysparseMatrixFromShape(rows=3, cols=3, nonZerosPerRow=1, exactNonZeros=True) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) # doctest: +IGNORE_EXCEPTION_DETAIL + + This is also true if multiple values are accumulated into the + same matrix entry. + + >>> m = _PysparseMatrixFromShape(rows=3, cols=3, nonZerosPerRow=1, exactNonZeros=True) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,0], [2,1,1,1]) # doctest: +IGNORE_EXCEPTION_DETAIL + + Preallocation can be specified row-by-row + + >>> m = _PysparseMatrixFromShape(rows=3, cols=3, + ... nonZerosPerRow=[2, 1, 1]) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) + + Preallocating on the wrong rows is not an error + + >>> m = _PysparseMatrixFromShape(rows=3, cols=3, + ... nonZerosPerRow=[1, 2, 1]) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) + + even when `exactNonZeros` is specified. + + >>> m = _PysparseMatrixFromShape(rows=3, cols=3, + ... nonZerosPerRow=[1, 2, 1], + ... exactNonZeros=True) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) """ pass @@ -689,7 +751,7 @@ def __init__(self, size): --- 1.000000 --- --- --- 1.000000 """ - _PysparseMatrixFromShape.__init__(self, rows=size, cols=size, bandwidth=1) + _PysparseMatrixFromShape.__init__(self, rows=size, cols=size, nonZerosPerRow=1) ids = numerix.arange(size) self.put(numerix.ones(size, 'd'), ids, ids) diff --git a/fipy/matrices/scipyMatrix.py b/fipy/matrices/scipyMatrix.py index 708cf8e139..d0b0f2a90f 100644 --- a/fipy/matrices/scipyMatrix.py +++ b/fipy/matrices/scipyMatrix.py @@ -293,7 +293,7 @@ def CSR(self): Examples -------- - >>> L = _ScipyMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _ScipyMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> L.addAt([1.73,2.2,8.4,3.9,1.23], [1,2,0,0,1], [2,2,0,0,2]) >>> ptrs, cols, data = L.CSR @@ -321,7 +321,7 @@ def LIL(self): Examples -------- - >>> L = _ScipyMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _ScipyMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> L.addAt([1.73,2.2,8.4,3.9,1.23], [1,2,0,0,1], [2,2,0,0,2]) >>> rows, data = L.LIL @@ -377,17 +377,19 @@ def T(self): class _ScipyMatrixFromShape(_ScipyMatrix): - def __init__(self, rows, cols, bandwidth=0, sizeHint=None, matrix=None, storeZeros=True): + def __init__(self, rows, cols, + nonZerosPerRow=0, exactNonZeros=False, + matrix=None, storeZeros=True): """Instantiates and wraps a scipy sparse matrix Parameters ---------- mesh : ~fipy.meshes.mesh.Mesh The `Mesh` to assemble the matrix for. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + *ignored* + exactNonZeros : bool + *ignored* matrix : ~scipy.sparse.csr_matrix Pre-assembled SciPy matrix to use for storage. storeZeros : bool @@ -399,7 +401,8 @@ def __init__(self, rows, cols, bandwidth=0, sizeHint=None, matrix=None, storeZer super(_ScipyMatrixFromShape, self).__init__(matrix=matrix) class _ScipyBaseMeshMatrix(_ScipyMatrixFromShape): - def __init__(self, mesh, rows, cols, bandwidth=0, sizeHint=None, + def __init__(self, mesh, rows, cols, + nonZerosPerRow=0, exactNonZeros=False, matrix=None, storeZeros=True): """Creates a `_ScipyMatrixFromShape` associated with a `Mesh`. @@ -411,10 +414,10 @@ def __init__(self, mesh, rows, cols, bandwidth=0, sizeHint=None, The number of local matrix rows. cols : int The number of local matrix columns. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros. + nonZerosPerRow : int or array_like of int + *ignored* + exactNonZeros : bool + *ignored* matrix : ~scipy.sparse.csr_matrix Pre-assembled SciPy matrix to use for storage. storeZeros : bool @@ -424,8 +427,8 @@ def __init__(self, mesh, rows, cols, bandwidth=0, sizeHint=None, super(_ScipyBaseMeshMatrix, self).__init__(rows=rows, cols=cols, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, storeZeros=storeZeros) @@ -442,8 +445,9 @@ def _getGhostedValues(self, var): return var.value class _ScipyRowMeshMatrix(_ScipyBaseMeshMatrix): - def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, - sizeHint=None, matrix=None, storeZeros=True): + def __init__(self, mesh, cols, numberOfEquations=1, + nonZerosPerRow=0, exactNonZeros=False, + matrix=None, storeZeros=True): """Creates a `_ScipyBaseMeshMatrix` with rows associated with equations. Parameters @@ -455,10 +459,10 @@ def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, numberOfEquations : int The rows of the matrix are determined by `numberOfEquations * mesh.numberOfCells`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + *ignored* + exactNonZeros : bool + *ignored* matrix : ~scipy.sparse.csr_matrix Pre-assembled SciPy matrix to use for storage. storeZeros : bool @@ -469,14 +473,15 @@ def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, super(_ScipyRowMeshMatrix, self).__init__(mesh=mesh, rows=numberOfEquations * mesh.numberOfCells, cols=cols, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, storeZeros=storeZeros) class _ScipyColMeshMatrix(_ScipyBaseMeshMatrix): - def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, - sizeHint=None, matrix=None, storeZeros=True): + def __init__(self, mesh, rows, numberOfVariables=1, + nonZerosPerRow=0, exactNonZeros=False, + matrix=None, storeZeros=True): """Creates a `_ScipyBaseMeshMatrix` with columns associated with solution variables. Parameters @@ -488,10 +493,10 @@ def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, numberOfVariables : int The columns of the matrix are determined by `numberOfVariables * mesh.globalNumberOfCells`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros. + nonZerosPerRow : int or array_like of int + *ignored* + exactNonZeros : bool + *ignored* matrix : ~scipy.sparse.csr_matrix Pre-assembled SciPy matrix to use for storage. storeZeros : bool @@ -502,14 +507,14 @@ def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, super(_ScipyColMeshMatrix, self).__init__(mesh=mesh, rows=rows, cols=numberOfVariables * mesh.numberOfCells, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, storeZeros=storeZeros) class _ScipyMeshMatrix(_ScipyRowMeshMatrix): def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, - bandwidth=0, sizeHint=None, matrix=None, storeZeros=True): + nonZerosPerRow=0, exactNonZeros=False, matrix=None, storeZeros=True): """Creates a `_ScipyBaseMeshMatrix` associated with equations and variables. Parameters @@ -522,10 +527,10 @@ def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, numberOfEquations : int The rows of the matrix are determined by `numberOfEquations * mesh.numberOfCells`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + *ignored* + exactNonZeros : bool + *ignored* matrix : ~scipy.sparse.csr_matrix Pre-assembled SciPy matrix to use for storage. storeZeros : bool @@ -535,8 +540,8 @@ def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, super(_ScipyMeshMatrix, self).__init__(mesh=mesh, cols=numberOfVariables * mesh.numberOfCells, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, numberOfEquations=numberOfEquations, storeZeros=storeZeros) @@ -602,7 +607,7 @@ def __init__(self, size): --- 1.000000 --- --- --- 1.000000 """ - _ScipyMatrixFromShape.__init__(self, rows=size, cols=size, bandwidth=1) + _ScipyMatrixFromShape.__init__(self, rows=size, cols=size, nonZerosPerRow=1) ids = numerix.arange(size) self.put(numerix.ones(size, 'd'), ids, ids) diff --git a/fipy/matrices/trilinosMatrix.py b/fipy/matrices/trilinosMatrix.py index ddee99cffa..e487c06d48 100644 --- a/fipy/matrices/trilinosMatrix.py +++ b/fipy/matrices/trilinosMatrix.py @@ -38,23 +38,25 @@ class _TrilinosMatrix(_SparseMatrix): Allows basic python operations __add__, __sub__ etc. Facilitate matrix populating in an easy way. """ - def __init__(self, matrix, bandwidth=None): + def __init__(self, matrix, nonZerosPerRow=None): """ Parameters ---------- matrix : Epetra.CrsMatrix The internal Trilinos matrix - bandwidth : int - The proposed band width of the matrix. + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). """ self.matrix = matrix self.comm = matrix.Comm() - if bandwidth is None: - self.bandwidth = ((matrix.NumGlobalNonzeros() + matrix.NumGlobalRows() - 1) + if nonZerosPerRow is None: + self.nonZerosPerRow = ((matrix.NumGlobalNonzeros() + matrix.NumGlobalRows() - 1) // matrix.NumGlobalRows()) else: - self.bandwidth = bandwidth + self.nonZerosPerRow = nonZerosPerRow super(_TrilinosMatrix, self).__init__() @@ -523,7 +525,7 @@ def _getDistributedMatrix(self): DistributedMap = Epetra.Map(totalElements, 0, self.comm) RootToDist = Epetra.Import(DistributedMap, self.rangeMap) - DistMatrix = Epetra.CrsMatrix(Epetra.Copy, DistributedMap, (self.bandwidth*3) // 2) + DistMatrix = Epetra.CrsMatrix(Epetra.Copy, DistributedMap, (self.nonZerosPerRow*3) // 2) DistMatrix.Import(self.matrix, RootToDist, Epetra.Insert) @@ -554,7 +556,7 @@ def CSR(self): Examples -------- - >>> L = _TrilinosMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _TrilinosMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> L.addAt([1.73,2.2,8.4,3.9,1.23], [1,2,0,0,1], [2,2,0,0,2]) >>> ptrs, cols, data = L.CSR @@ -589,7 +591,7 @@ def LIL(self): Examples -------- - >>> L = _TrilinosMatrixFromShape(rows=3, cols=3, bandwidth=3) + >>> L = _TrilinosMatrixFromShape(rows=3, cols=3, nonZerosPerRow=3) >>> L.put([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) >>> L.addAt([1.73,2.2,8.4,3.9,1.23], [1,2,0,0,1], [2,2,0,0,2]) >>> rows, data = L.LIL @@ -711,7 +713,8 @@ def T(self): return _TrilinosMatrix(matrix=A_T_bis) class _TrilinosMatrixFromShape(_TrilinosMatrix): - def __init__(self, rows, cols, bandwidth=1, sizeHint=None, matrix=None): + def __init__(self, rows, cols, + nonZerosPerRow=1, exactNonZeros=False, matrix=None): """Instantiates and wraps an `Epetra.CrsMatrix` Parameters @@ -720,22 +723,21 @@ def __init__(self, rows, cols, bandwidth=1, sizeHint=None, matrix=None): The number of matrix rows cols : int The number of matrix columns - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 1). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~PyTrilinos.Epetra.CrsMatrix Pre-assembled Trilinos matrix to use for storage. """ self.rows = rows self.cols = cols - size = max(rows, cols) - if sizeHint is not None and bandwidth == 0: - bandwidth = (sizeHint + size - 1) // (size or 1) - else: - bandwidth = bandwidth - if matrix is None: # On Mar 26, 2011, at 11:04 AM, Williams, Alan B wrote: # @@ -743,15 +745,18 @@ def __init__(self, rows, cols, bandwidth=1, sizeHint=None, matrix=None): # matrices with only a row-map and allow the column-map to be # generated internally. Especially if you have a rectangular # matrix, which is even trickier to get correct. - matrix = Epetra.CrsMatrix(Epetra.Copy, self.rowMap, (bandwidth*3)//2) - # Leave extra bandwidth, to handle multiple insertions into the - # same spot. It's memory-inefficient, but it'll get cleaned up when - # FillComplete is called, and according to the Trilinos devs the - # performance boost will be worth it. + # Pre-allocate extra non-zeros, to handle multiple insertions into the + # same spot. It's memory-inefficient, but it'll get cleaned up when + # FillComplete is called, and according to the Trilinos devs the + # performance boost will be worth it. + nonZerosPerRow = numerix.asarray(nonZerosPerRow, dtype='int32') + StaticProfile = (nonZerosPerRow*3)//2 + + matrix = Epetra.CrsMatrix(Epetra.Copy, self.rowMap, StaticProfile, exactNonZeros) super(_TrilinosMatrixFromShape, self).__init__(matrix=matrix, - bandwidth=bandwidth) + nonZerosPerRow=nonZerosPerRow) @property def rowMap(self): @@ -772,7 +777,8 @@ def colMap(self): return Epetra.Map(self.cols, self.cols, 0, comm) class _TrilinosBaseMeshMatrix(_TrilinosMatrixFromShape): - def __init__(self, mesh, rows, cols, m2m, bandwidth=0, sizeHint=None, matrix=None): + def __init__(self, mesh, rows, cols, m2m, + nonZerosPerRow=0, exactNonZeros=False, matrix=None): """Creates a `_TrilinosMatrixFromShape` associated with a `Mesh`. Parameters @@ -785,10 +791,15 @@ def __init__(self, mesh, rows, cols, m2m, bandwidth=0, sizeHint=None, matrix=Non The number of local matrix columns. m2m : ~fipy.matrices.sparseMatrix._Mesh2Matrix Object to convert between mesh coordinates and matrix coordinates. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros. + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~PyTrilinos.Epetra.CrsMatrix Pre-assembled Trilinos matrix to use for storage. """ @@ -797,8 +808,8 @@ def __init__(self, mesh, rows, cols, m2m, bandwidth=0, sizeHint=None, matrix=Non super(_TrilinosBaseMeshMatrix, self).__init__(rows=rows, cols=cols, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix) @property @@ -821,7 +832,7 @@ def colMap(self): def copy(self): tmp = super(_TrilinosBaseMeshMatrix, self).copy() - copy = self.__class__(mesh=self.mesh, bandwidth=self.bandwidth) + copy = self.__class__(mesh=self.mesh, nonZerosPerRow=self.nonZerosPerRow) copy.matrix = tmp.matrix return copy @@ -936,8 +947,8 @@ def addAt(self, vector, id1, id2, overlapping=False): super(_TrilinosBaseMeshMatrix, self).addAt(vector=vector, id1=id1, id2=id2) class _TrilinosRowMeshMatrix(_TrilinosBaseMeshMatrix): - def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, - sizeHint=None, matrix=None, m2m=None): + def __init__(self, mesh, cols, numberOfEquations=1, + nonZerosPerRow=0, exactNonZeros=False, matrix=None, m2m=None): """Creates a `_TrilinosMatrixFromShape` with rows associated with equations. Parameters @@ -949,10 +960,15 @@ def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, numberOfEquations : int The local rows of the matrix are determined by `numberOfEquations * mesh._localNonOverlappingCellIDs`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros. + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 1). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~PyTrilinos.Epetra.CrsMatrix Pre-assembled Trilinos matrix to use for storage. m2m : ~fipy.matrices.sparseMatrix._RowMesh2Matrix @@ -967,8 +983,8 @@ def __init__(self, mesh, cols, numberOfEquations=1, bandwidth=0, rows=rows, cols=cols, m2m=m2m, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix) @property @@ -988,8 +1004,8 @@ def domainMap(self): return self.rowMap class _TrilinosColMeshMatrix(_TrilinosBaseMeshMatrix): - def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, - sizeHint=None, matrix=None): + def __init__(self, mesh, rows, numberOfVariables=1, + nonZerosPerRow=0, exactNonZeros=False, matrix=None): """Creates a `_TrilinosMatrixFromShape` with columns associated with solution variables. Parameters @@ -1001,10 +1017,15 @@ def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, numberOfVariables : int The local columns of the matrix are determined by `numberOfVariables * mesh.globalNumberOfCells`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros. + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~PyTrilinos.Epetra.CrsMatrix Pre-assembled Trilinos matrix to use for storage. """ @@ -1016,8 +1037,8 @@ def __init__(self, mesh, rows, numberOfVariables=1, bandwidth=0, rows=rows, cols=cols, m2m=m2m, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix) @property @@ -1041,7 +1062,7 @@ def domainMap(self): class _TrilinosMeshMatrix(_TrilinosRowMeshMatrix): def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, - bandwidth=0, sizeHint=None, matrix=None): + nonZerosPerRow=0, exactNonZeros=False, matrix=None): """Creates a `_TrilinosRowMeshMatrix` associated with equations and variables. Parameters @@ -1054,10 +1075,15 @@ def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, numberOfEquations : int The local rows of the matrix are determined by `numberOfEquations * len(mesh._localNonOverlappingCellIDs)`. - bandwidth : int - The proposed band width of the matrix. - sizeHint : int - Estimate of the number of non-zeros + nonZerosPerRow : int or array_like of int + The approximate number of sparse entries per row. Either a + typical number, or an iterable of values for each row + (default: 0). + exactNonZeros : bool + Whether `nonZerosPerRow` is exact or approximate. + Performance is improved preallocation is exact, but errors + can result if additional allocations are necessary. + (default: False). matrix : ~PyTrilinos.Epetra.CrsMatrix Pre-assembled Trilinos matrix to use for storage. """ @@ -1069,8 +1095,8 @@ def __init__(self, mesh, numberOfVariables=1, numberOfEquations=1, super(_TrilinosMeshMatrix, self).__init__(mesh=mesh, cols=cols, numberOfEquations=numberOfEquations, - bandwidth=bandwidth, - sizeHint=sizeHint, + nonZerosPerRow=nonZerosPerRow, + exactNonZeros=exactNonZeros, matrix=matrix, m2m=m2m) @@ -1091,7 +1117,7 @@ def flush(self): def _getMatrixProperty(self): if not hasattr(self, '_matrix'): self._matrix = _TrilinosMeshMatrix(self.mesh, - bandwidth=self.bandwidth, + nonZerosPerRow=self.nonZerosPerRow, numberOfVariables=self._m2m.numberOfVariables, numberOfEquations=self._m2m.numberOfEquations).matrix return super(_TrilinosMeshMatrix, self).matrix @@ -1196,6 +1222,53 @@ def __mul__(self, other): else: raise TypeError("%s: %s != (%d,)" % (self.__class__, str(shape), N)) + def _test(self): + """Tests + + >>> m = _TrilinosMatrixFromShape(rows=3, cols=3, nonZerosPerRow=1) + >>> m.addAt((1., 0., 2.), (0, 2, 1), (1, 2, 0)) + + Storing more than preallocated is an error when `exactNonZeros` is set + + >>> m = _TrilinosMatrixFromShape(rows=3, cols=3, nonZerosPerRow=1, exactNonZeros=True) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) # doctest: +IGNORE_EXCEPTION_DETAIL + Traceback (most recent call last): + ... + RuntimeError: Processor 0, error code -1 + + This is also true if multiple values are accumulated into the + same matrix entry. + + >>> m = _TrilinosMatrixFromShape(rows=3, cols=3, nonZerosPerRow=1, exactNonZeros=True) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,0], [2,1,1,1]) # doctest: +IGNORE_EXCEPTION_DETAIL + Traceback (most recent call last): + ... + RuntimeError: Processor 0, error code -1 + + Preallocation can be specified row-by-row + + >>> m = _TrilinosMatrixFromShape(rows=3, cols=3, + ... nonZerosPerRow=[2, 1, 1]) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) + + Preallocating on the wrong rows is not an error + + >>> m = _TrilinosMatrixFromShape(rows=3, cols=3, + ... nonZerosPerRow=[1, 2, 1]) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) + + but it is when `exactNonZeros` is specified. + + >>> m = _TrilinosMatrixFromShape(rows=3, cols=3, + ... nonZerosPerRow=[1, 2, 1], + ... exactNonZeros=True) + >>> m.addAt([3.,10.,numerix.pi,2.5], [0,0,1,2], [2,1,1,0]) # doctest: +IGNORE_EXCEPTION_DETAIL + Traceback (most recent call last): + ... + RuntimeError: Processor 0, error code -1 + """ + pass + class _TrilinosIdentityMatrix(_TrilinosMatrixFromShape): """ Represents a sparse identity matrix for Trilinos. @@ -1209,7 +1282,7 @@ def __init__(self, size): --- 1.000000 --- --- --- 1.000000 """ - _TrilinosMatrixFromShape.__init__(self, rows=size, cols=size, bandwidth=1) + _TrilinosMatrixFromShape.__init__(self, rows=size, cols=size, nonZerosPerRow=1) ids = numerix.arange(size) self.addAt(numerix.ones(size, 'l'), ids, ids) @@ -1225,7 +1298,7 @@ def __init__(self, mesh): --- 1.000000 --- --- --- 1.000000 """ - _TrilinosMeshMatrix.__init__(self, mesh=mesh, bandwidth=1) + _TrilinosMeshMatrix.__init__(self, mesh=mesh, nonZerosPerRow=1) size = mesh.numberOfCells ids = numerix.arange(size) self.addAt(numerix.ones(size, 'l'), ids, ids) diff --git a/fipy/meshes/gmshMesh.py b/fipy/meshes/gmshMesh.py index 9814c5ab62..56556b0b9c 100755 --- a/fipy/meshes/gmshMesh.py +++ b/fipy/meshes/gmshMesh.py @@ -2185,10 +2185,13 @@ def _test(self): ... os.close(ftmp) ... else: ... posFile = None + >>> parallelComm.Barrier() >>> posFile = parallelComm.bcast(posFile) + >>> parallelComm.Barrier() >>> f = openPOSFile(posFile, mode='w') # doctest: +GMSH >>> f.write(vol) # doctest: +GMSH >>> f.close() # doctest: +GMSH + >>> parallelComm.Barrier() >>> f = open(posFile, mode='r') # doctest: +GMSH >>> l = f.readlines() # doctest: +GMSH @@ -2244,6 +2247,7 @@ def _test(self): >>> print(numerix.allclose(a1, a2)) # doctest: +GMSH True + >>> parallelComm.Barrier() >>> if parallelComm.procID == 0: ... os.remove(posFile) diff --git a/fipy/terms/abstractDiffusionTerm.py b/fipy/terms/abstractDiffusionTerm.py index 8a49a752bf..b6241d9ffe 100644 --- a/fipy/terms/abstractDiffusionTerm.py +++ b/fipy/terms/abstractDiffusionTerm.py @@ -189,7 +189,8 @@ def __getCoefficientMatrix(self, SparseMatrix, var, coeff): ## print 'id1',id1 ## print 'id2',id2 - coefficientMatrix = SparseMatrix(mesh=mesh, bandwidth = mesh._maxFacesPerCell + 1) + facesPerCell = mesh._facesPerCell[..., mesh._localNonOverlappingCellIDs] + coefficientMatrix = SparseMatrix(mesh=mesh, nonZerosPerRow=facesPerCell + 1) interiorCoeff = numerix.take(coeff, interiorFaces, axis=-1).ravel() coefficientMatrix.addAt(interiorCoeff, id1.ravel(), id1.swapaxes(0, 1).ravel()) coefficientMatrix.addAt(-interiorCoeff, id1.ravel(), id2.swapaxes(0, 1).ravel()) @@ -208,7 +209,7 @@ def __getCoefficientMatrix(self, SparseMatrix, var, coeff): ## ## print interiorCoeff[:,:,0] ## ## print interiorCoeff[:,:,1] -## coefficientMatrix = SparseMatrix(mesh=mesh, bandwidth = mesh._maxFacesPerCell + 1) +## coefficientMatrix = SparseMatrix(mesh=mesh, nonZerosPerRow=mesh._facesPerCell + 1) ## ## print 'numerix.sum(interiorCoeff, -2)',numerix.sum(interiorCoeff, -2) ## ## print numerix.sum(interiorCoeff, -2).ravel() @@ -341,7 +342,7 @@ def __higherOrderbuildMatrix(self, var, SparseMatrix, boundaryConditions=(), dt= del lowerOrderBCs lowerOrderb = lowerOrderb / mesh.cellVolumes - volMatrix = SparseMatrix(mesh=var.mesh, bandwidth = 1) + volMatrix = SparseMatrix(mesh=var.mesh, nonZerosPerRow=1) volMatrix.addAtDiagonal(1. / mesh.cellVolumes) lowerOrderL = volMatrix * lowerOrderL @@ -416,7 +417,7 @@ def __higherOrderbuildMatrix(self, var, SparseMatrix, boundaryConditions=(), dt= else: - L = SparseMatrix(mesh=mesh) + L = SparseMatrix(mesh=mesh, nonZerosPerRow=1) L.addAtDiagonal(mesh.cellVolumes) b = numerix.zeros(len(var.ravel()), 'd') diff --git a/fipy/terms/cellTerm.py b/fipy/terms/cellTerm.py index d6394ab2c7..d746df6d5e 100644 --- a/fipy/terms/cellTerm.py +++ b/fipy/terms/cellTerm.py @@ -121,7 +121,7 @@ def _buildMatrixNoInline_(self, L, oldArray, b, dt, coeffVectors): def _buildMatrix(self, var, SparseMatrix, boundaryConditions=(), dt=None, transientGeomCoeff=None, diffusionGeomCoeff=None): b = numerix.zeros(var.shape, 'd').ravel() - L = SparseMatrix(mesh=var.mesh) + L = SparseMatrix(mesh=var.mesh, nonZerosPerRow=1) coeffVectors = self._getCoeffVectors_(var=var, transientGeomCoeff=transientGeomCoeff, diffusionGeomCoeff=diffusionGeomCoeff) diff --git a/fipy/terms/faceTerm.py b/fipy/terms/faceTerm.py index d1fdf66a09..351f92b2f8 100644 --- a/fipy/terms/faceTerm.py +++ b/fipy/terms/faceTerm.py @@ -140,7 +140,8 @@ def _buildMatrix(self, var, SparseMatrix, boundaryConditions=(), dt=None, transi id2 = numerix.take(id2, interiorFaces) b = numerix.zeros(var.shape, 'd').ravel() - L = SparseMatrix(mesh=mesh) + facesPerCell = mesh._facesPerCell[..., mesh._localNonOverlappingCellIDs] + L = SparseMatrix(mesh=mesh, nonZerosPerRow=facesPerCell + 1) weight = self._getWeight(var, transientGeomCoeff, diffusionGeomCoeff)