00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__vector_h
00013 #define __deal2__vector_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/logstream.h>
00018 #include <deal.II/base/exceptions.h>
00019 #include <deal.II/base/parallel.h>
00020 #include <deal.II/base/subscriptor.h>
00021 #include <boost/lambda/lambda.hpp>
00022 #include <boost/serialization/array.hpp>
00023 #include <boost/serialization/split_member.hpp>
00024
00025 #include <cstdio>
00026 #include <iostream>
00027 #include <cstring>
00028 #include <vector>
00029
00030 DEAL_II_NAMESPACE_OPEN
00031
00032
00033 #ifdef DEAL_II_USE_PETSC
00034 namespace PETScWrappers
00035 {
00036 class Vector;
00037 namespace MPI
00038 {
00039 class Vector;
00040 }
00041 }
00042 #endif
00043
00044 #ifdef DEAL_II_USE_TRILINOS
00045 namespace TrilinosWrappers
00046 {
00047 namespace MPI
00048 {
00049 class Vector;
00050 }
00051 class Vector;
00052 }
00053 #endif
00054
00055 template<typename number> class LAPACKFullMatrix;
00056
00057 template <typename> class BlockVector;
00058
00059 template <typename> class VectorView;
00060
00061
00062
00063
00088 template <typename Number>
00089 class Vector : public Subscriptor
00090 {
00091 public:
00098 typedef Number value_type;
00099 typedef value_type *pointer;
00100 typedef const value_type *const_pointer;
00101 typedef value_type *iterator;
00102 typedef const value_type *const_iterator;
00103 typedef value_type &reference;
00104 typedef const value_type &const_reference;
00105 typedef std::size_t size_type;
00106
00126 typedef typename numbers::NumberTraits<Number>::real_type real_type;
00127
00128 public:
00129
00138 Vector ();
00139
00149 Vector (const Vector<Number> &v);
00150
00151
00152 #ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG
00153
00173 template <typename OtherNumber>
00174 explicit
00175 Vector (const Vector<OtherNumber> &v);
00176 #endif
00177
00178 #ifdef DEAL_II_USE_PETSC
00179
00186 Vector (const PETScWrappers::Vector &v);
00187
00204 Vector (const PETScWrappers::MPI::Vector &v);
00205 #endif
00206
00207 #ifdef DEAL_II_USE_TRILINOS
00208
00229 Vector (const TrilinosWrappers::MPI::Vector &v);
00230
00240 Vector (const TrilinosWrappers::Vector &v);
00241 #endif
00242
00260 explicit Vector (const unsigned int n);
00261
00269 template <typename InputIterator>
00270 Vector (const InputIterator first,
00271 const InputIterator last);
00272
00279 virtual ~Vector ();
00280
00300 void compress () const;
00301
00327 virtual void reinit (const unsigned int N,
00328 const bool fast=false);
00329
00339 template <typename Number2>
00340 void reinit (const Vector<Number2> &V,
00341 const bool fast=false);
00342
00371 virtual void swap (Vector<Number> &v);
00372
00392 Vector<Number> & operator = (const Number s);
00393
00398 Vector<Number> & operator= (const Vector<Number> &c);
00399
00404 template <typename Number2>
00405 Vector<Number> & operator= (const Vector<Number2> &v);
00406
00412 Vector<Number> & operator= (const BlockVector<Number> &v);
00413
00414 #ifdef DEAL_II_USE_PETSC
00415
00423 Vector<Number> &
00424 operator = (const PETScWrappers::Vector &v);
00425
00444 Vector<Number> &
00445 operator = (const PETScWrappers::MPI::Vector &v);
00446 #endif
00447
00448
00449 #ifdef DEAL_II_USE_TRILINOS
00450
00470 Vector<Number> &
00471 operator = (const TrilinosWrappers::MPI::Vector &v);
00472
00482 Vector<Number> &
00483 operator = (const TrilinosWrappers::Vector &v);
00484 #endif
00485
00495 template <typename Number2>
00496 bool operator == (const Vector<Number2> &v) const;
00497
00506 template <typename Number2>
00507 bool operator != (const Vector<Number2> &v) const;
00508
00525 template <typename Number2>
00526 Number operator * (const Vector<Number2> &V) const;
00527
00531 real_type norm_sqr () const;
00532
00537 Number mean_value () const;
00538
00543 real_type l1_norm () const;
00544
00550 real_type l2_norm () const;
00551
00558 real_type lp_norm (const real_type p) const;
00559
00564 real_type linfty_norm () const;
00565
00569 unsigned int size () const;
00570
00579 bool all_zero () const;
00580
00595 bool is_non_negative () const;
00596
00604 iterator begin ();
00605
00610 const_iterator begin () const;
00611
00616 iterator end ();
00617
00622 const_iterator end () const;
00624
00625
00634 Number operator() (const unsigned int i) const;
00635
00640 Number& operator() (const unsigned int i);
00641
00648 Number operator[] (const unsigned int i) const;
00649
00656 Number& operator[] (const unsigned int i);
00658
00659
00664
00669 Vector<Number> & operator += (const Vector<Number> &V);
00670
00675 Vector<Number> & operator -= (const Vector<Number> &V);
00676
00685 template <typename OtherNumber>
00686 void add (const std::vector<unsigned int> &indices,
00687 const std::vector<OtherNumber> &values);
00688
00696 template <typename OtherNumber>
00697 void add (const std::vector<unsigned int> &indices,
00698 const Vector<OtherNumber> &values);
00699
00709 template <typename OtherNumber>
00710 void add (const unsigned int n_elements,
00711 const unsigned int *indices,
00712 const OtherNumber *values);
00713
00719 void add (const Number s);
00720
00725 void add (const Vector<Number> &V);
00726
00731 void add (const Number a, const Vector<Number> &V);
00732
00737 void add (const Number a, const Vector<Number> &V,
00738 const Number b, const Vector<Number> &W);
00739
00745 void sadd (const Number s,
00746 const Vector<Number> &V);
00747
00752 void sadd (const Number s,
00753 const Number a,
00754 const Vector<Number> &V);
00755
00759 void sadd (const Number s,
00760 const Number a,
00761 const Vector<Number> &V,
00762 const Number b,
00763 const Vector<Number> &W);
00764
00769 void sadd (const Number s,
00770 const Number a,
00771 const Vector<Number> &V,
00772 const Number b,
00773 const Vector<Number> &W,
00774 const Number c,
00775 const Vector<Number> &X);
00776
00787 void scale (const Number factor);
00788
00789
00795 Vector<Number> & operator *= (const Number factor);
00796
00802 Vector<Number> & operator /= (const Number factor);
00803
00813 void scale (const Vector<Number> &scaling_factors);
00814
00824 template <typename Number2>
00825 void scale (const Vector<Number2> &scaling_factors);
00826
00830 void equ (const Number a, const Vector<Number>& u);
00831
00835 template <typename Number2>
00836 void equ (const Number a, const Vector<Number2>& u);
00837
00841 void equ (const Number a, const Vector<Number>& u,
00842 const Number b, const Vector<Number>& v);
00843
00847 void equ (const Number a, const Vector<Number>& u,
00848 const Number b, const Vector<Number>& v,
00849 const Number c, const Vector<Number>& w);
00850
00867 void ratio (const Vector<Number> &a,
00868 const Vector<Number> &b);
00869
00885 void update_ghost_values () const;
00887
00888
00900 void print (const char* format = 0) const;
00901
00916 void print (std::ostream& out,
00917 const unsigned int precision = 3,
00918 const bool scientific = true,
00919 const bool across = true) const;
00920
00933 void print (LogStream& out,
00934 const unsigned int width = 6,
00935 const bool across = true) const;
00936
00946 void block_write (std::ostream &out) const;
00947
00967 void block_read (std::istream &in);
00968
00974 std::size_t memory_consumption () const;
00976
00981 template <class Archive>
00982 void save (Archive & ar, const unsigned int version) const;
00983
00988 template <class Archive>
00989 void load (Archive & ar, const unsigned int version);
00990
00991 BOOST_SERIALIZATION_SPLIT_MEMBER()
00992
00993 protected:
00994
01001 unsigned int vec_size;
01002
01015 unsigned int max_vec_size;
01016
01021 Number *val;
01022
01023
01024
01025
01026
01027 template <typename Number2> friend class Vector;
01028
01029
01030
01031
01032 friend class LAPACKFullMatrix<Number>;
01033
01034
01035
01036
01037 friend class VectorView<Number>;
01038 };
01039
01041
01042
01043
01044 #ifndef DOXYGEN
01045
01046
01047 template <typename Number>
01048 inline
01049 Vector<Number>::Vector ()
01050 :
01051 vec_size(0),
01052 max_vec_size(0),
01053 val(0)
01054 {}
01055
01056
01057
01058 template <typename Number>
01059 template <typename InputIterator>
01060 Vector<Number>::Vector (const InputIterator first, const InputIterator last)
01061 :
01062 vec_size (0),
01063 max_vec_size (0),
01064 val (0)
01065 {
01066
01067
01068
01069 reinit (std::distance (first, last), true);
01070 std::copy (first, last, begin());
01071 }
01072
01073
01074
01075 template <typename Number>
01076 inline
01077 Vector<Number>::Vector (const unsigned int n)
01078 :
01079 vec_size(0),
01080 max_vec_size(0),
01081 val(0)
01082 {
01083 reinit (n, false);
01084 }
01085
01086
01087
01088 template <typename Number>
01089 inline
01090 Vector<Number>::~Vector ()
01091 {
01092 if (val)
01093 {
01094 delete[] val;
01095 val=0;
01096 }
01097 }
01098
01099
01100
01101 template <typename Number>
01102 inline
01103 void Vector<Number>::reinit (const unsigned int n, const bool fast)
01104 {
01105 if (n==0)
01106 {
01107 if (val) delete[] val;
01108 val = 0;
01109 max_vec_size = vec_size = 0;
01110 return;
01111 };
01112
01113 if (n>max_vec_size)
01114 {
01115 if (val) delete[] val;
01116 val = new value_type[n];
01117 Assert (val != 0, ExcOutOfMemory());
01118 max_vec_size = n;
01119 };
01120 vec_size = n;
01121 if (fast == false)
01122 *this = static_cast<Number>(0);
01123 }
01124
01125
01126
01127 namespace internal
01128 {
01129 namespace Vector
01130 {
01131 template<typename T>
01132 void set_subrange (const T s,
01133 const unsigned int begin,
01134 const unsigned int end,
01135 ::Vector<T> &dst)
01136 {
01137 if (s == T())
01138 std::memset ((dst.begin()+begin),0,(end-begin)*sizeof(T));
01139 else
01140 std::fill (&*(dst.begin()+begin), &*(dst.begin()+end), s);
01141 }
01142
01143 template<typename T>
01144 void copy_subrange (const ::Vector<T>&src,
01145 const unsigned int begin,
01146 const unsigned int end,
01147 ::Vector<T> &dst)
01148 {
01149 memcpy(&*(dst.begin()+begin), &*(src.begin()+begin),
01150 (end-begin)*sizeof(T));
01151 }
01152
01153 template<typename T, typename U>
01154 void copy_subrange_ext (const ::Vector<T>&src,
01155 const unsigned int begin,
01156 const unsigned int end,
01157 ::Vector<U> &dst)
01158 {
01159 const T* q = src.begin()+begin;
01160 const T* const end_q = src.begin()+end;
01161 U* p = dst.begin()+begin;
01162 for (; q!=end_q; ++q, ++p)
01163 *p = *q;
01164 }
01165 }
01166 }
01167
01168
01169
01170 template <typename Number>
01171 inline
01172 Vector<Number> & Vector<Number>::operator = (const Number s)
01173 {
01174 Assert (numbers::is_finite(s), ExcNumberNotFinite());
01175 if (s != Number())
01176 Assert (vec_size!=0, ExcEmptyObject());
01177 if (vec_size>internal::Vector::minimum_parallel_grain_size)
01178 parallel::apply_to_subranges (0U, vec_size,
01179 std_cxx1x::bind(&internal::Vector::template
01180 set_subrange<Number>,
01181 s, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::ref(*this)),
01182 internal::Vector::minimum_parallel_grain_size);
01183 else if (vec_size > 0)
01184 internal::Vector::set_subrange<Number>(s, 0U, vec_size, *this);
01185
01186 return *this;
01187 }
01188
01189
01190
01191 #ifdef DEAL_II_BOOST_BIND_COMPILER_BUG
01192 template <>
01193 inline
01194 Vector<std::complex<float> > &
01195 Vector<std::complex<float> >::operator = (const std::complex<float> s)
01196 {
01197 Assert (numbers::is_finite(s), ExcNumberNotFinite());
01198 if (s != std::complex<float>())
01199 Assert (vec_size!=0, ExcEmptyObject());
01200 if (vec_size!=0)
01201 std::fill (begin(), end(), s);
01202
01203 return *this;
01204 }
01205 #endif
01206
01207
01208
01209 template <typename Number>
01210 inline
01211 Vector<Number> &
01212 Vector<Number>::operator = (const Vector<Number>& v)
01213 {
01214
01215
01216 if (PointerComparison::equal(this, &v) == true)
01217 return *this;
01218
01219 if (v.vec_size != vec_size)
01220 reinit (v.vec_size, true);
01221 if (vec_size>internal::Vector::minimum_parallel_grain_size)
01222 parallel::apply_to_subranges (0U, vec_size,
01223 std_cxx1x::bind(&internal::Vector::template
01224 copy_subrange<Number>,
01225 std_cxx1x::cref(v), std_cxx1x::_1, std_cxx1x::_2,
01226 std_cxx1x::ref(*this)),
01227 internal::Vector::minimum_parallel_grain_size);
01228 else if (vec_size > 0)
01229 internal::Vector::copy_subrange<Number>(v, 0U, vec_size, *this);
01230
01231 return *this;
01232 }
01233
01234
01235
01236 template <typename Number>
01237 template <typename Number2>
01238 inline
01239 Vector<Number> &
01240 Vector<Number>::operator = (const Vector<Number2>& v)
01241 {
01242 if (v.vec_size != vec_size)
01243 reinit (v.vec_size, true);
01244 if (vec_size>internal::Vector::minimum_parallel_grain_size)
01245 parallel::apply_to_subranges (0U, vec_size,
01246 std_cxx1x::bind(&internal::Vector::template
01247 copy_subrange_ext<Number2,Number>,
01248 std_cxx1x::cref(v), std_cxx1x::_1, std_cxx1x::_2,
01249 std_cxx1x::ref(*this)),
01250 internal::Vector::minimum_parallel_grain_size);
01251 else if (vec_size > 0)
01252 internal::Vector::copy_subrange_ext<Number2,Number>(v, 0U, vec_size, *this);
01253
01254 return *this;
01255 }
01256
01257
01258
01259 template <typename Number>
01260 inline
01261 unsigned int Vector<Number>::size () const
01262 {
01263 return vec_size;
01264 }
01265
01266
01267
01268 template <typename Number>
01269 inline
01270 typename Vector<Number>::iterator
01271 Vector<Number>::begin ()
01272 {
01273 return &val[0];
01274 }
01275
01276
01277
01278 template <typename Number>
01279 inline
01280 typename Vector<Number>::const_iterator
01281 Vector<Number>::begin () const
01282 {
01283 return &val[0];
01284 }
01285
01286
01287
01288 template <typename Number>
01289 inline
01290 typename Vector<Number>::iterator
01291 Vector<Number>::end ()
01292 {
01293 return &val[vec_size];
01294 }
01295
01296
01297
01298 template <typename Number>
01299 inline
01300 typename Vector<Number>::const_iterator
01301 Vector<Number>::end () const
01302 {
01303 return &val[vec_size];
01304 }
01305
01306
01307
01308 template <typename Number>
01309 inline
01310 Number Vector<Number>::operator() (const unsigned int i) const
01311 {
01312 Assert (i<vec_size, ExcIndexRange(i,0,vec_size));
01313 return val[i];
01314 }
01315
01316
01317
01318 template <typename Number>
01319 inline
01320 Number& Vector<Number>::operator() (const unsigned int i)
01321 {
01322 Assert (i<vec_size, ExcIndexRange(i,0,vec_size));
01323 return val[i];
01324 }
01325
01326
01327
01328 template <typename Number>
01329 inline
01330 Number Vector<Number>::operator[] (const unsigned int i) const
01331 {
01332 return operator()(i);
01333 }
01334
01335
01336
01337 template <typename Number>
01338 inline
01339 Number& Vector<Number>::operator[] (const unsigned int i)
01340 {
01341 return operator()(i);
01342 }
01343
01344
01345
01346 template <typename Number>
01347 inline
01348 Vector<Number> & Vector<Number>::operator *= (const Number factor)
01349 {
01350
01351 Assert (numbers::is_finite(factor),ExcNumberNotFinite());
01352
01353 scale (factor);
01354 return *this;
01355 }
01356
01357
01358
01359 template <typename Number>
01360 inline
01361 Vector<Number> &
01362 Vector<Number>::operator /= (const Number factor)
01363 {
01364 Assert (numbers::is_finite(factor),ExcNumberNotFinite());
01365 Assert (factor != Number(0.), ExcZero() );
01366
01367 this->operator *= (Number(1.)/factor);
01368 return *this;
01369 }
01370
01371
01372
01373 template <typename Number>
01374 inline
01375 void
01376 Vector<Number>::scale (const Number factor)
01377 {
01378 Assert (numbers::is_finite(factor),ExcNumberNotFinite());
01379
01380 Assert (vec_size!=0, ExcEmptyObject());
01381
01382 parallel::transform (val,
01383 val+vec_size,
01384 val,
01385 (factor*boost::lambda::_1),
01386 internal::Vector::minimum_parallel_grain_size);
01387 }
01388
01389
01390
01391 template <typename Number>
01392 template <typename OtherNumber>
01393 inline
01394 void
01395 Vector<Number>::add (const std::vector<unsigned int> &indices,
01396 const std::vector<OtherNumber> &values)
01397 {
01398 Assert (indices.size() == values.size(),
01399 ExcDimensionMismatch(indices.size(), values.size()));
01400 add (indices.size(), &indices[0], &values[0]);
01401 }
01402
01403
01404
01405 template <typename Number>
01406 template <typename OtherNumber>
01407 inline
01408 void
01409 Vector<Number>::add (const std::vector<unsigned int> &indices,
01410 const Vector<OtherNumber> &values)
01411 {
01412 Assert (indices.size() == values.size(),
01413 ExcDimensionMismatch(indices.size(), values.size()));
01414 add (indices.size(), &indices[0], values.val);
01415 }
01416
01417
01418
01419 template <typename Number>
01420 template <typename OtherNumber>
01421 inline
01422 void
01423 Vector<Number>::add (const unsigned int n_indices,
01424 const unsigned int *indices,
01425 const OtherNumber *values)
01426 {
01427 for (unsigned int i=0; i<n_indices; ++i)
01428 {
01429 Assert (indices[i] < vec_size, ExcIndexRange(indices[i],0,vec_size));
01430 Assert (numbers::is_finite(values[i]),
01431 ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
01432
01433 val[indices[i]] += values[i];
01434 }
01435 }
01436
01437
01438
01439 template <typename Number>
01440 inline
01441 void
01442 Vector<Number>::add (const Number a,
01443 const Vector<Number>& v)
01444 {
01445 Assert (numbers::is_finite(a),ExcNumberNotFinite());
01446
01447 Assert (vec_size!=0, ExcEmptyObject());
01448 Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
01449
01450 parallel::transform (val,
01451 val+vec_size,
01452 v.val,
01453 val,
01454 (boost::lambda::_1 + a*boost::lambda::_2),
01455 internal::Vector::minimum_parallel_grain_size);
01456 }
01457
01458
01459
01460 template <typename Number>
01461 inline
01462 void
01463 Vector<Number>::sadd (const Number x,
01464 const Number a,
01465 const Vector<Number>& v)
01466 {
01467 Assert (numbers::is_finite(x),ExcNumberNotFinite());
01468 Assert (numbers::is_finite(a),ExcNumberNotFinite());
01469
01470 Assert (vec_size!=0, ExcEmptyObject());
01471 Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
01472
01473 parallel::transform (val,
01474 val+vec_size,
01475 v.val,
01476 val,
01477 (x*boost::lambda::_1 + a*boost::lambda::_2),
01478 internal::Vector::minimum_parallel_grain_size);
01479 }
01480
01481
01482
01483
01484 template <typename Number>
01485 template <typename Number2>
01486 inline
01487 bool
01488 Vector<Number>::operator != (const Vector<Number2>& v) const
01489 {
01490 return ! (*this == v);
01491 }
01492
01493
01494
01495 template <typename Number>
01496 inline
01497 void
01498 Vector<Number>::compress () const
01499 {}
01500
01501
01502
01503 template <typename Number>
01504 inline
01505 void
01506 Vector<Number>::update_ghost_values () const
01507 {}
01508
01509
01510
01511
01512
01513
01514 template <typename Number>
01515 inline
01516 void
01517 Vector<Number>::swap (Vector<Number> &v)
01518 {
01519 std::swap (vec_size, v.vec_size);
01520 std::swap (max_vec_size, v.max_vec_size);
01521 std::swap (val, v.val);
01522 }
01523
01524
01525
01526 template <typename Number>
01527 template <class Archive>
01528 inline
01529 void
01530 Vector<Number>::save (Archive & ar, const unsigned int) const
01531 {
01532
01533
01534 ar & static_cast<const Subscriptor &>(*this);
01535
01536 ar & vec_size & max_vec_size ;
01537 ar & boost::serialization::make_array(val, max_vec_size);
01538 }
01539
01540
01541
01542 template <typename Number>
01543 template <class Archive>
01544 inline
01545 void
01546 Vector<Number>::load (Archive & ar, const unsigned int)
01547 {
01548
01549
01550 ar & static_cast<Subscriptor &>(*this);
01551
01552 ar & vec_size & max_vec_size ;
01553
01554 val = new Number[max_vec_size];
01555 ar & boost::serialization::make_array(val, max_vec_size);
01556 }
01557
01558
01559 #endif
01560
01561
01575 template <typename Number>
01576 inline
01577 void swap (Vector<Number> &u, Vector<Number> &v)
01578 {
01579 u.swap (v);
01580 }
01581
01582
01586 template <typename number>
01587 inline
01588 std::ostream&
01589 operator << (std::ostream& os, const Vector<number>& v)
01590 {
01591 v.print(os);
01592 return os;
01593 }
01594
01598 template <typename number>
01599 inline
01600 LogStream&
01601 operator << (LogStream& os, const Vector<number>& v)
01602 {
01603 v.print(os);
01604 return os;
01605 }
01606
01607
01610 DEAL_II_NAMESPACE_CLOSE
01611
01612 #endif
01613