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_matrix_base.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
25
27
28namespace PETScWrappers
29{
30 namespace MatrixIterators
31 {
32# ifndef DOXYGEN
33 void
34 MatrixBase::const_iterator::Accessor::visit_present_row()
35 {
36 // if we are asked to visit the past-the-end line (or a line that is not
37 // stored on the current processor), then simply release all our caches
38 // and go on with life
39 if (matrix->in_local_range(this->a_row) == false)
40 {
41 colnum_cache.reset();
42 value_cache.reset();
43
44 return;
45 }
46
47 // get a representation of the present row
48 PetscInt ncols;
49 const PetscInt * colnums;
50 const PetscScalar *values;
51
52 PetscErrorCode ierr =
53 MatGetRow(*matrix, this->a_row, &ncols, &colnums, &values);
54 AssertThrow(ierr == 0, ExcPETScError(ierr));
55
56 // copy it into our caches if the line
57 // isn't empty. if it is, then we've
58 // done something wrong, since we
59 // shouldn't have initialized an
60 // iterator for an empty line (what
61 // would it point to?)
62 Assert(ncols != 0, ExcInternalError());
63 colnum_cache =
64 std::make_shared<std::vector<size_type>>(colnums, colnums + ncols);
65 value_cache =
66 std::make_shared<std::vector<PetscScalar>>(values, values + ncols);
67
68 // and finally restore the matrix
69 ierr = MatRestoreRow(*matrix, this->a_row, &ncols, &colnums, &values);
70 AssertThrow(ierr == 0, ExcPETScError(ierr));
71 }
72# endif
73 } // namespace MatrixIterators
74
75
76
78 : matrix(nullptr)
79 , last_action(VectorOperation::unknown)
80 {}
81
82
84 : matrix(A)
85 , last_action(VectorOperation::unknown)
86 {
87 const PetscErrorCode ierr =
88 PetscObjectReference(reinterpret_cast<PetscObject>(matrix));
89 AssertThrow(ierr == 0, ExcPETScError(ierr));
90 }
91
92 void
94 {
96 ExcMessage("Cannot assign a new Mat."));
97 PetscErrorCode ierr =
98 PetscObjectReference(reinterpret_cast<PetscObject>(A));
99 AssertThrow(ierr == 0, ExcPETScError(ierr));
100 ierr = MatDestroy(&matrix);
101 AssertThrow(ierr == 0, ExcPETScError(ierr));
102 matrix = A;
103 }
104
106 {
107 PetscErrorCode ierr = MatDestroy(&matrix);
108 (void)ierr;
109 AssertNothrow(ierr == 0, ExcPETScError(ierr));
110 }
111
112 void
114 {
115 // destroy the matrix...
116 {
117 const PetscErrorCode ierr = MatDestroy(&matrix);
118 AssertThrow(ierr == 0, ExcPETScError(ierr));
119 }
120
121 // ...and replace it by an empty
122 // sequential matrix
123 const int m = 0, n = 0, n_nonzero_per_row = 0;
124 const PetscErrorCode ierr = MatCreateSeqAIJ(
125 PETSC_COMM_SELF, m, n, n_nonzero_per_row, nullptr, &matrix);
126 AssertThrow(ierr == 0, ExcPETScError(ierr));
127 }
128
129
130
131 MatrixBase &
133 {
134 (void)d;
136
138
139 const PetscErrorCode ierr = MatZeroEntries(matrix);
140 AssertThrow(ierr == 0, ExcPETScError(ierr));
141
142 return *this;
143 }
144
145
146
147 void
148 MatrixBase::clear_row(const size_type row, const PetscScalar new_diag_value)
149 {
150 std::vector<size_type> rows(1, row);
151 clear_rows(rows, new_diag_value);
152 }
153
154
155
156 void
157 MatrixBase::clear_rows(const std::vector<size_type> &rows,
158 const PetscScalar new_diag_value)
159 {
161
162 // now set all the entries of these rows
163 // to zero
164 const std::vector<PetscInt> petsc_rows(rows.begin(), rows.end());
165
166 // call the functions. note that we have
167 // to call them even if #rows is empty,
168 // since this is a collective operation
169 IS index_set;
170
171 ISCreateGeneral(get_mpi_communicator(),
172 rows.size(),
173 petsc_rows.data(),
174 PETSC_COPY_VALUES,
175 &index_set);
176
177 const PetscErrorCode ierr =
178 MatZeroRowsIS(matrix, index_set, new_diag_value, nullptr, nullptr);
179 AssertThrow(ierr == 0, ExcPETScError(ierr));
180 ISDestroy(&index_set);
181 }
182
183 void
184 MatrixBase::clear_rows_columns(const std::vector<size_type> &rows,
185 const PetscScalar new_diag_value)
186 {
188
189 // now set all the entries of these rows
190 // to zero
191 const std::vector<PetscInt> petsc_rows(rows.begin(), rows.end());
192
193 // call the functions. note that we have
194 // to call them even if #rows is empty,
195 // since this is a collective operation
196 IS index_set;
197
198 ISCreateGeneral(get_mpi_communicator(),
199 rows.size(),
200 petsc_rows.data(),
201 PETSC_COPY_VALUES,
202 &index_set);
203
204 const PetscErrorCode ierr =
205 MatZeroRowsColumnsIS(matrix, index_set, new_diag_value, nullptr, nullptr);
206 AssertThrow(ierr == 0, ExcPETScError(ierr));
207 ISDestroy(&index_set);
208 }
209
210
211
212 PetscScalar
213 MatrixBase::el(const size_type i, const size_type j) const
214 {
215 PetscInt petsc_i = i, petsc_j = j;
216
217 PetscScalar value;
218
219 const PetscErrorCode ierr =
220 MatGetValues(matrix, 1, &petsc_i, 1, &petsc_j, &value);
221 AssertThrow(ierr == 0, ExcPETScError(ierr));
222
223 return value;
224 }
225
226
227
228 PetscScalar
230 {
231 Assert(m() == n(), ExcNotQuadratic());
232
233 // this doesn't seem to work any
234 // different than any other element
235 return el(i, i);
236 }
237
238
239
240 void
242 {
243 {
244# ifdef DEBUG
245# ifdef DEAL_II_WITH_MPI
246 // Check that all processors agree that last_action is the same (or none!)
247
248 int my_int_last_action = last_action;
249 int all_int_last_action;
250
251 const int ierr = MPI_Allreduce(&my_int_last_action,
252 &all_int_last_action,
253 1,
254 MPI_INT,
255 MPI_BOR,
257 AssertThrowMPI(ierr);
258
259 AssertThrow(all_int_last_action !=
261 ExcMessage("Error: not all processors agree on the last "
262 "VectorOperation before this compress() call."));
263# endif
264# endif
265 }
266
270 "Missing compress() or calling with wrong VectorOperation argument."));
271
272 // flush buffers
273 PetscErrorCode ierr = MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY);
274 AssertThrow(ierr == 0, ExcPETScError(ierr));
275
276 ierr = MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY);
277 AssertThrow(ierr == 0, ExcPETScError(ierr));
278
280 }
281
282
283
286 {
287 PetscInt n_rows, n_cols;
288
289 const PetscErrorCode ierr = MatGetSize(matrix, &n_rows, &n_cols);
290 AssertThrow(ierr == 0, ExcPETScError(ierr));
291
292 return n_rows;
293 }
294
295
296
299 {
300 PetscInt n_rows, n_cols;
301
302 const PetscErrorCode ierr = MatGetSize(matrix, &n_rows, &n_cols);
303 AssertThrow(ierr == 0, ExcPETScError(ierr));
304
305 return n_cols;
306 }
307
308
309
312 {
313 PetscInt n_rows;
314
315 const PetscErrorCode ierr = MatGetLocalSize(matrix, &n_rows, nullptr);
316 AssertThrow(ierr == 0, ExcPETScError(ierr));
317
318 return n_rows;
319 }
320
321
322
323 std::pair<MatrixBase::size_type, MatrixBase::size_type>
325 {
326 PetscInt begin, end;
327
328 const PetscErrorCode ierr =
329 MatGetOwnershipRange(static_cast<const Mat &>(matrix), &begin, &end);
330 AssertThrow(ierr == 0, ExcPETScError(ierr));
331
332 return {begin, end};
333 }
334
335
336
339 {
340 PetscInt n_cols;
341
342 const PetscErrorCode ierr = MatGetLocalSize(matrix, nullptr, &n_cols);
343 AssertThrow(ierr == 0, ExcPETScError(ierr));
344
345 return n_cols;
346 }
347
348
349
350 std::pair<MatrixBase::size_type, MatrixBase::size_type>
352 {
353 PetscInt begin, end;
354
355 const PetscErrorCode ierr =
356 MatGetOwnershipRangeColumn(static_cast<const Mat &>(matrix),
357 &begin,
358 &end);
359 AssertThrow(ierr == 0, ExcPETScError(ierr));
360
361 return {begin, end};
362 }
363
364
365
366 std::uint64_t
368 {
369 MatInfo mat_info;
370 const PetscErrorCode ierr = MatGetInfo(matrix, MAT_GLOBAL_SUM, &mat_info);
371 AssertThrow(ierr == 0, ExcPETScError(ierr));
372
373 // MatInfo logs quantities as PetscLogDouble. So we need to cast it to match
374 // our interface.
375 return static_cast<std::uint64_t>(mat_info.nz_used);
376 }
377
378
379
382 {
383 // TODO: this function will probably only work if compress() was called on
384 // the matrix previously. however, we can't do this here, since it would
385 // impose global communication and one would have to make sure that this
386 // function is called the same number of times from all processors,
387 // something that is unreasonable. there should simply be a way in PETSc to
388 // query the number of entries in a row bypassing the call to compress(),
389 // but I can't find one
390 Assert(row < m(), ExcInternalError());
391
392 // get a representation of the present
393 // row
394 PetscInt ncols;
395 const PetscInt * colnums;
396 const PetscScalar *values;
397
398 // TODO: this is probably horribly inefficient; we should lobby for a way to
399 // query this information from PETSc
400 PetscErrorCode ierr = MatGetRow(*this, row, &ncols, &colnums, &values);
401 AssertThrow(ierr == 0, ExcPETScError(ierr));
402
403 // then restore the matrix and return the number of columns in this row as
404 // queried previously. Starting with PETSc 3.4, MatRestoreRow actually
405 // resets the last three arguments to nullptr, to avoid abuse of pointers
406 // now dangling. as a consequence, we need to save the size of the array
407 // and return the saved value.
408 const PetscInt ncols_saved = ncols;
409 ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values);
410 AssertThrow(ierr == 0, ExcPETScError(ierr));
411
412 return ncols_saved;
413 }
414
415
416 PetscReal
418 {
419 PetscReal result;
420
421 const PetscErrorCode ierr = MatNorm(matrix, NORM_1, &result);
422 AssertThrow(ierr == 0, ExcPETScError(ierr));
423
424 return result;
425 }
426
427
428
429 PetscReal
431 {
432 PetscReal result;
433
434 const PetscErrorCode ierr = MatNorm(matrix, NORM_INFINITY, &result);
435 AssertThrow(ierr == 0, ExcPETScError(ierr));
436
437 return result;
438 }
439
440
441
442 PetscReal
444 {
445 PetscReal result;
446
447 const PetscErrorCode ierr = MatNorm(matrix, NORM_FROBENIUS, &result);
448 AssertThrow(ierr == 0, ExcPETScError(ierr));
449
450 return result;
451 }
452
453
454 PetscScalar
456 {
457 VectorBase tmp(v);
458 vmult(tmp, v);
459 return tmp * v;
460 }
461
462
463 PetscScalar
465 const VectorBase &v) const
466 {
467 VectorBase tmp(u);
468 vmult(tmp, v);
469 return u * tmp;
470 }
471
472
473 PetscScalar
475 {
476 PetscScalar result;
477
478 const PetscErrorCode ierr = MatGetTrace(matrix, &result);
479 AssertThrow(ierr == 0, ExcPETScError(ierr));
480
481 return result;
482 }
483
484
485
486 MatrixBase &
487 MatrixBase::operator*=(const PetscScalar a)
488 {
489 const PetscErrorCode ierr = MatScale(matrix, a);
490 AssertThrow(ierr == 0, ExcPETScError(ierr));
491
492 return *this;
493 }
494
495
496
497 MatrixBase &
498 MatrixBase::operator/=(const PetscScalar a)
499 {
500 const PetscScalar factor = 1. / a;
501 const PetscErrorCode ierr = MatScale(matrix, factor);
502 AssertThrow(ierr == 0, ExcPETScError(ierr));
503
504 return *this;
505 }
506
507
508
509 MatrixBase &
510 MatrixBase::add(const PetscScalar factor, const MatrixBase &other)
511 {
512 const PetscErrorCode ierr =
513 MatAXPY(matrix, factor, other, DIFFERENT_NONZERO_PATTERN);
514 AssertThrow(ierr == 0, ExcPETScError(ierr));
515
516 return *this;
517 }
518
519
520 void
522 {
523 Assert(&src != &dst, ExcSourceEqualsDestination());
524
525 const PetscErrorCode ierr = MatMult(matrix, src, dst);
526 AssertThrow(ierr == 0, ExcPETScError(ierr));
527 }
528
529
530
531 void
533 {
534 Assert(&src != &dst, ExcSourceEqualsDestination());
535
536 const PetscErrorCode ierr = MatMultTranspose(matrix, src, dst);
537 AssertThrow(ierr == 0, ExcPETScError(ierr));
538 }
539
540
541
542 void
544 {
545 Assert(&src != &dst, ExcSourceEqualsDestination());
546
547 const PetscErrorCode ierr = MatMultAdd(matrix, src, dst, dst);
548 AssertThrow(ierr == 0, ExcPETScError(ierr));
549 }
550
551
552
553 void
555 {
556 Assert(&src != &dst, ExcSourceEqualsDestination());
557
558 const PetscErrorCode ierr = MatMultTransposeAdd(matrix, src, dst, dst);
559 AssertThrow(ierr == 0, ExcPETScError(ierr));
560 }
561
562
563 namespace internals
564 {
565 void
566 perform_mmult(const MatrixBase &inputleft,
567 const MatrixBase &inputright,
568 MatrixBase & result,
569 const VectorBase &V,
570 const bool transpose_left)
571 {
572 const bool use_vector = (V.size() == inputright.m() ? true : false);
573 if (transpose_left == false)
574 {
575 Assert(inputleft.n() == inputright.m(),
576 ExcDimensionMismatch(inputleft.n(), inputright.m()));
577 }
578 else
579 {
580 Assert(inputleft.m() == inputright.m(),
581 ExcDimensionMismatch(inputleft.m(), inputright.m()));
582 }
583
584 result.clear();
585
586 PetscErrorCode ierr;
587
588 if (use_vector == false)
589 {
590 if (transpose_left)
591 {
592 ierr = MatTransposeMatMult(inputleft,
593 inputright,
594 MAT_INITIAL_MATRIX,
595 PETSC_DEFAULT,
596 &result.petsc_matrix());
597 AssertThrow(ierr == 0, ExcPETScError(ierr));
598 }
599 else
600 {
601 ierr = MatMatMult(inputleft,
602 inputright,
603 MAT_INITIAL_MATRIX,
604 PETSC_DEFAULT,
605 &result.petsc_matrix());
606 AssertThrow(ierr == 0, ExcPETScError(ierr));
607 }
608 }
609 else
610 {
611 Mat tmp;
612 ierr = MatDuplicate(inputleft, MAT_COPY_VALUES, &tmp);
613 AssertThrow(ierr == 0, ExcPETScError(ierr));
614 if (transpose_left)
615 {
616# if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
617 ierr = MatTranspose(tmp, MAT_REUSE_MATRIX, &tmp);
618# else
619 ierr = MatTranspose(tmp, MAT_INPLACE_MATRIX, &tmp);
620# endif
621 AssertThrow(ierr == 0, ExcPETScError(ierr));
622 }
623 ierr = MatDiagonalScale(tmp, nullptr, V);
624 AssertThrow(ierr == 0, ExcPETScError(ierr));
625 ierr = MatMatMult(tmp,
626 inputright,
627 MAT_INITIAL_MATRIX,
628 PETSC_DEFAULT,
629 &result.petsc_matrix());
630 AssertThrow(ierr == 0, ExcPETScError(ierr));
631 ierr = MatDestroy(&tmp);
632 AssertThrow(ierr == 0, ExcPETScError(ierr));
633 }
634 }
635 } // namespace internals
636
637 void
639 const MatrixBase &B,
640 const VectorBase &V) const
641 {
642 internals::perform_mmult(*this, B, C, V, false);
643 }
644
645 void
647 const MatrixBase &B,
648 const VectorBase &V) const
649 {
650 internals::perform_mmult(*this, B, C, V, true);
651 }
652
653 PetscScalar
655 const VectorBase &x,
656 const VectorBase &b) const
657 {
658 // avoid the use of a temporary, and
659 // rather do one negation pass more than
660 // necessary
661 vmult(dst, x);
662 dst -= b;
663 dst *= -1;
664
665 return dst.l2_norm();
666 }
667
668
669
670 MatrixBase::operator Mat() const
671 {
672 return matrix;
673 }
674
675 Mat &
677 {
678 return matrix;
679 }
680
681 void
683 {
684# if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
685 const PetscErrorCode ierr = MatTranspose(matrix, MAT_REUSE_MATRIX, &matrix);
686# else
687 const PetscErrorCode ierr =
688 MatTranspose(matrix, MAT_INPLACE_MATRIX, &matrix);
689# endif
690 AssertThrow(ierr == 0, ExcPETScError(ierr));
691 }
692
693 PetscBool
694 MatrixBase::is_symmetric(const double tolerance)
695 {
696 PetscBool truth;
698 const PetscErrorCode ierr = MatIsSymmetric(matrix, tolerance, &truth);
699 AssertThrow(ierr == 0, ExcPETScError(ierr));
700 return truth;
701 }
702
703 PetscBool
704 MatrixBase::is_hermitian(const double tolerance)
705 {
706 PetscBool truth;
707
709 const PetscErrorCode ierr = MatIsHermitian(matrix, tolerance, &truth);
710 AssertThrow(ierr == 0, ExcPETScError(ierr));
711
712 return truth;
713 }
714
715 void
716 MatrixBase::write_ascii(const PetscViewerFormat format)
717 {
719 MPI_Comm comm = PetscObjectComm(reinterpret_cast<PetscObject>(matrix));
720
721 // Set options
722 PetscErrorCode ierr =
723 PetscViewerSetFormat(PETSC_VIEWER_STDOUT_(comm), format);
724 AssertThrow(ierr == 0, ExcPETScError(ierr));
725
726 // Write to screen
727 ierr = MatView(matrix, PETSC_VIEWER_STDOUT_(comm));
728 AssertThrow(ierr == 0, ExcPETScError(ierr));
729 }
730
731 void
732 MatrixBase::print(std::ostream &out, const bool /*alternative_output*/) const
733 {
734 PetscBool has;
735
736 PetscErrorCode ierr = MatHasOperation(matrix, MATOP_GET_ROW, &has);
737 AssertThrow(ierr == 0, ExcPETScError(ierr));
738
739 Mat vmatrix = matrix;
740 if (!has)
741 {
742 ierr = MatConvert(matrix, MATAIJ, MAT_INITIAL_MATRIX, &vmatrix);
743 AssertThrow(ierr == 0, ExcPETScError(ierr));
744 }
745
746 std::pair<MatrixBase::size_type, MatrixBase::size_type> loc_range =
747 local_range();
748
749 PetscInt ncols;
750 const PetscInt * colnums;
751 const PetscScalar *values;
752
754 for (row = loc_range.first; row < loc_range.second; ++row)
755 {
756 ierr = MatGetRow(vmatrix, row, &ncols, &colnums, &values);
757 AssertThrow(ierr == 0, ExcPETScError(ierr));
758
759 for (PetscInt col = 0; col < ncols; ++col)
760 {
761 out << "(" << row << "," << colnums[col] << ") " << values[col]
762 << std::endl;
763 }
764
765 ierr = MatRestoreRow(vmatrix, row, &ncols, &colnums, &values);
766 AssertThrow(ierr == 0, ExcPETScError(ierr));
767 }
768 if (vmatrix != matrix)
769 {
770 ierr = MatDestroy(&vmatrix);
771 AssertThrow(ierr == 0, ExcPETScError(ierr));
772 }
773 AssertThrow(out.fail() == false, ExcIO());
774 }
775
776
777
778 std::size_t
780 {
781 MatInfo info;
782 const PetscErrorCode ierr = MatGetInfo(matrix, MAT_LOCAL, &info);
783 AssertThrow(ierr == 0, ExcPETScError(ierr));
784
785 return sizeof(*this) + static_cast<size_type>(info.memory);
786 }
787
788} // namespace PETScWrappers
789
791
792#endif // DEAL_II_WITH_PETSC
void add(const size_type i, const size_type j, const PetscScalar value)
size_type row_length(const size_type row) const
std::size_t memory_consumption() const
VectorOperation::values last_action
void vmult(VectorBase &dst, const VectorBase &src) const
void clear_rows(const std::vector< size_type > &rows, const PetscScalar new_diag_value=0)
MPI_Comm get_mpi_communicator() const
PetscScalar diag_element(const size_type i) const
size_type local_domain_size() const
const_iterator begin() const
MatrixBase & operator/=(const PetscScalar factor)
void mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
PetscBool is_symmetric(const double tolerance=1.e-12)
std::pair< size_type, size_type > local_domain() const
const_iterator end() const
void Tvmult_add(VectorBase &dst, const VectorBase &src) const
void print(std::ostream &out, const bool alternative_output=false) const
PetscScalar el(const size_type i, const size_type j) const
MatrixBase & operator=(const MatrixBase &)=delete
PetscBool is_hermitian(const double tolerance=1.e-12)
PetscScalar matrix_scalar_product(const VectorBase &u, const VectorBase &v) const
void Tvmult(VectorBase &dst, const VectorBase &src) const
void clear_rows_columns(const std::vector< size_type > &row_and_column_indices, const PetscScalar new_diag_value=0)
MatrixBase & operator*=(const PetscScalar factor)
PetscScalar residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const
void Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
void vmult_add(VectorBase &dst, const VectorBase &src) const
std::pair< size_type, size_type > local_range() const
void compress(const VectorOperation::values operation)
PetscScalar matrix_norm_square(const VectorBase &v) const
std::uint64_t n_nonzero_elements() const
void clear_row(const size_type row, const PetscScalar new_diag_value=0)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
#define Assert(cond, exc)
#define AssertThrowMPI(error_code)
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ matrix
Contents is actually a matrix.
void perform_mmult(const MatrixBase &inputleft, const MatrixBase &inputright, MatrixBase &result, const VectorBase &V, const bool transpose_left)
const MPI_Comm comm