Reference documentation for deal.II version GIT 2f5445400b 2023-02-05 22:30:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
petsc_parallel_sparse_matrix.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2021 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.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 
18 #ifdef DEAL_II_WITH_PETSC
19 
20 # include <deal.II/base/mpi.h>
21 
23 # include <deal.II/lac/exceptions.h>
27 
29 
30 namespace PETScWrappers
31 {
32  namespace MPI
33  {
35  {
36  // just like for vectors: since we
37  // create an empty matrix, we can as
38  // well make it sequential
39  const int m = 0, n = 0, n_nonzero_per_row = 0;
40  const PetscErrorCode ierr = MatCreateSeqAIJ(
41  PETSC_COMM_SELF, m, n, n_nonzero_per_row, nullptr, &matrix);
42  AssertThrow(ierr == 0, ExcPETScError(ierr));
43  }
44 
46  : MatrixBase(A)
47  {}
48 
50  {
52  }
53 
54 
55 
56  template <typename SparsityPatternType>
58  const MPI_Comm & communicator,
59  const SparsityPatternType & sparsity_pattern,
60  const std::vector<size_type> &local_rows_per_process,
61  const std::vector<size_type> &local_columns_per_process,
62  const unsigned int this_process,
63  const bool preset_nonzero_locations)
64  {
65  do_reinit(communicator,
66  sparsity_pattern,
67  local_rows_per_process,
68  local_columns_per_process,
69  this_process,
70  preset_nonzero_locations);
71  }
72 
73 
74 
75  void
77  {
78  if (&other == this)
79  return;
80 
81  PetscErrorCode ierr = destroy_matrix(matrix);
82  AssertThrow(ierr == 0, ExcPETScError(ierr));
83 
84  ierr = MatDuplicate(other.matrix, MAT_DO_NOT_COPY_VALUES, &matrix);
85  AssertThrow(ierr == 0, ExcPETScError(ierr));
86  }
87 
88  template <typename SparsityPatternType>
89  void
90  SparseMatrix::reinit(const IndexSet & local_rows,
91  const IndexSet & local_active_rows,
92  const IndexSet & local_columns,
93  const IndexSet & local_active_columns,
94  const SparsityPatternType &sparsity_pattern,
95  const MPI_Comm & communicator)
96  {
97  // get rid of old matrix and generate a new one
98  const PetscErrorCode ierr = destroy_matrix(matrix);
99  AssertThrow(ierr == 0, ExcPETScError(ierr));
100 
101  do_reinit(communicator,
102  local_rows,
103  local_active_rows,
104  local_columns,
105  local_active_columns,
106  sparsity_pattern);
107  }
108 
109 
110  SparseMatrix &
112  {
114  return *this;
115  }
116 
117  void
119  {
120  if (&other == this)
121  return;
122 
123  const PetscErrorCode ierr =
124  MatCopy(other.matrix, matrix, SAME_NONZERO_PATTERN);
125  AssertThrow(ierr == 0, ExcPETScError(ierr));
126  }
127 
128 
129 
130  template <typename SparsityPatternType>
131  void
133  const MPI_Comm & communicator,
134  const SparsityPatternType & sparsity_pattern,
135  const std::vector<size_type> &local_rows_per_process,
136  const std::vector<size_type> &local_columns_per_process,
137  const unsigned int this_process,
138  const bool preset_nonzero_locations)
139  {
140  // get rid of old matrix and generate a new one
141  const PetscErrorCode ierr = destroy_matrix(matrix);
142  AssertThrow(ierr == 0, ExcPETScError(ierr));
143 
144 
145  do_reinit(communicator,
146  sparsity_pattern,
147  local_rows_per_process,
148  local_columns_per_process,
149  this_process,
150  preset_nonzero_locations);
151  }
152 
153 
154 
155  template <typename SparsityPatternType>
156  void
157  SparseMatrix::reinit(const IndexSet & local_rows,
158  const SparsityPatternType &sparsity_pattern,
159  const MPI_Comm & communicator)
160  {
161  do_reinit(communicator, local_rows, local_rows, sparsity_pattern);
162  }
163 
164  template <typename SparsityPatternType>
165  void
166  SparseMatrix::reinit(const IndexSet & local_rows,
167  const IndexSet & local_columns,
168  const SparsityPatternType &sparsity_pattern,
169  const MPI_Comm & communicator)
170  {
171  // get rid of old matrix and generate a new one
172  const PetscErrorCode ierr = destroy_matrix(matrix);
173  AssertThrow(ierr == 0, ExcPETScError(ierr));
174 
175  do_reinit(communicator, local_rows, local_columns, sparsity_pattern);
176  }
177 
178 
179 
180  template <typename SparsityPatternType>
181  void
182  SparseMatrix::do_reinit(const MPI_Comm & communicator,
183  const IndexSet & local_rows,
184  const IndexSet & local_columns,
185  const SparsityPatternType &sparsity_pattern)
186  {
187  Assert(sparsity_pattern.n_rows() == local_rows.size(),
188  ExcMessage(
189  "SparsityPattern and IndexSet have different number of rows"));
190  Assert(
191  sparsity_pattern.n_cols() == local_columns.size(),
192  ExcMessage(
193  "SparsityPattern and IndexSet have different number of columns"));
194  Assert(local_rows.is_contiguous() && local_columns.is_contiguous(),
195  ExcMessage("PETSc only supports contiguous row/column ranges"));
196  Assert(local_rows.is_ascending_and_one_to_one(communicator),
198 
199 # ifdef DEBUG
200  {
201  // check indexsets
202  types::global_dof_index row_owners =
203  Utilities::MPI::sum(local_rows.n_elements(), communicator);
204  types::global_dof_index col_owners =
205  Utilities::MPI::sum(local_columns.n_elements(), communicator);
206  Assert(row_owners == sparsity_pattern.n_rows(),
207  ExcMessage(
208  std::string(
209  "Each row has to be owned by exactly one owner (n_rows()=") +
210  std::to_string(sparsity_pattern.n_rows()) +
211  " but sum(local_rows.n_elements())=" +
212  std::to_string(row_owners) + ")"));
213  Assert(
214  col_owners == sparsity_pattern.n_cols(),
215  ExcMessage(
216  std::string(
217  "Each column has to be owned by exactly one owner (n_cols()=") +
218  std::to_string(sparsity_pattern.n_cols()) +
219  " but sum(local_columns.n_elements())=" +
220  std::to_string(col_owners) + ")"));
221  }
222 # endif
223 
224 
225  // create the matrix. We do not set row length but set the
226  // correct SparsityPattern later.
227  PetscErrorCode ierr = MatCreate(communicator, &matrix);
228  AssertThrow(ierr == 0, ExcPETScError(ierr));
229 
230  ierr = MatSetSizes(matrix,
231  local_rows.n_elements(),
232  local_columns.n_elements(),
233  sparsity_pattern.n_rows(),
234  sparsity_pattern.n_cols());
235  AssertThrow(ierr == 0, ExcPETScError(ierr));
236 
237  // Use MATAIJ which dispatches to SEQAIJ
238  // if the size of the communicator is 1,
239  // and to MPIAIJ otherwise.
240  ierr = MatSetType(matrix, MATAIJ);
241  AssertThrow(ierr == 0, ExcPETScError(ierr));
242 
243 
244  // next preset the exact given matrix
245  // entries with zeros. this doesn't avoid any
246  // memory allocations, but it at least
247  // avoids some searches later on. the
248  // key here is that we can use the
249  // matrix set routines that set an
250  // entire row at once, not a single
251  // entry at a time
252  //
253  // for the usefulness of this option
254  // read the documentation of this
255  // class.
256  // if (preset_nonzero_locations == true)
257  if (local_rows.n_elements() > 0)
258  {
259  // MatXXXAIJSetPreallocationCSR
260  // can be used to allocate the sparsity
261  // pattern of a matrix
262 
263  const PetscInt local_row_start = local_rows.nth_index_in_set(0);
264  const PetscInt local_row_end =
265  local_row_start + local_rows.n_elements();
266 
267 
268  // first set up the column number
269  // array for the rows to be stored
270  // on the local processor. have one
271  // dummy entry at the end to make
272  // sure petsc doesn't read past the
273  // end
274  std::vector<PetscInt>
275 
276  rowstart_in_window(local_row_end - local_row_start + 1, 0),
277  colnums_in_window;
278  {
279  unsigned int n_cols = 0;
280  for (PetscInt i = local_row_start; i < local_row_end; ++i)
281  {
282  const PetscInt row_length = sparsity_pattern.row_length(i);
283  rowstart_in_window[i + 1 - local_row_start] =
284  rowstart_in_window[i - local_row_start] + row_length;
285  n_cols += row_length;
286  }
287  colnums_in_window.resize(n_cols + 1, -1);
288  }
289 
290  // now copy over the information
291  // from the sparsity pattern.
292  {
293  PetscInt *ptr = colnums_in_window.data();
294  for (PetscInt i = local_row_start; i < local_row_end; ++i)
295  for (typename SparsityPatternType::iterator p =
296  sparsity_pattern.begin(i);
297  p != sparsity_pattern.end(i);
298  ++p, ++ptr)
299  *ptr = p->column();
300  }
301 
302 
303  // then call the petsc functions
304  // that summarily allocates these
305  // entries.
306  // Here we both call the specific API since this is how
307  // PETSc polymorphism works. If the matrix is of type MPIAIJ,
308  // the second call is dummy. If the matrix is of type SEQAIJ,
309  // the first call is dummy.
310  ierr = MatMPIAIJSetPreallocationCSR(matrix,
311  rowstart_in_window.data(),
312  colnums_in_window.data(),
313  nullptr);
314  ierr = MatSeqAIJSetPreallocationCSR(matrix,
315  rowstart_in_window.data(),
316  colnums_in_window.data(),
317  nullptr);
318  AssertThrow(ierr == 0, ExcPETScError(ierr));
319  }
320  else
321  {
322  PetscInt i = 0;
323 
324  ierr = MatSeqAIJSetPreallocationCSR(matrix, &i, &i, nullptr);
325  AssertThrow(ierr == 0, ExcPETScError(ierr));
326  ierr = MatMPIAIJSetPreallocationCSR(matrix, &i, &i, nullptr);
327  AssertThrow(ierr == 0, ExcPETScError(ierr));
328  }
330 
331  {
334  }
335  }
336 
337 
338  template <typename SparsityPatternType>
339  void
341  const MPI_Comm & communicator,
342  const SparsityPatternType & sparsity_pattern,
343  const std::vector<size_type> &local_rows_per_process,
344  const std::vector<size_type> &local_columns_per_process,
345  const unsigned int this_process,
346  const bool preset_nonzero_locations)
347  {
348  Assert(local_rows_per_process.size() == local_columns_per_process.size(),
349  ExcDimensionMismatch(local_rows_per_process.size(),
350  local_columns_per_process.size()));
351  Assert(this_process < local_rows_per_process.size(), ExcInternalError());
353 
354  // for each row that we own locally, we
355  // have to count how many of the
356  // entries in the sparsity pattern lie
357  // in the column area we have locally,
358  // and how many aren't. for this, we
359  // first have to know which areas are
360  // ours
361  size_type local_row_start = 0;
362  for (unsigned int p = 0; p < this_process; ++p)
363  local_row_start += local_rows_per_process[p];
364  const size_type local_row_end =
365  local_row_start + local_rows_per_process[this_process];
366 
367  // create the matrix. We
368  // do not set row length but set the
369  // correct SparsityPattern later.
370  PetscErrorCode ierr = MatCreate(communicator, &matrix);
371  AssertThrow(ierr == 0, ExcPETScError(ierr));
372 
373  ierr = MatSetSizes(matrix,
374  local_rows_per_process[this_process],
375  local_columns_per_process[this_process],
376  sparsity_pattern.n_rows(),
377  sparsity_pattern.n_cols());
378  AssertThrow(ierr == 0, ExcPETScError(ierr));
379 
380  // Use MATAIJ which dispatches to SEQAIJ
381  // if the size of the communicator is 1,
382  // and to MPIAIJ otherwise.
383  ierr = MatSetType(matrix, MATAIJ);
384  AssertThrow(ierr == 0, ExcPETScError(ierr));
385 
386  // next preset the exact given matrix
387  // entries with zeros, if the user
388  // requested so. this doesn't avoid any
389  // memory allocations, but it at least
390  // avoids some searches later on. the
391  // key here is that we can use the
392  // matrix set routines that set an
393  // entire row at once, not a single
394  // entry at a time
395  //
396  // for the usefulness of this option
397  // read the documentation of this
398  // class.
399  if (preset_nonzero_locations == true)
400  {
401  // MatXXXAIJSetPreallocationCSR
402  // can be used to allocate the sparsity
403  // pattern of a matrix if it is already
404  // available:
405 
406  // first set up the column number
407  // array for the rows to be stored
408  // on the local processor. have one
409  // dummy entry at the end to make
410  // sure petsc doesn't read past the
411  // end
412  std::vector<PetscInt>
413 
414  rowstart_in_window(local_row_end - local_row_start + 1, 0),
415  colnums_in_window;
416  {
417  size_type n_cols = 0;
418  for (size_type i = local_row_start; i < local_row_end; ++i)
419  {
420  const size_type row_length = sparsity_pattern.row_length(i);
421  rowstart_in_window[i + 1 - local_row_start] =
422  rowstart_in_window[i - local_row_start] + row_length;
423  n_cols += row_length;
424  }
425  colnums_in_window.resize(n_cols + 1, -1);
426  }
427 
428  // now copy over the information
429  // from the sparsity pattern.
430  {
431  PetscInt *ptr = colnums_in_window.data();
432  for (size_type i = local_row_start; i < local_row_end; ++i)
433  for (typename SparsityPatternType::iterator p =
434  sparsity_pattern.begin(i);
435  p != sparsity_pattern.end(i);
436  ++p, ++ptr)
437  *ptr = p->column();
438  }
439 
440 
441  // then call the petsc function
442  // that summarily allocates these
443  // entries.
444  // Here we both call the specific API since this is how
445  // PETSc polymorphism works. If the matrix is of type MPIAIJ,
446  // the second call is dummy. If the matrix is of type SEQAIJ,
447  // the first call is dummy.
448  ierr = MatSeqAIJSetPreallocationCSR(matrix,
449  rowstart_in_window.data(),
450  colnums_in_window.data(),
451  nullptr);
452  ierr = MatMPIAIJSetPreallocationCSR(matrix,
453  rowstart_in_window.data(),
454  colnums_in_window.data(),
455  nullptr);
456  AssertThrow(ierr == 0, ExcPETScError(ierr));
457 
460  }
461  }
462 
463  // BDDC
464  template <typename SparsityPatternType>
465  void
466  SparseMatrix::do_reinit(const MPI_Comm & communicator,
467  const IndexSet & local_rows,
468  const IndexSet & local_active_rows,
469  const IndexSet & local_columns,
470  const IndexSet & local_active_columns,
471  const SparsityPatternType &sparsity_pattern)
472  {
473 # if DEAL_II_PETSC_VERSION_GTE(3, 10, 0)
474  Assert(sparsity_pattern.n_rows() == local_rows.size(),
475  ExcMessage(
476  "SparsityPattern and IndexSet have different number of rows."));
477  Assert(
478  sparsity_pattern.n_cols() == local_columns.size(),
479  ExcMessage(
480  "SparsityPattern and IndexSet have different number of columns"));
481  Assert(local_rows.is_contiguous() && local_columns.is_contiguous(),
482  ExcMessage("PETSc only supports contiguous row/column ranges"));
483  Assert(local_rows.is_ascending_and_one_to_one(communicator),
485 
486 # ifdef DEBUG
487  {
488  // check indexsets
489  types::global_dof_index row_owners =
490  Utilities::MPI::sum(local_rows.n_elements(), communicator);
491  types::global_dof_index col_owners =
492  Utilities::MPI::sum(local_columns.n_elements(), communicator);
493  Assert(row_owners == sparsity_pattern.n_rows(),
494  ExcMessage(
495  std::string(
496  "Each row has to be owned by exactly one owner (n_rows()=") +
497  std::to_string(sparsity_pattern.n_rows()) +
498  " but sum(local_rows.n_elements())=" +
499  std::to_string(row_owners) + ")"));
500  Assert(
501  col_owners == sparsity_pattern.n_cols(),
502  ExcMessage(
503  std::string(
504  "Each column has to be owned by exactly one owner (n_cols()=") +
505  std::to_string(sparsity_pattern.n_cols()) +
506  " but sum(local_columns.n_elements())=" +
507  std::to_string(col_owners) + ")"));
508  }
509 # endif
510  PetscErrorCode ierr;
511 
512  // create the local to global mappings as arrays.
513  IndexSet::size_type n_local_active_rows = local_active_rows.n_elements();
514  IndexSet::size_type n_local_active_cols =
515  local_active_columns.n_elements();
516  std::vector<PetscInt> idx_glob_row(n_local_active_rows);
517  std::vector<PetscInt> idx_glob_col(n_local_active_cols);
518  for (IndexSet::size_type k = 0; k < n_local_active_rows; ++k)
519  {
520  idx_glob_row[k] = local_active_rows.nth_index_in_set(k);
521  }
522  for (IndexSet::size_type k = 0; k < n_local_active_cols; ++k)
523  {
524  idx_glob_col[k] = local_active_columns.nth_index_in_set(k);
525  }
526 
527 
528  IS is_glob_row, is_glob_col;
529  // Create row index set
530  ISLocalToGlobalMapping l2gmap_row;
531  ierr = ISCreateGeneral(communicator,
532  n_local_active_rows,
533  idx_glob_row.data(),
534  PETSC_COPY_VALUES,
535  &is_glob_row);
536  AssertThrow(ierr == 0, ExcPETScError(ierr));
537  ierr = ISLocalToGlobalMappingCreateIS(is_glob_row, &l2gmap_row);
538  AssertThrow(ierr == 0, ExcPETScError(ierr));
539  ierr = ISDestroy(&is_glob_row);
540  AssertThrow(ierr == 0, ExcPETScError(ierr));
541  ierr =
542  ISLocalToGlobalMappingViewFromOptions(l2gmap_row, nullptr, "-view_map");
543  AssertThrow(ierr == 0, ExcPETScError(ierr));
544 
545  // Create column index set
546  ISLocalToGlobalMapping l2gmap_col;
547  ierr = ISCreateGeneral(communicator,
548  n_local_active_cols,
549  idx_glob_col.data(),
550  PETSC_COPY_VALUES,
551  &is_glob_col);
552  AssertThrow(ierr == 0, ExcPETScError(ierr));
553  ierr = ISLocalToGlobalMappingCreateIS(is_glob_col, &l2gmap_col);
554  AssertThrow(ierr == 0, ExcPETScError(ierr));
555  ierr = ISDestroy(&is_glob_col);
556  AssertThrow(ierr == 0, ExcPETScError(ierr));
557  ierr =
558  ISLocalToGlobalMappingViewFromOptions(l2gmap_col, nullptr, "-view_map");
559  AssertThrow(ierr == 0, ExcPETScError(ierr));
560 
561  // create the matrix with the IS constructor.
562  ierr = MatCreateIS(communicator,
563  1,
564  local_rows.n_elements(),
565  local_columns.n_elements(),
566  sparsity_pattern.n_rows(),
567  sparsity_pattern.n_cols(),
568  l2gmap_row,
569  l2gmap_col,
570  &matrix);
571  AssertThrow(ierr == 0, ExcPETScError(ierr));
572  ierr = ISLocalToGlobalMappingDestroy(&l2gmap_row);
573  AssertThrow(ierr == 0, ExcPETScError(ierr));
574  ierr = ISLocalToGlobalMappingDestroy(&l2gmap_col);
575  AssertThrow(ierr == 0, ExcPETScError(ierr));
576 
577  // next preset the exact given matrix
578  // entries with zeros. This doesn't avoid any
579  // memory allocations, but it at least
580  // avoids some searches later on. the
581  // key here is that we can use the
582  // matrix set routines that set an
583  // entire row at once, not a single
584  // entry at a time.
585  //
586  // for the usefulness of this option
587  // read the documentation of this
588  // class.
589 
590  Mat local_matrix; // In the MATIS case, we use the local matrix instead
591  ierr = MatISGetLocalMat(matrix, &local_matrix);
592  AssertThrow(ierr == 0, ExcPETScError(ierr));
593  ierr = MatSetType(local_matrix,
594  MATSEQAIJ); // SEQ as it is local! TODO: Allow for
595  // OpenMP parallelization in local node.
596  AssertThrow(ierr == 0, ExcPETScError(ierr));
597  if (local_rows.n_elements() > 0)
598  {
599  // MatSEQAIJSetPreallocationCSR
600  // can be used to allocate the sparsity
601  // pattern of a matrix. Local matrices start from 0 (MATIS).
602  const PetscInt local_row_start = 0;
603  const PetscInt local_row_end = local_active_rows.n_elements();
604 
605  // first set up the column number
606  // array for the rows to be stored
607  // on the local processor.
608  std::vector<PetscInt> rowstart_in_window(local_row_end -
609  local_row_start + 1,
610  0),
611  colnums_in_window;
612  unsigned int global_row_index = 0;
613  {
614  unsigned int n_cols = 0;
615  unsigned int global_row_index = 0;
616  for (PetscInt i = local_row_start; i < local_row_end; ++i)
617  {
618  global_row_index = local_active_rows.nth_index_in_set(i);
619  const PetscInt row_length =
620  sparsity_pattern.row_length(global_row_index);
621  rowstart_in_window[i + 1 - local_row_start] =
622  rowstart_in_window[i - local_row_start] + row_length;
623  n_cols += row_length;
624  }
625  colnums_in_window.resize(n_cols + 1, -1);
626  }
627 
628 
629  // now copy over the information
630  // from the sparsity pattern. For this we first invert the column
631  // index set.
632  std::map<unsigned int, unsigned int> loc_act_cols_inv;
633  for (unsigned int i = 0; i < local_active_columns.n_elements(); ++i)
634  {
635  loc_act_cols_inv[local_active_columns.nth_index_in_set(i)] = i;
636  }
637 
638  {
639  PetscInt *ptr = colnums_in_window.data();
640  for (PetscInt i = local_row_start; i < local_row_end; ++i)
641  {
642  global_row_index = local_active_rows.nth_index_in_set(i);
643  for (typename SparsityPatternType::iterator p =
644  sparsity_pattern.begin(global_row_index);
645  p != sparsity_pattern.end(global_row_index);
646  ++p, ++ptr)
647  *ptr = loc_act_cols_inv[p->column()];
648  }
649  }
650 
651  // then call the petsc function
652  // that summarily allocates these
653  // entries:
654  ierr = MatSeqAIJSetPreallocationCSR(local_matrix,
655  rowstart_in_window.data(),
656  colnums_in_window.data(),
657  nullptr);
658  AssertThrow(ierr == 0, ExcPETScError(ierr));
659  }
660  else
661  {
662  PetscInt i = 0;
663  ierr = MatSeqAIJSetPreallocationCSR(local_matrix, &i, &i, nullptr);
664  AssertThrow(ierr == 0, ExcPETScError(ierr));
665  }
667 
668  {
669  close_matrix(local_matrix);
670  set_keep_zero_rows(local_matrix);
671  }
672  ierr = MatISRestoreLocalMat(matrix, &local_matrix);
673  AssertThrow(ierr == 0, ExcPETScError(ierr));
674 # else
675  {
676  // Use this to avoid unused variables warning
677  (void)local_rows;
678  (void)local_active_rows;
679  (void)local_columns;
680  (void)local_active_columns;
681  (void)sparsity_pattern;
682  AssertThrow(false,
683  ExcMessage(
684  "BDDC preconditioner requires PETSc 3.10.0 or newer"));
685  }
686 # endif
687  }
688 
689 # ifndef DOXYGEN
690  // explicit instantiations
691  //
692  template SparseMatrix::SparseMatrix(const MPI_Comm &,
693  const SparsityPattern &,
694  const std::vector<size_type> &,
695  const std::vector<size_type> &,
696  const unsigned int,
697  const bool);
698  template SparseMatrix::SparseMatrix(const MPI_Comm &,
699  const DynamicSparsityPattern &,
700  const std::vector<size_type> &,
701  const std::vector<size_type> &,
702  const unsigned int,
703  const bool);
704 
705  template void
706  SparseMatrix::reinit(const MPI_Comm &,
707  const SparsityPattern &,
708  const std::vector<size_type> &,
709  const std::vector<size_type> &,
710  const unsigned int,
711  const bool);
712  template void
713  SparseMatrix::reinit(const MPI_Comm &,
714  const DynamicSparsityPattern &,
715  const std::vector<size_type> &,
716  const std::vector<size_type> &,
717  const unsigned int,
718  const bool);
719 
720  template void
722  const SparsityPattern &,
723  const MPI_Comm &);
724 
725  template void
727  const IndexSet &,
728  const SparsityPattern &,
729  const MPI_Comm &);
730 
731  template void
733  const DynamicSparsityPattern &,
734  const MPI_Comm &);
735 
736  template void
738  const IndexSet &,
739  const DynamicSparsityPattern &,
740  const MPI_Comm &);
741 
742  template void
743  SparseMatrix::do_reinit(const MPI_Comm &,
744  const SparsityPattern &,
745  const std::vector<size_type> &,
746  const std::vector<size_type> &,
747  const unsigned int,
748  const bool);
749  template void
750  SparseMatrix::do_reinit(const MPI_Comm &,
751  const DynamicSparsityPattern &,
752  const std::vector<size_type> &,
753  const std::vector<size_type> &,
754  const unsigned int,
755  const bool);
756 
757  template void
758  SparseMatrix::do_reinit(const MPI_Comm &,
759  const IndexSet &,
760  const IndexSet &,
761  const SparsityPattern &);
762 
763  template void
764  SparseMatrix::do_reinit(const MPI_Comm &,
765  const IndexSet &,
766  const IndexSet &,
767  const DynamicSparsityPattern &);
768 
769  template void
771  const IndexSet &,
772  const IndexSet &,
773  const IndexSet &,
774  const SparsityPattern &,
775  const MPI_Comm &);
776  template void
778  const IndexSet &,
779  const IndexSet &,
780  const IndexSet &,
781  const DynamicSparsityPattern &,
782  const MPI_Comm &);
783 
784  template void
785  SparseMatrix::do_reinit(const MPI_Comm &,
786  const IndexSet &,
787  const IndexSet &,
788  const IndexSet &,
789  const IndexSet &,
790  const SparsityPattern &);
791  template void
792  SparseMatrix::do_reinit(const MPI_Comm &,
793  const IndexSet &,
794  const IndexSet &,
795  const IndexSet &,
796  const IndexSet &,
797  const DynamicSparsityPattern &);
798 # endif
799 
800 
801  PetscScalar
803  {
804  Vector tmp(v);
805  vmult(tmp, v);
806  // note, that v*tmp returns sum_i conjugate(v)_i * tmp_i
807  return v * tmp;
808  }
809 
810  PetscScalar
812  {
813  Vector tmp(v);
814  vmult(tmp, v);
815  // note, that v*tmp returns sum_i conjugate(v)_i * tmp_i
816  return u * tmp;
817  }
818 
819  IndexSet
821  {
822  PetscInt n_rows, n_cols, n_loc_rows, n_loc_cols, min, max;
823  PetscErrorCode ierr;
824 
825  ierr = MatGetSize(matrix, &n_rows, &n_cols);
826  AssertThrow(ierr == 0, ExcPETScError(ierr));
827 
828  ierr = MatGetLocalSize(matrix, &n_loc_rows, &n_loc_cols);
829  AssertThrow(ierr == 0, ExcPETScError(ierr));
830 
831  ierr = MatGetOwnershipRangeColumn(matrix, &min, &max);
832  AssertThrow(ierr == 0, ExcPETScError(ierr));
833 
834  Assert(n_loc_cols == max - min,
835  ExcMessage(
836  "PETSc is requiring non contiguous memory allocation."));
837 
838  IndexSet indices(n_cols);
839  indices.add_range(min, max);
840  indices.compress();
841 
842  return indices;
843  }
844 
845  IndexSet
847  {
848  PetscInt n_rows, n_cols, n_loc_rows, n_loc_cols, min, max;
849  PetscErrorCode ierr;
850 
851  ierr = MatGetSize(matrix, &n_rows, &n_cols);
852  AssertThrow(ierr == 0, ExcPETScError(ierr));
853 
854  ierr = MatGetLocalSize(matrix, &n_loc_rows, &n_loc_cols);
855  AssertThrow(ierr == 0, ExcPETScError(ierr));
856 
857  ierr = MatGetOwnershipRange(matrix, &min, &max);
858  AssertThrow(ierr == 0, ExcPETScError(ierr));
859 
860  Assert(n_loc_rows == max - min,
861  ExcMessage(
862  "PETSc is requiring non contiguous memory allocation."));
863 
864  IndexSet indices(n_rows);
865  indices.add_range(min, max);
866  indices.compress();
867 
868  return indices;
869  }
870 
871  void
873  const SparseMatrix &B,
874  const MPI::Vector & V) const
875  {
876  // Simply forward to the protected member function of the base class
877  // that takes abstract matrix and vector arguments (to which the compiler
878  // automatically casts the arguments).
879  MatrixBase::mmult(C, B, V);
880  }
881 
882  void
884  const SparseMatrix &B,
885  const MPI::Vector & V) const
886  {
887  // Simply forward to the protected member function of the base class
888  // that takes abstract matrix and vector arguments (to which the compiler
889  // automatically casts the arguments).
890  MatrixBase::Tmmult(C, B, V);
891  }
892 
893  } // namespace MPI
894 } // namespace PETScWrappers
895 
896 
898 
899 #endif // DEAL_II_WITH_PETSC
bool is_contiguous() const
Definition: index_set.h:1772
size_type size() const
Definition: index_set.h:1641
size_type n_elements() const
Definition: index_set.h:1789
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1677
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1837
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:855
void compress() const
Definition: index_set.h:1649
types::global_dof_index size_type
Definition: index_set.h:75
SparseMatrix & operator=(const value_type d)
void copy_from(const SparseMatrix &other)
void mmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
void reinit(const MPI_Comm &communicator, const SparsityPatternType &sparsity_pattern, const std::vector< size_type > &local_rows_per_process, const std::vector< size_type > &local_columns_per_process, const unsigned int this_process, const bool preset_nonzero_locations=true)
PetscScalar matrix_scalar_product(const Vector &u, const Vector &v) const
void Tmmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
PetscScalar matrix_norm_square(const Vector &v) const
void do_reinit(const MPI_Comm &comm, const SparsityPatternType &sparsity_pattern, const std::vector< size_type > &local_rows_per_process, const std::vector< size_type > &local_columns_per_process, const unsigned int this_process, const bool preset_nonzero_locations)
size_type row_length(const size_type row) const
void vmult(VectorBase &dst, const VectorBase &src) const
void mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
types::global_dof_index size_type
MatrixBase & operator=(const MatrixBase &)=delete
void Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
void compress(const VectorOperation::values operation)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1583
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1672
static const char A
static const char V
void set_keep_zero_rows(Mat &matrix)
PetscErrorCode destroy_matrix(Mat &matrix)
void close_matrix(Mat &matrix)
std::string to_string(const T &t)
Definition: patterns.h:2394
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TrilinosWrappers::types::int_type global_row_index(const Epetra_CrsMatrix &matrix, const ::types::global_dof_index i)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:81