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_vector_base.h> 17 18 #ifdef DEAL_II_WITH_PETSC 19 20 # include <deal.II/base/memory_consumption.h> 21 # include <deal.II/base/multithread_info.h> 22 23 # include <deal.II/lac/exceptions.h> 24 # include <deal.II/lac/petsc_compatibility.h> 25 # include <deal.II/lac/petsc_vector.h> 26 27 # include <cmath> 28 29 DEAL_II_NAMESPACE_OPEN 30 31 namespace PETScWrappers 32 { 33 namespace internal 34 { 35 # ifndef DOXYGEN operator PetscScalar() const36 VectorReference::operator PetscScalar() const 37 { 38 AssertIndexRange(index, vector.size()); 39 40 // The vector may have ghost entries. In that case, we first need to 41 // figure out which elements we own locally, then get a pointer to the 42 // elements that are stored here (both the ones we own as well as the 43 // ghost elements). In this array, the locally owned elements come first 44 // followed by the ghost elements whose position we can get from an 45 // index set. 46 if (vector.ghosted) 47 { 48 PetscInt begin, end; 49 PetscErrorCode ierr = 50 VecGetOwnershipRange(vector.vector, &begin, &end); 51 AssertThrow(ierr == 0, ExcPETScError(ierr)); 52 53 Vec locally_stored_elements = PETSC_NULL; 54 ierr = VecGhostGetLocalForm(vector.vector, &locally_stored_elements); 55 AssertThrow(ierr == 0, ExcPETScError(ierr)); 56 57 PetscInt lsize; 58 ierr = VecGetSize(locally_stored_elements, &lsize); 59 AssertThrow(ierr == 0, ExcPETScError(ierr)); 60 61 PetscScalar *ptr; 62 ierr = VecGetArray(locally_stored_elements, &ptr); 63 AssertThrow(ierr == 0, ExcPETScError(ierr)); 64 65 PetscScalar value; 66 67 if (index >= static_cast<size_type>(begin) && 68 index < static_cast<size_type>(end)) 69 { 70 // local entry 71 value = *(ptr + index - begin); 72 } 73 else 74 { 75 // ghost entry 76 const size_type ghostidx = 77 vector.ghost_indices.index_within_set(index); 78 79 AssertIndexRange(ghostidx + end - begin, lsize); 80 value = *(ptr + ghostidx + end - begin); 81 } 82 83 ierr = VecRestoreArray(locally_stored_elements, &ptr); 84 AssertThrow(ierr == 0, ExcPETScError(ierr)); 85 86 ierr = 87 VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements); 88 AssertThrow(ierr == 0, ExcPETScError(ierr)); 89 90 return value; 91 } 92 93 94 // first verify that the requested 95 // element is actually locally 96 // available 97 PetscInt begin, end; 98 99 PetscErrorCode ierr = VecGetOwnershipRange(vector.vector, &begin, &end); 100 AssertThrow(ierr == 0, ExcPETScError(ierr)); 101 102 AssertThrow((index >= static_cast<size_type>(begin)) && 103 (index < static_cast<size_type>(end)), 104 ExcAccessToNonlocalElement(index, begin, end - 1)); 105 106 PetscInt idx = index; 107 PetscScalar value; 108 ierr = VecGetValues(vector.vector, 1, &idx, &value); 109 AssertThrow(ierr == 0, ExcPETScError(ierr)); 110 111 return value; 112 } 113 # endif 114 } // namespace internal 115 VectorBase()116 VectorBase::VectorBase() 117 : vector(nullptr) 118 , ghosted(false) 119 , last_action(::dealii::VectorOperation::unknown) 120 , obtained_ownership(true) 121 { 122 Assert(MultithreadInfo::is_running_single_threaded(), 123 ExcMessage("PETSc does not support multi-threaded access, set " 124 "the thread limit to 1 in MPI_InitFinalize().")); 125 } 126 127 128 VectorBase(const VectorBase & v)129 VectorBase::VectorBase(const VectorBase &v) 130 : Subscriptor() 131 , ghosted(v.ghosted) 132 , ghost_indices(v.ghost_indices) 133 , last_action(::dealii::VectorOperation::unknown) 134 , obtained_ownership(true) 135 { 136 Assert(MultithreadInfo::is_running_single_threaded(), 137 ExcMessage("PETSc does not support multi-threaded access, set " 138 "the thread limit to 1 in MPI_InitFinalize().")); 139 140 PetscErrorCode ierr = VecDuplicate(v.vector, &vector); 141 AssertThrow(ierr == 0, ExcPETScError(ierr)); 142 143 ierr = VecCopy(v.vector, vector); 144 AssertThrow(ierr == 0, ExcPETScError(ierr)); 145 } 146 147 148 VectorBase(const Vec & v)149 VectorBase::VectorBase(const Vec &v) 150 : Subscriptor() 151 , vector(v) 152 , ghosted(false) 153 , last_action(::dealii::VectorOperation::unknown) 154 , obtained_ownership(false) 155 { 156 Assert(MultithreadInfo::is_running_single_threaded(), 157 ExcMessage("PETSc does not support multi-threaded access, set " 158 "the thread limit to 1 in MPI_InitFinalize().")); 159 } 160 161 162 ~VectorBase()163 VectorBase::~VectorBase() 164 { 165 if (obtained_ownership) 166 { 167 const PetscErrorCode ierr = VecDestroy(&vector); 168 AssertNothrow(ierr == 0, ExcPETScError(ierr)); 169 (void)ierr; 170 } 171 } 172 173 174 175 void clear()176 VectorBase::clear() 177 { 178 if (obtained_ownership) 179 { 180 const PetscErrorCode ierr = VecDestroy(&vector); 181 AssertThrow(ierr == 0, ExcPETScError(ierr)); 182 } 183 184 ghosted = false; 185 ghost_indices.clear(); 186 last_action = ::dealii::VectorOperation::unknown; 187 obtained_ownership = true; 188 } 189 190 191 192 VectorBase & operator =(const PetscScalar s)193 VectorBase::operator=(const PetscScalar s) 194 { 195 AssertIsFinite(s); 196 197 // TODO[TH]: assert(is_compressed()) 198 199 PetscErrorCode ierr = VecSet(vector, s); 200 AssertThrow(ierr == 0, ExcPETScError(ierr)); 201 202 if (has_ghost_elements()) 203 { 204 Vec ghost = PETSC_NULL; 205 ierr = VecGhostGetLocalForm(vector, &ghost); 206 AssertThrow(ierr == 0, ExcPETScError(ierr)); 207 208 ierr = VecSet(ghost, s); 209 AssertThrow(ierr == 0, ExcPETScError(ierr)); 210 211 ierr = VecGhostRestoreLocalForm(vector, &ghost); 212 AssertThrow(ierr == 0, ExcPETScError(ierr)); 213 } 214 215 return *this; 216 } 217 218 219 220 bool operator ==(const VectorBase & v) const221 VectorBase::operator==(const VectorBase &v) const 222 { 223 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size())); 224 225 PetscBool flag; 226 const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag); 227 AssertThrow(ierr == 0, ExcPETScError(ierr)); 228 229 return (flag == PETSC_TRUE); 230 } 231 232 233 234 bool operator !=(const VectorBase & v) const235 VectorBase::operator!=(const VectorBase &v) const 236 { 237 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size())); 238 239 PetscBool flag; 240 const PetscErrorCode ierr = VecEqual(vector, v.vector, &flag); 241 AssertThrow(ierr == 0, ExcPETScError(ierr)); 242 243 return (flag == PETSC_FALSE); 244 } 245 246 247 248 VectorBase::size_type size() const249 VectorBase::size() const 250 { 251 PetscInt sz; 252 const PetscErrorCode ierr = VecGetSize(vector, &sz); 253 AssertThrow(ierr == 0, ExcPETScError(ierr)); 254 255 return sz; 256 } 257 258 259 260 VectorBase::size_type local_size() const261 VectorBase::local_size() const 262 { 263 PetscInt sz; 264 const PetscErrorCode ierr = VecGetLocalSize(vector, &sz); 265 AssertThrow(ierr == 0, ExcPETScError(ierr)); 266 267 return sz; 268 } 269 270 271 272 std::pair<VectorBase::size_type, VectorBase::size_type> local_range() const273 VectorBase::local_range() const 274 { 275 PetscInt begin, end; 276 const PetscErrorCode ierr = 277 VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end); 278 AssertThrow(ierr == 0, ExcPETScError(ierr)); 279 280 return std::make_pair(begin, end); 281 } 282 283 284 285 void set(const std::vector<size_type> & indices,const std::vector<PetscScalar> & values)286 VectorBase::set(const std::vector<size_type> & indices, 287 const std::vector<PetscScalar> &values) 288 { 289 Assert(indices.size() == values.size(), 290 ExcMessage("Function called with arguments of different sizes")); 291 do_set_add_operation(indices.size(), indices.data(), values.data(), false); 292 } 293 294 295 296 void add(const std::vector<size_type> & indices,const std::vector<PetscScalar> & values)297 VectorBase::add(const std::vector<size_type> & indices, 298 const std::vector<PetscScalar> &values) 299 { 300 Assert(indices.size() == values.size(), 301 ExcMessage("Function called with arguments of different sizes")); 302 do_set_add_operation(indices.size(), indices.data(), values.data(), true); 303 } 304 305 306 307 void add(const std::vector<size_type> & indices,const::dealii::Vector<PetscScalar> & values)308 VectorBase::add(const std::vector<size_type> & indices, 309 const ::dealii::Vector<PetscScalar> &values) 310 { 311 Assert(indices.size() == values.size(), 312 ExcMessage("Function called with arguments of different sizes")); 313 do_set_add_operation(indices.size(), indices.data(), values.begin(), true); 314 } 315 316 317 318 void add(const size_type n_elements,const size_type * indices,const PetscScalar * values)319 VectorBase::add(const size_type n_elements, 320 const size_type * indices, 321 const PetscScalar *values) 322 { 323 do_set_add_operation(n_elements, indices, values, true); 324 } 325 326 327 operator *(const VectorBase & vec) const328 PetscScalar VectorBase::operator*(const VectorBase &vec) const 329 { 330 Assert(size() == vec.size(), ExcDimensionMismatch(size(), vec.size())); 331 332 PetscScalar result; 333 334 // For complex vectors, VecDot() computes 335 // val = (x,y) = y^H x, 336 // where y^H denotes the conjugate transpose of y. 337 // Note that this corresponds to the usual "mathematicians'" 338 // complex inner product where the SECOND argument gets the 339 // complex conjugate, which is also how we document this 340 // operation. 341 const PetscErrorCode ierr = VecDot(vec.vector, vector, &result); 342 AssertThrow(ierr == 0, ExcPETScError(ierr)); 343 344 return result; 345 } 346 347 348 349 PetscScalar add_and_dot(const PetscScalar a,const VectorBase & V,const VectorBase & W)350 VectorBase::add_and_dot(const PetscScalar a, 351 const VectorBase &V, 352 const VectorBase &W) 353 { 354 this->add(a, V); 355 return *this * W; 356 } 357 358 359 360 void compress(const VectorOperation::values operation)361 VectorBase::compress(const VectorOperation::values operation) 362 { 363 { 364 # ifdef DEBUG 365 # ifdef DEAL_II_WITH_MPI 366 // Check that all processors agree that last_action is the same (or none!) 367 368 int my_int_last_action = last_action; 369 int all_int_last_action; 370 371 const int ierr = MPI_Allreduce(&my_int_last_action, 372 &all_int_last_action, 373 1, 374 MPI_INT, 375 MPI_BOR, 376 get_mpi_communicator()); 377 AssertThrowMPI(ierr); 378 379 AssertThrow(all_int_last_action != (::dealii::VectorOperation::add | 380 ::dealii::VectorOperation::insert), 381 ExcMessage("Error: not all processors agree on the last " 382 "VectorOperation before this compress() call.")); 383 # endif 384 # endif 385 } 386 387 AssertThrow( 388 last_action == ::dealii::VectorOperation::unknown || 389 last_action == operation, 390 ExcMessage( 391 "Missing compress() or calling with wrong VectorOperation argument.")); 392 393 // note that one may think that 394 // we only need to do something 395 // if in fact the state is 396 // anything but 397 // last_action::unknown. but 398 // that's not true: one 399 // frequently gets into 400 // situations where only one 401 // processor (or a subset of 402 // processors) actually writes 403 // something into a vector, but 404 // we still need to call 405 // VecAssemblyBegin/End on all 406 // processors. 407 PetscErrorCode ierr = VecAssemblyBegin(vector); 408 AssertThrow(ierr == 0, ExcPETScError(ierr)); 409 ierr = VecAssemblyEnd(vector); 410 AssertThrow(ierr == 0, ExcPETScError(ierr)); 411 412 // reset the last action field to 413 // indicate that we're back to a 414 // pristine state 415 last_action = ::dealii::VectorOperation::unknown; 416 } 417 418 419 420 VectorBase::real_type norm_sqr() const421 VectorBase::norm_sqr() const 422 { 423 const real_type d = l2_norm(); 424 return d * d; 425 } 426 427 428 429 PetscScalar mean_value() const430 VectorBase::mean_value() const 431 { 432 // We can only use our more efficient 433 // routine in the serial case. 434 if (dynamic_cast<const PETScWrappers::MPI::Vector *>(this) != nullptr) 435 { 436 PetscScalar sum; 437 const PetscErrorCode ierr = VecSum(vector, &sum); 438 AssertThrow(ierr == 0, ExcPETScError(ierr)); 439 return sum / static_cast<PetscReal>(size()); 440 } 441 442 // get a representation of the vector and 443 // loop over all the elements 444 PetscScalar * start_ptr; 445 PetscErrorCode ierr = VecGetArray(vector, &start_ptr); 446 AssertThrow(ierr == 0, ExcPETScError(ierr)); 447 448 PetscScalar mean = 0; 449 { 450 PetscScalar sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0; 451 452 // use modern processors better by 453 // allowing pipelined commands to be 454 // executed in parallel 455 const PetscScalar *ptr = start_ptr; 456 const PetscScalar *eptr = ptr + (size() / 4) * 4; 457 while (ptr != eptr) 458 { 459 sum0 += *ptr++; 460 sum1 += *ptr++; 461 sum2 += *ptr++; 462 sum3 += *ptr++; 463 } 464 // add up remaining elements 465 while (ptr != start_ptr + size()) 466 sum0 += *ptr++; 467 468 mean = (sum0 + sum1 + sum2 + sum3) / static_cast<PetscReal>(size()); 469 } 470 471 // restore the representation of the 472 // vector 473 ierr = VecRestoreArray(vector, &start_ptr); 474 AssertThrow(ierr == 0, ExcPETScError(ierr)); 475 476 return mean; 477 } 478 479 480 VectorBase::real_type l1_norm() const481 VectorBase::l1_norm() const 482 { 483 real_type d; 484 485 const PetscErrorCode ierr = VecNorm(vector, NORM_1, &d); 486 AssertThrow(ierr == 0, ExcPETScError(ierr)); 487 488 return d; 489 } 490 491 492 493 VectorBase::real_type l2_norm() const494 VectorBase::l2_norm() const 495 { 496 real_type d; 497 498 const PetscErrorCode ierr = VecNorm(vector, NORM_2, &d); 499 AssertThrow(ierr == 0, ExcPETScError(ierr)); 500 501 return d; 502 } 503 504 505 506 VectorBase::real_type lp_norm(const real_type p) const507 VectorBase::lp_norm(const real_type p) const 508 { 509 // get a representation of the vector and 510 // loop over all the elements 511 PetscScalar * start_ptr; 512 PetscErrorCode ierr = VecGetArray(vector, &start_ptr); 513 AssertThrow(ierr == 0, ExcPETScError(ierr)); 514 515 real_type norm = 0; 516 { 517 real_type sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0; 518 519 // use modern processors better by 520 // allowing pipelined commands to be 521 // executed in parallel 522 const PetscScalar *ptr = start_ptr; 523 const PetscScalar *eptr = ptr + (size() / 4) * 4; 524 while (ptr != eptr) 525 { 526 sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p); 527 sum1 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p); 528 sum2 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p); 529 sum3 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p); 530 } 531 // add up remaining elements 532 while (ptr != start_ptr + size()) 533 sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p); 534 535 norm = std::pow(sum0 + sum1 + sum2 + sum3, 1. / p); 536 } 537 538 // restore the representation of the 539 // vector 540 ierr = VecRestoreArray(vector, &start_ptr); 541 AssertThrow(ierr == 0, ExcPETScError(ierr)); 542 543 return norm; 544 } 545 546 547 548 VectorBase::real_type linfty_norm() const549 VectorBase::linfty_norm() const 550 { 551 real_type d; 552 553 const PetscErrorCode ierr = VecNorm(vector, NORM_INFINITY, &d); 554 AssertThrow(ierr == 0, ExcPETScError(ierr)); 555 556 return d; 557 } 558 559 560 561 VectorBase::real_type min() const562 VectorBase::min() const 563 { 564 PetscInt p; 565 real_type d; 566 567 const PetscErrorCode ierr = VecMin(vector, &p, &d); 568 AssertThrow(ierr == 0, ExcPETScError(ierr)); 569 570 return d; 571 } 572 573 574 VectorBase::real_type max() const575 VectorBase::max() const 576 { 577 PetscInt p; 578 real_type d; 579 580 const PetscErrorCode ierr = VecMax(vector, &p, &d); 581 AssertThrow(ierr == 0, ExcPETScError(ierr)); 582 583 return d; 584 } 585 586 587 588 bool all_zero() const589 VectorBase::all_zero() const 590 { 591 // get a representation of the vector and 592 // loop over all the elements 593 PetscScalar * start_ptr; 594 PetscErrorCode ierr = VecGetArray(vector, &start_ptr); 595 AssertThrow(ierr == 0, ExcPETScError(ierr)); 596 597 const PetscScalar *ptr = start_ptr, *eptr = start_ptr + local_size(); 598 bool flag = true; 599 while (ptr != eptr) 600 { 601 if (*ptr != value_type()) 602 { 603 flag = false; 604 break; 605 } 606 ++ptr; 607 } 608 609 // restore the representation of the 610 // vector 611 ierr = VecRestoreArray(vector, &start_ptr); 612 AssertThrow(ierr == 0, ExcPETScError(ierr)); 613 614 return flag; 615 } 616 617 618 namespace internal 619 { 620 template <typename T> 621 bool is_non_negative(const T & t)622 is_non_negative(const T &t) 623 { 624 return t >= 0; 625 } 626 627 628 629 template <typename T> 630 bool is_non_negative(const std::complex<T> &)631 is_non_negative(const std::complex<T> &) 632 { 633 Assert(false, 634 ExcMessage("You can't ask a complex value " 635 "whether it is non-negative.")) return true; 636 } 637 } // namespace internal 638 639 640 641 bool is_non_negative() const642 VectorBase::is_non_negative() const 643 { 644 // get a representation of the vector and 645 // loop over all the elements 646 PetscScalar * start_ptr; 647 PetscErrorCode ierr = VecGetArray(vector, &start_ptr); 648 AssertThrow(ierr == 0, ExcPETScError(ierr)); 649 650 const PetscScalar *ptr = start_ptr, *eptr = start_ptr + local_size(); 651 bool flag = true; 652 while (ptr != eptr) 653 { 654 if (!internal::is_non_negative(*ptr)) 655 { 656 flag = false; 657 break; 658 } 659 ++ptr; 660 } 661 662 // restore the representation of the 663 // vector 664 ierr = VecRestoreArray(vector, &start_ptr); 665 AssertThrow(ierr == 0, ExcPETScError(ierr)); 666 667 return flag; 668 } 669 670 671 672 VectorBase & operator *=(const PetscScalar a)673 VectorBase::operator*=(const PetscScalar a) 674 { 675 Assert(!has_ghost_elements(), ExcGhostsPresent()); 676 AssertIsFinite(a); 677 678 const PetscErrorCode ierr = VecScale(vector, a); 679 AssertThrow(ierr == 0, ExcPETScError(ierr)); 680 681 return *this; 682 } 683 684 685 686 VectorBase & operator /=(const PetscScalar a)687 VectorBase::operator/=(const PetscScalar a) 688 { 689 Assert(!has_ghost_elements(), ExcGhostsPresent()); 690 AssertIsFinite(a); 691 692 const PetscScalar factor = 1. / a; 693 AssertIsFinite(factor); 694 695 const PetscErrorCode ierr = VecScale(vector, factor); 696 AssertThrow(ierr == 0, ExcPETScError(ierr)); 697 698 return *this; 699 } 700 701 702 703 VectorBase & operator +=(const VectorBase & v)704 VectorBase::operator+=(const VectorBase &v) 705 { 706 Assert(!has_ghost_elements(), ExcGhostsPresent()); 707 const PetscErrorCode ierr = VecAXPY(vector, 1, v); 708 AssertThrow(ierr == 0, ExcPETScError(ierr)); 709 710 return *this; 711 } 712 713 714 715 VectorBase & operator -=(const VectorBase & v)716 VectorBase::operator-=(const VectorBase &v) 717 { 718 Assert(!has_ghost_elements(), ExcGhostsPresent()); 719 const PetscErrorCode ierr = VecAXPY(vector, -1, v); 720 AssertThrow(ierr == 0, ExcPETScError(ierr)); 721 722 return *this; 723 } 724 725 726 727 void add(const PetscScalar s)728 VectorBase::add(const PetscScalar s) 729 { 730 Assert(!has_ghost_elements(), ExcGhostsPresent()); 731 AssertIsFinite(s); 732 733 const PetscErrorCode ierr = VecShift(vector, s); 734 AssertThrow(ierr == 0, ExcPETScError(ierr)); 735 } 736 737 738 739 void add(const PetscScalar a,const VectorBase & v)740 VectorBase::add(const PetscScalar a, const VectorBase &v) 741 { 742 Assert(!has_ghost_elements(), ExcGhostsPresent()); 743 AssertIsFinite(a); 744 745 const PetscErrorCode ierr = VecAXPY(vector, a, v); 746 AssertThrow(ierr == 0, ExcPETScError(ierr)); 747 } 748 749 750 751 void add(const PetscScalar a,const VectorBase & v,const PetscScalar b,const VectorBase & w)752 VectorBase::add(const PetscScalar a, 753 const VectorBase &v, 754 const PetscScalar b, 755 const VectorBase &w) 756 { 757 Assert(!has_ghost_elements(), ExcGhostsPresent()); 758 AssertIsFinite(a); 759 AssertIsFinite(b); 760 761 const PetscScalar weights[2] = {a, b}; 762 Vec addends[2] = {v.vector, w.vector}; 763 764 const PetscErrorCode ierr = VecMAXPY(vector, 2, weights, addends); 765 AssertThrow(ierr == 0, ExcPETScError(ierr)); 766 } 767 768 769 770 void sadd(const PetscScalar s,const VectorBase & v)771 VectorBase::sadd(const PetscScalar s, const VectorBase &v) 772 { 773 Assert(!has_ghost_elements(), ExcGhostsPresent()); 774 AssertIsFinite(s); 775 776 const PetscErrorCode ierr = VecAYPX(vector, s, v); 777 AssertThrow(ierr == 0, ExcPETScError(ierr)); 778 } 779 780 781 782 void sadd(const PetscScalar s,const PetscScalar a,const VectorBase & v)783 VectorBase::sadd(const PetscScalar s, 784 const PetscScalar a, 785 const VectorBase &v) 786 { 787 Assert(!has_ghost_elements(), ExcGhostsPresent()); 788 AssertIsFinite(s); 789 AssertIsFinite(a); 790 791 // there is nothing like a AXPAY 792 // operation in Petsc, so do it in two 793 // steps 794 *this *= s; 795 add(a, v); 796 } 797 798 799 800 void scale(const VectorBase & factors)801 VectorBase::scale(const VectorBase &factors) 802 { 803 Assert(!has_ghost_elements(), ExcGhostsPresent()); 804 const PetscErrorCode ierr = VecPointwiseMult(vector, factors, vector); 805 AssertThrow(ierr == 0, ExcPETScError(ierr)); 806 } 807 808 809 810 void equ(const PetscScalar a,const VectorBase & v)811 VectorBase::equ(const PetscScalar a, const VectorBase &v) 812 { 813 Assert(!has_ghost_elements(), ExcGhostsPresent()); 814 AssertIsFinite(a); 815 816 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size())); 817 818 // there is no simple operation for this 819 // in PETSc. there are multiple ways to 820 // emulate it, we choose this one: 821 const PetscErrorCode ierr = VecCopy(v.vector, vector); 822 AssertThrow(ierr == 0, ExcPETScError(ierr)); 823 824 *this *= a; 825 } 826 827 828 829 void write_ascii(const PetscViewerFormat format)830 VectorBase::write_ascii(const PetscViewerFormat format) 831 { 832 // TODO[TH]:assert(is_compressed()) 833 834 // Set options 835 PetscErrorCode ierr = 836 PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, format); 837 AssertThrow(ierr == 0, ExcPETScError(ierr)); 838 839 // Write to screen 840 ierr = VecView(vector, PETSC_VIEWER_STDOUT_WORLD); 841 AssertThrow(ierr == 0, ExcPETScError(ierr)); 842 } 843 844 845 846 void print(std::ostream & out,const unsigned int precision,const bool scientific,const bool across) const847 VectorBase::print(std::ostream & out, 848 const unsigned int precision, 849 const bool scientific, 850 const bool across) const 851 { 852 AssertThrow(out, ExcIO()); 853 854 // get a representation of the vector and 855 // loop over all the elements 856 PetscScalar * val; 857 PetscErrorCode ierr = VecGetArray(vector, &val); 858 859 AssertThrow(ierr == 0, ExcPETScError(ierr)); 860 861 // save the state of out stream 862 const std::ios::fmtflags old_flags = out.flags(); 863 const unsigned int old_precision = out.precision(precision); 864 865 out.precision(precision); 866 if (scientific) 867 out.setf(std::ios::scientific, std::ios::floatfield); 868 else 869 out.setf(std::ios::fixed, std::ios::floatfield); 870 871 if (across) 872 for (size_type i = 0; i < local_size(); ++i) 873 out << val[i] << ' '; 874 else 875 for (size_type i = 0; i < local_size(); ++i) 876 out << val[i] << std::endl; 877 out << std::endl; 878 879 // reset output format 880 out.flags(old_flags); 881 out.precision(old_precision); 882 883 // restore the representation of the 884 // vector 885 ierr = VecRestoreArray(vector, &val); 886 AssertThrow(ierr == 0, ExcPETScError(ierr)); 887 888 AssertThrow(out, ExcIO()); 889 } 890 891 892 893 void swap(VectorBase & v)894 VectorBase::swap(VectorBase &v) 895 { 896 const PetscErrorCode ierr = VecSwap(vector, v.vector); 897 AssertThrow(ierr == 0, ExcPETScError(ierr)); 898 } 899 900 901 operator const Vec&() const902 VectorBase::operator const Vec &() const 903 { 904 return vector; 905 } 906 907 908 std::size_t memory_consumption() const909 VectorBase::memory_consumption() const 910 { 911 std::size_t mem = sizeof(Vec) + sizeof(last_action) + 912 MemoryConsumption::memory_consumption(ghosted) + 913 MemoryConsumption::memory_consumption(ghost_indices); 914 915 // TH: I am relatively sure that PETSc is 916 // storing the local data in a contiguous 917 // block without indices: 918 mem += local_size() * sizeof(PetscScalar); 919 // assume that PETSc is storing one index 920 // and one double per ghost element 921 if (ghosted) 922 mem += ghost_indices.n_elements() * (sizeof(PetscScalar) + sizeof(int)); 923 924 // TODO[TH]: size of constant memory for PETSc? 925 return mem; 926 } 927 928 929 930 void do_set_add_operation(const size_type n_elements,const size_type * indices,const PetscScalar * values,const bool add_values)931 VectorBase::do_set_add_operation(const size_type n_elements, 932 const size_type * indices, 933 const PetscScalar *values, 934 const bool add_values) 935 { 936 ::dealii::VectorOperation::values action = 937 (add_values ? ::dealii::VectorOperation::add : 938 ::dealii::VectorOperation::insert); 939 Assert((last_action == action) || 940 (last_action == ::dealii::VectorOperation::unknown), 941 internal::VectorReference::ExcWrongMode(action, last_action)); 942 Assert(!has_ghost_elements(), ExcGhostsPresent()); 943 // VecSetValues complains if we 944 // come with an empty 945 // vector. however, it is not a 946 // collective operation, so we 947 // can skip the call if necessary 948 // (unlike the above calls) 949 if (n_elements != 0) 950 { 951 const PetscInt *petsc_indices = 952 reinterpret_cast<const PetscInt *>(indices); 953 954 const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES); 955 const PetscErrorCode ierr = 956 VecSetValues(vector, n_elements, petsc_indices, values, mode); 957 AssertThrow(ierr == 0, ExcPETScError(ierr)); 958 } 959 960 // set the mode here, independent of whether we have actually 961 // written elements or whether the list was empty 962 last_action = action; 963 } 964 965 } // namespace PETScWrappers 966 967 DEAL_II_NAMESPACE_CLOSE 968 969 #endif // DEAL_II_WITH_PETSC 970