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