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_templates_h
17 #define dealii_vector_templates_h
18
19
20 #include <deal.II/base/config.h>
21
22 #include <deal.II/base/numbers.h>
23 #include <deal.II/base/template_constraints.h>
24
25 #include <deal.II/lac/block_vector.h>
26 #include <deal.II/lac/exceptions.h>
27 #include <deal.II/lac/vector.h>
28 #include <deal.II/lac/vector_operations_internal.h>
29
30 #ifdef DEAL_II_WITH_PETSC
31 # include <deal.II/lac/petsc_vector_base.h>
32 #endif
33
34 #ifdef DEAL_II_WITH_TRILINOS
35 # include <deal.II/lac/trilinos_vector.h>
36 #endif
37
38
39 #include <algorithm>
40 #include <cmath>
41 #include <cstdio>
42 #include <cstring>
43 #include <iomanip>
44 #include <iostream>
45 #include <memory>
46
47 DEAL_II_NAMESPACE_OPEN
48
49 template <typename Number>
Vector(const Vector<Number> & v)50 Vector<Number>::Vector(const Vector<Number> &v)
51 : Subscriptor()
52 {
53 *this = v;
54 }
55
56
57
58 template <typename Number>
59 void
apply_givens_rotation(const std::array<Number,3> & csr,const size_type i,const size_type k)60 Vector<Number>::apply_givens_rotation(const std::array<Number, 3> &csr,
61 const size_type i,
62 const size_type k)
63 {
64 auto & V = *this;
65 const Number t = V(i);
66 V(i) = csr[0] * V(i) + csr[1] * V(k);
67 V(k) = -csr[1] * t + csr[0] * V(k);
68 }
69
70
71
72 template <typename Number>
73 template <typename OtherNumber>
Vector(const Vector<OtherNumber> & v)74 Vector<Number>::Vector(const Vector<OtherNumber> &v)
75 {
76 *this = v;
77 }
78
79
80
81 #ifdef DEAL_II_WITH_PETSC
82 namespace internal
83 {
84 template <typename Number>
85 void
copy_petsc_vector(const PETScWrappers::VectorBase & v,::dealii::Vector<Number> & out)86 copy_petsc_vector(const PETScWrappers::VectorBase &v,
87 ::dealii::Vector<Number> & out)
88 {
89 if (v.size() == 0)
90 {
91 out.reinit(0);
92 return;
93 }
94 // Create a sequential PETSc vector and then copy over the entries into
95 // the deal.II vector.
96 Vec sequential_vector;
97 VecScatter scatter_context;
98
99 PetscErrorCode ierr =
100 VecScatterCreateToAll(v, &scatter_context, &sequential_vector);
101 AssertThrow(ierr == 0, ExcPETScError(ierr));
102
103 ierr = VecScatterBegin(
104 scatter_context, v, sequential_vector, INSERT_VALUES, SCATTER_FORWARD);
105 AssertThrow(ierr == 0, ExcPETScError(ierr));
106 ierr = VecScatterEnd(
107 scatter_context, v, sequential_vector, INSERT_VALUES, SCATTER_FORWARD);
108 AssertThrow(ierr == 0, ExcPETScError(ierr));
109
110 PetscScalar *start_ptr;
111 ierr = VecGetArray(sequential_vector, &start_ptr);
112 AssertThrow(ierr == 0, ExcPETScError(ierr));
113
114 const PETScWrappers::VectorBase::size_type v_size = v.size();
115 if (out.size() != v_size)
116 out.reinit(v_size, true);
117
118 internal::VectorOperations::copy(start_ptr,
119 start_ptr + out.size(),
120 out.begin());
121 ierr = VecRestoreArray(sequential_vector, &start_ptr);
122 AssertThrow(ierr == 0, ExcPETScError(ierr));
123
124 ierr = VecScatterDestroy(&scatter_context);
125 AssertNothrow(ierr == 0, ExcPETScError(ierr));
126 ierr = VecDestroy(&sequential_vector);
127 AssertNothrow(ierr == 0, ExcPETScError(ierr));
128 }
129 } // namespace internal
130
131
132
133 template <typename Number>
Vector(const PETScWrappers::VectorBase & v)134 Vector<Number>::Vector(const PETScWrappers::VectorBase &v)
135 {
136 internal::copy_petsc_vector(v, *this);
137 }
138 #endif
139
140
141 #ifdef DEAL_II_WITH_TRILINOS
142
143 template <typename Number>
Vector(const TrilinosWrappers::MPI::Vector & v)144 Vector<Number>::Vector(const TrilinosWrappers::MPI::Vector &v)
145 : values(v.size())
146 {
147 if (size() != 0)
148 {
149 // Copy the distributed vector to
150 // a local one at all processors
151 // that know about the original vector.
152 // TODO: There could
153 // be a better solution than
154 // this, but it has not yet been
155 // found.
156 TrilinosWrappers::MPI::Vector localized_vector;
157 localized_vector.reinit(complete_index_set(size()),
158 v.get_mpi_communicator());
159 localized_vector.reinit(v, false, true);
160
161 Assert(localized_vector.size() == size(),
162 ExcDimensionMismatch(localized_vector.size(), size()));
163
164 // get a representation of the vector
165 // and copy it
166 TrilinosScalar **start_ptr;
167
168 int ierr = localized_vector.trilinos_vector().ExtractView(&start_ptr);
169 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
170
171 std::copy(start_ptr[0], start_ptr[0] + size(), begin());
172
173 maybe_reset_thread_partitioner();
174 }
175 }
176
177 #endif
178
179
180 template <typename Number>
181 inline Vector<Number> &
182 Vector<Number>::operator=(const Vector<Number> &v)
183 {
184 if (PointerComparison::equal(this, &v))
185 return *this;
186
187 if (size() != v.size())
188 reinit(v, true);
189
190 if (0 < size())
191 {
192 dealii::internal::VectorOperations::Vector_copy<Number, Number> copier(
193 v.begin(), begin());
194 internal::VectorOperations::parallel_for(copier,
195 0,
196 size(),
197 thread_loop_partitioner);
198 }
199
200 return *this;
201 }
202
203
204
205 template <typename Number>
206 template <typename Number2>
207 inline Vector<Number> &
208 Vector<Number>::operator=(const Vector<Number2> &v)
209 {
210 if (size() != v.size())
211 reinit(v, true);
212
213 dealii::internal::VectorOperations::Vector_copy<Number, Number2> copier(
214 v.begin(), begin());
215 internal::VectorOperations::parallel_for(copier,
216 0,
217 size(),
218 thread_loop_partitioner);
219
220 return *this;
221 }
222
223
224
225 template <typename Number>
226 inline void
reinit(const size_type n,const bool omit_zeroing_entries)227 Vector<Number>::reinit(const size_type n, const bool omit_zeroing_entries)
228 {
229 do_reinit(n, omit_zeroing_entries, true);
230 }
231
232
233
234 template <typename Number>
235 inline void
grow_or_shrink(const size_type n)236 Vector<Number>::grow_or_shrink(const size_type n)
237 {
238 values.resize(n);
239
240 maybe_reset_thread_partitioner();
241 }
242
243
244
245 template <typename Number>
246 template <typename Number2>
247 void
reinit(const Vector<Number2> & v,const bool omit_zeroing_entries)248 Vector<Number>::reinit(const Vector<Number2> &v,
249 const bool omit_zeroing_entries)
250 {
251 do_reinit(v.size(), omit_zeroing_entries, false);
252 thread_loop_partitioner = v.thread_loop_partitioner;
253 }
254
255
256
257 template <typename Number>
258 bool
all_zero()259 Vector<Number>::all_zero() const
260 {
261 Assert(size() != 0, ExcEmptyObject());
262
263 for (size_type i = 0; i < size(); ++i)
264 if (values[i] != Number())
265 return false;
266 return true;
267 }
268
269
270
271 template <typename Number>
272 bool
is_non_negative()273 Vector<Number>::is_non_negative() const
274 {
275 Assert(size() != 0, ExcEmptyObject());
276
277 for (size_type i = 0; i < size(); ++i)
278 if (!internal::VectorOperations::is_non_negative(values[i]))
279 return false;
280
281 return true;
282 }
283
284
285
286 template <typename Number>
287 Vector<Number> &
288 Vector<Number>::operator=(const Number s)
289 {
290 AssertIsFinite(s);
291 if (s != Number())
292 Assert(size() != 0, ExcEmptyObject());
293
294 if (size() > 0)
295 {
296 internal::VectorOperations::Vector_set<Number> setter(s, values.begin());
297 internal::VectorOperations::parallel_for(setter,
298 0,
299 size(),
300 thread_loop_partitioner);
301 }
302
303 return *this;
304 }
305
306
307
308 template <typename Number>
309 Vector<Number> &
310 Vector<Number>::operator*=(const Number factor)
311 {
312 AssertIsFinite(factor);
313
314 Assert(size() != 0, ExcEmptyObject());
315
316 internal::VectorOperations::Vectorization_multiply_factor<Number>
317 vector_multiply(values.begin(), factor);
318
319 internal::VectorOperations::parallel_for(vector_multiply,
320 0,
321 size(),
322 thread_loop_partitioner);
323
324 return *this;
325 }
326
327
328
329 template <typename Number>
330 void
add(const Number a,const Vector<Number> & v)331 Vector<Number>::add(const Number a, const Vector<Number> &v)
332 {
333 AssertIsFinite(a);
334
335 Assert(size() != 0, ExcEmptyObject());
336 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
337
338 internal::VectorOperations::Vectorization_add_av<Number> vector_add_av(
339 values.begin(), v.values.begin(), a);
340 internal::VectorOperations::parallel_for(vector_add_av,
341 0,
342 size(),
343 thread_loop_partitioner);
344 }
345
346
347
348 template <typename Number>
349 void
sadd(const Number x,const Number a,const Vector<Number> & v)350 Vector<Number>::sadd(const Number x, const Number a, const Vector<Number> &v)
351 {
352 AssertIsFinite(x);
353 AssertIsFinite(a);
354
355 Assert(size() != 0, ExcEmptyObject());
356 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
357
358 internal::VectorOperations::Vectorization_sadd_xav<Number> vector_sadd_xav(
359 values.begin(), v.values.begin(), a, x);
360 internal::VectorOperations::parallel_for(vector_sadd_xav,
361 0,
362 size(),
363 thread_loop_partitioner);
364 }
365
366
367
368 template <typename Number>
369 template <typename Number2>
370 Number Vector<Number>::operator*(const Vector<Number2> &v) const
371 {
372 Assert(size() != 0, ExcEmptyObject());
373
374 if (PointerComparison::equal(this, &v))
375 return norm_sqr();
376
377 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
378
379 Number sum;
380 internal::VectorOperations::Dot<Number, Number2> dot(values.begin(),
381 v.values.begin());
382 internal::VectorOperations::parallel_reduce(
383 dot, 0, size(), sum, thread_loop_partitioner);
384 AssertIsFinite(sum);
385
386 return sum;
387 }
388
389
390
391 template <typename Number>
392 typename Vector<Number>::real_type
norm_sqr()393 Vector<Number>::norm_sqr() const
394 {
395 Assert(size() != 0, ExcEmptyObject());
396
397 real_type sum;
398 internal::VectorOperations::Norm2<Number, real_type> norm2(values.begin());
399 internal::VectorOperations::parallel_reduce(
400 norm2, 0, size(), sum, thread_loop_partitioner);
401
402 AssertIsFinite(sum);
403
404 return sum;
405 }
406
407
408
409 template <typename Number>
410 Number
mean_value()411 Vector<Number>::mean_value() const
412 {
413 Assert(size() != 0, ExcEmptyObject());
414
415 Number sum;
416 internal::VectorOperations::MeanValue<Number> mean(values.begin());
417 internal::VectorOperations::parallel_reduce(
418 mean, 0, size(), sum, thread_loop_partitioner);
419
420 return sum / real_type(size());
421 }
422
423
424
425 template <typename Number>
426 typename Vector<Number>::real_type
l1_norm()427 Vector<Number>::l1_norm() const
428 {
429 Assert(size() != 0, ExcEmptyObject());
430
431 real_type sum;
432 internal::VectorOperations::Norm1<Number, real_type> norm1(values.begin());
433 internal::VectorOperations::parallel_reduce(
434 norm1, 0, size(), sum, thread_loop_partitioner);
435
436 return sum;
437 }
438
439
440
441 template <typename Number>
442 typename Vector<Number>::real_type
l2_norm()443 Vector<Number>::l2_norm() const
444 {
445 // if l2_norm()^2 is finite and non-zero, the answer is computed as
446 // std::sqrt(norm_sqr()). If norm_sqr() is infinite or zero, the l2 norm
447 // might still be finite. In that case, recompute it (this is a rare case,
448 // so working on the vector twice is uncritical and paid off by the extended
449 // precision) using the BLAS approach with a weight, see e.g. dnrm2.f.
450 Assert(size() != 0, ExcEmptyObject());
451
452 real_type norm_square;
453 internal::VectorOperations::Norm2<Number, real_type> norm2(values.begin());
454 internal::VectorOperations::parallel_reduce(
455 norm2, 0, size(), norm_square, thread_loop_partitioner);
456 if (numbers::is_finite(norm_square) &&
457 norm_square >= std::numeric_limits<real_type>::min())
458 return static_cast<typename Vector<Number>::real_type>(
459 std::sqrt(norm_square));
460 else
461 {
462 real_type scale = 0.;
463 real_type sum = 1.;
464 for (size_type i = 0; i < size(); ++i)
465 {
466 if (values[i] != Number())
467 {
468 const real_type abs_x =
469 numbers::NumberTraits<Number>::abs(values[i]);
470 if (scale < abs_x)
471 {
472 sum = 1 + sum * (scale / abs_x) * (scale / abs_x);
473 scale = abs_x;
474 }
475 else
476 sum += (abs_x / scale) * (abs_x / scale);
477 }
478 }
479 AssertIsFinite(scale * std::sqrt(sum));
480 return static_cast<typename Vector<Number>::real_type>(scale *
481 std::sqrt(sum));
482 }
483 }
484
485
486
487 template <typename Number>
488 typename Vector<Number>::real_type
lp_norm(const real_type p)489 Vector<Number>::lp_norm(const real_type p) const
490 {
491 Assert(size() != 0, ExcEmptyObject());
492
493 if (p == 1.)
494 return l1_norm();
495 else if (p == 2.)
496 return l2_norm();
497
498 real_type sum;
499 internal::VectorOperations::NormP<Number, real_type> normp(values.begin(), p);
500 internal::VectorOperations::parallel_reduce(
501 normp, 0, size(), sum, thread_loop_partitioner);
502
503 if (numbers::is_finite(sum) && sum >= std::numeric_limits<real_type>::min())
504 return std::pow(sum, static_cast<real_type>(1. / p));
505 else
506 {
507 real_type scale = 0.;
508 real_type sum = 1.;
509 for (size_type i = 0; i < size(); ++i)
510 {
511 if (values[i] != Number())
512 {
513 const real_type abs_x =
514 numbers::NumberTraits<Number>::abs(values[i]);
515 if (scale < abs_x)
516 {
517 sum = 1. + sum * std::pow(scale / abs_x, p);
518 scale = abs_x;
519 }
520 else
521 sum += std::pow(abs_x / scale, p);
522 }
523 }
524 return scale * std::pow(sum, static_cast<real_type>(1. / p));
525 }
526 }
527
528
529
530 template <typename Number>
531 typename Vector<Number>::real_type
linfty_norm()532 Vector<Number>::linfty_norm() const
533 {
534 Assert(size() != 0, ExcEmptyObject());
535
536 real_type max = 0.;
537
538 for (size_type i = 0; i < size(); ++i)
539 max = std::max(numbers::NumberTraits<Number>::abs(values[i]), max);
540
541 return max;
542 }
543
544
545
546 template <typename Number>
547 Number
add_and_dot(const Number a,const Vector<Number> & V,const Vector<Number> & W)548 Vector<Number>::add_and_dot(const Number a,
549 const Vector<Number> &V,
550 const Vector<Number> &W)
551 {
552 Assert(size() != 0, ExcEmptyObject());
553 AssertDimension(size(), V.size());
554 AssertDimension(size(), W.size());
555
556 Number sum;
557 internal::VectorOperations::AddAndDot<Number> adder(values.begin(),
558 V.values.begin(),
559 W.values.begin(),
560 a);
561 internal::VectorOperations::parallel_reduce(
562 adder, 0, size(), sum, thread_loop_partitioner);
563 AssertIsFinite(sum);
564
565 return sum;
566 }
567
568
569
570 template <typename Number>
571 Vector<Number> &
572 Vector<Number>::operator+=(const Vector<Number> &v)
573 {
574 Assert(size() != 0, ExcEmptyObject());
575 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
576
577 internal::VectorOperations::Vectorization_add_v<Number> vector_add(
578 values.begin(), v.values.begin());
579 internal::VectorOperations::parallel_for(vector_add,
580 0,
581 size(),
582 thread_loop_partitioner);
583 return *this;
584 }
585
586
587
588 template <typename Number>
589 Vector<Number> &
590 Vector<Number>::operator-=(const Vector<Number> &v)
591 {
592 Assert(size() != 0, ExcEmptyObject());
593 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
594
595 internal::VectorOperations::Vectorization_subtract_v<Number> vector_subtract(
596 values.begin(), v.values.begin());
597 internal::VectorOperations::parallel_for(vector_subtract,
598 0,
599 size(),
600 thread_loop_partitioner);
601
602 return *this;
603 }
604
605
606
607 template <typename Number>
608 void
add(const Number v)609 Vector<Number>::add(const Number v)
610 {
611 Assert(size() != 0, ExcEmptyObject());
612
613 internal::VectorOperations::Vectorization_add_factor<Number> vector_add(
614 values.begin(), v);
615 internal::VectorOperations::parallel_for(vector_add,
616 0,
617 size(),
618 thread_loop_partitioner);
619 }
620
621
622
623 template <typename Number>
624 void
add(const Number a,const Vector<Number> & v,const Number b,const Vector<Number> & w)625 Vector<Number>::add(const Number a,
626 const Vector<Number> &v,
627 const Number b,
628 const Vector<Number> &w)
629 {
630 AssertIsFinite(a);
631 AssertIsFinite(b);
632
633 Assert(size() != 0, ExcEmptyObject());
634 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
635 Assert(size() == w.size(), ExcDimensionMismatch(size(), w.size()));
636
637 internal::VectorOperations::Vectorization_add_avpbw<Number> vector_add(
638 values.begin(), v.values.begin(), w.values.begin(), a, b);
639 internal::VectorOperations::parallel_for(vector_add,
640 0,
641 size(),
642 thread_loop_partitioner);
643 }
644
645
646
647 template <typename Number>
648 void
sadd(const Number x,const Vector<Number> & v)649 Vector<Number>::sadd(const Number x, const Vector<Number> &v)
650 {
651 AssertIsFinite(x);
652
653 Assert(size() != 0, ExcEmptyObject());
654 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
655
656 internal::VectorOperations::Vectorization_sadd_xv<Number> vector_sadd(
657 values.begin(), v.values.begin(), x);
658 internal::VectorOperations::parallel_for(vector_sadd,
659 0,
660 size(),
661 thread_loop_partitioner);
662 }
663
664
665
666 template <typename Number>
667 void
scale(const Vector<Number> & s)668 Vector<Number>::scale(const Vector<Number> &s)
669 {
670 Assert(size() != 0, ExcEmptyObject());
671 Assert(size() == s.size(), ExcDimensionMismatch(size(), s.size()));
672
673 internal::VectorOperations::Vectorization_scale<Number> vector_scale(
674 values.begin(), s.values.begin());
675 internal::VectorOperations::parallel_for(vector_scale,
676 0,
677 size(),
678 thread_loop_partitioner);
679 }
680
681
682
683 template <typename Number>
684 template <typename Number2>
685 void
scale(const Vector<Number2> & s)686 Vector<Number>::scale(const Vector<Number2> &s)
687 {
688 Assert(size() != 0, ExcEmptyObject());
689 Assert(size() == s.size(), ExcDimensionMismatch(size(), s.size()));
690
691 for (size_type i = 0; i < size(); ++i)
692 values[i] *= Number(s.values[i]);
693 }
694
695
696
697 template <typename Number>
698 void
equ(const Number a,const Vector<Number> & u)699 Vector<Number>::equ(const Number a, const Vector<Number> &u)
700 {
701 AssertIsFinite(a);
702
703 Assert(size() != 0, ExcEmptyObject());
704 Assert(size() == u.size(), ExcDimensionMismatch(size(), u.size()));
705
706 internal::VectorOperations::Vectorization_equ_au<Number> vector_equ(
707 values.begin(), u.values.begin(), a);
708 internal::VectorOperations::parallel_for(vector_equ,
709 0,
710 size(),
711 thread_loop_partitioner);
712 }
713
714
715
716 template <typename Number>
717 template <typename Number2>
718 void
equ(const Number a,const Vector<Number2> & u)719 Vector<Number>::equ(const Number a, const Vector<Number2> &u)
720 {
721 AssertIsFinite(a);
722
723 Assert(size() != 0, ExcEmptyObject());
724 Assert(size() == u.size(), ExcDimensionMismatch(size(), u.size()));
725
726 // set the result vector to a*u. we have to
727 // convert the elements of u to the type of
728 // the result vector. this is necessary
729 // because
730 // operator*(complex<float>,complex<double>)
731 // is not defined by default
732 for (size_type i = 0; i < size(); ++i)
733 values[i] = a * Number(u.values[i]);
734 }
735
736
737
738 template <typename Number>
739 Vector<Number> &
740 Vector<Number>::operator=(const BlockVector<Number> &v)
741 {
742 if (v.size() != size())
743 reinit(v.size(), true);
744
745 size_type this_index = 0;
746 for (size_type b = 0; b < v.n_blocks(); ++b)
747 for (size_type i = 0; i < v.block(b).size(); ++i, ++this_index)
748 values[this_index] = v.block(b)(i);
749
750 return *this;
751 }
752
753
754
755 #ifdef DEAL_II_WITH_PETSC
756 template <typename Number>
757 Vector<Number> &
758 Vector<Number>::operator=(const PETScWrappers::VectorBase &v)
759 {
760 internal::copy_petsc_vector(v, *this);
761 return *this;
762 }
763 #endif
764
765
766 #ifdef DEAL_II_WITH_TRILINOS
767
768 template <typename Number>
769 Vector<Number> &
770 Vector<Number>::operator=(const TrilinosWrappers::MPI::Vector &v)
771 {
772 if (v.size() != size())
773 reinit(v.size(), true);
774 if (size() != 0)
775 {
776 // Copy the distributed vector to
777 // a local one at all processors
778 // that know about the original vector.
779 // TODO: There could
780 // be a better solution than
781 // this, but it has not yet been
782 // found.
783 TrilinosWrappers::MPI::Vector localized_vector;
784 localized_vector.reinit(complete_index_set(size()),
785 v.get_mpi_communicator());
786 localized_vector.reinit(v, false, true);
787
788 Assert(localized_vector.size() == size(),
789 ExcDimensionMismatch(localized_vector.size(), size()));
790
791 // get a representation of the vector
792 // and copy it
793 TrilinosScalar **start_ptr;
794
795 int ierr = localized_vector.trilinos_vector().ExtractView(&start_ptr);
796 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
797
798 std::copy(start_ptr[0], start_ptr[0] + size(), begin());
799 }
800
801 return *this;
802 }
803
804 #endif
805
806 template <typename Number>
807 template <typename Number2>
808 bool
809 Vector<Number>::operator==(const Vector<Number2> &v) const
810 {
811 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
812
813 // compare the two vector. we have to
814 // convert the elements of v to the type of
815 // the result vector. this is necessary
816 // because
817 // operator==(complex<float>,complex<double>)
818 // is not defined by default
819 for (size_type i = 0; i < size(); ++i)
820 if (values[i] != Number(v.values[i]))
821 return false;
822
823 return true;
824 }
825
826
827
828 template <typename Number>
829 void
print(std::ostream & out,const unsigned int precision,const bool scientific,const bool across)830 Vector<Number>::print(std::ostream & out,
831 const unsigned int precision,
832 const bool scientific,
833 const bool across) const
834 {
835 Assert(size() != 0, ExcEmptyObject());
836 AssertThrow(out, ExcIO());
837
838 std::ios::fmtflags old_flags = out.flags();
839 unsigned int old_precision = out.precision(precision);
840
841 out.precision(precision);
842 if (scientific)
843 out.setf(std::ios::scientific, std::ios::floatfield);
844 else
845 out.setf(std::ios::fixed, std::ios::floatfield);
846
847 if (across)
848 for (size_type i = 0; i < size(); ++i)
849 out << values[i] << ' ';
850 else
851 for (size_type i = 0; i < size(); ++i)
852 out << values[i] << std::endl;
853 out << std::endl;
854
855 AssertThrow(out, ExcIO());
856 // reset output format
857 out.flags(old_flags);
858 out.precision(old_precision);
859 }
860
861
862
863 template <typename Number>
864 void
block_write(std::ostream & out)865 Vector<Number>::block_write(std::ostream &out) const
866 {
867 AssertThrow(out, ExcIO());
868
869 // other version of the following
870 // out << size() << std::endl << '[';
871
872 // reason: operator<< seems to use some resources that lead to problems in a
873 // multithreaded environment. We convert the size index to
874 // unsigned long long int that is at least 64 bits to be able to output it on
875 // all platforms, since std::uint64_t is not in C.
876 const unsigned long long int sz = size();
877 char buf[16];
878
879 std::sprintf(buf, "%llu", sz);
880 std::strcat(buf, "\n[");
881
882 out.write(buf, std::strlen(buf));
883 out.write(reinterpret_cast<const char *>(begin()),
884 reinterpret_cast<const char *>(end()) -
885 reinterpret_cast<const char *>(begin()));
886
887 // out << ']';
888 const char outro = ']';
889 out.write(&outro, 1);
890
891 AssertThrow(out, ExcIO());
892 }
893
894
895
896 template <typename Number>
897 void
block_read(std::istream & in)898 Vector<Number>::block_read(std::istream &in)
899 {
900 AssertThrow(in, ExcIO());
901
902 size_type sz;
903
904 char buf[16];
905
906
907 in.getline(buf, 16, '\n');
908 sz = std::atoi(buf);
909
910 // fast initialization, since the
911 // data elements are overwritten anyway
912 reinit(sz, true);
913
914 char c;
915 // in >> c;
916 in.read(&c, 1);
917 AssertThrow(c == '[', ExcIO());
918
919 in.read(reinterpret_cast<char *>(begin()),
920 reinterpret_cast<const char *>(end()) -
921 reinterpret_cast<const char *>(begin()));
922
923 // in >> c;
924 in.read(&c, 1);
925 AssertThrow(c == ']', ExcIO());
926 }
927
928
929
930 template <typename Number>
931 IndexSet
locally_owned_elements()932 Vector<Number>::locally_owned_elements() const
933 {
934 return complete_index_set(size());
935 }
936
937
938
939 template <typename Number>
940 std::size_t
memory_consumption()941 Vector<Number>::memory_consumption() const
942 {
943 return sizeof(*this) + values.memory_consumption() - sizeof(values);
944 }
945
946
947
948 template <typename Number>
949 void
maybe_reset_thread_partitioner()950 Vector<Number>::maybe_reset_thread_partitioner()
951 {
952 if (size() >= 4 * internal::VectorImplementation::minimum_parallel_grain_size)
953 {
954 if (thread_loop_partitioner == nullptr)
955 thread_loop_partitioner =
956 std::make_shared<parallel::internal::TBBPartitioner>();
957 }
958 else
959 thread_loop_partitioner.reset();
960 }
961
962
963
964 template <typename Number>
965 void
do_reinit(const size_type new_size,const bool omit_zeroing_entries,const bool reset_partitioner)966 Vector<Number>::do_reinit(const size_type new_size,
967 const bool omit_zeroing_entries,
968 const bool reset_partitioner)
969 {
970 if (new_size <= size())
971 {
972 if (new_size == 0)
973 {
974 values.clear();
975 }
976 else
977 {
978 values.resize_fast(new_size);
979 if (!omit_zeroing_entries)
980 values.fill();
981 }
982 }
983 else
984 {
985 // otherwise size() < new_size and we must allocate
986 AlignedVector<Number> new_values;
987 new_values.resize_fast(new_size);
988 if (!omit_zeroing_entries)
989 new_values.fill();
990 new_values.swap(values);
991 }
992
993 if (reset_partitioner)
994 maybe_reset_thread_partitioner();
995 }
996
997
998 DEAL_II_NAMESPACE_CLOSE
999
1000 #endif
1001