Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
vector.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: vector.templates.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__vector_templates_h
18 #define __deal2__vector_templates_h
19 
20 
21 #include <deal.II/base/template_constraints.h>
22 #include <deal.II/base/numbers.h>
23 #include <deal.II/base/parallel.h>
24 #include <deal.II/base/thread_management.h>
25 #include <deal.II/lac/vector.h>
26 #include <deal.II/lac/block_vector.h>
27 
28 #ifdef DEAL_II_WITH_PETSC
29 # include <deal.II/lac/petsc_vector.h>
30 # include <deal.II/lac/petsc_parallel_vector.h>
31 #endif
32 
33 #ifdef DEAL_II_WITH_TRILINOS
34 # include <deal.II/lac/trilinos_vector.h>
35 #endif
36 
37 #include <boost/lambda/lambda.hpp>
38 #include <cmath>
39 #include <cstring>
40 #include <algorithm>
41 #include <iostream>
42 #include <iomanip>
43 
45 
46 
47 namespace internal
48 {
49  template <typename T>
50  bool is_non_negative (const T &t)
51  {
52  return t >= 0;
53  }
54 
55 
56  template <typename T>
57  bool is_non_negative (const std::complex<T> &)
58  {
59  Assert (false,
60  ExcMessage ("Complex numbers do not have an ordering."));
61 
62  return false;
63  }
64 
65 
66  template <typename T>
67  void print (const T &t,
68  const char *format)
69  {
70  if (format != 0)
71  std::printf (format, t);
72  else
73  std::printf (" %5.2f", double(t));
74  }
75 
76 
77 
78  template <typename T>
79  void print (const std::complex<T> &t,
80  const char *format)
81  {
82  if (format != 0)
83  std::printf (format, t.real(), t.imag());
84  else
85  std::printf (" %5.2f+%5.2fi",
86  double(t.real()), double(t.imag()));
87  }
88 
89  // call std::copy, except for in
90  // the case where we want to copy
91  // from std::complex to a
92  // non-complex type
93  template <typename T, typename U>
94  void copy (const T *begin,
95  const T *end,
96  U *dest)
97  {
98  std::copy (begin, end, dest);
99  }
100 
101  template <typename T, typename U>
102  void copy (const std::complex<T> *begin,
103  const std::complex<T> *end,
104  std::complex<U> *dest)
105  {
106  std::copy (begin, end, dest);
107  }
108 
109  template <typename T, typename U>
110  void copy (const std::complex<T> *,
111  const std::complex<T> *,
112  U *)
113  {
114  Assert (false, ExcMessage ("Can't convert a vector of complex numbers "
115  "into a vector of reals/doubles"));
116  }
117 }
118 
119 
120 
121 
122 template <typename Number>
124  :
125  Subscriptor(),
126  vec_size(v.size()),
127  max_vec_size(v.size()),
128  val(0)
129 {
130  if (vec_size != 0)
131  {
132  val = new Number[max_vec_size];
133  Assert (val != 0, ExcOutOfMemory());
134  *this = v;
135  }
136 }
137 
138 
139 #ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG
140 
141 template <typename Number>
142 template <typename OtherNumber>
144  :
145  Subscriptor(),
146  vec_size(v.size()),
147  max_vec_size(v.size()),
148  val(0)
149 {
150  if (vec_size != 0)
151  {
152  val = new Number[max_vec_size];
153  Assert (val != 0, ExcOutOfMemory());
154  std::copy (v.begin(), v.end(), begin());
155  }
156 }
157 
158 #endif
159 
160 #ifdef DEAL_II_WITH_PETSC
161 
162 
163 template <typename Number>
165  :
166  Subscriptor(),
167  vec_size(v.size()),
168  max_vec_size(v.size()),
169  val(0)
170 {
171  if (vec_size != 0)
172  {
173  val = new Number[max_vec_size];
174  Assert (val != 0, ExcOutOfMemory());
175 
176  // get a representation of the vector
177  // and copy it
178  PetscScalar *start_ptr;
179  int ierr = VecGetArray (static_cast<const Vec &>(v), &start_ptr);
180  AssertThrow (ierr == 0, ExcPETScError(ierr));
181 
182  internal::copy (start_ptr, start_ptr+vec_size, begin());
183 
184  // restore the representation of the
185  // vector
186  ierr = VecRestoreArray (static_cast<const Vec &>(v), &start_ptr);
187  AssertThrow (ierr == 0, ExcPETScError(ierr));
188  }
189 }
190 
191 
192 
193 template <typename Number>
195  :
196  Subscriptor(),
197  vec_size(0),
198  max_vec_size(0),
199  val(0)
200 {
201  if (v.size() != 0)
202  {
203  // do this in a two-stage process:
204  // first convert to a sequential petsc
205  // vector, then copy that
206  PETScWrappers::Vector seq (v);
207  *this = seq;
208  }
209 }
210 
211 #endif
212 
213 #ifdef DEAL_II_WITH_TRILINOS
214 
215 template <typename Number>
217  :
218  Subscriptor(),
219  vec_size(v.size()),
220  max_vec_size(v.size()),
221  val(0)
222 {
223  if (vec_size != 0)
224  {
225  val = new Number[max_vec_size];
226  Assert (val != 0, ExcOutOfMemory());
227 
228  // Copy the distributed vector to
229  // a local one at all
230  // processors. TODO: There could
231  // be a better solution than
232  // this, but it has not yet been
233  // found.
234  TrilinosWrappers::Vector localized_vector (v);
235 
236  // get a representation of the vector
237  // and copy it
238  TrilinosScalar **start_ptr;
239 
240  int ierr = localized_vector.trilinos_vector().ExtractView (&start_ptr);
241  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
242 
243  std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());
244  }
245 }
246 
247 
248 
249 template <typename Number>
251  :
252  Subscriptor(),
253  vec_size(v.size()),
254  max_vec_size(v.size()),
255  val(0)
256 {
257  if (vec_size != 0)
258  {
259  val = new Number[max_vec_size];
260  Assert (val != 0, ExcOutOfMemory());
261 
262  // get a representation of the vector
263  // and copy it
264  TrilinosScalar **start_ptr;
265 
266  int ierr = v.trilinos_vector().ExtractView (&start_ptr);
267  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
268 
269  std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());
270  }
271 }
272 
273 #endif
274 
275 template <typename Number>
276 template <typename Number2>
278  const bool fast)
279 {
280  reinit (v.size(), fast);
281 }
282 
283 // Moved to vector.h as an inline function by Luca Heltai on
284 // 2009/04/12 to prevent strange compiling errors, after making swap
285 // virtual.
286 // template <typename Number>
287 // void
288 // Vector<Number>::swap (Vector<Number> &v)
289 // {
290 // std::swap (vec_size, v.vec_size);
291 // std::swap (max_vec_size, v.max_vec_size);
292 // std::swap (val, v.val);
293 // }
294 
295 
296 
297 template <typename Number>
298 bool
300 {
302 
303  for (size_type i=0; i<vec_size; ++i)
304  if (val[i] != Number(0))
305  return false;
306  return true;
307 }
308 
309 
310 
311 template <typename Number>
312 bool
314 {
316 
317  for (size_type i=0; i<vec_size; ++i)
318  if ( ! internal::is_non_negative (val[i]))
319  return false;
320 
321  return true;
322 }
323 
324 
325 
326 namespace internal
327 {
328  namespace Vector
329  {
330  template <typename T>
331  void set_subrange (const T s,
332  const typename ::Vector<T>::size_type begin,
333  const typename ::Vector<T>::size_type end,
334  ::Vector<T> &dst)
335  {
336  if (s == T())
337  std::memset ((dst.begin()+begin),0,(end-begin)*sizeof(T));
338  else
339  std::fill (&*(dst.begin()+begin), &*(dst.begin()+end), s);
340  }
341 
342 
343  template <typename T>
344  void copy_subrange (const typename ::Vector<T>::size_type begin,
345  const typename ::Vector<T>::size_type end,
346  const ::Vector<T> &src,
347  ::Vector<T> &dst)
348  {
349  memcpy(&*(dst.begin()+begin), &*(src.begin()+begin),
350  (end-begin)*sizeof(T));
351  }
352 
353 
354  template <typename T, typename U>
355  void copy_subrange (const typename ::Vector<T>::size_type begin,
356  const typename ::Vector<T>::size_type end,
357  const ::Vector<T> &src,
358  ::Vector<U> &dst)
359  {
360  const T *q = src.begin()+begin;
361  const T *const end_q = src.begin()+end;
362  U *p = dst.begin()+begin;
363  for (; q!=end_q; ++q, ++p)
364  *p = *q;
365  }
366 
367 
368  template<typename T, typename U>
369  void copy_subrange_wrap (const typename ::Vector<T>::size_type begin,
370  const typename ::Vector<T>::size_type end,
371  const ::Vector<T> &src,
372  ::Vector<U> &dst)
373  {
374  copy_subrange (begin, end, src, dst);
375  }
376 
377 
378  template <typename T, typename U>
379  void copy_vector (const ::Vector<T> &src,
380  ::Vector<U> &dst)
381  {
382  if (PointerComparison::equal(&src, &dst))
383  return;
384 
385  const typename ::Vector<T>::size_type vec_size = src.size();
386  const typename ::Vector<U>::size_type dst_size = dst.size();
387  if (dst_size != vec_size)
388  dst.reinit (vec_size, true);
389  if (vec_size>internal::Vector::minimum_parallel_grain_size)
390  parallel::apply_to_subranges (0U, vec_size,
391  std_cxx1x::bind(&internal::Vector::template
392  copy_subrange_wrap <T,U>,
393  std_cxx1x::_1,
394  std_cxx1x::_2,
395  std_cxx1x::cref(src),
396  std_cxx1x::ref(dst)),
397  internal::Vector::minimum_parallel_grain_size);
398  else if (vec_size > 0)
399  copy_subrange (0U, vec_size, src, dst);
400  }
401  }
402 }
403 
404 
405 
406 template <typename Number>
409 {
411  if (s != Number())
412  Assert (vec_size!=0, ExcEmptyObject());
413  if (vec_size>internal::Vector::minimum_parallel_grain_size)
414  parallel::apply_to_subranges (0U, vec_size,
415  std_cxx1x::bind(&internal::Vector::template
416  set_subrange<Number>,
417  s, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::ref(*this)),
418  internal::Vector::minimum_parallel_grain_size);
419  else if (vec_size > 0)
420  internal::Vector::set_subrange<Number>(s, 0U, vec_size, *this);
421 
422  return *this;
423 }
424 
425 
426 
427 #ifdef DEAL_II_BOOST_BIND_COMPILER_BUG
428 template <>
430 Vector<std::complex<float> >::operator = (const std::complex<float> s)
431 {
433  if (s != std::complex<float>())
434  Assert (vec_size!=0, ExcEmptyObject());
435  if (vec_size!=0)
436  std::fill (begin(), end(), s);
437 
438  return *this;
439 }
440 #endif
441 
442 
443 
444 template <typename Number>
446 {
448 
449  Assert (vec_size!=0, ExcEmptyObject());
450 
452  val+vec_size,
453  val,
454  (factor*boost::lambda::_1),
455  internal::Vector::minimum_parallel_grain_size);
456 
457  return *this;
458 }
459 
460 
461 
462 template <typename Number>
463 void
464 Vector<Number>::add (const Number a,
465  const Vector<Number> &v)
466 {
468 
469  Assert (vec_size!=0, ExcEmptyObject());
470  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
471 
473  val+vec_size,
474  v.val,
475  val,
476  (boost::lambda::_1 + a*boost::lambda::_2),
477  internal::Vector::minimum_parallel_grain_size);
478 }
479 
480 
481 
482 template <typename Number>
483 void
484 Vector<Number>::sadd (const Number x,
485  const Number a,
486  const Vector<Number> &v)
487 {
490 
491  Assert (vec_size!=0, ExcEmptyObject());
492  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
493 
495  val+vec_size,
496  v.val,
497  val,
498  (x*boost::lambda::_1 + a*boost::lambda::_2),
499  internal::Vector::minimum_parallel_grain_size);
500 }
501 
502 
503 
504 namespace internal
505 {
506  namespace Vector
507  {
508  // All sums over all the vector entries (l2-norm, inner product, etc.) are
509  // performed with the same code, using a templated operation defined here
510  template <typename Number, typename Number2>
511  struct InnerProd
512  {
513  Number
514  operator() (const Number *&X, const Number2 *&Y, const Number &) const
515  {
516  return *X++ * Number(numbers::NumberTraits<Number2>::conjugate(*Y++));
517  }
518  };
519 
520  template <typename Number, typename RealType>
521  struct Norm2
522  {
523  RealType
524  operator() (const Number *&X, const Number *&, const RealType &) const
525  {
527  }
528  };
529 
530  template <typename Number, typename RealType>
531  struct Norm1
532  {
533  RealType
534  operator() (const Number *&X, const Number *&, const RealType &) const
535  {
537  }
538  };
539 
540  template <typename Number, typename RealType>
541  struct NormP
542  {
543  RealType
544  operator() (const Number *&X, const Number *&, const RealType &p) const
545  {
546  return std::pow(numbers::NumberTraits<Number>::abs(*X++), p);
547  }
548  };
549 
550  template <typename Number>
551  struct MeanValue
552  {
553  Number
554  operator() (const Number *&X, const Number *&, const Number &) const
555  {
556  return *X++;
557  }
558  };
559 
560  // this is the main working loop for all vector sums using the templated
561  // operation above. it accumulates the sums using a block-wise summation
562  // algorithm with post-update. this blocked algorithm has been proposed in
563  // a similar form by Castaldo, Whaley and Chronopoulos (SIAM
564  // J. Sci. Comput. 31, 1156-1174, 2008) and we use the smallest possible
565  // block size, 2. Sometimes it is referred to as pairwise summation. The
566  // worst case error made by this algorithm is on the order O(eps *
567  // log2(vec_size)), whereas a naive summation is O(eps * vec_size). Even
568  // though the Kahan summation is even more accurate with an error O(eps)
569  // by carrying along remainders not captured by the main sum, that involves
570  // additional costs which are not worthwhile. See the Wikipedia article on
571  // the Kahan summation algorithm.
572 
573  // The algorithm implemented here has the additional benefit that it is
574  // easily parallelized without changing the order of how the elements are
575  // added (floating point addition is not associative). For the same vector
576  // size and minimum_parallel_grainsize, the blocks are always the
577  // same and added pairwise. At the innermost level, eight values are added
578  // consecutively in order to better balance multiplications and additions.
579 
580  // The code returns the result as the last argument in order to make
581  // spawning tasks simpler and use automatic template deduction.
582  template <typename Operation, typename Number, typename Number2,
583  typename ResultType, typename size_type>
584  void accumulate (const Operation &op,
585  const Number *X,
586  const Number2 *Y,
587  const ResultType power,
588  const size_type vec_size,
589  ResultType &result)
590  {
591  if (vec_size <= 4096)
592  {
593  // the vector is short enough so we perform the summation. first
594  // work on the regular part. The innermost 32 values are expanded in
595  // order to obtain known loop bounds for most of the work.
596  const Number *X_original = X;
597  ResultType outer_results [128];
598  size_type n_chunks = vec_size / 32;
599  const size_type remainder = vec_size % 32;
600  Assert (remainder == 0 || n_chunks < 128, ExcInternalError());
601 
602  for (size_type i=0; i<n_chunks; ++i)
603  {
604  ResultType r0 = op(X, Y, power);
605  for (size_type j=1; j<8; ++j)
606  r0 += op(X, Y, power);
607  ResultType r1 = op(X, Y, power);
608  for (size_type j=1; j<8; ++j)
609  r1 += op(X, Y, power);
610  r0 += r1;
611  r1 = op(X, Y, power);
612  for (size_type j=1; j<8; ++j)
613  r1 += op(X, Y, power);
614  ResultType r2 = op(X, Y, power);
615  for (size_type j=1; j<8; ++j)
616  r2 += op(X, Y, power);
617  r1 += r2;
618  r0 += r1;
619  outer_results[i] = r0;
620  }
621 
622  // now work on the remainder, i.e., the last
623  // up to 32 values. Use switch statement with
624  // fall-through to work on these values.
625  if (remainder > 0)
626  {
627  const size_type inner_chunks = remainder / 8;
628  Assert (inner_chunks <= 3, ExcInternalError());
629  const size_type remainder_inner = remainder % 8;
630  ResultType r0 = ResultType(), r1 = ResultType(),
631  r2 = ResultType();
632  switch (inner_chunks)
633  {
634  case 3:
635  r2 = op(X, Y, power);
636  for (size_type j=1; j<8; ++j)
637  r2 += op(X, Y, power);
638  // no break
639  case 2:
640  r1 = op(X, Y, power);
641  for (size_type j=1; j<8; ++j)
642  r1 += op(X, Y, power);
643  r1 += r2;
644  // no break
645  case 1:
646  r2 = op(X, Y, power);
647  for (size_type j=1; j<8; ++j)
648  r2 += op(X, Y, power);
649  // no break
650  default:
651  for (size_type j=0; j<remainder_inner; ++j)
652  r0 += op(X, Y, power);
653  r0 += r2;
654  r0 += r1;
655  outer_results[n_chunks] = r0;
656  break;
657  }
658  n_chunks++;
659  }
660  AssertDimension(static_cast<size_type> (X - X_original), vec_size);
661 
662  // now sum the results from the chunks
663  // recursively
664  while (n_chunks > 1)
665  {
666  if (n_chunks % 2 == 1)
667  outer_results[n_chunks++] = ResultType();
668  for (size_type i=0; i<n_chunks; i+=2)
669  outer_results[i/2] = outer_results[i] + outer_results[i+1];
670  n_chunks /= 2;
671  }
672  result = outer_results[0];
673  }
674 #ifdef DEAL_II_WITH_THREADS
675  else if (vec_size > 4 * internal::Vector::minimum_parallel_grain_size)
676  {
677  // split the vector into smaller pieces to be
678  // worked on recursively and create tasks for
679  // them. Make pieces divisible by 1024.
680  const size_type new_size = (vec_size / 4096) * 1024;
681  ResultType r0, r1, r2, r3;
682  Threads::TaskGroup<> task_group;
683  task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
684  ResultType,size_type>,
685  op, X, Y, power, new_size, r0);
686  task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
687  ResultType,size_type>,
688  op, X+new_size, Y+new_size, power,
689  new_size, r1);
690  task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
691  ResultType,size_type>,
692  op, X+2*new_size, Y+2*new_size, power,
693  new_size, r2);
694  task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
695  ResultType,size_type>,
696  op, X+3*new_size, Y+3*new_size, power,
697  vec_size-3*new_size, r3);
698  task_group.join_all();
699  r0 += r1;
700  r2 += r3;
701  result = r0 + r2;
702  }
703 #endif
704  else
705  {
706  // split vector into four pieces and work on
707  // the pieces recursively. Make pieces (except last)
708  // divisible by 1024.
709  const size_type new_size = (vec_size / 4096) * 1024;
710  ResultType r0, r1, r2, r3;
711  accumulate (op, X, Y, power, new_size, r0);
712  accumulate (op, X+new_size, Y+new_size, power, new_size, r1);
713  accumulate (op, X+2*new_size, Y+2*new_size, power, new_size, r2);
714  accumulate (op, X+3*new_size, Y+3*new_size, power, vec_size-3*new_size, r3);
715  r0 += r1;
716  r2 += r3;
717  result = r0 + r2;
718  }
719  }
720  }
721 }
722 
723 
724 
725 template <typename Number>
726 template <typename Number2>
728 {
729  Assert (vec_size!=0, ExcEmptyObject());
730 
731  if (PointerComparison::equal (this, &v))
732  return norm_sqr();
733 
734  Assert (vec_size == v.size(),
735  ExcDimensionMismatch(vec_size, v.size()));
736 
737  Number sum;
738  internal::Vector::accumulate (internal::Vector::InnerProd<Number,Number2>(),
739  val, v.val, Number(), vec_size, sum);
741 
742  return sum;
743 }
744 
745 
746 template <typename Number>
749 {
750  Assert (vec_size!=0, ExcEmptyObject());
751 
752  real_type sum;
753  internal::Vector::accumulate (internal::Vector::Norm2<Number,real_type>(),
754  val, val, real_type(), vec_size, sum);
755 
757 
758  return sum;
759 }
760 
761 
762 template <typename Number>
764 {
765  Assert (vec_size!=0, ExcEmptyObject());
766 
767  Number sum;
768  internal::Vector::accumulate (internal::Vector::MeanValue<Number>(),
769  val, val, Number(), vec_size, sum);
770 
771  return sum / real_type(size());
772 }
773 
774 
775 
776 template <typename Number>
779 {
780  Assert (vec_size!=0, ExcEmptyObject());
781 
782  real_type sum;
783  internal::Vector::accumulate (internal::Vector::Norm1<Number,real_type>(),
784  val, val, real_type(), vec_size, sum);
785 
786  return sum;
787 }
788 
789 
790 
791 template <typename Number>
794 {
795  // if l2_norm()^2 is finite and non-zero, the answer is computed as
796  // std::sqrt(norm_sqr()). If norm_sqr() is infinite or zero, the l2 norm
797  // might still be finite. In that case, recompute it (this is a rare case,
798  // so working on the vector twice is uncritical and paid off by the extended
799  // precision) using the BLAS approach with a weight, see e.g. dnrm2.f.
800  Assert (vec_size!=0, ExcEmptyObject());
801 
802  real_type norm_square;
803  internal::Vector::accumulate (internal::Vector::Norm2<Number,real_type>(),
804  val, val, real_type(), vec_size, norm_square);
805  if (numbers::is_finite(norm_square) &&
806  norm_square >= std::numeric_limits<real_type>::min())
807  return std::sqrt(norm_square);
808  else
809  {
810  real_type scale = 0.;
811  real_type sum = 1.;
812  for (size_type i=0; i<vec_size; ++i)
813  {
814  if (val[i] != Number())
815  {
816  const real_type abs_x =
818  if (scale < abs_x)
819  {
820  sum = 1. + sum * (scale/abs_x) * (scale/abs_x);
821  scale = abs_x;
822  }
823  else
824  sum += (abs_x/scale) * (abs_x/scale);
825  }
826  }
827  Assert(numbers::is_finite(scale)*std::sqrt(sum), ExcNumberNotFinite());
828  return scale * std::sqrt(sum);
829  }
830 }
831 
832 
833 
834 template <typename Number>
837 {
838  Assert (vec_size!=0, ExcEmptyObject());
839 
840  if (p == 1.)
841  return l1_norm();
842  else if (p == 2.)
843  return l2_norm();
844 
845  real_type sum;
846  internal::Vector::accumulate (internal::Vector::NormP<Number,real_type>(),
847  val, val, p, vec_size, sum);
848 
849  if (numbers::is_finite(sum) && sum >= std::numeric_limits<real_type>::min())
850  return std::pow(sum, static_cast<real_type>(1./p));
851  else
852  {
853  real_type scale = 0.;
854  real_type sum = 1.;
855  for (size_type i=0; i<vec_size; ++i)
856  {
857  if (val[i] != Number())
858  {
859  const real_type abs_x =
861  if (scale < abs_x)
862  {
863  sum = 1. + sum * std::pow(scale/abs_x, p);
864  scale = abs_x;
865  }
866  else
867  sum += std::pow(abs_x/scale, p);
868  }
869  }
870  return scale * std::pow(sum, static_cast<real_type>(1./p));
871  }
872 }
873 
874 
875 
876 template <typename Number>
879 {
880  Assert (vec_size!=0, ExcEmptyObject());
881 
882  real_type max = 0.;
883 
884  for (size_type i=0; i<vec_size; ++i)
885  max = std::max (numbers::NumberTraits<Number>::abs(val[i]), max);
886 
887  return max;
888 }
889 
890 
891 template <typename Number>
893 {
894  Assert (vec_size!=0, ExcEmptyObject());
895 
896  add (v);
897  return *this;
898 }
899 
900 
901 template <typename Number>
903 {
904  Assert (vec_size!=0, ExcEmptyObject());
905  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
906 
907  if (vec_size>internal::Vector::minimum_parallel_grain_size)
909  val+vec_size,
910  v.val,
911  val,
912  (boost::lambda::_1 - boost::lambda::_2),
913  internal::Vector::minimum_parallel_grain_size);
914  else if (vec_size > 0)
915  for (size_type i=0; i<vec_size; ++i)
916  val[i] -= v.val[i];
917 
918  return *this;
919 }
920 
921 
922 template <typename Number>
923 void Vector<Number>::add (const Number v)
924 {
925  Assert (vec_size!=0, ExcEmptyObject());
926 
927  if (vec_size>internal::Vector::minimum_parallel_grain_size)
929  val+vec_size,
930  val,
931  (boost::lambda::_1 + v),
932  internal::Vector::minimum_parallel_grain_size);
933  else if (vec_size > 0)
934  for (size_type i=0; i<vec_size; ++i)
935  val[i] += v;
936 }
937 
938 
939 template <typename Number>
941 {
942  Assert (vec_size!=0, ExcEmptyObject());
943  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
944 
945  if (vec_size>internal::Vector::minimum_parallel_grain_size)
947  val+vec_size,
948  v.val,
949  val,
950  (boost::lambda::_1 + boost::lambda::_2),
951  internal::Vector::minimum_parallel_grain_size);
952  else if (vec_size > 0)
953  for (size_type i=0; i<vec_size; ++i)
954  val[i] += v.val[i];
955 }
956 
957 
958 template <typename Number>
959 void Vector<Number>::add (const Number a, const Vector<Number> &v,
960  const Number b, const Vector<Number> &w)
961 {
964 
965  Assert (vec_size!=0, ExcEmptyObject());
966  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
967  Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size));
968 
969  if (vec_size>internal::Vector::minimum_parallel_grain_size)
971  val+vec_size,
972  v.val,
973  w.val,
974  val,
975  (boost::lambda::_1 + a*boost::lambda::_2 + b*boost::lambda::_3),
976  internal::Vector::minimum_parallel_grain_size);
977  else if (vec_size > 0)
978  for (size_type i=0; i<vec_size; ++i)
979  val[i] += a * v.val[i] + b * w.val[i];
980 }
981 
982 
983 template <typename Number>
984 void Vector<Number>::sadd (const Number x,
985  const Vector<Number> &v)
986 {
988 
989  Assert (vec_size!=0, ExcEmptyObject());
990  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
991 
992  if (vec_size>internal::Vector::minimum_parallel_grain_size)
994  val+vec_size,
995  v.val,
996  val,
997  (x*boost::lambda::_1 + boost::lambda::_2),
998  internal::Vector::minimum_parallel_grain_size);
999  else if (vec_size > 0)
1000  for (size_type i=0; i<vec_size; ++i)
1001  val[i] = x * val[i] + v.val[i];
1002 }
1003 
1004 
1005 
1006 template <typename Number>
1007 void Vector<Number>::sadd (const Number x, const Number a,
1008  const Vector<Number> &v, const Number b,
1009  const Vector<Number> &w)
1010 {
1014 
1015  Assert (vec_size!=0, ExcEmptyObject());
1016  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
1017  Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size));
1018 
1019  if (vec_size>internal::Vector::minimum_parallel_grain_size)
1021  val+vec_size,
1022  v.val,
1023  w.val,
1024  val,
1025  (x*boost::lambda::_1 + a*boost::lambda::_2 + b*boost::lambda::_3),
1026  internal::Vector::minimum_parallel_grain_size);
1027  else if (vec_size > 0)
1028  for (size_type i=0; i<vec_size; ++i)
1029  val[i] = x*val[i] + a * v.val[i] + b * w.val[i];
1030 }
1031 
1032 
1033 template <typename Number>
1034 void Vector<Number>::sadd (const Number x, const Number a,
1035  const Vector<Number> &v, const Number b,
1036  const Vector<Number> &w, const Number c,
1037  const Vector<Number> &y)
1038 {
1039  sadd (x, a, v, b, w);
1040  add (c, y);
1041 }
1042 
1043 
1044 
1045 template <typename Number>
1047 {
1048  Assert (vec_size!=0, ExcEmptyObject());
1049  Assert (vec_size == s.vec_size, ExcDimensionMismatch(vec_size, s.vec_size));
1050 
1051  if (vec_size>internal::Vector::minimum_parallel_grain_size)
1053  val+vec_size,
1054  s.val,
1055  val,
1056  (boost::lambda::_1*boost::lambda::_2),
1057  internal::Vector::minimum_parallel_grain_size);
1058  else if (vec_size > 0)
1059  for (size_type i=0; i<vec_size; ++i)
1060  val[i] *= s.val[i];
1061 }
1062 
1063 
1064 
1065 template <typename Number>
1066 template <typename Number2>
1068 {
1069  Assert (vec_size!=0, ExcEmptyObject());
1070  Assert (vec_size == s.vec_size, ExcDimensionMismatch(vec_size, s.vec_size));
1071 
1072  for (size_type i=0; i<vec_size; ++i)
1073  val[i] *= Number(s.val[i]);
1074 }
1075 
1076 
1077 
1078 template <typename Number>
1079 void Vector<Number>::equ (const Number a,
1080  const Vector<Number> &u)
1081 {
1083 
1084  Assert (vec_size!=0, ExcEmptyObject());
1085  Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size));
1086 
1087  if (vec_size>internal::Vector::minimum_parallel_grain_size)
1088  parallel::transform (u.val,
1089  u.val+u.vec_size,
1090  val,
1091  (a*boost::lambda::_1),
1092  internal::Vector::minimum_parallel_grain_size);
1093  else if (vec_size > 0)
1094  for (size_type i=0; i<vec_size; ++i)
1095  val[i] = a * u.val[i];
1096 }
1097 
1098 
1099 
1100 template <typename Number>
1101 template <typename Number2>
1102 void Vector<Number>::equ (const Number a,
1103  const Vector<Number2> &u)
1104 {
1106 
1107  Assert (vec_size!=0, ExcEmptyObject());
1108  Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size));
1109 
1110  // set the result vector to a*u. we have to
1111  // convert the elements of u to the type of
1112  // the result vector. this is necessary
1113  // because
1114  // operator*(complex<float>,complex<double>)
1115  // is not defined by default
1116  for (size_type i=0; i<vec_size; ++i)
1117  val[i] = a * Number(u.val[i]);
1118 }
1119 
1120 
1121 
1122 template <typename Number>
1123 void Vector<Number>::equ (const Number a, const Vector<Number> &u,
1124  const Number b, const Vector<Number> &v)
1125 {
1128 
1129  Assert (vec_size!=0, ExcEmptyObject());
1130  Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size));
1131  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
1132 
1133  if (vec_size>internal::Vector::minimum_parallel_grain_size)
1134  parallel::transform (u.val,
1135  u.val+u.vec_size,
1136  v.val,
1137  val,
1138  (a*boost::lambda::_1 + b*boost::lambda::_2),
1139  internal::Vector::minimum_parallel_grain_size);
1140  else if (vec_size > 0)
1141  for (size_type i=0; i<vec_size; ++i)
1142  val[i] = a * u.val[i] + b * v.val[i];
1143 }
1144 
1145 
1146 template <typename Number>
1147 void Vector<Number>::equ (const Number a, const Vector<Number> &u,
1148  const Number b, const Vector<Number> &v,
1149  const Number c, const Vector<Number> &w)
1150 {
1151  Assert (vec_size!=0, ExcEmptyObject());
1152  Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size));
1153  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
1154  Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size));
1155 
1156  if (vec_size>internal::Vector::minimum_parallel_grain_size)
1157  parallel::transform (u.val,
1158  u.val+u.vec_size,
1159  v.val,
1160  w.val,
1161  val,
1162  (a*boost::lambda::_1 + b*boost::lambda::_2 + c*boost::lambda::_3),
1163  internal::Vector::minimum_parallel_grain_size);
1164  else if (vec_size > 0)
1165  for (size_type i=0; i<vec_size; ++i)
1166  val[i] = a * u.val[i] + b * v.val[i] + c * w.val[i];
1167 }
1168 
1169 
1170 template <typename Number>
1172  const Vector<Number> &b)
1173 {
1174  Assert (vec_size!=0, ExcEmptyObject());
1175  Assert (a.vec_size == b.vec_size,
1176  ExcDimensionMismatch (a.vec_size, b.vec_size));
1177 
1178  // no need to reinit with zeros, since
1179  // we overwrite them anyway
1180  reinit (a.size(), true);
1181 
1182  if (vec_size>internal::Vector::minimum_parallel_grain_size)
1183  parallel::transform (a.val,
1184  a.val+a.vec_size,
1185  b.val,
1186  val,
1187  (boost::lambda::_1 / boost::lambda::_2),
1188  internal::Vector::minimum_parallel_grain_size);
1189  else if (vec_size > 0)
1190  for (size_type i=0; i<vec_size; ++i)
1191  val[i] = a.val[i]/b.val[i];
1192 }
1193 
1194 
1195 
1196 template <typename Number>
1199 {
1200  if (v.size() != vec_size)
1201  reinit (v.size(), true);
1202 
1203  size_type this_index = 0;
1204  for (size_type b=0; b<v.n_blocks(); ++b)
1205  for (size_type i=0; i<v.block(b).size(); ++i, ++this_index)
1206  val[this_index] = v.block(b)(i);
1207 
1208  return *this;
1209 }
1210 
1211 
1212 
1213 #ifdef DEAL_II_WITH_PETSC
1214 
1215 template <typename Number>
1218 {
1219  if (v.size() != vec_size)
1220  reinit (v.size(), true);
1221  if (vec_size != 0)
1222  {
1223  // get a representation of the vector
1224  // and copy it
1225  PetscScalar *start_ptr;
1226  int ierr = VecGetArray (static_cast<const Vec &>(v), &start_ptr);
1227  AssertThrow (ierr == 0, ExcPETScError(ierr));
1228 
1229  internal::copy (start_ptr, start_ptr+vec_size, begin());
1230 
1231  // restore the representation of the
1232  // vector
1233  ierr = VecRestoreArray (static_cast<const Vec &>(v), &start_ptr);
1234  AssertThrow (ierr == 0, ExcPETScError(ierr));
1235  }
1236 
1237  return *this;
1238 }
1239 
1240 
1241 
1242 template <typename Number>
1245 {
1246  // do this in a two-stage process:
1247  // first convert to a sequential petsc
1248  // vector, then copy that
1249  PETScWrappers::Vector seq (v);
1250  *this = seq;
1251 
1252  return *this;
1253 }
1254 
1255 #endif
1256 
1257 
1258 #ifdef DEAL_II_WITH_TRILINOS
1259 
1260 template <typename Number>
1263 {
1264  // Generate a localized version
1265  // of the Trilinos vectors and
1266  // then call the other =
1267  // operator.
1268  TrilinosWrappers::Vector localized_vector (v);
1269  *this = localized_vector;
1270  return *this;
1271 }
1272 
1273 
1274 
1275 template <typename Number>
1278 {
1279  if (v.size() != vec_size)
1280  reinit (v.size(), true);
1281  if (vec_size != 0)
1282  {
1283  // get a representation of the vector
1284  // and copy it
1285  TrilinosScalar **start_ptr;
1286  int ierr = v.trilinos_vector().ExtractView (&start_ptr);
1287  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1288 
1289  std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());
1290  }
1291 
1292  return *this;
1293 }
1294 
1295 #endif
1296 
1297 template <typename Number>
1298 template <typename Number2>
1299 bool
1301 {
1302  Assert (vec_size!=0, ExcEmptyObject());
1303  Assert (vec_size == v.size(), ExcDimensionMismatch(vec_size, v.size()));
1304 
1305  // compare the two vector. we have to
1306  // convert the elements of v to the type of
1307  // the result vector. this is necessary
1308  // because
1309  // operator==(complex<float>,complex<double>)
1310  // is not defined by default
1311  for (size_type i=0; i<vec_size; ++i)
1312  if (val[i] != Number(v.val[i]))
1313  return false;
1314 
1315  return true;
1316 }
1317 
1318 
1319 
1320 template <typename Number>
1321 void Vector<Number>::print (const char *format) const
1322 {
1323  Assert (vec_size!=0, ExcEmptyObject());
1324 
1325  for (size_type j=0; j<size(); ++j)
1326  internal::print (val[j], format);
1327  std::printf ("\n");
1328 }
1329 
1330 
1331 
1332 template <typename Number>
1333 void Vector<Number>::print (std::ostream &out,
1334  const unsigned int precision,
1335  const bool scientific,
1336  const bool across) const
1337 {
1338  Assert (vec_size!=0, ExcEmptyObject());
1339  AssertThrow (out, ExcIO());
1340 
1341  std::ios::fmtflags old_flags = out.flags();
1342  unsigned int old_precision = out.precision (precision);
1343 
1344  out.precision (precision);
1345  if (scientific)
1346  out.setf (std::ios::scientific, std::ios::floatfield);
1347  else
1348  out.setf (std::ios::fixed, std::ios::floatfield);
1349 
1350  if (across)
1351  for (size_type i=0; i<size(); ++i)
1352  out << val[i] << ' ';
1353  else
1354  for (size_type i=0; i<size(); ++i)
1355  out << val[i] << std::endl;
1356  out << std::endl;
1357 
1358  AssertThrow (out, ExcIO());
1359  // reset output format
1360  out.flags (old_flags);
1361  out.precision(old_precision);
1362 }
1363 
1364 
1365 
1366 template <typename Number>
1367 void
1368 Vector<Number>::print (LogStream &out, const unsigned int width, const bool across) const
1369 {
1370  Assert (vec_size!=0, ExcEmptyObject());
1371 
1372  if (across)
1373  for (size_type i=0; i<size(); ++i)
1374  out << std::setw(width) << val[i] << ' ';
1375  else
1376  for (size_type i=0; i<size(); ++i)
1377  out << val[i] << std::endl;
1378  out << std::endl;
1379 }
1380 
1381 
1382 template <typename Number>
1383 void Vector<Number>::block_write (std::ostream &out) const
1384 {
1385  AssertThrow (out, ExcIO());
1386 
1387  // other version of the following
1388  // out << size() << std::endl << '[';
1389  // reason: operator<< seems to use
1390  // some resources that lead to
1391  // problems in a multithreaded
1392  // environment
1393  const size_type sz = size();
1394  char buf[16];
1395 
1396 #ifdef DEAL_II_USE_LARGE_INDEX_TYPE
1397  std::sprintf(buf, "%llu", sz);
1398 #else
1399  std::sprintf(buf, "%u", sz);
1400 #endif
1401  std::strcat(buf, "\n[");
1402 
1403  out.write(buf, std::strlen(buf));
1404  out.write (reinterpret_cast<const char *>(begin()),
1405  reinterpret_cast<const char *>(end())
1406  - reinterpret_cast<const char *>(begin()));
1407 
1408  // out << ']';
1409  const char outro = ']';
1410  out.write (&outro, 1);
1411 
1412  AssertThrow (out, ExcIO());
1413 }
1414 
1415 
1416 
1417 template <typename Number>
1418 void Vector<Number>::block_read (std::istream &in)
1419 {
1420  AssertThrow (in, ExcIO());
1421 
1422  size_type sz;
1423 
1424  char buf[16];
1425 
1426 
1427  in.getline(buf,16,'\n');
1428  sz=std::atoi(buf);
1429 
1430  // fast initialization, since the
1431  // data elements are overwritten anyway
1432  reinit (sz, true);
1433 
1434  char c;
1435  // in >> c;
1436  in.read (&c, 1);
1437  AssertThrow (c=='[', ExcIO());
1438 
1439  in.read (reinterpret_cast<char *>(begin()),
1440  reinterpret_cast<const char *>(end())
1441  - reinterpret_cast<const char *>(begin()));
1442 
1443  // in >> c;
1444  in.read (&c, 1);
1445  AssertThrow (c==']', ExcIO());
1446 }
1447 
1448 
1449 
1450 template <typename Number>
1451 IndexSet
1453 {
1454  return complete_index_set(size());
1455 }
1456 
1457 
1458 
1459 template <typename Number>
1460 std::size_t
1462 {
1463  return sizeof(*this) + (max_vec_size * sizeof(Number));
1464 }
1465 
1466 
1467 DEAL_II_NAMESPACE_CLOSE
1468 
1469 #endif
Vector< Number > & operator=(const Number s)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
virtual void reinit(const size_type N, const bool fast=false)
Number operator*(const Vector< Number2 > &V) const
::ExceptionBase & ExcOutOfMemory()
::ExceptionBase & ExcMessage(std::string arg1)
bool is_non_negative() const
void block_read(std::istream &in)
static real_type abs(const number &x)
Definition: numbers.h:316
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
numbers::NumberTraits< Number >::real_type real_type
Definition: vector.h:147
void apply_to_subranges(const RangeType &begin, const typename identity< RangeType >::type &end, const Function &f, const unsigned int grainsize)
Definition: parallel.h:464
bool is_finite(const double x)
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
Vector< Number > & operator+=(const Vector< Number > &V)
Vector< Number > & operator-=(const Vector< Number > &V)
::ExceptionBase & ExcIO()
real_type norm_sqr() const
void ratio(const Vector< Number > &a, const Vector< Number > &b)
real_type l2_norm() const
IndexSet locally_owned_elements() const
size_type size() const
size_type max_vec_size
Definition: vector.h:1104
static real_type abs_square(const number &x)
Definition: numbers.h:307
#define Assert(cond, exc)
Definition: exceptions.h:299
::ExceptionBase & ExcEmptyObject()
bool all_zero() const
real_type linfty_norm() const
static bool equal(const T *p1, const T *p2)
BlockType & block(const unsigned int i)
iterator begin()
Number * val
Definition: vector.h:1110
Vector< Number > & operator*=(const Number factor)
iterator end()
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:178
void print(const char *format=0) const
size_type vec_size
Definition: vector.h:1090
::ExceptionBase & ExcNumberNotFinite()
void equ(const Number a, const Vector< Number > &u)
void scale(const Number factor) DEAL_II_DEPRECATED
Definition: vector.h:873
void block_write(std::ostream &out) const
bool operator==(const Vector< Number2 > &v) const
real_type l1_norm() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void sadd(const Number s, const Vector< Number > &V)
const Epetra_MultiVector & trilinos_vector() const
Task< RT > new_task(const std_cxx1x::function< RT()> &function)
::ExceptionBase & ExcInternalError()
std::size_t memory_consumption() const
Number mean_value() const
size_type size() const
std::size_t size() const
real_type lp_norm(const real_type p) const