1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2008 - 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 #ifndef dealii_trilinos_vector_h 17 #define dealii_trilinos_vector_h 18 19 20 #include <deal.II/base/config.h> 21 22 #ifdef DEAL_II_WITH_TRILINOS 23 # include <deal.II/base/index_set.h> 24 # include <deal.II/base/mpi.h> 25 # include <deal.II/base/subscriptor.h> 26 # include <deal.II/base/utilities.h> 27 28 # include <deal.II/lac/exceptions.h> 29 # include <deal.II/lac/vector.h> 30 # include <deal.II/lac/vector_operation.h> 31 # include <deal.II/lac/vector_type_traits.h> 32 33 # include <Epetra_ConfigDefs.h> 34 35 # include <memory> 36 # include <utility> 37 # include <vector> 38 # ifdef DEAL_II_WITH_MPI // only if MPI is installed 39 # include <Epetra_MpiComm.h> 40 # include <mpi.h> 41 # else 42 # include <Epetra_SerialComm.h> 43 # endif 44 # include <Epetra_FEVector.h> 45 # include <Epetra_LocalMap.h> 46 # include <Epetra_Map.h> 47 48 DEAL_II_NAMESPACE_OPEN 49 50 // Forward declarations 51 # ifndef DOXYGEN 52 namespace LinearAlgebra 53 { 54 // Forward declaration 55 template <typename Number> 56 class ReadWriteVector; 57 } // namespace LinearAlgebra 58 # endif 59 60 /** 61 * @addtogroup TrilinosWrappers 62 * @{ 63 */ 64 65 /** 66 * A namespace in which wrapper classes for Trilinos objects reside. 67 * 68 * @ingroup TrilinosWrappers 69 */ 70 namespace TrilinosWrappers 71 { 72 class SparseMatrix; 73 74 /** 75 * @cond internal 76 */ 77 78 /** 79 * A namespace for internal implementation details of the TrilinosWrapper 80 * members. 81 * 82 * @ingroup TrilinosWrappers 83 */ 84 namespace internal 85 { 86 /** 87 * Declare type for container size. 88 */ 89 using size_type = dealii::types::global_dof_index; 90 91 /** 92 * This class implements a wrapper for accessing the Trilinos vector in 93 * the same way as we access deal.II objects: it is initialized with a 94 * vector and an element within it, and has a conversion operator to 95 * extract the scalar value of this element. It also has a variety of 96 * assignment operator for writing to this one element. 97 * 98 * @ingroup TrilinosWrappers 99 */ 100 class VectorReference 101 { 102 private: 103 /** 104 * Constructor. It is made private so as to only allow the actual vector 105 * class to create it. 106 */ 107 VectorReference(MPI::Vector &vector, const size_type index); 108 109 public: 110 /** 111 * Copy constructor. 112 */ 113 VectorReference(const VectorReference &) = default; 114 115 /** 116 * This looks like a copy operator, but does something different than 117 * usual. In particular, it does not copy the member variables of this 118 * reference. Rather, it handles the situation where we have two vectors 119 * @p v and @p w, and assign elements like in <tt>v(i)=w(i)</tt>. Here, 120 * both left and right hand side of the assignment have data type 121 * VectorReference, but what we really mean is to assign the vector 122 * elements represented by the two references. This operator implements 123 * this operation. Note also that this allows us to make the assignment 124 * operator const. 125 */ 126 const VectorReference & 127 operator=(const VectorReference &r) const; 128 129 /** 130 * Same as above but for non-const reference objects. 131 */ 132 VectorReference & 133 operator=(const VectorReference &r); 134 135 /** 136 * Set the referenced element of the vector to <tt>s</tt>. 137 */ 138 const VectorReference & 139 operator=(const TrilinosScalar &s) const; 140 141 /** 142 * Add <tt>s</tt> to the referenced element of the vector-> 143 */ 144 const VectorReference & 145 operator+=(const TrilinosScalar &s) const; 146 147 /** 148 * Subtract <tt>s</tt> from the referenced element of the vector-> 149 */ 150 const VectorReference & 151 operator-=(const TrilinosScalar &s) const; 152 153 /** 154 * Multiply the referenced element of the vector by <tt>s</tt>. 155 */ 156 const VectorReference & 157 operator*=(const TrilinosScalar &s) const; 158 159 /** 160 * Divide the referenced element of the vector by <tt>s</tt>. 161 */ 162 const VectorReference & 163 operator/=(const TrilinosScalar &s) const; 164 165 /** 166 * Convert the reference to an actual value, i.e. return the value of 167 * the referenced element of the vector. 168 */ 169 operator TrilinosScalar() const; 170 171 /** 172 * Exception 173 */ 174 DeclException1(ExcTrilinosError, 175 int, 176 << "An error with error number " << arg1 177 << " occurred while calling a Trilinos function"); 178 179 private: 180 /** 181 * Point to the vector we are referencing. 182 */ 183 MPI::Vector &vector; 184 185 /** 186 * Index of the referenced element of the vector. 187 */ 188 const size_type index; 189 190 // Make the vector class a friend, so that it can create objects of the 191 // present type. 192 friend class ::dealii::TrilinosWrappers::MPI::Vector; 193 }; 194 } // namespace internal 195 /** 196 * @endcond 197 */ 198 199 # ifndef DEAL_II_WITH_64BIT_INDICES 200 // define a helper function that queries the global ID of local ID of 201 // an Epetra_BlockMap object by calling either the 32- or 64-bit 202 // function necessary. 203 inline int gid(const Epetra_BlockMap & map,int i)204 gid(const Epetra_BlockMap &map, int i) 205 { 206 return map.GID(i); 207 } 208 # else 209 // define a helper function that queries the global ID of local ID of 210 // an Epetra_BlockMap object by calling either the 32- or 64-bit 211 // function necessary. 212 inline long long int gid(const Epetra_BlockMap & map,int i)213 gid(const Epetra_BlockMap &map, int i) 214 { 215 return map.GID64(i); 216 } 217 # endif 218 219 /** 220 * Namespace for Trilinos vector classes that work in parallel over MPI. 221 * 222 * @ingroup TrilinosWrappers 223 */ 224 namespace MPI 225 { 226 class BlockVector; 227 228 /** 229 * This class implements a wrapper to use the Trilinos distributed vector 230 * class Epetra_FEVector, the (parallel) partitioning of which 231 * is governed by an Epetra_Map. 232 * The Epetra_FEVector is precisely the kind of vector 233 * we deal with all the time - we probably get it from some assembly 234 * process, where also entries not locally owned might need to written and 235 * hence need to be forwarded to the owner. 236 * 237 * The interface of this class is modeled after the existing Vector class in 238 * deal.II. It has almost the same member functions, and is often 239 * exchangeable. However, since Trilinos only supports a single scalar type 240 * (double), it is not templated, and only works with that type. 241 * 242 * Note that Trilinos only guarantees that operations do what you expect 243 * if the function @p GlobalAssemble has been called after vector assembly 244 * in order to distribute the data. This is necessary since some processes 245 * might have accumulated data of elements that are not owned by 246 * themselves, but must be sent to the owning process. In order to avoid 247 * using the wrong data, you need to call Vector::compress() before you 248 * actually use the vectors. 249 * 250 * <h3>Parallel communication model</h3> 251 * 252 * The parallel functionality of Trilinos is built on top of the Message 253 * Passing Interface (MPI). MPI's communication model is built on 254 * collective communications: if one process wants something from another, 255 * that other process has to be willing to accept this communication. A 256 * process cannot query data from another process by calling a remote 257 * function, without that other process expecting such a transaction. The 258 * consequence is that most of the operations in the base class of this 259 * class have to be called collectively. For example, if you want to 260 * compute the l2 norm of a parallel vector, @em all processes across 261 * which this vector is shared have to call the @p l2_norm function. If 262 * you don't do this, but instead only call the @p l2_norm function on one 263 * process, then the following happens: This one process will call one of 264 * the collective MPI functions and wait for all the other processes to 265 * join in on this. Since the other processes don't call this function, 266 * you will either get a time-out on the first process, or, worse, by the 267 * time the next a call to a Trilinos function generates an MPI message on 268 * the other processes, you will get a cryptic message that only a subset 269 * of processes attempted a communication. These bugs can be very hard to 270 * figure out, unless you are well-acquainted with the communication model 271 * of MPI, and know which functions may generate MPI messages. 272 * 273 * One particular case, where an MPI message may be generated unexpectedly 274 * is discussed below. 275 * 276 * 277 * <h3>Accessing individual elements of a vector</h3> 278 * 279 * Trilinos does of course allow read access to individual elements of a 280 * vector, but in the distributed case only to elements that are stored 281 * locally. We implement this through calls like <tt>d=vec(i)</tt>. 282 * However, if you access an element outside the locally stored range, an 283 * exception is generated. 284 * 285 * In contrast to read access, Trilinos (and the respective deal.II 286 * wrapper classes) allow to write (or add) to individual elements of 287 * vectors, even if they are stored on a different process. You can do 288 * this by writing into or adding to elements using the syntax 289 * <tt>vec(i)=d</tt> or <tt>vec(i)+=d</tt>, or similar operations. There 290 * is one catch, however, that may lead to very confusing error messages: 291 * Trilinos requires application programs to call the compress() function 292 * when they switch from performing a set of operations that add to 293 * elements, to performing a set of operations that write to elements. The 294 * reasoning is that all processes might accumulate addition operations to 295 * elements, even if multiple processes write to the same elements. By the 296 * time we call compress() the next time, all these additions are 297 * executed. However, if one process adds to an element, and another 298 * overwrites to it, the order of execution would yield non-deterministic 299 * behavior if we don't make sure that a synchronization with compress() 300 * happens in between. 301 * 302 * In order to make sure these calls to compress() happen at the 303 * appropriate time, the deal.II wrappers keep a state variable that store 304 * which is the presently allowed operation: additions or writes. If it 305 * encounters an operation of the opposite kind, it calls compress() and 306 * flips the state. This can sometimes lead to very confusing behavior, in 307 * code that may for example look like this: 308 * 309 * @code 310 * TrilinosWrappers::MPI::Vector vector; 311 * // do some write operations on the vector 312 * for (size_type i=0; i<vector->size(); ++i) 313 * vector(i) = i; 314 * 315 * // do some additions to vector elements, but 316 * // only for some elements 317 * for (size_type i=0; i<vector->size(); ++i) 318 * if (some_condition(i) == true) 319 * vector(i) += 1; 320 * 321 * // do another collective operation 322 * const double norm = vector->l2_norm(); 323 * @endcode 324 * 325 * This code can run into trouble: by the time we see the first addition 326 * operation, we need to flush the overwrite buffers for the vector, and 327 * the deal.II library will do so by calling compress(). However, it will 328 * only do so for all processes that actually do an addition -- if the 329 * condition is never true for one of the processes, then this one will 330 * not get to the actual compress() call, whereas all the other ones do. 331 * This gets us into trouble, since all the other processes hang in the 332 * call to flush the write buffers, while the one other process advances 333 * to the call to compute the l2 norm. At this time, you will get an error 334 * that some operation was attempted by only a subset of processes. This 335 * behavior may seem surprising, unless you know that write/addition 336 * operations on single elements may trigger this behavior. 337 * 338 * The problem described here may be avoided by placing additional calls 339 * to compress(), or making sure that all processes do the same type of 340 * operations at the same time, for example by placing zero additions if 341 * necessary. 342 * 343 * 344 * <h3>Ghost elements of vectors</h3> 345 * 346 * Parallel vectors come in two kinds: without and with ghost elements. 347 * Vectors without ghost elements uniquely partition the vector elements 348 * between processors: each vector entry has exactly one processor that 349 * owns it. For such vectors, you can read those elements that the 350 * processor you are currently on owns, and you can write into any element 351 * whether you own it or not: if you don't own it, the value written or 352 * added to a vector element will be shipped to the processor that owns 353 * this vector element the next time you call compress(), as described 354 * above. 355 * 356 * What we call a 'ghosted' vector (see 357 * @ref GlossGhostedVector "vectors with ghost elements" 358 * ) is simply a view of the parallel vector where the element 359 * distributions overlap. The 'ghosted' Trilinos vector in itself has no 360 * idea of which entries are ghosted and which are locally owned. In fact, 361 * a ghosted vector may not even store all of the elements a non- ghosted 362 * vector would store on the current processor. Consequently, for 363 * Trilinos vectors, there is no notion of an 'owner' of vector elements 364 * in the way we have it in the non-ghost case view. 365 * 366 * This explains why we do not allow writing into ghosted vectors on the 367 * Trilinos side: Who would be responsible for taking care of the 368 * duplicated entries, given that there is not such information as locally 369 * owned indices? In other words, since a processor doesn't know which 370 * other processors own an element, who would it send a value to if one 371 * were to write to it? The only possibility would be to send this 372 * information to <i>all</i> other processors, but that is clearly not 373 * practical. Thus, we only allow reading from ghosted vectors, which 374 * however we do very often. 375 * 376 * So how do you fill a ghosted vector if you can't write to it? This only 377 * happens through the assignment with a non-ghosted vector. It can go 378 * both ways (non-ghosted is assigned to a ghosted vector, and a ghosted 379 * vector is assigned to a non-ghosted one; the latter one typically only 380 * requires taking out the locally owned part as most often ghosted 381 * vectors store a superset of elements of non-ghosted ones). In general, 382 * you send data around with that operation and it all depends on the 383 * different views of the two vectors. Trilinos also allows you to get 384 * subvectors out of a big vector that way. 385 * 386 * 387 * <h3>Thread safety of Trilinos vectors</h3> 388 * 389 * When writing into Trilinos vectors from several threads in shared 390 * memory, several things must be kept in mind as there is no built-in 391 * locks in this class to prevent data races. Simultaneous access to the 392 * same vector entry at the same time results in data races and must be 393 * explicitly avoided by the user. However, it is possible to access 394 * <b>different</b> entries of the vector from several threads 395 * simultaneously when only one MPI process is present or the vector has 396 * been constructed with an additional index set for ghost entries in 397 * write mode. 398 * 399 * @ingroup TrilinosWrappers 400 * @ingroup Vectors 401 * 2008, 2009, 2017 402 */ 403 class Vector : public Subscriptor 404 { 405 public: 406 /** 407 * Declare some of the standard types used in all containers. These types 408 * parallel those in the <tt>C</tt> standard libraries 409 * <tt>vector<...></tt> class. 410 */ 411 using value_type = TrilinosScalar; 412 using real_type = TrilinosScalar; 413 using size_type = dealii::types::global_dof_index; 414 using iterator = value_type *; 415 using const_iterator = const value_type *; 416 using reference = internal::VectorReference; 417 using const_reference = const internal::VectorReference; 418 419 /** 420 * @name 1: Basic Object-handling 421 */ 422 //@{ 423 /** 424 * Default constructor that generates an empty (zero size) vector. The 425 * function <tt>reinit()</tt> will have to give the vector the correct 426 * size and distribution among processes in case of an MPI run. 427 */ 428 Vector(); 429 430 /** 431 * Copy constructor using the given vector. 432 */ 433 Vector(const Vector &v); 434 435 /** 436 * This constructor takes an IndexSet that defines how to distribute the 437 * individual components among the MPI processors. Since it also 438 * includes information about the size of the vector, this is all we 439 * need to generate a %parallel vector. 440 * 441 * Depending on whether the @p parallel_partitioning argument uniquely 442 * subdivides elements among processors or not, the resulting vector may 443 * or may not have ghost elements. See the general documentation of this 444 * class for more information. 445 * 446 * In case the provided IndexSet forms an overlapping partitioning, 447 * it is not clear which elements are owned by which process and 448 * locally_owned_elements() will return an IndexSet of size zero. 449 * 450 * @see 451 * @ref GlossGhostedVector "vectors with ghost elements" 452 */ 453 explicit Vector(const IndexSet ¶llel_partitioning, 454 const MPI_Comm &communicator = MPI_COMM_WORLD); 455 456 /** 457 * Creates a ghosted parallel vector. 458 * 459 * Depending on whether the @p ghost argument uniquely subdivides 460 * elements among processors or not, the resulting vector may or may not 461 * have ghost elements. See the general documentation of this class for 462 * more information. 463 * 464 * @see 465 * @ref GlossGhostedVector "vectors with ghost elements" 466 */ 467 Vector(const IndexSet &local, 468 const IndexSet &ghost, 469 const MPI_Comm &communicator = MPI_COMM_WORLD); 470 471 /** 472 * Copy constructor from the TrilinosWrappers vector class. Since a 473 * vector of this class does not necessarily need to be distributed 474 * among processes, the user needs to supply us with an IndexSet and an 475 * MPI communicator that set the partitioning details. 476 * 477 * Depending on whether the @p parallel_partitioning argument uniquely 478 * subdivides elements among processors or not, the resulting vector may 479 * or may not have ghost elements. See the general documentation of this 480 * class for more information. 481 * 482 * @see 483 * @ref GlossGhostedVector "vectors with ghost elements" 484 */ 485 Vector(const IndexSet ¶llel_partitioning, 486 const Vector & v, 487 const MPI_Comm &communicator = MPI_COMM_WORLD); 488 489 /** 490 * Copy-constructor from deal.II vectors. Sets the dimension to that of 491 * the given vector, and copies all the elements. 492 * 493 * Depending on whether the @p parallel_partitioning argument uniquely 494 * subdivides elements among processors or not, the resulting vector may 495 * or may not have ghost elements. See the general documentation of this 496 * class for more information. 497 * 498 * @see 499 * @ref GlossGhostedVector "vectors with ghost elements" 500 */ 501 template <typename Number> 502 Vector(const IndexSet & parallel_partitioning, 503 const dealii::Vector<Number> &v, 504 const MPI_Comm & communicator = MPI_COMM_WORLD); 505 506 /** 507 * Move constructor. Creates a new vector by stealing the internal data 508 * of the vector @p v. 509 */ 510 Vector(Vector &&v) noexcept; 511 512 /** 513 * Destructor. 514 */ 515 ~Vector() override = default; 516 517 /** 518 * Release all memory and return to a state just like after having called 519 * the default constructor. 520 */ 521 void 522 clear(); 523 524 /** 525 * Reinit functionality. This function sets the calling vector to the 526 * dimension and the parallel distribution of the input vector, but does 527 * not copy the elements in <tt>v</tt>. If <tt>omit_zeroing_entries</tt> 528 * is not <tt>true</tt>, the elements in the vector are initialized with 529 * zero. If it is set to <tt>true</tt>, the vector entries are in an 530 * unspecified state and the user has to set all elements. In the 531 * current implementation, this method does not touch the vector entries 532 * in case the vector layout is unchanged from before, otherwise entries 533 * are set to zero. Note that this behavior might change between 534 * releases without notification. 535 * 536 * This function has a third argument, <tt>allow_different_maps</tt>, 537 * that allows for an exchange of data between two equal-sized vectors 538 * (but being distributed differently among the processors). A trivial 539 * application of this function is to generate a replication of a whole 540 * vector on each machine, when the calling vector is built with a map 541 * consisting of all indices on each process, and <tt>v</tt> 542 * is a distributed vector. In this case, the variable 543 * <tt>omit_zeroing_entries</tt> needs to be set to <tt>false</tt>, 544 * since it does not make sense to exchange data between differently 545 * parallelized vectors without touching the elements. 546 */ 547 void 548 reinit(const Vector &v, 549 const bool omit_zeroing_entries = false, 550 const bool allow_different_maps = false); 551 552 /** 553 * Reinit functionality. This function destroys the old vector content 554 * and generates a new one based on the input partitioning. The flag 555 * <tt>omit_zeroing_entries</tt> determines whether the vector should be 556 * filled with zero (false). If the flag is set to <tt>true</tt>, the 557 * vector entries are in an unspecified state and the user has to set 558 * all elements. In the current implementation, this method still sets 559 * the entries to zero, but this might change between releases without 560 * notification. 561 * 562 * Depending on whether the @p parallel_partitioning argument uniquely 563 * subdivides elements among processors or not, the resulting vector may 564 * or may not have ghost elements. See the general documentation of this 565 * class for more information. 566 * 567 * In case @p parallel_partitioning is overlapping, it is not clear which 568 * process should own which elements. Hence, locally_owned_elements() 569 * returns an empty IndexSet in this case. 570 * 571 * @see 572 * @ref GlossGhostedVector "vectors with ghost elements" 573 */ 574 void 575 reinit(const IndexSet ¶llel_partitioning, 576 const MPI_Comm &communicator = MPI_COMM_WORLD, 577 const bool omit_zeroing_entries = false); 578 579 /** 580 * Reinit functionality. This function destroys the old vector content 581 * and generates a new one based on the input partitioning. In addition 582 * to just specifying one index set as in all the other methods above, 583 * this method allows to supply an additional set of ghost entries. 584 * There are two different versions of a vector that can be created. If 585 * the flag @p vector_writable is set to @p false, the vector only 586 * allows read access to the joint set of @p parallel_partitioning and 587 * @p ghost_entries. The effect of the reinit method is then equivalent 588 * to calling the other reinit method with an index set containing both 589 * the locally owned entries and the ghost entries. 590 * 591 * If the flag @p vector_writable is set to true, this creates an 592 * alternative storage scheme for ghost elements that allows multiple 593 * threads to write into the vector (for the other reinit methods, only 594 * one thread is allowed to write into the ghost entries at a time). 595 * 596 * Depending on whether the @p ghost_entries argument uniquely 597 * subdivides elements among processors or not, the resulting vector may 598 * or may not have ghost elements. See the general documentation of this 599 * class for more information. 600 * 601 * @see 602 * @ref GlossGhostedVector "vectors with ghost elements" 603 */ 604 void 605 reinit(const IndexSet &locally_owned_entries, 606 const IndexSet &ghost_entries, 607 const MPI_Comm &communicator = MPI_COMM_WORLD, 608 const bool vector_writable = false); 609 610 /** 611 * Create vector by merging components from a block vector. 612 */ 613 void 614 reinit(const BlockVector &v, const bool import_data = false); 615 616 /** 617 * Compress the underlying representation of the Trilinos object, i.e. 618 * flush the buffers of the vector object if it has any. This function is 619 * necessary after writing into a vector element-by-element and before 620 * anything else can be done on it. 621 * 622 * The (defaulted) argument can be used to specify the compress mode 623 * (<code>Add</code> or <code>Insert</code>) in case the vector has not 624 * been written to since the last time this function was called. The 625 * argument is ignored if the vector has been added or written to since 626 * the last time compress() was called. 627 * 628 * See 629 * @ref GlossCompress "Compressing distributed objects" 630 * for more information. 631 */ 632 void 633 compress(::dealii::VectorOperation::values operation); 634 635 /** 636 * Set all components of the vector to the given number @p s. Simply 637 * pass this down to the base class, but we still need to declare this 638 * function to make the example given in the discussion about making the 639 * constructor explicit work. 640 * the constructor explicit work. 641 * 642 * Since the semantics of assigning a scalar to a vector are not 643 * immediately clear, this operator can only be used if you want 644 * to set the entire vector to zero. This allows the intuitive notation 645 * <tt>v=0</tt>. 646 */ 647 Vector & 648 operator=(const TrilinosScalar s); 649 650 /** 651 * Copy the given vector. Resize the present vector if necessary. In 652 * this case, also the Epetra_Map that designs the parallel partitioning 653 * is taken from the input vector. 654 */ 655 Vector & 656 operator=(const Vector &v); 657 658 /** 659 * Move the given vector. This operator replaces the present vector with 660 * @p v by efficiently swapping the internal data structures. 661 */ 662 Vector & 663 operator=(Vector &&v) noexcept; 664 665 /** 666 * Another copy function. This one takes a deal.II vector and copies it 667 * into a TrilinosWrapper vector. Note that since we do not provide any 668 * Epetra_map that tells about the partitioning of the vector among the 669 * MPI processes, the size of the TrilinosWrapper vector has to be the 670 * same as the size of the input vector. 671 */ 672 template <typename Number> 673 Vector & 674 operator=(const ::dealii::Vector<Number> &v); 675 676 /** 677 * This reinit function is meant to be used for parallel calculations 678 * where some non-local data has to be used. The typical situation where 679 * one needs this function is the call of the 680 * FEValues<dim>::get_function_values function (or of some derivatives) 681 * in parallel. Since it is usually faster to retrieve the data in 682 * advance, this function can be called before the assembly forks out to 683 * the different processors. What this function does is the following: 684 * It takes the information in the columns of the given matrix and looks 685 * which data couples between the different processors. That data is 686 * then queried from the input vector. Note that you should not write to 687 * the resulting vector any more, since the some data can be stored 688 * several times on different processors, leading to unpredictable 689 * results. In particular, such a vector cannot be used for matrix- 690 * vector products as for example done during the solution of linear 691 * systems. 692 */ 693 void 694 import_nonlocal_data_for_fe( 695 const dealii::TrilinosWrappers::SparseMatrix &matrix, 696 const Vector & vector); 697 698 /** 699 * Imports all the elements present in the vector's IndexSet from the 700 * input vector @p rwv. VectorOperation::values @p operation is used to decide if 701 * the elements in @p rwv should be added to the current vector or replace the 702 * current elements. 703 */ 704 void 705 import(const LinearAlgebra::ReadWriteVector<double> &rwv, 706 const VectorOperation::values operation); 707 708 709 /** 710 * Test for equality. This function assumes that the present vector and 711 * the one to compare with have the same size already, since comparing 712 * vectors of different sizes makes not much sense anyway. 713 */ 714 bool 715 operator==(const Vector &v) const; 716 717 /** 718 * Test for inequality. This function assumes that the present vector and 719 * the one to compare with have the same size already, since comparing 720 * vectors of different sizes makes not much sense anyway. 721 */ 722 bool 723 operator!=(const Vector &v) const; 724 725 /** 726 * Return the global dimension of the vector. 727 */ 728 size_type 729 size() const; 730 731 /** 732 * Return the local dimension of the vector, i.e. the number of elements 733 * stored on the present MPI process. For sequential vectors, this number 734 * is the same as size(), but for parallel vectors it may be smaller. 735 * 736 * To figure out which elements exactly are stored locally, use 737 * local_range(). 738 * 739 * If the vector contains ghost elements, they are included in this 740 * number. 741 */ 742 size_type 743 local_size() const; 744 745 /** 746 * Return a pair of indices indicating which elements of this vector are 747 * stored locally. The first number is the index of the first element 748 * stored, the second the index of the one past the last one that is 749 * stored locally. If this is a sequential vector, then the result will be 750 * the pair <code>(0,N)</code>, otherwise it will be a pair 751 * <code>(i,i+n)</code>, where <code>n=local_size()</code> and 752 * <code>i</code> is the first element of the vector stored on this 753 * processor, corresponding to the half open interval $[i,i+n)$ 754 * 755 * @note The description above is true most of the time, but not always. 756 * In particular, Trilinos vectors need not store contiguous ranges of 757 * elements such as $[i,i+n)$. Rather, it can store vectors where the 758 * elements are distributed in an arbitrary way across all processors and 759 * each processor simply stores a particular subset, not necessarily 760 * contiguous. In this case, this function clearly makes no sense since it 761 * could, at best, return a range that includes all elements that are 762 * stored locally. Thus, the function only succeeds if the locally stored 763 * range is indeed contiguous. It will trigger an assertion if the local 764 * portion of the vector is not contiguous. 765 */ 766 std::pair<size_type, size_type> 767 local_range() const; 768 769 /** 770 * Return whether @p index is in the local range or not, see also 771 * local_range(). 772 * 773 * @note The same limitation for the applicability of this function 774 * applies as listed in the documentation of local_range(). 775 */ 776 bool 777 in_local_range(const size_type index) const; 778 779 /** 780 * Return an index set that describes which elements of this vector are 781 * owned by the current processor. Note that this index set does not 782 * include elements this vector may store locally as ghost elements but 783 * that are in fact owned by another processor. As a consequence, the 784 * index sets returned on different processors if this is a distributed 785 * vector will form disjoint sets that add up to the complete index set. 786 * Obviously, if a vector is created on only one processor, then the 787 * result would satisfy 788 * @code 789 * vec.locally_owned_elements() == complete_index_set (vec.size()) 790 * @endcode 791 */ 792 IndexSet 793 locally_owned_elements() const; 794 795 /** 796 * Return if the vector contains ghost elements. This answer is true if 797 * there are ghost elements on at least one process. 798 * 799 * @see 800 * @ref GlossGhostedVector "vectors with ghost elements" 801 */ 802 bool 803 has_ghost_elements() const; 804 805 /** 806 * This function only exists for compatibility with the @p 807 * LinearAlgebra::distributed::Vector class and does nothing: this class 808 * implements ghost value updates in a different way that is a better fit 809 * with the underlying Trilinos vector object. 810 */ 811 void 812 update_ghost_values() const; 813 814 /** 815 * Return the scalar (inner) product of two vectors. The vectors must have 816 * the same size. 817 */ 818 TrilinosScalar operator*(const Vector &vec) const; 819 820 /** 821 * Return the square of the $l_2$-norm. 822 */ 823 real_type 824 norm_sqr() const; 825 826 /** 827 * Mean value of the elements of this vector. 828 */ 829 TrilinosScalar 830 mean_value() const; 831 832 /** 833 * Compute the minimal value of the elements of this vector. 834 */ 835 TrilinosScalar 836 min() const; 837 838 /** 839 * Compute the maximal value of the elements of this vector. 840 */ 841 TrilinosScalar 842 max() const; 843 844 /** 845 * $l_1$-norm of the vector. The sum of the absolute values. 846 */ 847 real_type 848 l1_norm() const; 849 850 /** 851 * $l_2$-norm of the vector. The square root of the sum of the squares of 852 * the elements. 853 */ 854 real_type 855 l2_norm() const; 856 857 /** 858 * $l_p$-norm of the vector. The <i>p</i>th root of the sum of the 859 * <i>p</i>th powers of the absolute values of the elements. 860 */ 861 real_type 862 lp_norm(const TrilinosScalar p) const; 863 864 /** 865 * Maximum absolute value of the elements. 866 */ 867 real_type 868 linfty_norm() const; 869 870 /** 871 * Performs a combined operation of a vector addition and a subsequent 872 * inner product, returning the value of the inner product. In other 873 * words, the result of this function is the same as if the user called 874 * @code 875 * this->add(a, V); 876 * return_value = *this * W; 877 * @endcode 878 * 879 * The reason this function exists is for compatibility with deal.II's own 880 * vector classes which can implement this functionality with less memory 881 * transfer. However, for Trilinos vectors such a combined operation is 882 * not natively supported and thus the cost is completely equivalent as 883 * calling the two methods separately. 884 * 885 * For complex-valued vectors, the scalar product in the second step is 886 * implemented as 887 * $\left<v,w\right>=\sum_i v_i \bar{w_i}$. 888 */ 889 TrilinosScalar 890 add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W); 891 892 /** 893 * Return whether the vector contains only elements with value zero. This 894 * is a collective operation. This function is expensive, because 895 * potentially all elements have to be checked. 896 */ 897 bool 898 all_zero() const; 899 900 /** 901 * Return @p true if the vector has no negative entries, i.e. all entries 902 * are zero or positive. This function is used, for example, to check 903 * whether refinement indicators are really all positive (or zero). 904 */ 905 bool 906 is_non_negative() const; 907 //@} 908 909 910 /** 911 * @name 2: Data-Access 912 */ 913 //@{ 914 915 /** 916 * Provide access to a given element, both read and write. 917 * 918 * When using a vector distributed with MPI, this operation only makes 919 * sense for elements that are actually present on the calling processor. 920 * Otherwise, an exception is thrown. 921 */ 922 reference 923 operator()(const size_type index); 924 925 /** 926 * Provide read-only access to an element. 927 * 928 * When using a vector distributed with MPI, this operation only makes 929 * sense for elements that are actually present on the calling processor. 930 * Otherwise, an exception is thrown. 931 */ 932 TrilinosScalar 933 operator()(const size_type index) const; 934 935 /** 936 * Provide access to a given element, both read and write. 937 * 938 * Exactly the same as operator(). 939 */ 940 reference operator[](const size_type index); 941 942 /** 943 * Provide read-only access to an element. 944 * 945 * Exactly the same as operator(). 946 */ 947 TrilinosScalar operator[](const size_type index) const; 948 949 /** 950 * Instead of getting individual elements of a vector via operator(), 951 * this function allows getting a whole set of elements at once. The 952 * indices of the elements to be read are stated in the first argument, 953 * the corresponding values are returned in the second. 954 * 955 * If the current vector is called @p v, then this function is the equivalent 956 * to the code 957 * @code 958 * for (unsigned int i=0; i<indices.size(); ++i) 959 * values[i] = v[indices[i]]; 960 * @endcode 961 * 962 * @pre The sizes of the @p indices and @p values arrays must be identical. 963 */ 964 void 965 extract_subvector_to(const std::vector<size_type> &indices, 966 std::vector<TrilinosScalar> & values) const; 967 968 /** 969 * Instead of getting individual elements of a vector via operator(), 970 * this function allows getting a whole set of elements at once. In 971 * contrast to the previous function, this function obtains the 972 * indices of the elements by dereferencing all elements of the iterator 973 * range provided by the first two arguments, and puts the vector 974 * values into memory locations obtained by dereferencing a range 975 * of iterators starting at the location pointed to by the third 976 * argument. 977 * 978 * If the current vector is called @p v, then this function is the equivalent 979 * to the code 980 * @code 981 * ForwardIterator indices_p = indices_begin; 982 * OutputIterator values_p = values_begin; 983 * while (indices_p != indices_end) 984 * { 985 * *values_p = v[*indices_p]; 986 * ++indices_p; 987 * ++values_p; 988 * } 989 * @endcode 990 * 991 * @pre It must be possible to write into as many memory locations 992 * starting at @p values_begin as there are iterators between 993 * @p indices_begin and @p indices_end. 994 */ 995 template <typename ForwardIterator, typename OutputIterator> 996 void 997 extract_subvector_to(ForwardIterator indices_begin, 998 const ForwardIterator indices_end, 999 OutputIterator values_begin) const; 1000 1001 /** 1002 * Make the Vector class a bit like the <tt>vector<></tt> class of the C++ 1003 * standard library by returning iterators to the start and end of the 1004 * locally owned elements of this vector. The ordering of local elements 1005 * corresponds to the one given by the global indices in case the vector 1006 * is constructed from an IndexSet or other methods in deal.II (note that 1007 * an Epetra_Map can contain elements in arbitrary orders, though). 1008 * 1009 * It holds that end() - begin() == local_size(). 1010 */ 1011 iterator 1012 begin(); 1013 1014 /** 1015 * Return constant iterator to the start of the locally owned elements of 1016 * the vector. 1017 */ 1018 const_iterator 1019 begin() const; 1020 1021 /** 1022 * Return an iterator pointing to the element past the end of the array of 1023 * locally owned entries. 1024 */ 1025 iterator 1026 end(); 1027 1028 /** 1029 * Return a constant iterator pointing to the element past the end of the 1030 * array of the locally owned entries. 1031 */ 1032 const_iterator 1033 end() const; 1034 1035 //@} 1036 1037 1038 /** 1039 * @name 3: Modification of vectors 1040 */ 1041 //@{ 1042 1043 /** 1044 * A collective set operation: instead of setting individual elements of a 1045 * vector, this function allows to set a whole set of elements at once. 1046 * The indices of the elements to be set are stated in the first argument, 1047 * the corresponding values in the second. 1048 */ 1049 void 1050 set(const std::vector<size_type> & indices, 1051 const std::vector<TrilinosScalar> &values); 1052 1053 /** 1054 * This is a second collective set operation. As a difference, this 1055 * function takes a deal.II vector of values. 1056 */ 1057 void 1058 set(const std::vector<size_type> & indices, 1059 const ::dealii::Vector<TrilinosScalar> &values); 1060 1061 /** 1062 * This collective set operation is of lower level and can handle anything 1063 * else — the only thing you have to provide is an address where all 1064 * the indices are stored and the number of elements to be set. 1065 */ 1066 void 1067 set(const size_type n_elements, 1068 const size_type * indices, 1069 const TrilinosScalar *values); 1070 1071 /** 1072 * A collective add operation: This function adds a whole set of values 1073 * stored in @p values to the vector components specified by @p indices. 1074 */ 1075 void 1076 add(const std::vector<size_type> & indices, 1077 const std::vector<TrilinosScalar> &values); 1078 1079 /** 1080 * This is a second collective add operation. As a difference, this 1081 * function takes a deal.II vector of values. 1082 */ 1083 void 1084 add(const std::vector<size_type> & indices, 1085 const ::dealii::Vector<TrilinosScalar> &values); 1086 1087 /** 1088 * Take an address where <tt>n_elements</tt> are stored contiguously and 1089 * add them into the vector. Handles all cases which are not covered by 1090 * the other two <tt>add()</tt> functions above. 1091 */ 1092 void 1093 add(const size_type n_elements, 1094 const size_type * indices, 1095 const TrilinosScalar *values); 1096 1097 /** 1098 * Multiply the entire vector by a fixed factor. 1099 */ 1100 Vector & 1101 operator*=(const TrilinosScalar factor); 1102 1103 /** 1104 * Divide the entire vector by a fixed factor. 1105 */ 1106 Vector & 1107 operator/=(const TrilinosScalar factor); 1108 1109 /** 1110 * Add the given vector to the present one. 1111 */ 1112 Vector & 1113 operator+=(const Vector &V); 1114 1115 /** 1116 * Subtract the given vector from the present one. 1117 */ 1118 Vector & 1119 operator-=(const Vector &V); 1120 1121 /** 1122 * Addition of @p s to all components. Note that @p s is a scalar and not 1123 * a vector. 1124 */ 1125 void 1126 add(const TrilinosScalar s); 1127 1128 /** 1129 * Simple vector addition, equal to the <tt>operator+=</tt>. 1130 * 1131 * Though, if the second argument <tt>allow_different_maps</tt> is set, 1132 * then it is possible to add data from a vector that uses a different 1133 * map, i.e., a vector whose elements are split across processors 1134 * differently. This may include vectors with ghost elements, for example. 1135 * In general, however, adding vectors with a different element-to- 1136 * processor map requires communicating data among processors and, 1137 * consequently, is a slower operation than when using vectors using the 1138 * same map. 1139 */ 1140 void 1141 add(const Vector &V, const bool allow_different_maps = false); 1142 1143 /** 1144 * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>. 1145 */ 1146 void 1147 add(const TrilinosScalar a, const Vector &V); 1148 1149 /** 1150 * Multiple addition of scaled vectors, i.e. <tt>*this += a*V + b*W</tt>. 1151 */ 1152 void 1153 add(const TrilinosScalar a, 1154 const Vector & V, 1155 const TrilinosScalar b, 1156 const Vector & W); 1157 1158 /** 1159 * Scaling and simple vector addition, i.e. <tt>*this = s*(*this) + 1160 * V</tt>. 1161 */ 1162 void 1163 sadd(const TrilinosScalar s, const Vector &V); 1164 1165 /** 1166 * Scaling and simple addition, i.e. <tt>*this = s*(*this) + a*V</tt>. 1167 */ 1168 void 1169 sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V); 1170 1171 /** 1172 * Scale each element of this vector by the corresponding element in the 1173 * argument. This function is mostly meant to simulate multiplication (and 1174 * immediate re-assignment) by a diagonal scaling matrix. 1175 */ 1176 void 1177 scale(const Vector &scaling_factors); 1178 1179 /** 1180 * Assignment <tt>*this = a*V</tt>. 1181 */ 1182 void 1183 equ(const TrilinosScalar a, const Vector &V); 1184 //@} 1185 1186 /** 1187 * @name 4: Mixed stuff 1188 */ 1189 //@{ 1190 1191 /** 1192 * Return a const reference to the underlying Trilinos Epetra_MultiVector 1193 * class. 1194 */ 1195 const Epetra_MultiVector & 1196 trilinos_vector() const; 1197 1198 /** 1199 * Return a (modifiable) reference to the underlying Trilinos 1200 * Epetra_FEVector class. 1201 */ 1202 Epetra_FEVector & 1203 trilinos_vector(); 1204 1205 /** 1206 * Return a const reference to the underlying Trilinos Epetra_BlockMap 1207 * that sets the parallel partitioning of the vector. 1208 */ 1209 const Epetra_BlockMap & 1210 trilinos_partitioner() const; 1211 1212 /** 1213 * Print to a stream. @p precision denotes the desired precision with 1214 * which values shall be printed, @p scientific whether scientific 1215 * notation shall be used. If @p across is @p true then the vector is 1216 * printed in a line, while if @p false then the elements are printed on a 1217 * separate line each. 1218 */ 1219 void 1220 print(std::ostream & out, 1221 const unsigned int precision = 3, 1222 const bool scientific = true, 1223 const bool across = true) const; 1224 1225 /** 1226 * Swap the contents of this vector and the other vector @p v. One could 1227 * do this operation with a temporary variable and copying over the data 1228 * elements, but this function is significantly more efficient since it 1229 * only swaps the pointers to the data of the two vectors and therefore 1230 * does not need to allocate temporary storage and move data around. Note 1231 * that the vectors need to be of the same size and base on the same map. 1232 * 1233 * This function is analogous to the @p swap function of all C++ 1234 * standard containers. Also, there is a global function 1235 * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in 1236 * analogy to standard functions. 1237 */ 1238 void 1239 swap(Vector &v); 1240 1241 /** 1242 * Estimate for the memory consumption in bytes. 1243 */ 1244 std::size_t 1245 memory_consumption() const; 1246 1247 /** 1248 * Return a reference to the MPI communicator object in use with this 1249 * object. 1250 */ 1251 const MPI_Comm & 1252 get_mpi_communicator() const; 1253 //@} 1254 1255 /** 1256 * Exception 1257 */ 1258 DeclException0(ExcDifferentParallelPartitioning); 1259 1260 /** 1261 * Exception 1262 */ 1263 DeclException1(ExcTrilinosError, 1264 int, 1265 << "An error with error number " << arg1 1266 << " occurred while calling a Trilinos function"); 1267 1268 /** 1269 * Exception 1270 */ 1271 DeclException4( 1272 ExcAccessToNonLocalElement, 1273 size_type, 1274 size_type, 1275 size_type, 1276 size_type, 1277 << "You tried to access element " << arg1 1278 << " of a distributed vector, but this element is not stored " 1279 << "on the current processor. Note: There are " << arg2 1280 << " elements stored " 1281 << "on the current processor from within the range " << arg3 1282 << " through " << arg4 1283 << " but Trilinos vectors need not store contiguous " 1284 << "ranges on each processor, and not every element in " 1285 << "this range may in fact be stored locally."); 1286 1287 private: 1288 /** 1289 * Trilinos doesn't allow to mix additions to matrix entries and 1290 * overwriting them (to make synchronization of parallel computations 1291 * simpler). The way we do it is to, for each access operation, store 1292 * whether it is an insertion or an addition. If the previous one was of 1293 * different type, then we first have to flush the Trilinos buffers; 1294 * otherwise, we can simply go on. Luckily, Trilinos has an object for 1295 * this which does already all the parallel communications in such a case, 1296 * so we simply use their model, which stores whether the last operation 1297 * was an addition or an insertion. 1298 */ 1299 Epetra_CombineMode last_action; 1300 1301 /** 1302 * A boolean variable to hold information on whether the vector is 1303 * compressed or not. 1304 */ 1305 bool compressed; 1306 1307 /** 1308 * Whether this vector has ghost elements. This is true on all processors 1309 * even if only one of them has any ghost elements. 1310 */ 1311 bool has_ghosts; 1312 1313 /** 1314 * Pointer to the actual Epetra vector object. This may represent a vector 1315 * that is in fact distributed among multiple processors. The object 1316 * requires an existing Epetra_Map for storing data when setting it up. 1317 */ 1318 std::unique_ptr<Epetra_FEVector> vector; 1319 1320 /** 1321 * A vector object in Trilinos to be used for collecting the non-local 1322 * elements if the vector was constructed with an additional IndexSet 1323 * describing ghost elements. 1324 */ 1325 std::unique_ptr<Epetra_MultiVector> nonlocal_vector; 1326 1327 /** 1328 * An IndexSet storing the indices this vector owns exclusively. 1329 */ 1330 IndexSet owned_elements; 1331 1332 // Make the reference class a friend. 1333 friend class internal::VectorReference; 1334 }; 1335 1336 1337 1338 // ------------------- inline and template functions -------------- 1339 1340 1341 /** 1342 * Global function @p swap which overloads the default implementation of 1343 * the C++ standard library which uses a temporary object. The function 1344 * simply exchanges the data of the two vectors. 1345 * 1346 * @relatesalso TrilinosWrappers::MPI::Vector 1347 */ 1348 inline void swap(Vector & u,Vector & v)1349 swap(Vector &u, Vector &v) 1350 { 1351 u.swap(v); 1352 } 1353 } // namespace MPI 1354 1355 # ifndef DOXYGEN 1356 1357 namespace internal 1358 { VectorReference(MPI::Vector & vector,const size_type index)1359 inline VectorReference::VectorReference(MPI::Vector & vector, 1360 const size_type index) 1361 : vector(vector) 1362 , index(index) 1363 {} 1364 1365 1366 inline const VectorReference & 1367 VectorReference::operator=(const VectorReference &r) const 1368 { 1369 // as explained in the class 1370 // documentation, this is not the copy 1371 // operator. so simply pass on to the 1372 // "correct" assignment operator 1373 *this = static_cast<TrilinosScalar>(r); 1374 1375 return *this; 1376 } 1377 1378 1379 1380 inline VectorReference & 1381 VectorReference::operator=(const VectorReference &r) 1382 { 1383 // as above 1384 *this = static_cast<TrilinosScalar>(r); 1385 1386 return *this; 1387 } 1388 1389 1390 inline const VectorReference & 1391 VectorReference::operator=(const TrilinosScalar &value) const 1392 { 1393 vector.set(1, &index, &value); 1394 return *this; 1395 } 1396 1397 1398 1399 inline const VectorReference & 1400 VectorReference::operator+=(const TrilinosScalar &value) const 1401 { 1402 vector.add(1, &index, &value); 1403 return *this; 1404 } 1405 1406 1407 1408 inline const VectorReference & 1409 VectorReference::operator-=(const TrilinosScalar &value) const 1410 { 1411 TrilinosScalar new_value = -value; 1412 vector.add(1, &index, &new_value); 1413 return *this; 1414 } 1415 1416 1417 1418 inline const VectorReference & 1419 VectorReference::operator*=(const TrilinosScalar &value) const 1420 { 1421 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value; 1422 vector.set(1, &index, &new_value); 1423 return *this; 1424 } 1425 1426 1427 1428 inline const VectorReference & 1429 VectorReference::operator/=(const TrilinosScalar &value) const 1430 { 1431 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value; 1432 vector.set(1, &index, &new_value); 1433 return *this; 1434 } 1435 } // namespace internal 1436 1437 namespace MPI 1438 { 1439 inline bool in_local_range(const size_type index)1440 Vector::in_local_range(const size_type index) const 1441 { 1442 std::pair<size_type, size_type> range = local_range(); 1443 1444 return ((index >= range.first) && (index < range.second)); 1445 } 1446 1447 1448 1449 inline IndexSet locally_owned_elements()1450 Vector::locally_owned_elements() const 1451 { 1452 Assert(owned_elements.size() == size(), 1453 ExcMessage( 1454 "The locally owned elements have not been properly initialized!" 1455 " This happens for example if this object has been initialized" 1456 " with exactly one overlapping IndexSet.")); 1457 return owned_elements; 1458 } 1459 1460 1461 1462 inline bool has_ghost_elements()1463 Vector::has_ghost_elements() const 1464 { 1465 return has_ghosts; 1466 } 1467 1468 1469 1470 inline void update_ghost_values()1471 Vector::update_ghost_values() const 1472 {} 1473 1474 1475 1476 inline internal::VectorReference operator()1477 Vector::operator()(const size_type index) 1478 { 1479 return internal::VectorReference(*this, index); 1480 } 1481 1482 1483 1484 inline internal::VectorReference Vector::operator[](const size_type index) 1485 { 1486 return operator()(index); 1487 } 1488 1489 1490 1491 inline TrilinosScalar Vector::operator[](const size_type index) const 1492 { 1493 return operator()(index); 1494 } 1495 1496 1497 1498 inline void extract_subvector_to(const std::vector<size_type> & indices,std::vector<TrilinosScalar> & values)1499 Vector::extract_subvector_to(const std::vector<size_type> &indices, 1500 std::vector<TrilinosScalar> & values) const 1501 { 1502 for (size_type i = 0; i < indices.size(); ++i) 1503 values[i] = operator()(indices[i]); 1504 } 1505 1506 1507 1508 template <typename ForwardIterator, typename OutputIterator> 1509 inline void extract_subvector_to(ForwardIterator indices_begin,const ForwardIterator indices_end,OutputIterator values_begin)1510 Vector::extract_subvector_to(ForwardIterator indices_begin, 1511 const ForwardIterator indices_end, 1512 OutputIterator values_begin) const 1513 { 1514 while (indices_begin != indices_end) 1515 { 1516 *values_begin = operator()(*indices_begin); 1517 indices_begin++; 1518 values_begin++; 1519 } 1520 } 1521 1522 1523 1524 inline Vector::iterator begin()1525 Vector::begin() 1526 { 1527 return (*vector)[0]; 1528 } 1529 1530 1531 1532 inline Vector::iterator end()1533 Vector::end() 1534 { 1535 return (*vector)[0] + local_size(); 1536 } 1537 1538 1539 1540 inline Vector::const_iterator begin()1541 Vector::begin() const 1542 { 1543 return (*vector)[0]; 1544 } 1545 1546 1547 1548 inline Vector::const_iterator end()1549 Vector::end() const 1550 { 1551 return (*vector)[0] + local_size(); 1552 } 1553 1554 1555 1556 inline void set(const std::vector<size_type> & indices,const std::vector<TrilinosScalar> & values)1557 Vector::set(const std::vector<size_type> & indices, 1558 const std::vector<TrilinosScalar> &values) 1559 { 1560 // if we have ghost values, do not allow 1561 // writing to this vector at all. 1562 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1563 1564 Assert(indices.size() == values.size(), 1565 ExcDimensionMismatch(indices.size(), values.size())); 1566 1567 set(indices.size(), indices.data(), values.data()); 1568 } 1569 1570 1571 1572 inline void set(const std::vector<size_type> & indices,const::dealii::Vector<TrilinosScalar> & values)1573 Vector::set(const std::vector<size_type> & indices, 1574 const ::dealii::Vector<TrilinosScalar> &values) 1575 { 1576 // if we have ghost values, do not allow 1577 // writing to this vector at all. 1578 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1579 1580 Assert(indices.size() == values.size(), 1581 ExcDimensionMismatch(indices.size(), values.size())); 1582 1583 set(indices.size(), indices.data(), values.begin()); 1584 } 1585 1586 1587 1588 inline void set(const size_type n_elements,const size_type * indices,const TrilinosScalar * values)1589 Vector::set(const size_type n_elements, 1590 const size_type * indices, 1591 const TrilinosScalar *values) 1592 { 1593 // if we have ghost values, do not allow 1594 // writing to this vector at all. 1595 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1596 1597 if (last_action == Add) 1598 { 1599 const int ierr = vector->GlobalAssemble(Add); 1600 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1601 } 1602 1603 if (last_action != Insert) 1604 last_action = Insert; 1605 1606 for (size_type i = 0; i < n_elements; ++i) 1607 { 1608 const TrilinosWrappers::types::int_type row = indices[i]; 1609 const TrilinosWrappers::types::int_type local_row = 1610 vector->Map().LID(row); 1611 if (local_row != -1) 1612 (*vector)[0][local_row] = values[i]; 1613 else 1614 { 1615 const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]); 1616 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1617 compressed = false; 1618 } 1619 // in set operation, do not use the pre-allocated vector for nonlocal 1620 // entries even if it exists. This is to ensure that we really only 1621 // set the elements touched by the set() method and not all contained 1622 // in the nonlocal entries vector (there is no way to distinguish them 1623 // on the receiving processor) 1624 } 1625 } 1626 1627 1628 1629 inline void add(const std::vector<size_type> & indices,const std::vector<TrilinosScalar> & values)1630 Vector::add(const std::vector<size_type> & indices, 1631 const std::vector<TrilinosScalar> &values) 1632 { 1633 // if we have ghost values, do not allow 1634 // writing to this vector at all. 1635 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1636 Assert(indices.size() == values.size(), 1637 ExcDimensionMismatch(indices.size(), values.size())); 1638 1639 add(indices.size(), indices.data(), values.data()); 1640 } 1641 1642 1643 1644 inline void add(const std::vector<size_type> & indices,const::dealii::Vector<TrilinosScalar> & values)1645 Vector::add(const std::vector<size_type> & indices, 1646 const ::dealii::Vector<TrilinosScalar> &values) 1647 { 1648 // if we have ghost values, do not allow 1649 // writing to this vector at all. 1650 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1651 Assert(indices.size() == values.size(), 1652 ExcDimensionMismatch(indices.size(), values.size())); 1653 1654 add(indices.size(), indices.data(), values.begin()); 1655 } 1656 1657 1658 1659 inline void add(const size_type n_elements,const size_type * indices,const TrilinosScalar * values)1660 Vector::add(const size_type n_elements, 1661 const size_type * indices, 1662 const TrilinosScalar *values) 1663 { 1664 // if we have ghost values, do not allow 1665 // writing to this vector at all. 1666 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1667 1668 if (last_action != Add) 1669 { 1670 if (last_action == Insert) 1671 { 1672 const int ierr = vector->GlobalAssemble(Insert); 1673 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1674 } 1675 last_action = Add; 1676 } 1677 1678 for (size_type i = 0; i < n_elements; ++i) 1679 { 1680 const size_type row = indices[i]; 1681 const TrilinosWrappers::types::int_type local_row = vector->Map().LID( 1682 static_cast<TrilinosWrappers::types::int_type>(row)); 1683 if (local_row != -1) 1684 (*vector)[0][local_row] += values[i]; 1685 else if (nonlocal_vector.get() == nullptr) 1686 { 1687 const int ierr = vector->SumIntoGlobalValues( 1688 1, 1689 reinterpret_cast<const TrilinosWrappers::types::int_type *>( 1690 &row), 1691 &values[i]); 1692 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1693 compressed = false; 1694 } 1695 else 1696 { 1697 // use pre-allocated vector for non-local entries if it exists for 1698 // addition operation 1699 const TrilinosWrappers::types::int_type my_row = 1700 nonlocal_vector->Map().LID( 1701 static_cast<TrilinosWrappers::types::int_type>(row)); 1702 Assert(my_row != -1, 1703 ExcMessage( 1704 "Attempted to write into off-processor vector entry " 1705 "that has not be specified as being writable upon " 1706 "initialization")); 1707 (*nonlocal_vector)[0][my_row] += values[i]; 1708 compressed = false; 1709 } 1710 } 1711 } 1712 1713 1714 1715 inline Vector::size_type size()1716 Vector::size() const 1717 { 1718 # ifndef DEAL_II_WITH_64BIT_INDICES 1719 return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID(); 1720 # else 1721 return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64(); 1722 # endif 1723 } 1724 1725 1726 1727 inline Vector::size_type local_size()1728 Vector::local_size() const 1729 { 1730 return vector->Map().NumMyElements(); 1731 } 1732 1733 1734 1735 inline std::pair<Vector::size_type, Vector::size_type> local_range()1736 Vector::local_range() const 1737 { 1738 # ifndef DEAL_II_WITH_64BIT_INDICES 1739 const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID(); 1740 const TrilinosWrappers::types::int_type end = 1741 vector->Map().MaxMyGID() + 1; 1742 # else 1743 const TrilinosWrappers::types::int_type begin = 1744 vector->Map().MinMyGID64(); 1745 const TrilinosWrappers::types::int_type end = 1746 vector->Map().MaxMyGID64() + 1; 1747 # endif 1748 1749 Assert( 1750 end - begin == vector->Map().NumMyElements(), 1751 ExcMessage( 1752 "This function only makes sense if the elements that this " 1753 "vector stores on the current processor form a contiguous range. " 1754 "This does not appear to be the case for the current vector.")); 1755 1756 return std::make_pair(begin, end); 1757 } 1758 1759 1760 1761 inline TrilinosScalar Vector::operator*(const Vector &vec) const 1762 { 1763 Assert(vector->Map().SameAs(vec.vector->Map()), 1764 ExcDifferentParallelPartitioning()); 1765 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1766 1767 TrilinosScalar result; 1768 1769 const int ierr = vector->Dot(*(vec.vector), &result); 1770 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1771 1772 return result; 1773 } 1774 1775 1776 1777 inline Vector::real_type norm_sqr()1778 Vector::norm_sqr() const 1779 { 1780 const TrilinosScalar d = l2_norm(); 1781 return d * d; 1782 } 1783 1784 1785 1786 inline TrilinosScalar mean_value()1787 Vector::mean_value() const 1788 { 1789 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1790 1791 TrilinosScalar mean; 1792 const int ierr = vector->MeanValue(&mean); 1793 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1794 1795 return mean; 1796 } 1797 1798 1799 1800 inline TrilinosScalar min()1801 Vector::min() const 1802 { 1803 TrilinosScalar min_value; 1804 const int ierr = vector->MinValue(&min_value); 1805 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1806 1807 return min_value; 1808 } 1809 1810 1811 1812 inline TrilinosScalar max()1813 Vector::max() const 1814 { 1815 TrilinosScalar max_value; 1816 const int ierr = vector->MaxValue(&max_value); 1817 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1818 1819 return max_value; 1820 } 1821 1822 1823 1824 inline Vector::real_type l1_norm()1825 Vector::l1_norm() const 1826 { 1827 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1828 1829 TrilinosScalar d; 1830 const int ierr = vector->Norm1(&d); 1831 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1832 1833 return d; 1834 } 1835 1836 1837 1838 inline Vector::real_type l2_norm()1839 Vector::l2_norm() const 1840 { 1841 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1842 1843 TrilinosScalar d; 1844 const int ierr = vector->Norm2(&d); 1845 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1846 1847 return d; 1848 } 1849 1850 1851 1852 inline Vector::real_type lp_norm(const TrilinosScalar p)1853 Vector::lp_norm(const TrilinosScalar p) const 1854 { 1855 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1856 1857 TrilinosScalar norm = 0; 1858 TrilinosScalar sum = 0; 1859 const size_type n_local = local_size(); 1860 1861 // loop over all the elements because 1862 // Trilinos does not support lp norms 1863 for (size_type i = 0; i < n_local; ++i) 1864 sum += std::pow(std::fabs((*vector)[0][i]), p); 1865 1866 norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p)); 1867 1868 return norm; 1869 } 1870 1871 1872 1873 inline Vector::real_type linfty_norm()1874 Vector::linfty_norm() const 1875 { 1876 // while we disallow the other 1877 // norm operations on ghosted 1878 // vectors, this particular norm 1879 // is safe to run even in the 1880 // presence of ghost elements 1881 TrilinosScalar d; 1882 const int ierr = vector->NormInf(&d); 1883 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1884 1885 return d; 1886 } 1887 1888 1889 1890 inline TrilinosScalar add_and_dot(const TrilinosScalar a,const Vector & V,const Vector & W)1891 Vector::add_and_dot(const TrilinosScalar a, 1892 const Vector & V, 1893 const Vector & W) 1894 { 1895 this->add(a, V); 1896 return *this * W; 1897 } 1898 1899 1900 1901 // inline also scalar products, vector 1902 // additions etc. since they are all 1903 // representable by a single Trilinos 1904 // call. This reduces the overhead of the 1905 // wrapper class. 1906 inline Vector & 1907 Vector::operator*=(const TrilinosScalar a) 1908 { 1909 AssertIsFinite(a); 1910 1911 const int ierr = vector->Scale(a); 1912 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1913 1914 return *this; 1915 } 1916 1917 1918 1919 inline Vector & 1920 Vector::operator/=(const TrilinosScalar a) 1921 { 1922 AssertIsFinite(a); 1923 1924 const TrilinosScalar factor = 1. / a; 1925 1926 AssertIsFinite(factor); 1927 1928 const int ierr = vector->Scale(factor); 1929 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1930 1931 return *this; 1932 } 1933 1934 1935 1936 inline Vector & 1937 Vector::operator+=(const Vector &v) 1938 { 1939 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size())); 1940 Assert(vector->Map().SameAs(v.vector->Map()), 1941 ExcDifferentParallelPartitioning()); 1942 1943 const int ierr = vector->Update(1.0, *(v.vector), 1.0); 1944 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1945 1946 return *this; 1947 } 1948 1949 1950 1951 inline Vector & 1952 Vector::operator-=(const Vector &v) 1953 { 1954 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size())); 1955 Assert(vector->Map().SameAs(v.vector->Map()), 1956 ExcDifferentParallelPartitioning()); 1957 1958 const int ierr = vector->Update(-1.0, *(v.vector), 1.0); 1959 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1960 1961 return *this; 1962 } 1963 1964 1965 1966 inline void add(const TrilinosScalar s)1967 Vector::add(const TrilinosScalar s) 1968 { 1969 // if we have ghost values, do not allow 1970 // writing to this vector at all. 1971 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1972 AssertIsFinite(s); 1973 1974 size_type n_local = local_size(); 1975 for (size_type i = 0; i < n_local; ++i) 1976 (*vector)[0][i] += s; 1977 } 1978 1979 1980 1981 inline void add(const TrilinosScalar a,const Vector & v)1982 Vector::add(const TrilinosScalar a, const Vector &v) 1983 { 1984 // if we have ghost values, do not allow 1985 // writing to this vector at all. 1986 Assert(!has_ghost_elements(), ExcGhostsPresent()); 1987 Assert(local_size() == v.local_size(), 1988 ExcDimensionMismatch(local_size(), v.local_size())); 1989 1990 AssertIsFinite(a); 1991 1992 const int ierr = vector->Update(a, *(v.vector), 1.); 1993 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 1994 } 1995 1996 1997 1998 inline void add(const TrilinosScalar a,const Vector & v,const TrilinosScalar b,const Vector & w)1999 Vector::add(const TrilinosScalar a, 2000 const Vector & v, 2001 const TrilinosScalar b, 2002 const Vector & w) 2003 { 2004 // if we have ghost values, do not allow 2005 // writing to this vector at all. 2006 Assert(!has_ghost_elements(), ExcGhostsPresent()); 2007 Assert(local_size() == v.local_size(), 2008 ExcDimensionMismatch(local_size(), v.local_size())); 2009 Assert(local_size() == w.local_size(), 2010 ExcDimensionMismatch(local_size(), w.local_size())); 2011 2012 AssertIsFinite(a); 2013 AssertIsFinite(b); 2014 2015 const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.); 2016 2017 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 2018 } 2019 2020 2021 2022 inline void sadd(const TrilinosScalar s,const Vector & v)2023 Vector::sadd(const TrilinosScalar s, const Vector &v) 2024 { 2025 // if we have ghost values, do not allow 2026 // writing to this vector at all. 2027 Assert(!has_ghost_elements(), ExcGhostsPresent()); 2028 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size())); 2029 2030 AssertIsFinite(s); 2031 2032 // We assume that the vectors have the same Map 2033 // if the local size is the same and if the vectors are not ghosted 2034 if (local_size() == v.local_size() && !v.has_ghost_elements()) 2035 { 2036 Assert(this->vector->Map().SameAs(v.vector->Map()) == true, 2037 ExcDifferentParallelPartitioning()); 2038 const int ierr = vector->Update(1., *(v.vector), s); 2039 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 2040 } 2041 else 2042 { 2043 (*this) *= s; 2044 this->add(v, true); 2045 } 2046 } 2047 2048 2049 2050 inline void sadd(const TrilinosScalar s,const TrilinosScalar a,const Vector & v)2051 Vector::sadd(const TrilinosScalar s, 2052 const TrilinosScalar a, 2053 const Vector & v) 2054 { 2055 // if we have ghost values, do not allow 2056 // writing to this vector at all. 2057 Assert(!has_ghost_elements(), ExcGhostsPresent()); 2058 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size())); 2059 AssertIsFinite(s); 2060 AssertIsFinite(a); 2061 2062 // We assume that the vectors have the same Map 2063 // if the local size is the same and if the vectors are not ghosted 2064 if (local_size() == v.local_size() && !v.has_ghost_elements()) 2065 { 2066 Assert(this->vector->Map().SameAs(v.vector->Map()) == true, 2067 ExcDifferentParallelPartitioning()); 2068 const int ierr = vector->Update(a, *(v.vector), s); 2069 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 2070 } 2071 else 2072 { 2073 (*this) *= s; 2074 Vector tmp = v; 2075 tmp *= a; 2076 this->add(tmp, true); 2077 } 2078 } 2079 2080 2081 2082 inline void scale(const Vector & factors)2083 Vector::scale(const Vector &factors) 2084 { 2085 // if we have ghost values, do not allow 2086 // writing to this vector at all. 2087 Assert(!has_ghost_elements(), ExcGhostsPresent()); 2088 Assert(local_size() == factors.local_size(), 2089 ExcDimensionMismatch(local_size(), factors.local_size())); 2090 2091 const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0); 2092 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 2093 } 2094 2095 2096 2097 inline void equ(const TrilinosScalar a,const Vector & v)2098 Vector::equ(const TrilinosScalar a, const Vector &v) 2099 { 2100 // if we have ghost values, do not allow 2101 // writing to this vector at all. 2102 Assert(!has_ghost_elements(), ExcGhostsPresent()); 2103 AssertIsFinite(a); 2104 2105 // If we don't have the same map, copy. 2106 if (vector->Map().SameAs(v.vector->Map()) == false) 2107 { 2108 this->sadd(0., a, v); 2109 } 2110 else 2111 { 2112 // Otherwise, just update 2113 int ierr = vector->Update(a, *v.vector, 0.0); 2114 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 2115 2116 last_action = Zero; 2117 } 2118 } 2119 2120 2121 2122 inline const Epetra_MultiVector & trilinos_vector()2123 Vector::trilinos_vector() const 2124 { 2125 return static_cast<const Epetra_MultiVector &>(*vector); 2126 } 2127 2128 2129 2130 inline Epetra_FEVector & trilinos_vector()2131 Vector::trilinos_vector() 2132 { 2133 return *vector; 2134 } 2135 2136 2137 2138 inline const Epetra_BlockMap & trilinos_partitioner()2139 Vector::trilinos_partitioner() const 2140 { 2141 return vector->Map(); 2142 } 2143 2144 2145 2146 inline const MPI_Comm & get_mpi_communicator()2147 Vector::get_mpi_communicator() const 2148 { 2149 static MPI_Comm comm; 2150 2151 # ifdef DEAL_II_WITH_MPI 2152 2153 const Epetra_MpiComm *mpi_comm = 2154 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm()); 2155 comm = mpi_comm->Comm(); 2156 2157 # else 2158 2159 comm = MPI_COMM_SELF; 2160 2161 # endif 2162 2163 return comm; 2164 } 2165 2166 template <typename number> Vector(const IndexSet & parallel_partitioner,const dealii::Vector<number> & v,const MPI_Comm & communicator)2167 Vector::Vector(const IndexSet & parallel_partitioner, 2168 const dealii::Vector<number> &v, 2169 const MPI_Comm & communicator) 2170 { 2171 *this = 2172 Vector(parallel_partitioner.make_trilinos_map(communicator, true), v); 2173 owned_elements = parallel_partitioner; 2174 } 2175 2176 2177 2178 inline Vector & 2179 Vector::operator=(const TrilinosScalar s) 2180 { 2181 AssertIsFinite(s); 2182 2183 int ierr = vector->PutScalar(s); 2184 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 2185 2186 if (nonlocal_vector.get() != nullptr) 2187 { 2188 ierr = nonlocal_vector->PutScalar(0.); 2189 AssertThrow(ierr == 0, ExcTrilinosError(ierr)); 2190 } 2191 2192 return *this; 2193 } 2194 } /* end of namespace MPI */ 2195 2196 # endif /* DOXYGEN */ 2197 2198 } /* end of namespace TrilinosWrappers */ 2199 2200 /*@}*/ 2201 2202 2203 namespace internal 2204 { 2205 namespace LinearOperatorImplementation 2206 { 2207 template <typename> 2208 class ReinitHelper; 2209 2210 /** 2211 * A helper class used internally in linear_operator.h. Specialization for 2212 * TrilinosWrappers::MPI::Vector. 2213 */ 2214 template <> 2215 class ReinitHelper<TrilinosWrappers::MPI::Vector> 2216 { 2217 public: 2218 template <typename Matrix> 2219 static void reinit_range_vector(const Matrix & matrix,TrilinosWrappers::MPI::Vector & v,bool omit_zeroing_entries)2220 reinit_range_vector(const Matrix & matrix, 2221 TrilinosWrappers::MPI::Vector &v, 2222 bool omit_zeroing_entries) 2223 { 2224 v.reinit(matrix.locally_owned_range_indices(), 2225 matrix.get_mpi_communicator(), 2226 omit_zeroing_entries); 2227 } 2228 2229 template <typename Matrix> 2230 static void reinit_domain_vector(const Matrix & matrix,TrilinosWrappers::MPI::Vector & v,bool omit_zeroing_entries)2231 reinit_domain_vector(const Matrix & matrix, 2232 TrilinosWrappers::MPI::Vector &v, 2233 bool omit_zeroing_entries) 2234 { 2235 v.reinit(matrix.locally_owned_domain_indices(), 2236 matrix.get_mpi_communicator(), 2237 omit_zeroing_entries); 2238 } 2239 }; 2240 2241 } // namespace LinearOperatorImplementation 2242 } /* namespace internal */ 2243 2244 2245 2246 /** 2247 * Declare dealii::TrilinosWrappers::MPI::Vector as distributed vector. 2248 */ 2249 template <> 2250 struct is_serial_vector<TrilinosWrappers::MPI::Vector> : std::false_type 2251 {}; 2252 2253 2254 DEAL_II_NAMESPACE_CLOSE 2255 2256 #endif // DEAL_II_WITH_TRILINOS 2257 2258 /*---------------------------- trilinos_vector.h ---------------------------*/ 2259 2260 #endif // dealii_trilinos_vector_h 2261