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