1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 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_la_vector_templates_h
17 #define dealii_la_vector_templates_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/lac/la_vector.h>
22 #include <deal.II/lac/vector_operation.h>
23 #include <deal.II/lac/vector_operations_internal.h>
24 
25 #include <boost/io/ios_state.hpp>
26 
27 #include <iomanip>
28 #include <iostream>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
32 namespace LinearAlgebra
33 {
34   template <typename Number>
35   void
reinit(const size_type size,const bool omit_zeroing_entries)36   Vector<Number>::reinit(const size_type size, const bool omit_zeroing_entries)
37   {
38     ReadWriteVector<Number>::reinit(size, omit_zeroing_entries);
39   }
40 
41 
42 
43   template <typename Number>
44   template <typename Number2>
45   void
reinit(const ReadWriteVector<Number2> & in_vector,const bool omit_zeroing_entries)46   Vector<Number>::reinit(const ReadWriteVector<Number2> &in_vector,
47                          const bool                      omit_zeroing_entries)
48   {
49     ReadWriteVector<Number>::reinit(in_vector, omit_zeroing_entries);
50   }
51 
52 
53 
54   template <typename Number>
55   void
reinit(const IndexSet & locally_stored_indices,const bool omit_zeroing_entries)56   Vector<Number>::reinit(const IndexSet &locally_stored_indices,
57                          const bool      omit_zeroing_entries)
58   {
59     ReadWriteVector<Number>::reinit(locally_stored_indices,
60                                     omit_zeroing_entries);
61   }
62 
63 
64 
65   template <typename Number>
66   void
reinit(const VectorSpaceVector<Number> & V,const bool omit_zeroing_entries)67   Vector<Number>::reinit(const VectorSpaceVector<Number> &V,
68                          const bool                       omit_zeroing_entries)
69   {
70     // Check that casting will work.
71     Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
72            ExcVectorTypeNotCompatible());
73 
74     // Downcast V. If fails, throws an exception.
75     const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
76     Assert(down_V.size() == this->size(),
77            ExcMessage(
78              "Cannot add two vectors with different numbers of elements"));
79 
80     ReadWriteVector<Number>::reinit(down_V, omit_zeroing_entries);
81   }
82 
83 
84 
85   template <typename Number>
86   Vector<Number> &
87   Vector<Number>::operator=(const Vector<Number> &in_vector)
88   {
89     if (PointerComparison::equal(this, &in_vector))
90       return *this;
91 
92     this->thread_loop_partitioner = in_vector.thread_loop_partitioner;
93     if (this->size() != in_vector.size())
94       this->reinit(in_vector, true);
95 
96     dealii::internal::VectorOperations::Vector_copy<Number, Number> copier(
97       in_vector.values.get(), this->values.get());
98     dealii::internal::VectorOperations::parallel_for(
99       copier,
100       static_cast<size_type>(0),
101       this->size(),
102       this->thread_loop_partitioner);
103 
104     return *this;
105   }
106 
107 
108 
109   template <typename Number>
110   template <typename Number2>
111   Vector<Number> &
112   Vector<Number>::operator=(const Vector<Number2> &in_vector)
113   {
114     this->thread_loop_partitioner = in_vector.thread_loop_partitioner;
115     if (this->size() != in_vector.size())
116       ReadWriteVector<Number>::reinit(in_vector, true);
117 
118     dealii::internal::VectorOperations::Vector_copy<Number, Number2> copier(
119       in_vector.values.get(), this->values.get());
120     dealii::internal::VectorOperations::parallel_for(
121       copier,
122       static_cast<size_type>(0),
123       this->size(),
124       this->thread_loop_partitioner);
125 
126     return *this;
127   }
128 
129 
130 
131   template <typename Number>
132   Vector<Number> &
133   Vector<Number>::operator=(const Number s)
134   {
135     Assert(s == static_cast<Number>(0),
136            ExcMessage("Only 0 can be assigned to a vector."));
137     (void)s;
138 
139     dealii::internal::VectorOperations::Vector_set<Number> setter(
140       Number(), this->values.get());
141     dealii::internal::VectorOperations::parallel_for(
142       setter, 0, this->size(), this->thread_loop_partitioner);
143 
144     return *this;
145   }
146 
147 
148 
149   template <typename Number>
150   Vector<Number> &
151   Vector<Number>::operator*=(const Number factor)
152   {
153     AssertIsFinite(factor);
154 
155     dealii::internal::VectorOperations::Vectorization_multiply_factor<Number>
156       vector_multiply(this->values.get(), factor);
157     dealii::internal::VectorOperations::parallel_for(
158       vector_multiply,
159       static_cast<size_type>(0),
160       this->size(),
161       this->thread_loop_partitioner);
162 
163     return *this;
164   }
165 
166 
167 
168   template <typename Number>
169   Vector<Number> &
170   Vector<Number>::operator/=(const Number factor)
171   {
172     AssertIsFinite(factor);
173     Assert(factor != Number(0.), ExcZero());
174     this->operator*=(Number(1.) / factor);
175 
176     return *this;
177   }
178 
179 
180 
181   template <typename Number>
182   Vector<Number> &
183   Vector<Number>::operator+=(const VectorSpaceVector<Number> &V)
184   {
185     // Check that casting will work.
186     Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
187            ExcVectorTypeNotCompatible());
188 
189     // Downcast V. If fails, throws an exception.
190     const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
191     Assert(down_V.size() == this->size(),
192            ExcMessage(
193              "Cannot add two vectors with different numbers of elements"));
194 
195     dealii::internal::VectorOperations::Vectorization_add_v<Number> vector_add(
196       this->values.get(), down_V.values.get());
197     dealii::internal::VectorOperations::parallel_for(
198       vector_add, 0, this->size(), this->thread_loop_partitioner);
199 
200     return *this;
201   }
202 
203 
204 
205   template <typename Number>
206   Vector<Number> &
207   Vector<Number>::operator-=(const VectorSpaceVector<Number> &V)
208   {
209     // Check that casting will work.
210     Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
211            ExcVectorTypeNotCompatible());
212 
213     // Downcast V. If fails, throws an exception.
214     const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
215     Assert(down_V.size() == this->size(),
216            ExcMessage(
217              "Cannot subtract two vectors with different numbers of elements"));
218     dealii::internal::VectorOperations::Vectorization_subtract_v<Number>
219       vector_subtract(this->values.get(), down_V.values.get());
220     dealii::internal::VectorOperations::parallel_for(
221       vector_subtract, 0, this->size(), this->thread_loop_partitioner);
222 
223     return *this;
224   }
225 
226 
227 
228   template <typename Number>
229   Number Vector<Number>::operator*(const VectorSpaceVector<Number> &V) const
230   {
231     // Check that casting will work.
232     Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
233            ExcVectorTypeNotCompatible());
234 
235     // Downcast V. If fails, throws an exception.
236     const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
237     Assert(down_V.size() == this->size(),
238            ExcMessage("Cannot compute the scalar product "
239                       "of two vectors with different numbers of elements"));
240     Number                                                  sum;
241     dealii::internal::VectorOperations::Dot<Number, Number> dot(
242       this->values.get(), down_V.values.get());
243     dealii::internal::VectorOperations::parallel_reduce(
244       dot, 0, this->size(), sum, this->thread_loop_partitioner);
245 
246     return sum;
247   }
248 
249 
250 
251   template <typename Number>
252   void
import(const ReadWriteVector<Number> &,VectorOperation::values,std::shared_ptr<const CommunicationPatternBase>)253   Vector<Number>::import(const ReadWriteVector<Number> &,
254                          VectorOperation::values,
255                          std::shared_ptr<const CommunicationPatternBase>)
256   {
257     AssertThrow(false, ExcMessage("This function is not implemented."));
258   }
259 
260 
261 
262   template <typename Number>
263   inline void
add(const Number a)264   Vector<Number>::add(const Number a)
265   {
266     AssertIsFinite(a);
267 
268     dealii::internal::VectorOperations::Vectorization_add_factor<Number>
269       vector_add(this->values.get(), a);
270     dealii::internal::VectorOperations::parallel_for(
271       vector_add, 0, this->size(), this->thread_loop_partitioner);
272   }
273 
274 
275 
276   template <typename Number>
277   void
add(const Number a,const VectorSpaceVector<Number> & V)278   Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V)
279   {
280     // Check that casting will work.
281     Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
282            ExcVectorTypeNotCompatible());
283 
284     // Downcast V. If fails, throws an exception.
285     const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
286     AssertIsFinite(a);
287     Assert(down_V.size() == this->size(),
288            ExcMessage(
289              "Cannot add two vectors with different numbers of elements"));
290 
291     dealii::internal::VectorOperations::Vectorization_add_av<Number>
292       vector_add_av(this->values.get(), down_V.values.get(), a);
293     dealii::internal::VectorOperations::parallel_for(
294       vector_add_av, 0, this->size(), this->thread_loop_partitioner);
295   }
296 
297 
298 
299   template <typename Number>
300   void
add(const Number a,const VectorSpaceVector<Number> & V,const Number b,const VectorSpaceVector<Number> & W)301   Vector<Number>::add(const Number                     a,
302                       const VectorSpaceVector<Number> &V,
303                       const Number                     b,
304                       const VectorSpaceVector<Number> &W)
305   {
306     // Check that casting will work.
307     Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
308            ExcVectorTypeNotCompatible());
309     // Check that casting will work.
310     Assert(dynamic_cast<const Vector<Number> *>(&W) != nullptr,
311            ExcVectorTypeNotCompatible());
312 
313     // Downcast V. If fails, throws an exception.
314     const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
315     // Downcast W. If fails, throws an exception.
316     const Vector<Number> &down_W = dynamic_cast<const Vector<Number> &>(W);
317     AssertIsFinite(a);
318     Assert(down_V.size() == this->size(),
319            ExcMessage(
320              "Cannot add two vectors with different numbers of elements"));
321     AssertIsFinite(b);
322     Assert(down_W.size() == this->size(),
323            ExcMessage(
324              "Cannot add two vectors with different numbers of elements"));
325 
326     dealii::internal::VectorOperations::Vectorization_add_avpbw<Number>
327       vector_add(
328         this->values.get(), down_V.values.get(), down_W.values.get(), a, b);
329     dealii::internal::VectorOperations::parallel_for(
330       vector_add, 0, this->size(), this->thread_loop_partitioner);
331   }
332 
333 
334 
335   template <typename Number>
336   void
sadd(const Number s,const Number a,const VectorSpaceVector<Number> & V)337   Vector<Number>::sadd(const Number                     s,
338                        const Number                     a,
339                        const VectorSpaceVector<Number> &V)
340   {
341     AssertIsFinite(s);
342     AssertIsFinite(a);
343 
344     // Check that casting will work.
345     Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
346            ExcVectorTypeNotCompatible());
347 
348     // Downcast V. It fails, throws an exception.
349     const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
350     dealii::internal::VectorOperations::Vectorization_sadd_xav<Number>
351       vector_sadd_xav(this->values.get(), down_V.values.get(), a, s);
352     dealii::internal::VectorOperations::parallel_for(
353       vector_sadd_xav, 0, this->size(), this->thread_loop_partitioner);
354   }
355 
356 
357 
358   template <typename Number>
359   void
scale(const VectorSpaceVector<Number> & scaling_factors)360   Vector<Number>::scale(const VectorSpaceVector<Number> &scaling_factors)
361   {
362     // Check that casting will work.
363     Assert(dynamic_cast<const Vector<Number> *>(&scaling_factors) != nullptr,
364            ExcVectorTypeNotCompatible());
365 
366     // Downcast scaling_factors. If fails, throws an exception.
367     const Vector<Number> &down_scaling_factors =
368       dynamic_cast<const Vector<Number> &>(scaling_factors);
369     Assert(down_scaling_factors.size() == this->size(),
370            ExcMessage(
371              "Cannot add two vectors with different numbers of elements"));
372 
373     dealii::internal::VectorOperations::Vectorization_scale<Number>
374       vector_scale(this->values.get(), down_scaling_factors.values.get());
375     dealii::internal::VectorOperations::parallel_for(
376       vector_scale, 0, this->size(), this->thread_loop_partitioner);
377   }
378 
379 
380 
381   template <typename Number>
382   void
equ(const Number a,const VectorSpaceVector<Number> & V)383   Vector<Number>::equ(const Number a, const VectorSpaceVector<Number> &V)
384   {
385     AssertIsFinite(a);
386 
387     // Check that casting will work.
388     Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
389            ExcVectorTypeNotCompatible());
390 
391     // Downcast V. If fails, throws an exception.
392     const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
393     dealii::internal::VectorOperations::Vectorization_equ_au<Number> vector_equ(
394       this->values.get(), down_V.values.get(), a);
395     dealii::internal::VectorOperations::parallel_for(
396       vector_equ, 0, this->size(), this->thread_loop_partitioner);
397   }
398 
399 
400 
401   template <typename Number>
402   bool
all_zero()403   Vector<Number>::all_zero() const
404   {
405     Assert(this->size(), ExcEmptyObject());
406 
407     const size_type size = this->size();
408     for (size_type i = 0; i < size; ++i)
409       if (this->values[i] != Number())
410         return false;
411 
412     return true;
413   }
414 
415 
416 
417   template <typename Number>
418   typename Vector<Number>::value_type
mean_value()419   Vector<Number>::mean_value() const
420   {
421     Assert(this->size(), ExcEmptyObject());
422 
423     using real_type = typename VectorSpaceVector<Number>::real_type;
424     value_type                                            sum;
425     dealii::internal::VectorOperations::MeanValue<Number> mean_value(
426       this->values.get());
427     dealii::internal::VectorOperations::parallel_reduce(
428       mean_value, 0, this->size(), sum, this->thread_loop_partitioner);
429 
430     return sum / static_cast<real_type>(this->size());
431   }
432 
433 
434 
435   template <typename Number>
436   typename VectorSpaceVector<Number>::real_type
l1_norm()437   Vector<Number>::l1_norm() const
438   {
439     Assert(this->size(), ExcEmptyObject());
440 
441     using real_type = typename VectorSpaceVector<Number>::real_type;
442     real_type                                                    sum;
443     dealii::internal::VectorOperations::Norm1<Number, real_type> norm1(
444       this->values.get());
445     dealii::internal::VectorOperations::parallel_reduce(
446       norm1, 0, this->size(), sum, this->thread_loop_partitioner);
447 
448     return sum;
449   }
450 
451 
452 
453   template <typename Number>
454   typename VectorSpaceVector<Number>::real_type
l2_norm()455   Vector<Number>::l2_norm() const
456   {
457     Assert(this->size(), ExcEmptyObject());
458 
459     // if l2_norm()^2 is finite and non-zero, the answer is computed as
460     // std::sqrt(norm_sqr()). If norm_sqrt() is infinite or zero, the l2 norm
461     // might still be finite. In that case, recompute ig (this is a rare case,
462     // so working on the vector twice is uncritical and paid off by the extended
463     // precision) using the BLAS approach with a weight, see e.g. dnrm2.f.
464     using real_type = typename VectorSpaceVector<Number>::real_type;
465     real_type                                                    norm_square;
466     dealii::internal::VectorOperations::Norm2<Number, real_type> norm2(
467       this->values.get());
468     dealii::internal::VectorOperations::parallel_reduce(
469       norm2, 0, this->size(), norm_square, this->thread_loop_partitioner);
470     if (numbers::is_finite(norm_square) &&
471         norm_square >= std::numeric_limits<real_type>::min())
472       return std::sqrt(norm_square);
473     else
474       {
475         real_type       scale = 0.;
476         real_type       sum   = 1.;
477         const size_type size  = this->size();
478         for (size_type i = 0; i < size; ++i)
479           {
480             if (this->values[i] != Number())
481               {
482                 const real_type abs_x =
483                   numbers::NumberTraits<Number>::abs(this->values[i]);
484                 if (scale < abs_x)
485                   {
486                     sum   = 1. + sum * (scale / abs_x) * (scale / abs_x);
487                     scale = abs_x;
488                   }
489                 else
490                   sum += (abs_x / scale) * (abs_x / scale);
491               }
492           }
493         AssertIsFinite(scale * std::sqrt(sum));
494         return scale * std::sqrt(sum);
495       }
496   }
497 
498 
499 
500   template <typename Number>
501   typename VectorSpaceVector<Number>::real_type
linfty_norm()502   Vector<Number>::linfty_norm() const
503   {
504     typename ReadWriteVector<Number>::real_type norm = 0.;
505     const size_type                             size = this->size();
506     for (size_type i = 0; i < size; ++i)
507       norm = std::max(std::abs(this->values[i]), norm);
508 
509     return norm;
510   }
511 
512 
513 
514   template <typename Number>
515   Number
add_and_dot(const Number a,const VectorSpaceVector<Number> & V,const VectorSpaceVector<Number> & W)516   Vector<Number>::add_and_dot(const Number                     a,
517                               const VectorSpaceVector<Number> &V,
518                               const VectorSpaceVector<Number> &W)
519   {
520     // Check that casting will work.
521     Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
522            ExcVectorTypeNotCompatible());
523     // Check that casting will work.
524     Assert(dynamic_cast<const Vector<Number> *>(&W) != nullptr,
525            ExcVectorTypeNotCompatible());
526 
527     // Downcast V. If fails, throws an exception.
528     const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
529     // Downcast W. If fails, throws an exception.
530     const Vector<Number> &down_W = dynamic_cast<const Vector<Number> &>(W);
531     AssertIsFinite(a);
532     Assert(down_V.size() == this->size(),
533            ExcMessage(
534              "Cannot add two vectors with different numbers of elements"));
535     Assert(down_W.size() == this->size(),
536            ExcMessage(
537              "Cannot add two vectors with different numbers of elements"));
538 
539     Number                                                sum;
540     dealii::internal::VectorOperations::AddAndDot<Number> adder(
541       this->values.get(), down_V.values.get(), down_W.values.get(), a);
542     dealii::internal::VectorOperations::parallel_reduce(
543       adder, 0, this->size(), sum, this->thread_loop_partitioner);
544     AssertIsFinite(sum);
545 
546     return sum;
547   }
548 
549 
550 
551   template <typename Number>
552   void
print_as_numpy_array(std::ostream & out,const unsigned int precision)553   Vector<Number>::print_as_numpy_array(std::ostream &     out,
554                                        const unsigned int precision) const
555   {
556     AssertThrow(out, ExcIO());
557     boost::io::ios_flags_saver restore_flags(out);
558 
559     out.precision(precision);
560 
561     const unsigned int n_elements = this->n_elements();
562     for (unsigned int i = 0; i < n_elements; ++i)
563       out << this->values[i] << ' ';
564     out << '\n' << std::flush;
565 
566     AssertThrow(out, ExcIO());
567   }
568 
569 
570 
571   template <typename Number>
572   void
block_write(std::ostream & out)573   Vector<Number>::block_write(std::ostream &out) const
574   {
575     AssertThrow(out, ExcIO());
576 
577     // Other version of the following
578     //  out << size() << std::endl << '[';
579     // Reason: operator<< seems to use some resources that lead to problems in
580     // a multithreaded environment.  We convert the size index to
581     // unsigned long long int that is at least 64 bits to be able to output it
582     // on all platforms, since std::uint64_t is not in C.
583     const unsigned long long int sz = this->size();
584     char                         buf[16];
585 
586     std::sprintf(buf, "%llu", sz);
587     std::strcat(buf, "\n[");
588 
589     out.write(buf, std::strlen(buf));
590     out.write(reinterpret_cast<const char *>(this->begin()),
591               reinterpret_cast<const char *>(this->end()) -
592                 reinterpret_cast<const char *>(this->begin()));
593 
594     // out << ']';
595     const char outro = ']';
596     out.write(&outro, 1);
597 
598     AssertThrow(out, ExcIO());
599   }
600 
601 
602 
603   template <typename Number>
604   void
block_read(std::istream & in)605   Vector<Number>::block_read(std::istream &in)
606   {
607     AssertThrow(in, ExcIO());
608 
609     size_type sz;
610 
611     char buf[16];
612 
613     in.getline(buf, 16, '\n');
614     sz = std::atoi(buf);
615     ReadWriteVector<Number>::reinit(sz, true);
616 
617     char c;
618     // in >> c;
619     in.read(&c, 1);
620     AssertThrow(c == '[', ExcIO());
621 
622     in.read(reinterpret_cast<char *>(this->begin()),
623             reinterpret_cast<const char *>(this->end()) -
624               reinterpret_cast<const char *>(this->begin()));
625 
626     // in >> c;
627     in.read(&c, 1);
628     AssertThrow(c == ']', ExcIO());
629   }
630 } // namespace LinearAlgebra
631 
632 DEAL_II_NAMESPACE_CLOSE
633 
634 #endif
635