1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 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_vector_h
17 #define dealii_vector_h
18
19
20 #include <deal.II/base/config.h>
21
22 #include <deal.II/base/aligned_vector.h>
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/index_set.h>
25 #include <deal.II/base/subscriptor.h>
26
27 #include <deal.II/differentiation/ad/ad_number_traits.h>
28
29 #include <deal.II/lac/vector_operation.h>
30 #include <deal.II/lac/vector_type_traits.h>
31
32 #include <boost/serialization/split_member.hpp>
33
34 #include <algorithm>
35 #include <initializer_list>
36 #include <iosfwd>
37 #include <iterator>
38 #include <vector>
39
40 DEAL_II_NAMESPACE_OPEN
41
42
43 // Forward declarations
44 #ifndef DOXYGEN
45 # ifdef DEAL_II_WITH_PETSC
46 namespace PETScWrappers
47 {
48 class VectorBase;
49 }
50 # endif
51
52 # ifdef DEAL_II_WITH_TRILINOS
53 namespace TrilinosWrappers
54 {
55 namespace MPI
56 {
57 class Vector;
58 }
59 } // namespace TrilinosWrappers
60 # endif
61
62 template <typename number>
63 class LAPACKFullMatrix;
64
65 template <typename>
66 class BlockVector;
67
68 namespace parallel
69 {
70 namespace internal
71 {
72 class TBBPartitioner;
73 }
74 } // namespace parallel
75 #endif
76
77
78 /*! @addtogroup Vectors
79 *@{
80 */
81
82 /**
83 * A class that represents a vector of numerical elements. As for the
84 * other classes, in the
85 * @ref Vectors
86 * group, this class has a substantial
87 * number of member functions. These include:
88 * - functions that initialize the vector or change its size;
89 * - functions that compute properties of the vector, such as a variety of
90 * norms;
91 * - functions that allow reading from or writing to individual elements of the
92 * vector;
93 * - functions that implement algebraic operations for vectors, such as
94 * addition of vectors; and
95 * - functions that allow inputting and outputting the data stored by vectors.
96 *
97 * In contrast to the C++ standard library class `std::vector`, this class
98 * intends to implement not simply an array that allows access to its elements,
99 * but indeed a vector that is a member of the mathematical concept of a
100 * "vector space" suitable for numerical computations.
101 *
102 * @note Instantiations for this template are provided for <tt>@<float@>,
103 * @<double@>, @<std::complex@<float@>@>, @<std::complex@<double@>@></tt>;
104 * others can be generated in application programs (see the section on
105 * @ref Instantiations
106 * in the manual).
107 */
108 template <typename Number>
109 class Vector : public Subscriptor
110 {
111 public:
112 // The assertion in vector.templates.h for whether or not a number is
113 // finite is not compatible for AD number types.
114 static_assert(
115 !Differentiation::AD::is_ad_number<Number>::value,
116 "The Vector class does not support auto-differentiable numbers.");
117
118 /**
119 * Declare standard types used in all containers. These types parallel those
120 * in the <tt>C++</tt> standard libraries <tt>vector<...></tt> class.
121 */
122 using value_type = Number;
123 using pointer = value_type *;
124 using const_pointer = const value_type *;
125 using iterator = value_type *;
126 using const_iterator = const value_type *;
127 using reference = value_type &;
128 using const_reference = const value_type &;
129 using size_type = types::global_dof_index;
130
131 /**
132 * Declare a type that has holds real-valued numbers with the same precision
133 * as the template argument to this class. If the template argument of this
134 * class is a real data type, then real_type equals the template argument.
135 * If the template argument is a std::complex type then real_type equals the
136 * type underlying the complex numbers.
137 *
138 * This alias is used to represent the return type of norms.
139 */
140 using real_type = typename numbers::NumberTraits<Number>::real_type;
141
142 /**
143 * @name Basic object handling
144 */
145 //@{
146 /**
147 * Constructor. Create a vector of dimension zero.
148 */
149 Vector();
150
151 /**
152 * Copy constructor. Sets the dimension to that of the given vector, and
153 * copies all elements.
154 *
155 * We would like to make this constructor explicit, but standard containers
156 * insist on using it implicitly.
157 *
158 * @dealiiOperationIsMultithreaded
159 */
160 Vector(const Vector<Number> &v);
161
162 /**
163 * Move constructor. Creates a new vector by stealing the internal data of
164 * the vector @p v.
165 */
166 Vector(Vector<Number> &&v) noexcept = default;
167
168 /**
169 * Copy constructor taking a vector of another data type.
170 *
171 * This constructor will fail to compile if
172 * there is no conversion path from @p OtherNumber to @p Number. You may
173 * lose accuracy when copying to a vector with data elements with
174 * less accuracy.
175 */
176 template <typename OtherNumber>
177 explicit Vector(const Vector<OtherNumber> &v);
178
179 /**
180 * Copy constructor taking an object of type `std::initializer_list`. This
181 * constructor can be used to initialize a vector using a brace-enclosed
182 * list of numbers, such as in the following example:
183 * @code
184 * Vector<double> v({1,2,3});
185 * @endcode
186 * This creates a vector of size 3, whose (double precision) elements have
187 * values 1.0, 2.0, and 3.0.
188 *
189 * This constructor will fail to compile if
190 * there is no conversion path from @p OtherNumber to @p Number. You may
191 * lose accuracy when copying to a vector with data elements with
192 * less accuracy.
193 */
194 template <typename OtherNumber>
195 explicit Vector(const std::initializer_list<OtherNumber> &v);
196
197 #ifdef DEAL_II_WITH_PETSC
198 /**
199 * Another copy constructor: copy the values from a PETSc vector class. This
200 * copy constructor is only available if PETSc was detected during
201 * configuration time.
202 *
203 * Note that due to the communication model used in MPI, this operation can
204 * only succeed if all processes do it at the same time when <code>v</code>
205 * is a distributed vector: It is not possible for only one process to
206 * obtain a copy of a parallel vector while the other jobs do something
207 * else.
208 */
209 explicit Vector(const PETScWrappers::VectorBase &v);
210 #endif
211
212 #ifdef DEAL_II_WITH_TRILINOS
213 /**
214 * Another copy constructor: copy the values from a Trilinos wrapper vector.
215 * This copy constructor is only available if Trilinos was detected during
216 * configuration time.
217 *
218 * @note Due to the communication model used in MPI, this operation can
219 * only succeed if all processes that have knowledge of @p v
220 * (i.e. those given by <code>v.get_mpi_communicator()</code>) do it at
221 * the same time. This means that unless you use a split MPI communicator
222 * then it is not normally possible for only one or a subset of processes
223 * to obtain a copy of a parallel vector while the other jobs do something
224 * else. In other words, calling this function is a 'collective operation'
225 * that needs to be executed by all MPI processes that jointly share @p v.
226 */
227 explicit Vector(const TrilinosWrappers::MPI::Vector &v);
228 #endif
229
230 /**
231 * Constructor. Set dimension to @p n and initialize all elements with zero.
232 *
233 * The constructor is made explicit to avoid accidents like this:
234 * <tt>v=0;</tt>. Presumably, the user wants to set every element of the
235 * vector to zero, but instead, what happens is this call:
236 * <tt>v=Vector@<number@>(0);</tt>, i.e. the vector is replaced by one of
237 * length zero.
238 */
239 explicit Vector(const size_type n);
240
241 /**
242 * Initialize the vector with a given range of values pointed to by the
243 * iterators. This function is there in analogy to the @p std::vector class.
244 */
245 template <typename InputIterator>
246 Vector(const InputIterator first, const InputIterator last);
247
248 /**
249 * Destructor, deallocates memory. Made virtual to allow for derived classes
250 * to behave properly.
251 */
252 virtual ~Vector() override = default;
253
254 /**
255 * This function does nothing but exists for compatibility with the parallel
256 * vector classes.
257 *
258 * For the parallel vector wrapper class, this function compresses the
259 * underlying representation of the vector, i.e. flushes the buffers of the
260 * vector object if it has any. This function is necessary after writing
261 * into a vector element-by-element and before anything else can be done on
262 * it.
263 *
264 * However, for the implementation of this class, it is immaterial and thus
265 * an empty function.
266 */
267 void
268 compress(::dealii::VectorOperation::values operation =
269 ::dealii::VectorOperation::unknown) const;
270
271 /**
272 * Change the dimension of the vector to @p N. The reserved memory for this
273 * vector remains unchanged if possible, to make things faster; this may
274 * waste some memory, so keep this in mind. However, if <tt>N==0</tt> all
275 * memory is freed, i.e. if you want to resize the vector and release the
276 * memory not needed, you have to first call <tt>reinit(0)</tt> and then
277 * <tt>reinit(N)</tt>. This cited behaviour is analogous to that of the
278 * standard library containers.
279 *
280 * If @p omit_zeroing_entries is false, the vector is filled by zeros.
281 * Otherwise, the elements are left an unspecified state.
282 *
283 * This function is virtual in order to allow for derived classes to handle
284 * memory separately.
285 */
286 virtual void
287 reinit(const size_type N, const bool omit_zeroing_entries = false);
288
289 /**
290 * Same as above, but will preserve the values of vector upon resizing.
291 * If we new size is bigger, we will have
292 * \f[
293 * \mathbf V \rightarrow
294 * \left(
295 * \begin{array}{c}
296 * \mathbf V \\
297 * \mathbf 0
298 * \end{array}
299 * \right)
300 * \f]
301 * whereas if the desired size is smaller, then
302 * \f[
303 * \left(
304 * \begin{array}{c}
305 * \mathbf V_1 \\
306 * \mathbf V_2
307 * \end{array}
308 * \right)
309 * \rightarrow \mathbf V_1
310 * \f]
311 */
312 void
313 grow_or_shrink(const size_type N);
314
315 /**
316 * Apply <a href="https://en.wikipedia.org/wiki/Givens_rotation">Givens
317 * rotation</a>
318 * @p csr (a triplet of cosine, sine and radius, see
319 * Utilities::LinearAlgebra::givens_rotation())
320 * to the vector in the plane spanned by the @p i'th and @p k'th unit vectors.
321 */
322 void
323 apply_givens_rotation(const std::array<Number, 3> &csr,
324 const size_type i,
325 const size_type k);
326
327 /**
328 * Change the dimension to that of the vector @p V. The same applies as for
329 * the other @p reinit function.
330 *
331 * The elements of @p V are not copied, i.e. this function is the same as
332 * calling <tt>reinit (V.size(), omit_zeroing_entries)</tt>.
333 */
334 template <typename Number2>
335 void
336 reinit(const Vector<Number2> &V, const bool omit_zeroing_entries = false);
337
338 /**
339 * Swap the contents of this vector and the other vector @p v. One could do
340 * this operation with a temporary variable and copying over the data
341 * elements, but this function is significantly more efficient since it only
342 * swaps the pointers to the data of the two vectors and therefore does not
343 * need to allocate temporary storage and move data around.
344 *
345 * This function is analogous to the @p swap function of all C++
346 * standard containers. Also, there is a global function <tt>swap(u,v)</tt>
347 * that simply calls <tt>u.swap(v)</tt>, again in analogy to standard
348 * functions.
349 *
350 * This function is virtual in order to allow for derived classes to handle
351 * memory separately.
352 */
353 virtual void
354 swap(Vector<Number> &v);
355
356 /**
357 * Set all components of the vector to the given number @p s.
358 *
359 * Since the semantics of assigning a scalar to a vector are not immediately
360 * clear, this operator should really only be used if you want to set the
361 * entire vector to zero. This allows the intuitive notation <tt>v=0</tt>.
362 * Assigning other values is deprecated and may be disallowed in the future.
363 *
364 * @dealiiOperationIsMultithreaded
365 */
366 Vector<Number> &
367 operator=(const Number s);
368
369 /**
370 * Copy the given vector. Resize the present vector if necessary.
371 *
372 * @dealiiOperationIsMultithreaded
373 */
374 Vector<Number> &
375 operator=(const Vector<Number> &v);
376
377 /**
378 * Move the given vector. This operator replaces the present vector with
379 * the internal data of the vector @p v and resets @p v to the state it would
380 * have after being newly default-constructed.
381 */
382 Vector<Number> &
383 operator=(Vector<Number> &&v) noexcept = default;
384
385 /**
386 * Copy the given vector. Resize the present vector if necessary.
387 *
388 * @dealiiOperationIsMultithreaded
389 */
390 template <typename Number2>
391 Vector<Number> &
392 operator=(const Vector<Number2> &v);
393
394 /**
395 * Copy operator for assigning a block vector to a regular vector.
396 */
397 Vector<Number> &
398 operator=(const BlockVector<Number> &v);
399
400 #ifdef DEAL_II_WITH_PETSC
401 /**
402 * Another copy operator: copy the values from a PETSc wrapper vector
403 * class. This operator is only available if PETSc was detected during
404 * configuration time.
405 *
406 * Note that due to the communication model used in MPI, this operation can
407 * only succeed if all processes do it at the same time when <code>v</code>
408 * is a distributed vector: It is not possible for only one process to
409 * obtain a copy of a parallel vector while the other jobs do something
410 * else.
411 */
412 Vector<Number> &
413 operator=(const PETScWrappers::VectorBase &v);
414 #endif
415
416
417 #ifdef DEAL_II_WITH_TRILINOS
418 /**
419 * Another copy operator: copy the values from a (sequential or parallel,
420 * depending on the underlying compiler) Trilinos wrapper vector class. This
421 * operator is only available if Trilinos was detected during configuration
422 * time.
423 *
424 * @note Due to the communication model used in MPI, this operation can
425 * only succeed if all processes that have knowledge of @p v
426 * (i.e. those given by <code>v.get_mpi_communicator()</code>) do it at
427 * the same time. This means that unless you use a split MPI communicator
428 * then it is not normally possible for only one or a subset of processes
429 * to obtain a copy of a parallel vector while the other jobs do something
430 * else. In other words, calling this function is a 'collective operation'
431 * that needs to be executed by all MPI processes that jointly share @p v.
432 */
433 Vector<Number> &
434 operator=(const TrilinosWrappers::MPI::Vector &v);
435 #endif
436
437 /**
438 * Test for equality. This function assumes that the present vector and the
439 * one to compare with have the same size already, since comparing vectors
440 * of different sizes makes not much sense anyway.
441 */
442 template <typename Number2>
443 bool
444 operator==(const Vector<Number2> &v) const;
445
446 /**
447 * Test for inequality. This function assumes that the present vector and
448 * the one to compare with have the same size already, since comparing
449 * vectors of different sizes makes not much sense anyway.
450 */
451 template <typename Number2>
452 bool
453 operator!=(const Vector<Number2> &v) const;
454
455 //@}
456
457
458 /**
459 * @name Scalar products, norms and related operations
460 */
461 //@{
462
463 /**
464 * Return the scalar product of two vectors. The return type is the
465 * underlying type of @p this vector, so the return type and the accuracy
466 * with which it the result is computed depend on the order of the arguments
467 * of this vector.
468 *
469 * For complex vectors, the scalar product is implemented as
470 * $\left<v,w\right>=\sum_i v_i \bar{w_i}$.
471 *
472 * @dealiiOperationIsMultithreaded The algorithm uses pairwise summation
473 * with the same order of summation in every run, which gives fully
474 * repeatable results from one run to another.
475 */
476 template <typename Number2>
477 Number operator*(const Vector<Number2> &V) const;
478
479 /**
480 * Return the square of the $l_2$-norm.
481 *
482 * @dealiiOperationIsMultithreaded The algorithm uses pairwise summation
483 * with the same order of summation in every run, which gives fully
484 * repeatable results from one run to another.
485 */
486 real_type
487 norm_sqr() const;
488
489 /**
490 * Mean value of the elements of this vector.
491 *
492 * @dealiiOperationIsMultithreaded The algorithm uses pairwise summation
493 * with the same order of summation in every run, which gives fully
494 * repeatable results from one run to another.
495 */
496 Number
497 mean_value() const;
498
499 /**
500 * $l_1$-norm of the vector. The sum of the absolute values.
501 *
502 * @dealiiOperationIsMultithreaded The algorithm uses pairwise summation
503 * with the same order of summation in every run, which gives fully
504 * repeatable results from one run to another.
505 */
506 real_type
507 l1_norm() const;
508
509 /**
510 * $l_2$-norm of the vector. The square root of the sum of the squares of
511 * the elements.
512 *
513 * @dealiiOperationIsMultithreaded The algorithm uses pairwise summation
514 * with the same order of summation in every run, which gives fully
515 * repeatable results from one run to another.
516 */
517 real_type
518 l2_norm() const;
519
520 /**
521 * $l_p$-norm of the vector. The pth root of the sum of the pth powers of
522 * the absolute values of the elements.
523 *
524 * @dealiiOperationIsMultithreaded The algorithm uses pairwise summation
525 * with the same order of summation in every run, which gives fully
526 * repeatable results from one run to another.
527 */
528 real_type
529 lp_norm(const real_type p) const;
530
531 /**
532 * Maximum absolute value of the elements.
533 */
534 real_type
535 linfty_norm() const;
536
537 /**
538 * Performs a combined operation of a vector addition and a subsequent inner
539 * product, returning the value of the inner product. In other words, the
540 * result of this function is the same as if the user called
541 * @code
542 * this->add(a, V);
543 * return_value = *this * W;
544 * @endcode
545 *
546 * The reason this function exists is that this operation involves less
547 * memory transfer than calling the two functions separately. This method
548 * only needs to load three vectors, @p this, @p V, @p W, whereas calling
549 * separate methods means to load the calling vector @p this twice. Since
550 * most vector operations are memory transfer limited, this reduces the time
551 * by 25\% (or 50\% if @p W equals @p this).
552 *
553 * For complex-valued vectors, the scalar product in the second step is
554 * implemented as
555 * $\left<v,w\right>=\sum_i v_i \bar{w_i}$.
556 *
557 * @dealiiOperationIsMultithreaded The algorithm uses pairwise summation
558 * with the same order of summation in every run, which gives fully
559 * repeatable results from one run to another.
560 */
561 Number
562 add_and_dot(const Number a, const Vector<Number> &V, const Vector<Number> &W);
563
564 //@}
565
566
567 /**
568 * @name Data access
569 */
570 //@{
571
572 /**
573 * Return a pointer to the underlying data buffer.
574 */
575 pointer
576 data();
577
578 /**
579 * Return a const pointer to the underlying data buffer.
580 */
581 const_pointer
582 data() const;
583
584 /**
585 * Make the @p Vector class a bit like the <tt>vector<></tt> class of the
586 * C++ standard library by returning iterators to the start and end of the
587 * elements of this vector.
588 */
589 iterator
590 begin();
591
592 /**
593 * Return constant iterator to the start of the vectors.
594 */
595 const_iterator
596 begin() const;
597
598 /**
599 * Return an iterator pointing to the element past the end of the array.
600 */
601 iterator
602 end();
603
604 /**
605 * Return a constant iterator pointing to the element past the end of the
606 * array.
607 */
608 const_iterator
609 end() const;
610
611 /**
612 * Access the value of the @p ith component.
613 */
614 Number
615 operator()(const size_type i) const;
616
617 /**
618 * Access the @p ith component as a writeable reference.
619 */
620 Number &
621 operator()(const size_type i);
622
623 /**
624 * Access the value of the @p ith component.
625 *
626 * Exactly the same as operator().
627 */
628 Number operator[](const size_type i) const;
629
630 /**
631 * Access the @p ith component as a writeable reference.
632 *
633 * Exactly the same as operator().
634 */
635 Number &operator[](const size_type i);
636
637 /**
638 * Instead of getting individual elements of a vector via operator(),
639 * this function allows getting a whole set of elements at once. The
640 * indices of the elements to be read are stated in the first argument, the
641 * corresponding values are returned in the second.
642 *
643 * If the current vector is called @p v, then this function is the equivalent
644 * to the code
645 * @code
646 * for (unsigned int i = 0; i < indices.size(); ++i)
647 * values[i] = v[indices[i]];
648 * @endcode
649 *
650 * @pre The sizes of the @p indices and @p values arrays must be identical.
651 */
652 template <typename OtherNumber>
653 void
654 extract_subvector_to(const std::vector<size_type> &indices,
655 std::vector<OtherNumber> & values) const;
656
657 /**
658 * Instead of getting individual elements of a vector via operator(),
659 * this function allows getting a whole set of elements at once. In
660 * contrast to the previous function, this function obtains the
661 * indices of the elements by dereferencing all elements of the iterator
662 * range provided by the first two arguments, and puts the vector
663 * values into memory locations obtained by dereferencing a range
664 * of iterators starting at the location pointed to by the third
665 * argument.
666 *
667 * If the current vector is called @p v, then this function is the equivalent
668 * to the code
669 * @code
670 * ForwardIterator indices_p = indices_begin;
671 * OutputIterator values_p = values_begin;
672 * while (indices_p != indices_end)
673 * {
674 * *values_p = v[*indices_p];
675 * ++indices_p;
676 * ++values_p;
677 * }
678 * @endcode
679 *
680 * @pre It must be possible to write into as many memory locations
681 * starting at @p values_begin as there are iterators between
682 * @p indices_begin and @p indices_end.
683 */
684 template <typename ForwardIterator, typename OutputIterator>
685 void
686 extract_subvector_to(ForwardIterator indices_begin,
687 const ForwardIterator indices_end,
688 OutputIterator values_begin) const;
689 //@}
690
691
692 /**
693 * @name Modification of vectors
694 */
695 //@{
696
697 /**
698 * Add the given vector to the present one.
699 *
700 * @dealiiOperationIsMultithreaded
701 */
702 Vector<Number> &
703 operator+=(const Vector<Number> &V);
704
705 /**
706 * Subtract the given vector from the present one.
707 *
708 * @dealiiOperationIsMultithreaded
709 */
710 Vector<Number> &
711 operator-=(const Vector<Number> &V);
712
713 /**
714 * A collective add operation: This function adds a whole set of values
715 * stored in @p values to the vector components specified by @p indices.
716 */
717 template <typename OtherNumber>
718 void
719 add(const std::vector<size_type> & indices,
720 const std::vector<OtherNumber> &values);
721
722 /**
723 * This is a second collective add operation. As a difference, this function
724 * takes a deal.II vector of values.
725 */
726 template <typename OtherNumber>
727 void
728 add(const std::vector<size_type> &indices, const Vector<OtherNumber> &values);
729
730 /**
731 * Take an address where <tt>n_elements</tt> are stored contiguously and add
732 * them into the vector. Handles all cases which are not covered by the
733 * other two <tt>add()</tt> functions above.
734 */
735 template <typename OtherNumber>
736 void
737 add(const size_type n_elements,
738 const size_type * indices,
739 const OtherNumber *values);
740
741 /**
742 * Addition of @p s to all components. Note that @p s is a scalar and not a
743 * vector.
744 *
745 * @dealiiOperationIsMultithreaded
746 */
747 void
748 add(const Number s);
749
750 /**
751 * Multiple addition of scaled vectors, i.e. <tt>*this += a*V+b*W</tt>.
752 *
753 * @dealiiOperationIsMultithreaded
754 */
755 void
756 add(const Number a,
757 const Vector<Number> &V,
758 const Number b,
759 const Vector<Number> &W);
760
761 /**
762 * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
763 *
764 * @dealiiOperationIsMultithreaded
765 */
766 void
767 add(const Number a, const Vector<Number> &V);
768
769 /**
770 * Scaling and simple vector addition, i.e. <tt>*this = s*(*this)+V</tt>.
771 *
772 * @dealiiOperationIsMultithreaded
773 */
774 void
775 sadd(const Number s, const Vector<Number> &V);
776
777 /**
778 * Scaling and simple addition, i.e. <tt>*this = s*(*this)+a*V</tt>.
779 *
780 * @dealiiOperationIsMultithreaded
781 */
782 void
783 sadd(const Number s, const Number a, const Vector<Number> &V);
784
785 /**
786 * Scale each element of the vector by a constant value.
787 *
788 * @dealiiOperationIsMultithreaded
789 */
790 Vector<Number> &
791 operator*=(const Number factor);
792
793 /**
794 * Scale each element of the vector by the inverse of the given value.
795 *
796 * @dealiiOperationIsMultithreaded
797 */
798 Vector<Number> &
799 operator/=(const Number factor);
800
801 /**
802 * Scale each element of this vector by the corresponding element in the
803 * argument. This function is mostly meant to simulate multiplication (and
804 * immediate re-assignment) by a diagonal scaling matrix.
805 *
806 * @dealiiOperationIsMultithreaded
807 */
808 void
809 scale(const Vector<Number> &scaling_factors);
810
811 /**
812 * Scale each element of this vector by the corresponding element in the
813 * argument. This function is mostly meant to simulate multiplication (and
814 * immediate re-assignment) by a diagonal scaling matrix.
815 */
816 template <typename Number2>
817 void
818 scale(const Vector<Number2> &scaling_factors);
819
820 /**
821 * Assignment <tt>*this = a*u</tt>.
822 *
823 * @dealiiOperationIsMultithreaded
824 */
825 void
826 equ(const Number a, const Vector<Number> &u);
827
828 /**
829 * Assignment <tt>*this = a*u</tt>.
830 */
831 template <typename Number2>
832 void
833 equ(const Number a, const Vector<Number2> &u);
834
835 /**
836 * This function does nothing but exists for compatibility with the @p
837 * parallel vector classes (e.g., LinearAlgebra::distributed::Vector class).
838 */
839 void
840 update_ghost_values() const;
841 //@}
842
843
844 /**
845 * @name Input and output
846 */
847 //@{
848 /**
849 * Print to a stream. @p precision denotes the desired precision with which
850 * values shall be printed, @p scientific whether scientific notation shall
851 * be used. If @p across is @p true then the vector is printed in a line,
852 * while if @p false then the elements are printed on a separate line each.
853 */
854 void
855 print(std::ostream & out,
856 const unsigned int precision = 3,
857 const bool scientific = true,
858 const bool across = true) const;
859
860 /**
861 * Write the vector en bloc to a file. This is done in a binary mode, so the
862 * output is neither readable by humans nor (probably) by other computers
863 * using a different operating system or number format.
864 */
865 void
866 block_write(std::ostream &out) const;
867
868 /**
869 * Read a vector en block from a file. This is done using the inverse
870 * operations to the above function, so it is reasonably fast because the
871 * bitstream is not interpreted.
872 *
873 * The vector is resized if necessary.
874 *
875 * A primitive form of error checking is performed which will recognize the
876 * bluntest attempts to interpret some data as a vector stored bitwise to a
877 * file, but not more.
878 */
879 void
880 block_read(std::istream &in);
881
882 /**
883 * Write the data of this object to a stream for the purpose of
884 * serialization.
885 */
886 template <class Archive>
887 void
888 save(Archive &ar, const unsigned int version) const;
889
890 /**
891 * Read the data of this object from a stream for the purpose of
892 * serialization.
893 */
894 template <class Archive>
895 void
896 load(Archive &ar, const unsigned int version);
897
898 #ifdef DOXYGEN
899 /**
900 * Write and read the data of this object from a stream for the purpose
901 * of serialization.
902 */
903 template <class Archive>
904 void
905 serialize(Archive &archive, const unsigned int version);
906 #else
907 // This macro defines the serialize() method that is compatible with
908 // the templated save() and load() method that have been implemented.
909 BOOST_SERIALIZATION_SPLIT_MEMBER()
910 #endif
911
912 /**
913 * @}
914 */
915
916 /**
917 * @name Information about the object
918 */
919 //@{
920
921 /**
922 * Return true if the given global index is in the local range of this
923 * processor. Since this is not a distributed vector the method always
924 * returns true.
925 */
926 bool
927 in_local_range(const size_type global_index) const;
928
929 /**
930 * Return an index set that describes which elements of this vector are
931 * owned by the current processor. Note that this index set does not include
932 * elements this vector may store locally as ghost elements but that are in
933 * fact owned by another processor. As a consequence, the index sets
934 * returned on different processors if this is a distributed vector will
935 * form disjoint sets that add up to the complete index set. Obviously, if a
936 * vector is created on only one processor, then the result would satisfy
937 * @code
938 * vec.locally_owned_elements() == complete_index_set (vec.size())
939 * @endcode
940 *
941 * Since the current data type does not support parallel data storage across
942 * different processors, the returned index set is the complete index set.
943 */
944 IndexSet
945 locally_owned_elements() const;
946
947 /**
948 * Return dimension of the vector.
949 */
950 size_type
951 size() const;
952
953 /**
954 * Return whether the vector contains only elements with value zero. This
955 * function is mainly for internal consistency checks and should seldom be
956 * used when not in debug mode since it uses quite some time.
957 */
958 bool
959 all_zero() const;
960
961 /**
962 * Return @p true if the vector has no negative entries, i.e. all entries
963 * are zero or positive. This function is used, for example, to check
964 * whether refinement indicators are really all positive (or zero).
965 *
966 * The function obviously only makes sense if the template argument of this
967 * class is a real type. If it is a complex type, then an exception is
968 * thrown.
969 */
970 bool
971 is_non_negative() const;
972
973 /**
974 * Determine an estimate for the memory consumption (in bytes) of this
975 * object.
976 */
977 std::size_t
978 memory_consumption() const;
979
980 /**
981 * This function exists for compatibility with the @p
982 * parallel vector classes (e.g., LinearAlgebra::distributed::Vector class).
983 * Always returns false since this implementation is serial.
984 */
985 bool
986 has_ghost_elements() const;
987 //@}
988
989 private:
990 /**
991 * Array of elements owned by this vector.
992 */
993 AlignedVector<Number> values;
994
995 /**
996 * Convenience function used at the end of initialization or
997 * reinitialization. Resets (if necessary) the loop partitioner to the
998 * correct state, based on its current state and the length of the vector.
999 */
1000 void
1001 maybe_reset_thread_partitioner();
1002
1003 /**
1004 * Actual implementation of the reinit functions.
1005 */
1006 void
1007 do_reinit(const size_type new_size,
1008 const bool omit_zeroing_entries,
1009 const bool reset_partitioner);
1010
1011 /**
1012 * For parallel loops with TBB, this member variable stores the affinity
1013 * information of loops.
1014 */
1015 mutable std::shared_ptr<parallel::internal::TBBPartitioner>
1016 thread_loop_partitioner;
1017
1018 // Make all other vector types friends.
1019 template <typename Number2>
1020 friend class Vector;
1021 };
1022
1023 /*@}*/
1024 /*----------------------- Inline functions ----------------------------------*/
1025
1026
1027 #ifndef DOXYGEN
1028
1029
1030 //------------------------ declarations for explicit specializations
1031 template <>
1032 Vector<int>::real_type
1033 Vector<int>::lp_norm(const real_type) const;
1034
1035
1036 //------------------------ inline functions
1037
1038 template <typename Number>
Vector()1039 inline Vector<Number>::Vector()
1040 {
1041 // virtual functions called in constructors and destructors never use the
1042 // override in a derived class
1043 // for clarity be explicit on which function is called
1044 Vector<Number>::reinit(0);
1045 }
1046
1047
1048
1049 template <typename Number>
1050 template <typename OtherNumber>
Vector(const std::initializer_list<OtherNumber> & v)1051 Vector<Number>::Vector(const std::initializer_list<OtherNumber> &v)
1052 : Vector(v.begin(), v.end())
1053 {}
1054
1055
1056
1057 template <typename Number>
1058 template <typename InputIterator>
Vector(const InputIterator first,const InputIterator last)1059 Vector<Number>::Vector(const InputIterator first, const InputIterator last)
1060 {
1061 // allocate memory. do not initialize it, as we will copy over to it in a
1062 // second
1063 reinit(std::distance(first, last), true);
1064 std::copy(first, last, begin());
1065 }
1066
1067
1068
1069 template <typename Number>
Vector(const size_type n)1070 inline Vector<Number>::Vector(const size_type n)
1071 {
1072 // virtual functions called in constructors and destructors never use the
1073 // override in a derived class
1074 // for clarity be explicit on which function is called
1075 Vector<Number>::reinit(n, false);
1076 }
1077
1078
1079
1080 template <typename Number>
1081 inline typename Vector<Number>::size_type
size()1082 Vector<Number>::size() const
1083 {
1084 return values.size();
1085 }
1086
1087
1088 template <typename Number>
1089 inline bool
in_local_range(const size_type)1090 Vector<Number>::in_local_range(const size_type) const
1091 {
1092 return true;
1093 }
1094
1095
1096
1097 template <typename Number>
1098 inline typename Vector<Number>::pointer
data()1099 Vector<Number>::data()
1100 {
1101 return values.data();
1102 }
1103
1104
1105
1106 template <typename Number>
1107 inline typename Vector<Number>::const_pointer
data()1108 Vector<Number>::data() const
1109 {
1110 return values.data();
1111 }
1112
1113
1114
1115 template <typename Number>
1116 inline typename Vector<Number>::iterator
begin()1117 Vector<Number>::begin()
1118 {
1119 return values.begin();
1120 }
1121
1122
1123
1124 template <typename Number>
1125 inline typename Vector<Number>::const_iterator
begin()1126 Vector<Number>::begin() const
1127 {
1128 return values.begin();
1129 }
1130
1131
1132
1133 template <typename Number>
1134 inline typename Vector<Number>::iterator
end()1135 Vector<Number>::end()
1136 {
1137 return values.end();
1138 }
1139
1140
1141
1142 template <typename Number>
1143 inline typename Vector<Number>::const_iterator
end()1144 Vector<Number>::end() const
1145 {
1146 return values.end();
1147 }
1148
1149
1150
1151 template <typename Number>
1152 inline Number
operator()1153 Vector<Number>::operator()(const size_type i) const
1154 {
1155 AssertIndexRange(i, size());
1156 return values[i];
1157 }
1158
1159
1160
1161 template <typename Number>
1162 inline Number &
operator()1163 Vector<Number>::operator()(const size_type i)
1164 {
1165 AssertIndexRange(i, size());
1166 return values[i];
1167 }
1168
1169
1170
1171 template <typename Number>
1172 inline Number Vector<Number>::operator[](const size_type i) const
1173 {
1174 return operator()(i);
1175 }
1176
1177
1178
1179 template <typename Number>
1180 inline Number &Vector<Number>::operator[](const size_type i)
1181 {
1182 return operator()(i);
1183 }
1184
1185
1186
1187 template <typename Number>
1188 template <typename OtherNumber>
1189 inline void
extract_subvector_to(const std::vector<size_type> & indices,std::vector<OtherNumber> & values)1190 Vector<Number>::extract_subvector_to(const std::vector<size_type> &indices,
1191 std::vector<OtherNumber> & values) const
1192 {
1193 for (size_type i = 0; i < indices.size(); ++i)
1194 values[i] = operator()(indices[i]);
1195 }
1196
1197
1198
1199 template <typename Number>
1200 template <typename ForwardIterator, typename OutputIterator>
1201 inline void
extract_subvector_to(ForwardIterator indices_begin,const ForwardIterator indices_end,OutputIterator values_begin)1202 Vector<Number>::extract_subvector_to(ForwardIterator indices_begin,
1203 const ForwardIterator indices_end,
1204 OutputIterator values_begin) const
1205 {
1206 while (indices_begin != indices_end)
1207 {
1208 *values_begin = operator()(*indices_begin);
1209 indices_begin++;
1210 values_begin++;
1211 }
1212 }
1213
1214
1215
1216 template <typename Number>
1217 inline Vector<Number> &
1218 Vector<Number>::operator/=(const Number factor)
1219 {
1220 AssertIsFinite(factor);
1221 Assert(factor != Number(0.), ExcZero());
1222
1223 this->operator*=(Number(1.) / factor);
1224 return *this;
1225 }
1226
1227
1228
1229 template <typename Number>
1230 template <typename OtherNumber>
1231 inline void
add(const std::vector<size_type> & indices,const std::vector<OtherNumber> & values)1232 Vector<Number>::add(const std::vector<size_type> & indices,
1233 const std::vector<OtherNumber> &values)
1234 {
1235 Assert(indices.size() == values.size(),
1236 ExcDimensionMismatch(indices.size(), values.size()));
1237 add(indices.size(), indices.data(), values.data());
1238 }
1239
1240
1241
1242 template <typename Number>
1243 template <typename OtherNumber>
1244 inline void
add(const std::vector<size_type> & indices,const Vector<OtherNumber> & values)1245 Vector<Number>::add(const std::vector<size_type> &indices,
1246 const Vector<OtherNumber> & values)
1247 {
1248 Assert(indices.size() == values.size(),
1249 ExcDimensionMismatch(indices.size(), values.size()));
1250 add(indices.size(), indices.data(), values.values.begin());
1251 }
1252
1253
1254
1255 template <typename Number>
1256 template <typename OtherNumber>
1257 inline void
add(const size_type n_indices,const size_type * indices,const OtherNumber * values)1258 Vector<Number>::add(const size_type n_indices,
1259 const size_type * indices,
1260 const OtherNumber *values)
1261 {
1262 for (size_type i = 0; i < n_indices; ++i)
1263 {
1264 AssertIndexRange(indices[i], size());
1265 Assert(
1266 numbers::is_finite(values[i]),
1267 ExcMessage(
1268 "The given value is not finite but either infinite or Not A Number (NaN)"));
1269
1270 this->values[indices[i]] += values[i];
1271 }
1272 }
1273
1274
1275
1276 template <typename Number>
1277 template <typename Number2>
1278 inline bool
1279 Vector<Number>::operator!=(const Vector<Number2> &v) const
1280 {
1281 return !(*this == v);
1282 }
1283
1284
1285
1286 template <typename Number>
compress(::dealii::VectorOperation::values)1287 inline void Vector<Number>::compress(::dealii::VectorOperation::values) const
1288 {}
1289
1290
1291 template <typename Number>
1292 inline bool
has_ghost_elements()1293 Vector<Number>::has_ghost_elements() const
1294 {
1295 return false;
1296 }
1297
1298 template <typename Number>
1299 inline void
update_ghost_values()1300 Vector<Number>::update_ghost_values() const
1301 {}
1302
1303
1304
1305 // Moved from vector.templates.h as an inline function by Luca Heltai
1306 // on 2009/04/12 to prevent strange compiling errors, after making
1307 // swap virtual.
1308 template <typename Number>
1309 inline void
swap(Vector<Number> & v)1310 Vector<Number>::swap(Vector<Number> &v)
1311 {
1312 values.swap(v.values);
1313 std::swap(thread_loop_partitioner, v.thread_loop_partitioner);
1314 }
1315
1316
1317
1318 template <typename Number>
1319 template <class Archive>
1320 inline void
save(Archive & ar,const unsigned int)1321 Vector<Number>::save(Archive &ar, const unsigned int) const
1322 {
1323 // forward to serialization function in the base class.
1324 ar &static_cast<const Subscriptor &>(*this);
1325 ar &values;
1326 }
1327
1328
1329
1330 template <typename Number>
1331 template <class Archive>
1332 inline void
load(Archive & ar,const unsigned int)1333 Vector<Number>::load(Archive &ar, const unsigned int)
1334 {
1335 // the load stuff again from the archive
1336 ar &static_cast<Subscriptor &>(*this);
1337 ar &values;
1338 maybe_reset_thread_partitioner();
1339 }
1340
1341 #endif
1342
1343
1344 /*! @addtogroup Vectors
1345 *@{
1346 */
1347
1348
1349 /**
1350 * Global function @p swap which overloads the default implementation of the
1351 * C++ standard library which uses a temporary object. The function simply
1352 * exchanges the data of the two vectors.
1353 *
1354 * @relatesalso Vector
1355 */
1356 template <typename Number>
1357 inline void
swap(Vector<Number> & u,Vector<Number> & v)1358 swap(Vector<Number> &u, Vector<Number> &v)
1359 {
1360 u.swap(v);
1361 }
1362
1363
1364 /**
1365 * Output operator writing a vector to a stream. This operator outputs the
1366 * elements of the vector one by one, with a space between entries. Each entry
1367 * is formatted according to the flags set on the output stream.
1368 *
1369 * @relatesalso Vector
1370 */
1371 template <typename number>
1372 inline std::ostream &
1373 operator<<(std::ostream &out, const Vector<number> &v)
1374 {
1375 Assert(v.size() != 0, ExcEmptyObject());
1376 AssertThrow(out, ExcIO());
1377
1378 for (typename Vector<number>::size_type i = 0; i < v.size() - 1; ++i)
1379 out << v(i) << ' ';
1380 out << v(v.size() - 1);
1381
1382 AssertThrow(out, ExcIO());
1383
1384 return out;
1385 }
1386
1387 /*@}*/
1388
1389
1390 /**
1391 * Declare dealii::Vector< Number > as serial vector.
1392 *
1393 * @relatesalso Vector
1394 */
1395 template <typename Number>
1396 struct is_serial_vector<Vector<Number>> : std::true_type
1397 {};
1398
1399
1400 DEAL_II_NAMESPACE_CLOSE
1401
1402 #endif
1403