17 #ifndef __deal2__vector_templates_h
18 #define __deal2__vector_templates_h
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>
28 #ifdef DEAL_II_WITH_PETSC
29 # include <deal.II/lac/petsc_vector.h>
30 # include <deal.II/lac/petsc_parallel_vector.h>
33 #ifdef DEAL_II_WITH_TRILINOS
34 # include <deal.II/lac/trilinos_vector.h>
37 #include <boost/lambda/lambda.hpp>
50 bool is_non_negative (
const T &t)
57 bool is_non_negative (
const std::complex<T> &)
60 ExcMessage (
"Complex numbers do not have an ordering."));
67 void print (
const T &t,
71 std::printf (format, t);
73 std::printf (
" %5.2f",
double(t));
79 void print (
const std::complex<T> &t,
83 std::printf (format, t.real(), t.imag());
85 std::printf (
" %5.2f+%5.2fi",
86 double(t.real()),
double(t.imag()));
93 template <
typename T,
typename U>
94 void copy (
const T *begin,
98 std::copy (begin, end, dest);
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)
106 std::copy (begin, end, dest);
109 template <
typename T,
typename U>
110 void copy (
const std::complex<T> *,
111 const std::complex<T> *,
115 "into a vector of reals/doubles"));
122 template <
typename Number>
127 max_vec_size(v.size()),
139 #ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG
141 template <
typename Number>
142 template <
typename OtherNumber>
147 max_vec_size(v.size()),
160 #ifdef DEAL_II_WITH_PETSC
163 template <
typename Number>
168 max_vec_size(v.size()),
178 PetscScalar *start_ptr;
179 int ierr = VecGetArray (static_cast<const Vec &>(v), &start_ptr);
186 ierr = VecRestoreArray (static_cast<const Vec &>(v), &start_ptr);
193 template <
typename Number>
213 #ifdef DEAL_II_WITH_TRILINOS
215 template <
typename Number>
220 max_vec_size(v.size()),
238 TrilinosScalar **start_ptr;
240 int ierr = localized_vector.
trilinos_vector().ExtractView (&start_ptr);
249 template <
typename Number>
254 max_vec_size(v.size()),
264 TrilinosScalar **start_ptr;
275 template <
typename Number>
276 template <
typename Number2>
297 template <
typename Number>
303 for (size_type i=0; i<
vec_size; ++i)
304 if (
val[i] != Number(0))
311 template <
typename Number>
317 for (size_type i=0; i<
vec_size; ++i)
318 if ( ! internal::is_non_negative (
val[i]))
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,
337 std::memset ((dst.
begin()+begin),0,(end-begin)*
sizeof(T));
339 std::fill (&*(dst.
begin()+begin), &*(dst.
begin()+end), s);
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,
349 memcpy(&*(dst.
begin()+begin), &*(src.begin()+begin),
350 (end-begin)*
sizeof(T));
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,
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)
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,
374 copy_subrange (begin, end, src, dst);
378 template <
typename T,
typename U>
379 void copy_vector (const ::Vector<T> &src,
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)
391 std_cxx1x::bind(&internal::Vector::template
392 copy_subrange_wrap <T,U>,
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);
406 template <
typename Number>
413 if (vec_size>internal::Vector::minimum_parallel_grain_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);
427 #ifdef DEAL_II_BOOST_BIND_COMPILER_BUG
433 if (s != std::complex<float>())
436 std::fill (begin(), end(), s);
444 template <
typename Number>
454 (factor*boost::lambda::_1),
455 internal::Vector::minimum_parallel_grain_size);
462 template <
typename Number>
476 (boost::lambda::_1 + a*boost::lambda::_2),
477 internal::Vector::minimum_parallel_grain_size);
482 template <
typename Number>
498 (x*boost::lambda::_1 + a*boost::lambda::_2),
499 internal::Vector::minimum_parallel_grain_size);
510 template <
typename Number,
typename Number2>
514 operator() (
const Number *&X,
const Number2 *&Y,
const Number &)
const
520 template <
typename Number,
typename RealType>
524 operator() (
const Number *&X,
const Number *&,
const RealType &)
const
530 template <
typename Number,
typename RealType>
534 operator() (
const Number *&X,
const Number *&,
const RealType &)
const
540 template <
typename Number,
typename RealType>
544 operator() (
const Number *&X,
const Number *&,
const RealType &p)
const
550 template <
typename Number>
554 operator() (
const Number *&X,
const Number *&,
const Number &)
const
582 template <
typename Operation,
typename Number,
typename Number2,
584 void accumulate (
const Operation &op,
587 const ResultType power,
591 if (vec_size <= 4096)
596 const Number *X_original = X;
597 ResultType outer_results [128];
599 const size_type remainder = vec_size % 32;
604 ResultType r0 = op(X, Y, power);
606 r0 += op(X, Y, power);
607 ResultType r1 = op(X, Y, power);
609 r1 += op(X, Y, power);
611 r1 = op(X, Y, power);
613 r1 += op(X, Y, power);
614 ResultType r2 = op(X, Y, power);
616 r2 += op(X, Y, power);
619 outer_results[i] = r0;
627 const size_type inner_chunks = remainder / 8;
629 const size_type remainder_inner = remainder % 8;
630 ResultType r0 = ResultType(), r1 = ResultType(),
632 switch (inner_chunks)
635 r2 = op(X, Y, power);
637 r2 += op(X, Y, power);
640 r1 = op(X, Y, power);
642 r1 += op(X, Y, power);
646 r2 = op(X, Y, power);
648 r2 += op(X, Y, power);
651 for (
size_type j=0; j<remainder_inner; ++j)
652 r0 += op(X, Y, power);
655 outer_results[n_chunks] = r0;
666 if (n_chunks % 2 == 1)
667 outer_results[n_chunks++] = ResultType();
669 outer_results[i/2] = outer_results[i] + outer_results[i+1];
672 result = outer_results[0];
674 #ifdef DEAL_II_WITH_THREADS
675 else if (vec_size > 4 * internal::Vector::minimum_parallel_grain_size)
680 const size_type new_size = (vec_size / 4096) * 1024;
681 ResultType r0, r1, r2, r3;
685 op, X, Y, power, new_size, r0);
688 op, X+new_size, Y+new_size, power,
692 op, X+2*new_size, Y+2*new_size, power,
696 op, X+3*new_size, Y+3*new_size, power,
697 vec_size-3*new_size, r3);
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);
725 template <
typename Number>
726 template <
typename Number2>
739 val, v.
val, Number(), vec_size, sum);
746 template <
typename Number>
762 template <
typename Number>
769 val,
val, Number(), vec_size, sum);
776 template <
typename Number>
791 template <
typename Number>
806 norm_square >= std::numeric_limits<real_type>::min())
807 return std::sqrt(norm_square);
812 for (size_type i=0; i<
vec_size; ++i)
814 if (
val[i] != Number())
820 sum = 1. + sum * (scale/abs_x) * (scale/abs_x);
824 sum += (abs_x/
scale) * (abs_x/scale);
828 return scale * std::sqrt(sum);
834 template <
typename Number>
847 val,
val, p, vec_size, sum);
850 return std::pow(sum, static_cast<real_type>(1./p));
855 for (size_type i=0; i<
vec_size; ++i)
857 if (
val[i] != Number())
863 sum = 1. + sum * std::pow(scale/abs_x, p);
867 sum += std::pow(abs_x/scale, p);
870 return scale * std::pow(sum, static_cast<real_type>(1./p));
876 template <
typename Number>
884 for (size_type i=0; i<
vec_size; ++i)
891 template <
typename Number>
901 template <
typename Number>
907 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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)
922 template <
typename Number>
927 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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)
939 template <
typename Number>
945 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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)
958 template <
typename Number>
969 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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];
983 template <
typename Number>
992 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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];
1006 template <
typename Number>
1019 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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];
1033 template <
typename Number>
1039 sadd (x, a, v, b, w);
1045 template <
typename Number>
1051 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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)
1065 template <
typename Number>
1066 template <
typename Number2>
1072 for (size_type i=0; i<
vec_size; ++i)
1073 val[i] *= Number(s.
val[i]);
1078 template <
typename Number>
1087 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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];
1100 template <
typename Number>
1101 template <
typename Number2>
1116 for (size_type i=0; i<
vec_size; ++i)
1117 val[i] = a * Number(u.
val[i]);
1122 template <
typename Number>
1133 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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];
1146 template <
typename Number>
1156 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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];
1170 template <
typename Number>
1175 Assert (a.vec_size == b.vec_size,
1182 if (vec_size>internal::Vector::minimum_parallel_grain_size)
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];
1196 template <
typename Number>
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)
1213 #ifdef DEAL_II_WITH_PETSC
1215 template <
typename Number>
1225 PetscScalar *start_ptr;
1226 int ierr = VecGetArray (static_cast<const Vec &>(v), &start_ptr);
1229 internal::copy (start_ptr, start_ptr+vec_size,
begin());
1233 ierr = VecRestoreArray (static_cast<const Vec &>(v), &start_ptr);
1242 template <
typename Number>
1258 #ifdef DEAL_II_WITH_TRILINOS
1260 template <
typename Number>
1269 *
this = localized_vector;
1275 template <
typename Number>
1285 TrilinosScalar **start_ptr;
1289 std::copy (start_ptr[0], start_ptr[0]+vec_size,
begin());
1297 template <
typename Number>
1298 template <
typename Number2>
1311 for (size_type i=0; i<
vec_size; ++i)
1312 if (
val[i] != Number(v.
val[i]))
1320 template <
typename Number>
1325 for (size_type j=0; j<
size(); ++j)
1326 internal::print (
val[j], format);
1332 template <
typename Number>
1334 const unsigned int precision,
1335 const bool scientific,
1336 const bool across)
const
1341 std::ios::fmtflags old_flags = out.flags();
1342 unsigned int old_precision = out.precision (precision);
1344 out.precision (precision);
1346 out.setf (std::ios::scientific, std::ios::floatfield);
1348 out.setf (std::ios::fixed, std::ios::floatfield);
1351 for (size_type i=0; i<
size(); ++i)
1352 out <<
val[i] <<
' ';
1354 for (size_type i=0; i<
size(); ++i)
1355 out <<
val[i] << std::endl;
1360 out.flags (old_flags);
1361 out.precision(old_precision);
1366 template <
typename Number>
1373 for (size_type i=0; i<
size(); ++i)
1374 out << std::setw(width) <<
val[i] <<
' ';
1376 for (size_type i=0; i<
size(); ++i)
1377 out <<
val[i] << std::endl;
1382 template <
typename Number>
1393 const size_type sz =
size();
1396 #ifdef DEAL_II_USE_LARGE_INDEX_TYPE
1397 std::sprintf(buf,
"%llu", sz);
1399 std::sprintf(buf,
"%u", sz);
1401 std::strcat(buf,
"\n[");
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()));
1409 const char outro =
']';
1410 out.write (&outro, 1);
1417 template <
typename Number>
1427 in.getline(buf,16,
'\n');
1439 in.read (reinterpret_cast<char *>(
begin()),
1440 reinterpret_cast<const char *>(
end())
1441 - reinterpret_cast<const char *>(
begin()));
1450 template <
typename Number>
1454 return complete_index_set(
size());
1459 template <
typename Number>
1463 return sizeof(*this) + (
max_vec_size *
sizeof(Number));
1467 DEAL_II_NAMESPACE_CLOSE
Vector< Number > & operator=(const Number s)
#define AssertDimension(dim1, dim2)
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)
#define AssertThrow(cond, exc)
numbers::NumberTraits< Number >::real_type real_type
void apply_to_subranges(const RangeType &begin, const typename identity< RangeType >::type &end, const Function &f, const unsigned int grainsize)
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
static real_type abs_square(const number &x)
#define Assert(cond, exc)
::ExceptionBase & ExcEmptyObject()
real_type linfty_norm() const
static bool equal(const T *p1, const T *p2)
BlockType & block(const unsigned int i)
Vector< Number > & operator*=(const Number factor)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, Predicate &predicate, const unsigned int grainsize)
void print(const char *format=0) const
::ExceptionBase & ExcNumberNotFinite()
void equ(const Number a, const Vector< Number > &u)
unsigned int n_blocks() const
void scale(const Number factor) DEAL_II_DEPRECATED
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
real_type lp_norm(const real_type p) const