1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2004 - 2020 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 16 #include <deal.II/lac/petsc_matrix_base.h> 17 18 #ifdef DEAL_II_WITH_PETSC 19 20 # include <deal.II/lac/exceptions.h> 21 # include <deal.II/lac/petsc_compatibility.h> 22 # include <deal.II/lac/petsc_full_matrix.h> 23 # include <deal.II/lac/petsc_sparse_matrix.h> 24 # include <deal.II/lac/petsc_vector_base.h> 25 26 DEAL_II_NAMESPACE_OPEN 27 28 namespace PETScWrappers 29 { 30 namespace MatrixIterators 31 { 32 # ifndef DOXYGEN 33 void visit_present_row()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 MatrixBase()77 MatrixBase::MatrixBase() 78 : matrix(nullptr) 79 , last_action(VectorOperation::unknown) 80 {} 81 82 83 ~MatrixBase()84 MatrixBase::~MatrixBase() 85 { 86 destroy_matrix(matrix); 87 } 88 89 90 91 void clear()92 MatrixBase::clear() 93 { 94 // destroy the matrix... 95 { 96 const PetscErrorCode ierr = destroy_matrix(matrix); 97 AssertThrow(ierr == 0, ExcPETScError(ierr)); 98 } 99 100 // ...and replace it by an empty 101 // sequential matrix 102 const int m = 0, n = 0, n_nonzero_per_row = 0; 103 const PetscErrorCode ierr = MatCreateSeqAIJ( 104 PETSC_COMM_SELF, m, n, n_nonzero_per_row, nullptr, &matrix); 105 AssertThrow(ierr == 0, ExcPETScError(ierr)); 106 } 107 108 109 110 MatrixBase & operator =(const value_type d)111 MatrixBase::operator=(const value_type d) 112 { 113 (void)d; 114 Assert(d == value_type(), ExcScalarAssignmentOnlyForZeroValue()); 115 116 assert_is_compressed(); 117 118 const PetscErrorCode ierr = MatZeroEntries(matrix); 119 AssertThrow(ierr == 0, ExcPETScError(ierr)); 120 121 return *this; 122 } 123 124 125 126 void clear_row(const size_type row,const PetscScalar new_diag_value)127 MatrixBase::clear_row(const size_type row, const PetscScalar new_diag_value) 128 { 129 std::vector<size_type> rows(1, row); 130 clear_rows(rows, new_diag_value); 131 } 132 133 134 135 void clear_rows(const std::vector<size_type> & rows,const PetscScalar new_diag_value)136 MatrixBase::clear_rows(const std::vector<size_type> &rows, 137 const PetscScalar new_diag_value) 138 { 139 assert_is_compressed(); 140 141 // now set all the entries of these rows 142 // to zero 143 const std::vector<PetscInt> petsc_rows(rows.begin(), rows.end()); 144 145 // call the functions. note that we have 146 // to call them even if #rows is empty, 147 // since this is a collective operation 148 IS index_set; 149 150 ISCreateGeneral(get_mpi_communicator(), 151 rows.size(), 152 petsc_rows.data(), 153 PETSC_COPY_VALUES, 154 &index_set); 155 156 const PetscErrorCode ierr = 157 MatZeroRowsIS(matrix, index_set, new_diag_value, nullptr, nullptr); 158 AssertThrow(ierr == 0, ExcPETScError(ierr)); 159 ISDestroy(&index_set); 160 } 161 162 163 164 PetscScalar el(const size_type i,const size_type j) const165 MatrixBase::el(const size_type i, const size_type j) const 166 { 167 PetscInt petsc_i = i, petsc_j = j; 168 169 PetscScalar value; 170 171 const PetscErrorCode ierr = 172 MatGetValues(matrix, 1, &petsc_i, 1, &petsc_j, &value); 173 AssertThrow(ierr == 0, ExcPETScError(ierr)); 174 175 return value; 176 } 177 178 179 180 PetscScalar diag_element(const size_type i) const181 MatrixBase::diag_element(const size_type i) const 182 { 183 Assert(m() == n(), ExcNotQuadratic()); 184 185 // this doesn't seem to work any 186 // different than any other element 187 return el(i, i); 188 } 189 190 191 192 void compress(const VectorOperation::values operation)193 MatrixBase::compress(const VectorOperation::values operation) 194 { 195 { 196 # ifdef DEBUG 197 # ifdef DEAL_II_WITH_MPI 198 // Check that all processors agree that last_action is the same (or none!) 199 200 int my_int_last_action = last_action; 201 int all_int_last_action; 202 203 const int ierr = MPI_Allreduce(&my_int_last_action, 204 &all_int_last_action, 205 1, 206 MPI_INT, 207 MPI_BOR, 208 get_mpi_communicator()); 209 AssertThrowMPI(ierr); 210 211 AssertThrow(all_int_last_action != 212 (VectorOperation::add | VectorOperation::insert), 213 ExcMessage("Error: not all processors agree on the last " 214 "VectorOperation before this compress() call.")); 215 # endif 216 # endif 217 } 218 219 AssertThrow( 220 last_action == VectorOperation::unknown || last_action == operation, 221 ExcMessage( 222 "Missing compress() or calling with wrong VectorOperation argument.")); 223 224 // flush buffers 225 PetscErrorCode ierr = MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY); 226 AssertThrow(ierr == 0, ExcPETScError(ierr)); 227 228 ierr = MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY); 229 AssertThrow(ierr == 0, ExcPETScError(ierr)); 230 231 last_action = VectorOperation::unknown; 232 } 233 234 235 236 MatrixBase::size_type m() const237 MatrixBase::m() const 238 { 239 PetscInt n_rows, n_cols; 240 241 const PetscErrorCode ierr = MatGetSize(matrix, &n_rows, &n_cols); 242 AssertThrow(ierr == 0, ExcPETScError(ierr)); 243 244 return n_rows; 245 } 246 247 248 249 MatrixBase::size_type n() const250 MatrixBase::n() const 251 { 252 PetscInt n_rows, n_cols; 253 254 const PetscErrorCode ierr = MatGetSize(matrix, &n_rows, &n_cols); 255 AssertThrow(ierr == 0, ExcPETScError(ierr)); 256 257 return n_cols; 258 } 259 260 261 262 MatrixBase::size_type local_size() const263 MatrixBase::local_size() const 264 { 265 PetscInt n_rows, n_cols; 266 267 const PetscErrorCode ierr = MatGetLocalSize(matrix, &n_rows, &n_cols); 268 AssertThrow(ierr == 0, ExcPETScError(ierr)); 269 270 return n_rows; 271 } 272 273 274 275 std::pair<MatrixBase::size_type, MatrixBase::size_type> local_range() const276 MatrixBase::local_range() const 277 { 278 PetscInt begin, end; 279 280 const PetscErrorCode ierr = 281 MatGetOwnershipRange(static_cast<const Mat &>(matrix), &begin, &end); 282 AssertThrow(ierr == 0, ExcPETScError(ierr)); 283 284 return std::make_pair(begin, end); 285 } 286 287 288 289 MatrixBase::size_type n_nonzero_elements() const290 MatrixBase::n_nonzero_elements() const 291 { 292 MatInfo mat_info; 293 const PetscErrorCode ierr = MatGetInfo(matrix, MAT_GLOBAL_SUM, &mat_info); 294 AssertThrow(ierr == 0, ExcPETScError(ierr)); 295 296 return static_cast<size_type>(mat_info.nz_used); 297 } 298 299 300 301 MatrixBase::size_type row_length(const size_type row) const302 MatrixBase::row_length(const size_type row) const 303 { 304 // TODO: this function will probably only work if compress() was called on 305 // the matrix previously. however, we can't do this here, since it would 306 // impose global communication and one would have to make sure that this 307 // function is called the same number of times from all processors, 308 // something that is unreasonable. there should simply be a way in PETSc to 309 // query the number of entries in a row bypassing the call to compress(), 310 // but I can't find one 311 Assert(row < m(), ExcInternalError()); 312 313 // get a representation of the present 314 // row 315 PetscInt ncols; 316 const PetscInt * colnums; 317 const PetscScalar *values; 318 319 // TODO: this is probably horribly inefficient; we should lobby for a way to 320 // query this information from PETSc 321 PetscErrorCode ierr = MatGetRow(*this, row, &ncols, &colnums, &values); 322 AssertThrow(ierr == 0, ExcPETScError(ierr)); 323 324 // then restore the matrix and return the number of columns in this row as 325 // queried previously. Starting with PETSc 3.4, MatRestoreRow actually 326 // resets the last three arguments to nullptr, to avoid abuse of pointers 327 // now dangling. as a consequence, we need to save the size of the array 328 // and return the saved value. 329 const PetscInt ncols_saved = ncols; 330 ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values); 331 AssertThrow(ierr == 0, ExcPETScError(ierr)); 332 333 return ncols_saved; 334 } 335 336 337 PetscReal l1_norm() const338 MatrixBase::l1_norm() const 339 { 340 PetscReal result; 341 342 const PetscErrorCode ierr = MatNorm(matrix, NORM_1, &result); 343 AssertThrow(ierr == 0, ExcPETScError(ierr)); 344 345 return result; 346 } 347 348 349 350 PetscReal linfty_norm() const351 MatrixBase::linfty_norm() const 352 { 353 PetscReal result; 354 355 const PetscErrorCode ierr = MatNorm(matrix, NORM_INFINITY, &result); 356 AssertThrow(ierr == 0, ExcPETScError(ierr)); 357 358 return result; 359 } 360 361 362 363 PetscReal frobenius_norm() const364 MatrixBase::frobenius_norm() const 365 { 366 PetscReal result; 367 368 const PetscErrorCode ierr = MatNorm(matrix, NORM_FROBENIUS, &result); 369 AssertThrow(ierr == 0, ExcPETScError(ierr)); 370 371 return result; 372 } 373 374 375 PetscScalar matrix_norm_square(const VectorBase & v) const376 MatrixBase::matrix_norm_square(const VectorBase &v) const 377 { 378 VectorBase tmp(v); 379 vmult(tmp, v); 380 return tmp * v; 381 } 382 383 384 PetscScalar matrix_scalar_product(const VectorBase & u,const VectorBase & v) const385 MatrixBase::matrix_scalar_product(const VectorBase &u, 386 const VectorBase &v) const 387 { 388 VectorBase tmp(u); 389 vmult(tmp, v); 390 return u * tmp; 391 } 392 393 394 PetscScalar trace() const395 MatrixBase::trace() const 396 { 397 PetscScalar result; 398 399 const PetscErrorCode ierr = MatGetTrace(matrix, &result); 400 AssertThrow(ierr == 0, ExcPETScError(ierr)); 401 402 return result; 403 } 404 405 406 407 MatrixBase & operator *=(const PetscScalar a)408 MatrixBase::operator*=(const PetscScalar a) 409 { 410 const PetscErrorCode ierr = MatScale(matrix, a); 411 AssertThrow(ierr == 0, ExcPETScError(ierr)); 412 413 return *this; 414 } 415 416 417 418 MatrixBase & operator /=(const PetscScalar a)419 MatrixBase::operator/=(const PetscScalar a) 420 { 421 const PetscScalar factor = 1. / a; 422 const PetscErrorCode ierr = MatScale(matrix, factor); 423 AssertThrow(ierr == 0, ExcPETScError(ierr)); 424 425 return *this; 426 } 427 428 429 MatrixBase & add(const PetscScalar factor,const MatrixBase & other)430 MatrixBase::add(const PetscScalar factor, const MatrixBase &other) 431 { 432 const PetscErrorCode ierr = 433 MatAXPY(matrix, factor, other, DIFFERENT_NONZERO_PATTERN); 434 AssertThrow(ierr == 0, ExcPETScError(ierr)); 435 436 return *this; 437 } 438 439 440 441 MatrixBase & add(const MatrixBase & other,const PetscScalar factor)442 MatrixBase::add(const MatrixBase &other, const PetscScalar factor) 443 { 444 return add(factor, other); 445 } 446 447 448 void vmult(VectorBase & dst,const VectorBase & src) const449 MatrixBase::vmult(VectorBase &dst, const VectorBase &src) const 450 { 451 Assert(&src != &dst, ExcSourceEqualsDestination()); 452 453 const PetscErrorCode ierr = MatMult(matrix, src, dst); 454 AssertThrow(ierr == 0, ExcPETScError(ierr)); 455 } 456 457 458 459 void Tvmult(VectorBase & dst,const VectorBase & src) const460 MatrixBase::Tvmult(VectorBase &dst, const VectorBase &src) const 461 { 462 Assert(&src != &dst, ExcSourceEqualsDestination()); 463 464 const PetscErrorCode ierr = MatMultTranspose(matrix, src, dst); 465 AssertThrow(ierr == 0, ExcPETScError(ierr)); 466 } 467 468 469 470 void vmult_add(VectorBase & dst,const VectorBase & src) const471 MatrixBase::vmult_add(VectorBase &dst, const VectorBase &src) const 472 { 473 Assert(&src != &dst, ExcSourceEqualsDestination()); 474 475 const PetscErrorCode ierr = MatMultAdd(matrix, src, dst, dst); 476 AssertThrow(ierr == 0, ExcPETScError(ierr)); 477 } 478 479 480 481 void Tvmult_add(VectorBase & dst,const VectorBase & src) const482 MatrixBase::Tvmult_add(VectorBase &dst, const VectorBase &src) const 483 { 484 Assert(&src != &dst, ExcSourceEqualsDestination()); 485 486 const PetscErrorCode ierr = MatMultTransposeAdd(matrix, src, dst, dst); 487 AssertThrow(ierr == 0, ExcPETScError(ierr)); 488 } 489 490 491 namespace internals 492 { 493 void perform_mmult(const MatrixBase & inputleft,const MatrixBase & inputright,MatrixBase & result,const VectorBase & V,const bool transpose_left)494 perform_mmult(const MatrixBase &inputleft, 495 const MatrixBase &inputright, 496 MatrixBase & result, 497 const VectorBase &V, 498 const bool transpose_left) 499 { 500 const bool use_vector = (V.size() == inputright.m() ? true : false); 501 if (transpose_left == false) 502 { 503 Assert(inputleft.n() == inputright.m(), 504 ExcDimensionMismatch(inputleft.n(), inputright.m())); 505 } 506 else 507 { 508 Assert(inputleft.m() == inputright.m(), 509 ExcDimensionMismatch(inputleft.m(), inputright.m())); 510 } 511 512 result.clear(); 513 514 PetscErrorCode ierr; 515 516 if (use_vector == false) 517 { 518 if (transpose_left) 519 { 520 ierr = MatTransposeMatMult(inputleft, 521 inputright, 522 MAT_INITIAL_MATRIX, 523 PETSC_DEFAULT, 524 &result.petsc_matrix()); 525 AssertThrow(ierr == 0, ExcPETScError(ierr)); 526 } 527 else 528 { 529 ierr = MatMatMult(inputleft, 530 inputright, 531 MAT_INITIAL_MATRIX, 532 PETSC_DEFAULT, 533 &result.petsc_matrix()); 534 AssertThrow(ierr == 0, ExcPETScError(ierr)); 535 } 536 } 537 else 538 { 539 Mat tmp; 540 ierr = MatDuplicate(inputleft, MAT_COPY_VALUES, &tmp); 541 AssertThrow(ierr == 0, ExcPETScError(ierr)); 542 if (transpose_left) 543 { 544 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0) 545 ierr = MatTranspose(tmp, MAT_REUSE_MATRIX, &tmp); 546 # else 547 ierr = MatTranspose(tmp, MAT_INPLACE_MATRIX, &tmp); 548 # endif 549 AssertThrow(ierr == 0, ExcPETScError(ierr)); 550 } 551 ierr = MatDiagonalScale(tmp, nullptr, V); 552 AssertThrow(ierr == 0, ExcPETScError(ierr)); 553 ierr = MatMatMult(tmp, 554 inputright, 555 MAT_INITIAL_MATRIX, 556 PETSC_DEFAULT, 557 &result.petsc_matrix()); 558 AssertThrow(ierr == 0, ExcPETScError(ierr)); 559 ierr = PETScWrappers::destroy_matrix(tmp); 560 AssertThrow(ierr == 0, ExcPETScError(ierr)); 561 } 562 } 563 } // namespace internals 564 565 void mmult(MatrixBase & C,const MatrixBase & B,const VectorBase & V) const566 MatrixBase::mmult(MatrixBase & C, 567 const MatrixBase &B, 568 const VectorBase &V) const 569 { 570 internals::perform_mmult(*this, B, C, V, false); 571 } 572 573 void Tmmult(MatrixBase & C,const MatrixBase & B,const VectorBase & V) const574 MatrixBase::Tmmult(MatrixBase & C, 575 const MatrixBase &B, 576 const VectorBase &V) const 577 { 578 internals::perform_mmult(*this, B, C, V, true); 579 } 580 581 PetscScalar residual(VectorBase & dst,const VectorBase & x,const VectorBase & b) const582 MatrixBase::residual(VectorBase & dst, 583 const VectorBase &x, 584 const VectorBase &b) const 585 { 586 // avoid the use of a temporary, and 587 // rather do one negation pass more than 588 // necessary 589 vmult(dst, x); 590 dst -= b; 591 dst *= -1; 592 593 return dst.l2_norm(); 594 } 595 596 597 operator Mat() const598 MatrixBase::operator Mat() const 599 { 600 return matrix; 601 } 602 603 Mat & petsc_matrix()604 MatrixBase::petsc_matrix() 605 { 606 return matrix; 607 } 608 609 void transpose()610 MatrixBase::transpose() 611 { 612 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0) 613 const PetscErrorCode ierr = MatTranspose(matrix, MAT_REUSE_MATRIX, &matrix); 614 # else 615 const PetscErrorCode ierr = 616 MatTranspose(matrix, MAT_INPLACE_MATRIX, &matrix); 617 # endif 618 AssertThrow(ierr == 0, ExcPETScError(ierr)); 619 } 620 621 PetscBool is_symmetric(const double tolerance)622 MatrixBase::is_symmetric(const double tolerance) 623 { 624 PetscBool truth; 625 assert_is_compressed(); 626 const PetscErrorCode ierr = MatIsSymmetric(matrix, tolerance, &truth); 627 AssertThrow(ierr == 0, ExcPETScError(ierr)); 628 return truth; 629 } 630 631 PetscBool is_hermitian(const double tolerance)632 MatrixBase::is_hermitian(const double tolerance) 633 { 634 PetscBool truth; 635 636 assert_is_compressed(); 637 const PetscErrorCode ierr = MatIsHermitian(matrix, tolerance, &truth); 638 AssertThrow(ierr == 0, ExcPETScError(ierr)); 639 640 return truth; 641 } 642 643 void write_ascii(const PetscViewerFormat format)644 MatrixBase::write_ascii(const PetscViewerFormat format) 645 { 646 assert_is_compressed(); 647 648 // Set options 649 PetscErrorCode ierr = 650 PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, format); 651 AssertThrow(ierr == 0, ExcPETScError(ierr)); 652 653 // Write to screen 654 ierr = MatView(matrix, PETSC_VIEWER_STDOUT_WORLD); 655 AssertThrow(ierr == 0, ExcPETScError(ierr)); 656 } 657 658 void print(std::ostream & out,const bool) const659 MatrixBase::print(std::ostream &out, const bool /*alternative_output*/) const 660 { 661 std::pair<MatrixBase::size_type, MatrixBase::size_type> loc_range = 662 local_range(); 663 664 PetscInt ncols; 665 const PetscInt * colnums; 666 const PetscScalar *values; 667 668 MatrixBase::size_type row; 669 for (row = loc_range.first; row < loc_range.second; ++row) 670 { 671 PetscErrorCode ierr = MatGetRow(*this, row, &ncols, &colnums, &values); 672 AssertThrow(ierr == 0, ExcPETScError(ierr)); 673 674 for (PetscInt col = 0; col < ncols; ++col) 675 { 676 out << "(" << row << "," << colnums[col] << ") " << values[col] 677 << std::endl; 678 } 679 680 ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values); 681 AssertThrow(ierr == 0, ExcPETScError(ierr)); 682 } 683 684 AssertThrow(out, ExcIO()); 685 } 686 687 688 689 std::size_t memory_consumption() const690 MatrixBase::memory_consumption() const 691 { 692 MatInfo info; 693 const PetscErrorCode ierr = MatGetInfo(matrix, MAT_LOCAL, &info); 694 AssertThrow(ierr == 0, ExcPETScError(ierr)); 695 696 return sizeof(*this) + static_cast<size_type>(info.memory); 697 } 698 699 } // namespace PETScWrappers 700 701 DEAL_II_NAMESPACE_CLOSE 702 703 #endif // DEAL_II_WITH_PETSC 704