Chaste  Release::2018.1
PetscMatTools.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "PetscMatTools.hpp"
37 #include <algorithm>
38 #include <cassert>
39 
40 
42 // Implementation
44 
45 void PetscMatTools::SetElement(Mat matrix, PetscInt row, PetscInt col, double value)
46 {
47  PetscInt lo, hi;
48  GetOwnershipRange(matrix, lo, hi);
49 
50  if (row >= lo && row < hi)
51  {
52  MatSetValue(matrix, row, col, value, INSERT_VALUES);
53  }
54 }
55 
56 void PetscMatTools::AddToElement(Mat matrix, PetscInt row, PetscInt col, double value)
57 {
58  PetscInt lo, hi;
59  GetOwnershipRange(matrix, lo, hi);
60 
61  if (row >= lo && row < hi)
62  {
63  MatSetValue(matrix, row, col, value, ADD_VALUES);
64  }
65 }
66 
68 {
69  MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY);
70  MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY);
71 }
72 
74 {
75  MatAssemblyBegin(matrix, MAT_FLUSH_ASSEMBLY);
76  MatAssemblyEnd(matrix, MAT_FLUSH_ASSEMBLY);
77 }
78 
80 {
81  //Give full precision, scientific notation
82 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 7) // PETSc 3.7+
83  PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_MATLAB);
84 #else
85  PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_MATLAB);
86 #endif
87  MatView(matrix,PETSC_VIEWER_STDOUT_WORLD);
88 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 7) // PETSc 3.7+
89  PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
90 #endif
91 }
92 
93 void PetscMatTools::SetRow(Mat matrix, PetscInt row, double value)
94 {
95  PetscInt lo, hi;
96  GetOwnershipRange(matrix, lo, hi);
97 
98  if (row >= lo && row < hi)
99  {
100  PetscInt rows, cols;
101  MatGetSize(matrix, &rows, &cols);
102  for (PetscInt i=0; i<cols; i++)
103  {
104  SetElement(matrix, row, i, value);
105  }
106  }
107 }
108 
109 void PetscMatTools::ZeroRowsWithValueOnDiagonal(Mat matrix, std::vector<unsigned>& rRows, double diagonalValue)
110 {
111  Finalise(matrix);
112 
113  /*
114  * Important! PETSc by default will destroy the sparsity structure for this row and deallocate memory
115  * when the row is zeroed, and if there is a next timestep, the memory will have to reallocated when
116  * assembly to done again. This can kill performance. The following makes sure the zeroed rows are kept.
117  *
118  * Note: if the following lines are called, then once MatZeroRows() is called below, there will be an
119  * error if some rows have no entries at all. Hence for problems with dummy variables, Stokes flow for
120  * example, the identity block needs to be added before dirichlet BCs.
121  */
122 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 1) //PETSc 3.1 or later
123  MatSetOption(matrix, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
124 #elif (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 0) //PETSc 3.0
125  MatSetOption(matrix, MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
126 #else //PETSc 2.x.x
127  MatSetOption(matrix, MAT_KEEP_ZEROED_ROWS);
128 #endif
129 
130  PetscInt* rows = new PetscInt[rRows.size()];
131  for (unsigned i=0; i<rRows.size(); i++)
132  {
133  rows[i] = rRows[i];
134  }
135 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
136  IS is;
137  ISCreateGeneral(PETSC_COMM_WORLD, rRows.size(), rows, &is);
138  MatZeroRows(matrix, is, &diagonalValue);
139  ISDestroy(PETSC_DESTROY_PARAM(is));
140  /*
141  *
142 [2]PETSC ERROR: MatMissingDiagonal_SeqAIJ() line 1011 in src/mat/impls/aij/seq/aij.c
143 [2]PETSC ERROR: PETSc has generated inconsistent data!
144 [2]PETSC ERROR: Matrix is missing diagonal number 15!
145 [2]PETSC ERROR: MatILUFactorSymbolic_SeqAIJ() line 906 in src/mat/impls/aij/seq/aijfact.c
146  *
147  *
148  // There appears to be a problem with MatZeroRows not setting diagonals correctly
149  // While we are supporting PETSc 2.2, we have to do this the slow way
150 
151  AssembleFinal(matrix);
152  PetscInt lo, hi;
153  GetOwnershipRange(matrix, lo, hi);
154  PetscInt size=GetSize(matrix);
156  for (unsigned index=0; index<rRows.size(); index++)
157  {
158  PetscInt row = rRows[index];
159  if (row >= lo && row < hi)
160  {
161  std::vector<unsigned> non_zero_cols;
162  // This row is local, so zero it.
163  for (PetscInt column = 0; column < size; column++)
164  {
165  if (GetElement(matrix, row, column) != 0.0)
166  {
167  non_zero_cols.push_back(column);
168  }
169  }
170  // Separate "gets" from "sets" so that we don't have to keep going into "assembled" mode
171  for (unsigned i=0; i<non_zero_cols.size();i++)
172  {
173  SetElement(matrix, row, non_zero_cols[i], 0.0);
174  }
175  // Set the diagonal
176  SetElement(matrix, row, row, diagonalValue);
177  }
178  // Everyone communicate after row is finished
179  AssembleFinal(matrix);
180  }
181  */
182 #elif (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 2) //PETSc 3.2 or later
183  MatZeroRows(matrix, rRows.size(), rows, diagonalValue , nullptr, nullptr);
184 #else
185  MatZeroRows(matrix, rRows.size(), rows, diagonalValue);
186 #endif
187  delete [] rows;
188 }
189 
190 
191 void PetscMatTools::ZeroRowsAndColumnsWithValueOnDiagonal(Mat matrix, std::vector<unsigned> rowColIndices, double diagonalValue)
192 {
193  Finalise(matrix);
194 
195  // sort the vector as we will be repeatedly searching for entries in it
196  std::sort(rowColIndices.begin(), rowColIndices.end());
197 
198  PetscInt lo, hi;
199  GetOwnershipRange(matrix, lo, hi);
200  unsigned size = hi-lo;
201 
202  std::vector<unsigned>* cols_to_zero_per_row = new std::vector<unsigned>[size];
203 
204  // collect the columns to be zeroed, for each row. We don't zero yet as we would
205  // then have to repeatedly call Finalise before each MatGetRow
206  for (PetscInt row = lo; row < hi; row++)
207  {
208  // get all the non-zero cols for this row
209  PetscInt num_cols;
210  const PetscInt* cols;
211  MatGetRow(matrix, row, &num_cols, &cols, PETSC_NULL);
212 
213  // see which of these cols are in the list of cols to be zeroed
214  for (PetscInt i=0; i<num_cols; i++)
215  {
216  if (std::binary_search(rowColIndices.begin(), rowColIndices.end(), cols[i]))
217  {
218  cols_to_zero_per_row[row-lo].push_back(cols[i]);
219  }
220  }
221 
222  // this must be called for each MatGetRow
223  MatRestoreRow(matrix, row, &num_cols, &cols, PETSC_NULL);
224  }
225 
226  // Now zero columns of each row
227  for (PetscInt row = lo; row < hi; row++)
228  {
229  unsigned num_cols_to_zero_this_row = cols_to_zero_per_row[row-lo].size();
230 
231  if (num_cols_to_zero_this_row>0)
232  {
233  PetscInt* cols_to_zero = new PetscInt[num_cols_to_zero_this_row];
234  double* zeros = new double[num_cols_to_zero_this_row];
235  for (unsigned i=0; i<num_cols_to_zero_this_row; i++)
236  {
237  cols_to_zero[i] = cols_to_zero_per_row[row-lo][i];
238  zeros[i] = 0.0;
239  }
240 
241  PetscInt rows[1];
242  rows[0] = row;
243  MatSetValues(matrix, 1, rows, num_cols_to_zero_this_row, cols_to_zero, zeros, INSERT_VALUES);
244  delete [] cols_to_zero;
245  delete [] zeros;
246  }
247  }
248 
249  delete [] cols_to_zero_per_row;
250 
251  // Now zero the rows and add the diagonal entries
252  ZeroRowsWithValueOnDiagonal(matrix, rowColIndices, diagonalValue);
253 }
254 
256 {
257  Finalise(matrix);
258 
259  PetscInt lo, hi;
260  GetOwnershipRange(matrix, lo, hi);
261 
262  // Determine which rows in this column are non-zero (and therefore need to be zeroed)
263  std::vector<unsigned> nonzero_rows;
264  for (PetscInt row = lo; row < hi; row++)
265  {
266  if (GetElement(matrix, row, col) != 0.0)
267  {
268  nonzero_rows.push_back(row);
269  }
270  }
271 
272  // Set those rows to be zero by calling MatSetValues
273  unsigned size = nonzero_rows.size();
274  PetscInt* rows = new PetscInt[size];
275  PetscInt cols[1];
276  double* zeros = new double[size];
277 
278  cols[0] = col;
279 
280  for (unsigned i=0; i<size; i++)
281  {
282  rows[i] = nonzero_rows[i];
283  zeros[i] = 0.0;
284  }
285 
286  MatSetValues(matrix, size, rows, 1, cols, zeros, INSERT_VALUES);
287  delete [] rows;
288  delete [] zeros;
289 }
290 
292 {
293  MatZeroEntries(matrix);
294 }
295 
296 unsigned PetscMatTools::GetSize(Mat matrix)
297 {
298  PetscInt rows, cols;
299 
300  MatGetSize(matrix, &rows, &cols);
301  assert(rows == cols);
302  return (unsigned) rows;
303 }
304 
306 {
307  MatGetOwnershipRange(matrix, &lo, &hi);
308 }
309 
311 {
312  PetscInt lo, hi;
313  GetOwnershipRange(matrix, lo, hi);
314 
315  assert(lo <= row && row < hi);
316  PetscInt row_as_array[1];
317  row_as_array[0] = row;
318  PetscInt col_as_array[1];
319  col_as_array[0] = col;
320 
321  double ret_array[1];
322 
323  MatGetValues(matrix, 1, row_as_array, 1, col_as_array, ret_array);
324 
325  return ret_array[0];
326 }
327 
328 void PetscMatTools::SetOption(Mat matrix, MatOption option)
329 {
330 #if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
331  MatSetOption(matrix, option, PETSC_TRUE);
332 #else
333  MatSetOption(matrix, option);
334 #endif
335 }
336 
338 {
339  /*
340  * We need to make sure that lhs_ith_row doesn't ignore off processor entries when assembling,
341  * otherwise the VecSetValues call a few lines below will not work as expected.
342  */
343 
344  PetscInt lo, hi;
345  PetscMatTools::GetOwnershipRange(matrix, lo, hi);
346  unsigned size = PetscMatTools::GetSize(matrix);
347 
348  Vec mat_ith_row = PetscTools::CreateVec(size, hi-lo, false);
349 
350  PetscInt num_entries;
351  const PetscInt* column_indices;
352  const PetscScalar* values;
353 
354  bool am_row_owner = (PetscInt)rowIndex >= lo && (PetscInt)rowIndex < hi;
355 
356  /*
357  * Am I the owner of the row? If so get the non-zero entries and add them lhs_ith_row.
358  * In parallel, VecAssembly{Begin,End} will send values to the rest of processors.
359  */
360  if (am_row_owner)
361  {
362  MatGetRow(matrix, rowIndex, &num_entries, &column_indices, &values);
363  VecSetValues(mat_ith_row, num_entries, column_indices, values, INSERT_VALUES);
364  }
365 
366  VecAssemblyBegin(mat_ith_row);
367  VecAssemblyEnd(mat_ith_row);
368 
369  if (am_row_owner)
370  {
371  MatRestoreRow(matrix, rowIndex, &num_entries, &column_indices, &values);
372  }
373 
374  return mat_ith_row;
375 }
376 
377 
378 bool PetscMatTools::CheckEquality(const Mat mat1, const Mat mat2, double tol)
379 {
380  Mat y;
381  MatDuplicate(mat2, MAT_COPY_VALUES, &y);
382 
383  double minus_one = -1.0;
384 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
385  // MatAYPX(*a, X, Y) does Y = X + a*Y.
386  MatAYPX(&minus_one, mat1, y);
387 #elif (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 1) //PETSc 2.3.1
388  // MatAYPX( Y, a, X) does Y = a*Y + X.
389  MatAYPX(y, minus_one, mat1);
390 #else
391  // MatAYPX( Y, a, X, structure) does Y = a*Y + X.
392  MatAYPX(y, minus_one, mat1, DIFFERENT_NONZERO_PATTERN);
393 #endif
394  PetscReal norm;
395  MatNorm(y, NORM_INFINITY, &norm);
397 
398  return (norm < tol);
399 }
400 
401 bool PetscMatTools::CheckSymmetry(const Mat matrix, double tol)
402 {
403  Mat trans;
404 #if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
405  MatTranspose(matrix, MAT_INITIAL_MATRIX, &trans);
406 #else
407  MatTranspose(matrix, &trans);
408 #endif
409  bool is_symmetric = PetscMatTools::CheckEquality(matrix, trans, tol);
410  PetscTools::Destroy(trans);
411  return is_symmetric;
412 }
413 
415 {
416 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 3) //PETSc 3.3 or later
417  MatSetOption(matrix, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
418 #endif
419 }
420 
static void TurnOffVariableAllocationError(Mat matrix)
#define PETSC_DESTROY_PARAM(x)
Definition: PetscTools.hpp:70
static Vec CreateVec(int size, int localSize=PETSC_DECIDE, bool ignoreOffProcEntries=true)
Definition: PetscTools.cpp:214
static void GetOwnershipRange(Mat matrix, PetscInt &lo, PetscInt &hi)
static void AddToElement(Mat matrix, PetscInt row, PetscInt col, double value)
static Vec GetMatrixRowDistributed(Mat matrix, unsigned rowIndex)
static void Zero(Mat matrix)
static void Display(Mat matrix)
static double GetElement(Mat matrix, PetscInt row, PetscInt col)
static bool CheckEquality(const Mat mat1, const Mat mat2, double tol=1e-10)
static void SetElement(Mat matrix, PetscInt row, PetscInt col, double value)
static void ZeroRowsWithValueOnDiagonal(Mat matrix, std::vector< unsigned > &rRows, double diagonalValue)
static void ZeroRowsAndColumnsWithValueOnDiagonal(Mat matrix, std::vector< unsigned > rowColIndices, double diagonalValue)
static void SwitchWriteMode(Mat matrix)
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:352
static bool CheckSymmetry(const Mat matrix, double tol=1e-10)
static void SetOption(Mat matrix, MatOption option)
static void Finalise(Mat matrix)
static unsigned GetSize(Mat matrix)
static void SetRow(Mat matrix, PetscInt row, double value)
static void ZeroColumn(Mat matrix, PetscInt col)