Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
petsc_parallel_sparse_matrix.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2023 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
27
29
30namespace 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
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
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(),
191 "SparsityPattern and IndexSet have different number of rows"));
192 Assert(
193 sparsity_pattern.n_cols() == local_columns.size(),
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(),
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(),
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
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(),
478 "SparsityPattern and IndexSet have different number of rows."));
479 Assert(
480 sparsity_pattern.n_cols() == local_columns.size(),
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(),
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(),
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,
688 "BDDC preconditioner requires PETSc 3.10.0 or newer"));
689 }
690# endif
691 }
692
693# ifndef DOXYGEN
694 // explicit instantiations
695 //
697 const SparsityPattern &,
698 const std::vector<size_type> &,
699 const std::vector<size_type> &,
700 const unsigned int,
701 const bool);
704 const std::vector<size_type> &,
705 const std::vector<size_type> &,
706 const unsigned int,
707 const bool);
708
709 template void
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
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
738 const MPI_Comm);
739
740 template void
742 const IndexSet &,
744 const MPI_Comm);
745
746 template void
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
756 const std::vector<size_type> &,
757 const std::vector<size_type> &,
758 const unsigned int,
759 const bool);
760
761 template void
763 const IndexSet &,
764 const IndexSet &,
765 const SparsityPattern &);
766
767 template void
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 &,
786 const MPI_Comm);
787
788 template void
790 const IndexSet &,
791 const IndexSet &,
792 const IndexSet &,
793 const IndexSet &,
794 const SparsityPattern &);
795 template void
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
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,
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
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,
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
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
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void set_keep_zero_rows(Mat &matrix)
void close_matrix(Mat &matrix)
T sum(const T &t, const MPI_Comm mpi_communicator)