Reference documentation for deal.II version 8.4.1
petsc_sparse_matrix.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2015 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/lac/petsc_sparse_matrix.h>
17 
18 #ifdef DEAL_II_WITH_PETSC
19 
20 # include <deal.II/lac/petsc_vector.h>
21 # include <deal.II/lac/sparsity_pattern.h>
22 # include <deal.II/lac/dynamic_sparsity_pattern.h>
23 
24 DEAL_II_NAMESPACE_OPEN
25 
26 namespace PETScWrappers
27 {
28 
30  {
31  const int m=0, n=0, n_nonzero_per_row=0;
32  const int ierr
33  = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, n_nonzero_per_row,
34  0, &matrix);
35  AssertThrow (ierr == 0, ExcPETScError(ierr));
36  }
37 
38 
39 
41  const size_type n,
42  const size_type n_nonzero_per_row,
43  const bool is_symmetric)
44  {
45  do_reinit (m, n, n_nonzero_per_row, is_symmetric);
46  }
47 
48 
49 
51  const size_type n,
52  const std::vector<size_type> &row_lengths,
53  const bool is_symmetric)
54  {
55  do_reinit (m, n, row_lengths, is_symmetric);
56  }
57 
58 
59 
60  template <typename SparsityPatternType>
62  SparseMatrix (const SparsityPatternType &sparsity_pattern,
63  const bool preset_nonzero_locations)
64  {
65  do_reinit (sparsity_pattern, preset_nonzero_locations);
66  }
67 
68 
69 
70  SparseMatrix &
71  SparseMatrix::operator = (const double d)
72  {
74  return *this;
75  }
76 
77 
78 
79  void
81  const size_type n,
82  const size_type n_nonzero_per_row,
83  const bool is_symmetric)
84  {
85  // get rid of old matrix and generate a
86  // new one
87 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
88  const int ierr = MatDestroy (matrix);
89 #else
90  const int ierr = MatDestroy (&matrix);
91 #endif
92  AssertThrow (ierr == 0, ExcPETScError(ierr));
93 
94  do_reinit (m, n, n_nonzero_per_row, is_symmetric);
95  }
96 
97 
98 
99  void
101  const size_type n,
102  const std::vector<size_type> &row_lengths,
103  const bool is_symmetric)
104  {
105  // get rid of old matrix and generate a
106  // new one
107 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
108  const int ierr = MatDestroy (matrix);
109 #else
110  const int ierr = MatDestroy (&matrix);
111 #endif
112  AssertThrow (ierr == 0, ExcPETScError(ierr));
113 
114  do_reinit (m, n, row_lengths, is_symmetric);
115  }
116 
117 
118 
119  template <typename SparsityPatternType>
120  void
122  reinit (const SparsityPatternType &sparsity_pattern,
123  const bool preset_nonzero_locations)
124  {
125  // get rid of old matrix and generate a
126  // new one
127 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
128  const int ierr = MatDestroy (matrix);
129 #else
130  const int ierr = MatDestroy (&matrix);
131 #endif
132  AssertThrow (ierr == 0, ExcPETScError(ierr));
133 
134  do_reinit (sparsity_pattern, preset_nonzero_locations);
135  }
136 
137 
138 
139  const MPI_Comm &
141  {
142  static MPI_Comm comm;
143  PetscObjectGetComm((PetscObject)matrix, &comm);
144  return comm;
145  }
146 
147 
148 
149  void
151  const size_type n,
152  const size_type n_nonzero_per_row,
153  const bool is_symmetric)
154  {
155  // use the call sequence indicating only
156  // a maximal number of elements per row
157  // for all rows globally
158  const int ierr
159  = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, n_nonzero_per_row,
160  0, &matrix);
161  AssertThrow (ierr == 0, ExcPETScError(ierr));
162 
163  // set symmetric flag, if so requested
164  if (is_symmetric == true)
165  {
166 #if DEAL_II_PETSC_VERSION_LT(3,0,0)
167  const int ierr
168  = MatSetOption (matrix, MAT_SYMMETRIC);
169 #else
170  const int ierr
171  = MatSetOption (matrix, MAT_SYMMETRIC, PETSC_TRUE);
172 #endif
173 
174  AssertThrow (ierr == 0, ExcPETScError(ierr));
175  }
176  }
177 
178 
179 
180  void
182  const size_type n,
183  const std::vector<size_type> &row_lengths,
184  const bool is_symmetric)
185  {
186  Assert (row_lengths.size() == m,
187  ExcDimensionMismatch (row_lengths.size(), m));
188 
189  // use the call sequence indicating a
190  // maximal number of elements for each
191  // row individually. annoyingly, we
192  // always use unsigned ints for cases
193  // like this, while PETSc wants to see
194  // signed integers. so we have to
195  // convert, unless we want to play dirty
196  // tricks with conversions of pointers
197  const std::vector<PetscInt>
198  int_row_lengths (row_lengths.begin(), row_lengths.end());
199 
200  const int ierr
201  = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0,
202  &int_row_lengths[0], &matrix);
203  AssertThrow (ierr == 0, ExcPETScError(ierr));
204 
205  // set symmetric flag, if so requested
206  if (is_symmetric == true)
207  {
208 #if DEAL_II_PETSC_VERSION_LT(3,0,0)
209  const int ierr
210  = MatSetOption (matrix, MAT_SYMMETRIC);
211 #else
212  const int ierr
213  = MatSetOption (matrix, MAT_SYMMETRIC, PETSC_TRUE);
214 #endif
215 
216  AssertThrow (ierr == 0, ExcPETScError(ierr));
217  }
218  }
219 
220 
221 
222  template <typename SparsityPatternType>
223  void
224  SparseMatrix::do_reinit (const SparsityPatternType &sparsity_pattern,
225  const bool preset_nonzero_locations)
226  {
227  std::vector<size_type> row_lengths (sparsity_pattern.n_rows());
228  for (size_type i=0; i<sparsity_pattern.n_rows(); ++i)
229  row_lengths[i] = sparsity_pattern.row_length (i);
230 
231  do_reinit (sparsity_pattern.n_rows(),
232  sparsity_pattern.n_cols(),
233  row_lengths, false);
234 
235  // next preset the exact given matrix
236  // entries with zeros, if the user
237  // requested so. this doesn't avoid any
238  // memory allocations, but it at least
239  // avoids some searches later on. the
240  // key here is that we can use the
241  // matrix set routines that set an
242  // entire row at once, not a single
243  // entry at a time
244  //
245  // for the usefulness of this option
246  // read the documentation of this
247  // class.
248  if (preset_nonzero_locations == true)
249  {
250  std::vector<PetscInt> row_entries;
251  std::vector<PetscScalar> row_values;
252  for (size_type i=0; i<sparsity_pattern.n_rows(); ++i)
253  {
254  row_entries.resize (row_lengths[i]);
255  row_values.resize (row_lengths[i], 0.0);
256  for (size_type j=0; j<row_lengths[i]; ++j)
257  row_entries[j] = sparsity_pattern.column_number (i,j);
258 
259  const PetscInt int_row = i;
260  MatSetValues (matrix, 1, &int_row,
261  row_lengths[i], &row_entries[0],
262  &row_values[0], INSERT_VALUES);
263  }
264  compress (VectorOperation::insert);
265 
266 
267  // Tell PETSc that we are not
268  // planning on adding new entries
269  // to the matrix. Generate errors
270  // in debug mode.
271  int ierr;
272 #if DEAL_II_PETSC_VERSION_LT(3,0,0)
273 #ifdef DEBUG
274  ierr = MatSetOption (matrix, MAT_NEW_NONZERO_LOCATION_ERR);
275  AssertThrow (ierr == 0, ExcPETScError(ierr));
276 #else
277  ierr = MatSetOption (matrix, MAT_NO_NEW_NONZERO_LOCATIONS);
278  AssertThrow (ierr == 0, ExcPETScError(ierr));
279 #endif
280 #else
281 #ifdef DEBUG
282  ierr = MatSetOption (matrix, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
283  AssertThrow (ierr == 0, ExcPETScError(ierr));
284 #else
285  ierr = MatSetOption (matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
286  AssertThrow (ierr == 0, ExcPETScError(ierr));
287 #endif
288 #endif
289 
290  // Tell PETSc to keep the
291  // SparsityPattern entries even if
292  // we delete a row with
293  // clear_rows() which calls
294  // MatZeroRows(). Otherwise one can
295  // not write into that row
296  // afterwards.
297 #if DEAL_II_PETSC_VERSION_LT(3,0,0)
298  ierr = MatSetOption (matrix, MAT_KEEP_ZEROED_ROWS);
299  AssertThrow (ierr == 0, ExcPETScError(ierr));
300 #elif DEAL_II_PETSC_VERSION_LT(3,1,0)
301  ierr = MatSetOption (matrix, MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
302  AssertThrow (ierr == 0, ExcPETScError(ierr));
303 #else
304  ierr = MatSetOption (matrix, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
305  AssertThrow (ierr == 0, ExcPETScError(ierr));
306 #endif
307 
308  }
309  }
310 
311 
312  // Explicit instantiations
313  //
314  template
316  const bool);
317  template
319  const bool);
320 
321  template void
323  const bool);
324  template void
326  const bool);
327 
328  template void
330  const bool);
331  template void
333  const bool);
334 
335  PetscScalar
337  {
338  Vector tmp (v.size());
339  vmult (tmp, v);
340  return tmp*v;
341  }
342 
343  PetscScalar
345  const VectorBase &v) const
346  {
347  Vector tmp (v.size());
348  vmult (tmp, v);
349  return u*tmp;
350  }
351 }
352 
353 
354 DEAL_II_NAMESPACE_CLOSE
355 
356 #endif // DEAL_II_WITH_PETSC
#define AssertThrow(cond, exc)
Definition: exceptions.h:358
PetscScalar matrix_scalar_product(const VectorBase &u, const VectorBase &v) const
void compress(const VectorOperation::values operation)
SparseMatrix & operator=(const double d)
void reinit(const size_type m, const size_type n, const size_type n_nonzero_per_row, const bool is_symmetric=false)
#define Assert(cond, exc)
Definition: exceptions.h:294
MatrixBase & operator=(const value_type d)
types::global_dof_index size_type
virtual const MPI_Comm & get_mpi_communicator() const
void vmult(VectorBase &dst, const VectorBase &src) const
void do_reinit(const size_type m, const size_type n, const size_type n_nonzero_per_row, const bool is_symmetric=false)
PetscScalar matrix_norm_square(const VectorBase &v) const