00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__block_vector_base_h
00013 #define __deal2__block_vector_base_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/subscriptor.h>
00019 #include <deal.II/base/memory_consumption.h>
00020 #include <deal.II/lac/exceptions.h>
00021 #include <deal.II/lac/block_indices.h>
00022 #include <deal.II/lac/vector.h>
00023 #include <vector>
00024 #include <iterator>
00025 #include <cmath>
00026 #include <cstddef>
00027
00028 DEAL_II_NAMESPACE_OPEN
00029
00030
00035 template <typename> class BlockVectorBase;
00036
00037
00054 template <typename VectorType>
00055 struct IsBlockVector
00056 {
00057 private:
00058 struct yes_type { char c[1]; };
00059 struct no_type { char c[2]; };
00060
00066 template <typename T>
00067 static yes_type check_for_block_vector (const BlockVectorBase<T> *);
00068
00074 static no_type check_for_block_vector (...);
00075
00076 public:
00084 static const bool value = (sizeof(check_for_block_vector
00085 ((VectorType*)0))
00086 ==
00087 sizeof(yes_type));
00088 };
00089
00090
00091
00092 template <typename VectorType>
00093 const bool IsBlockVector<VectorType>::value;
00094
00095
00096
00097
00098 namespace internal
00099 {
00100
00106 namespace BlockVectorIterators
00107 {
00115 template <class BlockVectorType, bool constness>
00116 struct Types
00117 {
00118 };
00119
00120
00121
00132 template <class BlockVectorType>
00133 struct Types<BlockVectorType,false>
00134 {
00142 typedef typename BlockVectorType::BlockType Vector;
00143
00150 typedef BlockVectorType BlockVector;
00151
00157 typedef typename BlockVector::value_type value_type;
00158
00165 typedef typename Vector::reference dereference_type;
00166 };
00167
00168
00169
00180 template <class BlockVectorType>
00181 struct Types<BlockVectorType,true>
00182 {
00191 typedef const typename BlockVectorType::BlockType Vector;
00192
00200 typedef const BlockVectorType BlockVector;
00201
00208 typedef const typename BlockVector::value_type value_type;
00209
00219 typedef value_type dereference_type;
00220 };
00221
00222
00281 template <class BlockVectorType, bool constness>
00282 class Iterator :
00283 #ifdef HAVE_STD_ITERATOR_CLASS
00284 public std::iterator<std::random_access_iterator_tag,
00285 typename Types<BlockVectorType,constness>::value_type>
00286 #else
00287 random_access_iterator<typename Types<BlockVectorType,constness>::value_type,int>
00288 #endif
00289 {
00290 private:
00297 typedef Iterator<BlockVectorType,!constness> InverseConstnessIterator;
00298
00299 public:
00309 typedef
00310 typename Types<BlockVectorType,constness>::value_type
00311 value_type;
00312
00321 typedef std::random_access_iterator_tag iterator_type;
00322 typedef std::ptrdiff_t difference_type;
00323 typedef typename BlockVectorType::reference reference;
00324 typedef value_type *pointer;
00325
00326 typedef
00327 typename Types<BlockVectorType,constness>::dereference_type
00328 dereference_type;
00329
00337 typedef
00338 typename Types<BlockVectorType,constness>::BlockVector
00339 BlockVector;
00340
00355 Iterator (BlockVector &parent,
00356 const unsigned int global_index);
00357
00361 Iterator (const Iterator<BlockVectorType,constness> &c);
00362
00373 Iterator (const InverseConstnessIterator &c);
00374
00375 private:
00383 Iterator (BlockVector &parent,
00384 const unsigned int global_index,
00385 const unsigned int current_block,
00386 const unsigned int index_within_block,
00387 const unsigned int next_break_forward,
00388 const unsigned int next_break_backward);
00389
00390 public:
00391
00395 Iterator & operator = (const Iterator &c);
00396
00405 dereference_type operator * () const;
00406
00413 dereference_type operator [] (const difference_type d) const;
00414
00421 Iterator & operator ++ ();
00422
00431 Iterator operator ++ (int);
00432
00439 Iterator & operator -- ();
00440
00449 Iterator operator -- (int);
00450
00459 bool operator == (const Iterator &i) const;
00460
00466 bool operator == (const InverseConstnessIterator &i) const;
00467
00476 bool operator != (const Iterator &i) const;
00477
00483 bool operator != (const InverseConstnessIterator &i) const;
00484
00496 bool operator < (const Iterator &i) const;
00497
00503 bool operator < (const InverseConstnessIterator &i) const;
00504
00509 bool operator <= (const Iterator &i) const;
00510
00516 bool operator <= (const InverseConstnessIterator &i) const;
00517
00522 bool operator > (const Iterator &i) const;
00523
00529 bool operator > (const InverseConstnessIterator &i) const;
00530
00535 bool operator >= (const Iterator &i) const;
00536
00542 bool operator >= (const InverseConstnessIterator &i) const;
00543
00549 difference_type operator - (const Iterator &i) const;
00550
00555 difference_type operator - (const InverseConstnessIterator &i) const;
00556
00562 Iterator operator + (const difference_type &d) const;
00563
00569 Iterator operator - (const difference_type &d) const;
00570
00576 Iterator & operator += (const difference_type &d);
00577
00583 Iterator & operator -= (const difference_type &d);
00584
00591 DeclException0 (ExcPointerToDifferentVectors);
00595 DeclException0 (ExcCastingAwayConstness);
00597 private:
00608 BlockVector *parent;
00609
00615 unsigned int global_index;
00616
00623 unsigned int current_block;
00624 unsigned int index_within_block;
00625
00638 unsigned int next_break_forward;
00639 unsigned int next_break_backward;
00640
00644 void move_forward ();
00645
00649 void move_backward ();
00650
00651
00652 #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG
00653
00668 template <typename N, bool C>
00669 friend class Iterator;
00670 #else
00671 friend class InverseConstnessIterator;
00672 #endif
00673 };
00674 }
00675 }
00676
00677
00718 template <class VectorType>
00719 class BlockVectorBase : public Subscriptor
00720 {
00721 public:
00726 typedef VectorType BlockType;
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738 typedef typename BlockType::value_type value_type;
00739 typedef value_type *pointer;
00740 typedef const value_type *const_pointer;
00741 typedef internal::BlockVectorIterators::Iterator<BlockVectorBase,false> iterator;
00742 typedef internal::BlockVectorIterators::Iterator<BlockVectorBase,true> const_iterator;
00743 typedef typename BlockType::reference reference;
00744 typedef typename BlockType::const_reference const_reference;
00745
00746 typedef std::size_t size_type;
00747
00767 typedef typename BlockType::real_type real_type;
00768
00772 BlockVectorBase ();
00773
00785 void collect_sizes ();
00786
00791 void compress ();
00792
00796 BlockType &
00797 block (const unsigned int i);
00798
00802 const BlockType &
00803 block (const unsigned int i) const;
00804
00815 const BlockIndices &
00816 get_block_indices () const;
00817
00821 unsigned int n_blocks () const;
00822
00828 unsigned int size () const;
00829
00834 iterator begin ();
00835
00841 const_iterator begin () const;
00842
00847 iterator end ();
00848
00854 const_iterator end () const;
00855
00859 value_type operator() (const unsigned int i) const;
00860
00865 reference operator() (const unsigned int i);
00866
00872 value_type operator[] (const unsigned int i) const;
00873
00880 reference operator[] (const unsigned int i);
00881
00887 BlockVectorBase & operator = (const value_type s);
00888
00893 BlockVectorBase &
00894 operator= (const BlockVectorBase& V);
00895
00900 template <class VectorType2>
00901 BlockVectorBase &
00902 operator= (const BlockVectorBase<VectorType2> &V);
00903
00908 BlockVectorBase &
00909 operator = (const VectorType &v);
00910
00917 template <class VectorType2>
00918 bool
00919 operator == (const BlockVectorBase<VectorType2> &v) const;
00920
00924 value_type operator* (const BlockVectorBase& V) const;
00925
00929 real_type norm_sqr () const;
00930
00935 value_type mean_value () const;
00936
00941 real_type l1_norm () const;
00942
00948 real_type l2_norm () const;
00949
00955 real_type linfty_norm () const;
00956
00965 bool all_zero () const;
00966
00975 bool is_non_negative () const;
00976
00981 BlockVectorBase &
00982 operator += (const BlockVectorBase &V);
00983
00988 BlockVectorBase &
00989 operator -= (const BlockVectorBase &V);
00990
00991
01000 template <typename Number>
01001 void add (const std::vector<unsigned int> &indices,
01002 const std::vector<Number> &values);
01003
01011 template <typename Number>
01012 void add (const std::vector<unsigned int> &indices,
01013 const Vector<Number> &values);
01014
01024 template <typename Number>
01025 void add (const unsigned int n_elements,
01026 const unsigned int *indices,
01027 const Number *values);
01028
01035 void add (const value_type s);
01036
01042 void add (const BlockVectorBase& V);
01043
01048 void add (const value_type a, const BlockVectorBase& V);
01049
01054 void add (const value_type a, const BlockVectorBase& V,
01055 const value_type b, const BlockVectorBase& W);
01056
01061 void sadd (const value_type s, const BlockVectorBase& V);
01062
01067 void sadd (const value_type s, const value_type a, const BlockVectorBase& V);
01068
01073 void sadd (const value_type s, const value_type a,
01074 const BlockVectorBase& V,
01075 const value_type b, const BlockVectorBase& W);
01076
01081 void sadd (const value_type s, const value_type a,
01082 const BlockVectorBase& V,
01083 const value_type b, const BlockVectorBase& W,
01084 const value_type c, const BlockVectorBase& X);
01085
01091 BlockVectorBase & operator *= (const value_type factor);
01092
01098 BlockVectorBase & operator /= (const value_type factor);
01099
01105 template <class BlockVector2>
01106 void scale (const BlockVector2 &v);
01107
01111 template <class BlockVector2>
01112 void equ (const value_type a, const BlockVector2 &V);
01113
01118 void equ (const value_type a, const BlockVectorBase& V,
01119 const value_type b, const BlockVectorBase& W);
01120
01136 void update_ghost_values () const;
01137
01143 std::size_t memory_consumption () const;
01144
01145 protected:
01149 std::vector<VectorType> components;
01150
01157 BlockIndices block_indices;
01158
01164 #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG
01165 template <typename N, bool C>
01166 friend class internal::BlockVectorIterators::Iterator;
01167 #else
01168 friend class iterator;
01169 friend class const_iterator;
01170 #endif
01171
01172 template <typename> friend class BlockVectorBase;
01173 };
01174
01175
01178
01179
01180
01181 #ifndef DOXYGEN
01182 namespace internal
01183 {
01184 namespace BlockVectorIterators
01185 {
01186
01187 template <class BlockVectorType, bool constness>
01188 inline
01189 Iterator<BlockVectorType,constness>::
01190 Iterator (const Iterator<BlockVectorType,constness> &c)
01191 :
01192 parent (c.parent),
01193 global_index (c.global_index),
01194 current_block (c.current_block),
01195 index_within_block (c.index_within_block),
01196 next_break_forward (c.next_break_forward),
01197 next_break_backward (c.next_break_backward)
01198 {}
01199
01200
01201
01202 template <class BlockVectorType, bool constness>
01203 inline
01204 Iterator<BlockVectorType,constness>::
01205 Iterator (const InverseConstnessIterator &c)
01206 :
01207 parent (const_cast<BlockVectorType*>(c.parent)),
01208 global_index (c.global_index),
01209 current_block (c.current_block),
01210 index_within_block (c.index_within_block),
01211 next_break_forward (c.next_break_forward),
01212 next_break_backward (c.next_break_backward)
01213 {
01214
01215
01216
01217
01218
01219 Assert (constness==true, ExcCastingAwayConstness());
01220 }
01221
01222
01223
01224 template <class BlockVectorType, bool constness>
01225 inline
01226 Iterator<BlockVectorType,constness>::
01227 Iterator (BlockVector &parent,
01228 const unsigned int global_index,
01229 const unsigned int current_block,
01230 const unsigned int index_within_block,
01231 const unsigned int next_break_forward,
01232 const unsigned int next_break_backward)
01233 :
01234 parent (&parent),
01235 global_index (global_index),
01236 current_block (current_block),
01237 index_within_block (index_within_block),
01238 next_break_forward (next_break_forward),
01239 next_break_backward (next_break_backward)
01240 {
01241 }
01242
01243
01244
01245 template <class BlockVectorType, bool constness>
01246 inline
01247 Iterator<BlockVectorType,constness> &
01248 Iterator<BlockVectorType,constness>::
01249 operator = (const Iterator &c)
01250 {
01251 parent = c.parent;
01252 global_index = c.global_index;
01253 index_within_block = c.index_within_block;
01254 current_block = c.current_block;
01255 next_break_forward = c.next_break_forward;
01256 next_break_backward = c.next_break_backward;
01257
01258 return *this;
01259 }
01260
01261
01262
01263 template <class BlockVectorType, bool constness>
01264 inline
01265 typename Iterator<BlockVectorType,constness>::dereference_type
01266 Iterator<BlockVectorType,constness>::operator * () const
01267 {
01268 return parent->block(current_block)(index_within_block);
01269 }
01270
01271
01272
01273 template <class BlockVectorType, bool constness>
01274 inline
01275 typename Iterator<BlockVectorType,constness>::dereference_type
01276 Iterator<BlockVectorType,constness>::operator [] (const difference_type d) const
01277 {
01278
01279
01280
01281
01282
01283 if ((global_index+d >= next_break_backward) &&
01284 (global_index+d <= next_break_forward))
01285 return parent->block(current_block)(index_within_block + d);
01286
01287
01288
01289
01290
01291
01292
01293
01294 return (*parent)(global_index+d);
01295 }
01296
01297
01298
01299 template <class BlockVectorType, bool constness>
01300 inline
01301 Iterator<BlockVectorType,constness> &
01302 Iterator<BlockVectorType,constness>::operator ++ ()
01303 {
01304 move_forward ();
01305 return *this;
01306 }
01307
01308
01309
01310 template <class BlockVectorType, bool constness>
01311 inline
01312 Iterator<BlockVectorType,constness>
01313 Iterator<BlockVectorType,constness>::operator ++ (int)
01314 {
01315 const Iterator old_value = *this;
01316 move_forward ();
01317 return old_value;
01318 }
01319
01320
01321
01322 template <class BlockVectorType, bool constness>
01323 inline
01324 Iterator<BlockVectorType,constness> &
01325 Iterator<BlockVectorType,constness>::operator -- ()
01326 {
01327 move_backward ();
01328 return *this;
01329 }
01330
01331
01332
01333 template <class BlockVectorType, bool constness>
01334 inline
01335 Iterator<BlockVectorType,constness>
01336 Iterator<BlockVectorType,constness>::operator -- (int)
01337 {
01338 const Iterator old_value = *this;
01339 move_backward ();
01340 return old_value;
01341 }
01342
01343
01344
01345 template <class BlockVectorType, bool constness>
01346 inline
01347 bool
01348 Iterator<BlockVectorType,constness>::
01349 operator == (const Iterator &i) const
01350 {
01351 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01352
01353 return (global_index == i.global_index);
01354 }
01355
01356
01357
01358 template <class BlockVectorType, bool constness>
01359 inline
01360 bool
01361 Iterator<BlockVectorType,constness>::
01362 operator == (const InverseConstnessIterator &i) const
01363 {
01364 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01365
01366 return (global_index == i.global_index);
01367 }
01368
01369
01370
01371 template <class BlockVectorType, bool constness>
01372 inline
01373 bool
01374 Iterator<BlockVectorType,constness>::
01375 operator != (const Iterator &i) const
01376 {
01377 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01378
01379 return (global_index != i.global_index);
01380 }
01381
01382
01383
01384 template <class BlockVectorType, bool constness>
01385 inline
01386 bool
01387 Iterator<BlockVectorType,constness>::
01388 operator != (const InverseConstnessIterator &i) const
01389 {
01390 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01391
01392 return (global_index != i.global_index);
01393 }
01394
01395
01396
01397 template <class BlockVectorType, bool constness>
01398 inline
01399 bool
01400 Iterator<BlockVectorType,constness>::
01401 operator < (const Iterator &i) const
01402 {
01403 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01404
01405 return (global_index < i.global_index);
01406 }
01407
01408
01409
01410 template <class BlockVectorType, bool constness>
01411 inline
01412 bool
01413 Iterator<BlockVectorType,constness>::
01414 operator < (const InverseConstnessIterator &i) const
01415 {
01416 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01417
01418 return (global_index < i.global_index);
01419 }
01420
01421
01422
01423 template <class BlockVectorType, bool constness>
01424 inline
01425 bool
01426 Iterator<BlockVectorType,constness>::
01427 operator <= (const Iterator &i) const
01428 {
01429 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01430
01431 return (global_index <= i.global_index);
01432 }
01433
01434
01435
01436 template <class BlockVectorType, bool constness>
01437 inline
01438 bool
01439 Iterator<BlockVectorType,constness>::
01440 operator <= (const InverseConstnessIterator &i) const
01441 {
01442 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01443
01444 return (global_index <= i.global_index);
01445 }
01446
01447
01448
01449 template <class BlockVectorType, bool constness>
01450 inline
01451 bool
01452 Iterator<BlockVectorType,constness>::
01453 operator > (const Iterator &i) const
01454 {
01455 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01456
01457 return (global_index > i.global_index);
01458 }
01459
01460
01461
01462 template <class BlockVectorType, bool constness>
01463 inline
01464 bool
01465 Iterator<BlockVectorType,constness>::
01466 operator > (const InverseConstnessIterator &i) const
01467 {
01468 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01469
01470 return (global_index > i.global_index);
01471 }
01472
01473
01474
01475 template <class BlockVectorType, bool constness>
01476 inline
01477 bool
01478 Iterator<BlockVectorType,constness>::
01479 operator >= (const Iterator &i) const
01480 {
01481 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01482
01483 return (global_index >= i.global_index);
01484 }
01485
01486
01487
01488 template <class BlockVectorType, bool constness>
01489 inline
01490 bool
01491 Iterator<BlockVectorType,constness>::
01492 operator >= (const InverseConstnessIterator &i) const
01493 {
01494 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01495
01496 return (global_index >= i.global_index);
01497 }
01498
01499
01500
01501 template <class BlockVectorType, bool constness>
01502 inline
01503 typename Iterator<BlockVectorType,constness>::difference_type
01504 Iterator<BlockVectorType,constness>::
01505 operator - (const Iterator &i) const
01506 {
01507 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01508
01509 return (static_cast<signed int>(global_index) -
01510 static_cast<signed int>(i.global_index));
01511 }
01512
01513
01514
01515 template <class BlockVectorType, bool constness>
01516 inline
01517 typename Iterator<BlockVectorType,constness>::difference_type
01518 Iterator<BlockVectorType,constness>::
01519 operator - (const InverseConstnessIterator &i) const
01520 {
01521 Assert (parent == i.parent, ExcPointerToDifferentVectors());
01522
01523 return (static_cast<signed int>(global_index) -
01524 static_cast<signed int>(i.global_index));
01525 }
01526
01527
01528
01529 template <class BlockVectorType, bool constness>
01530 inline
01531 Iterator<BlockVectorType,constness>
01532 Iterator<BlockVectorType,constness>::
01533 operator + (const difference_type &d) const
01534 {
01535
01536
01537
01538
01539
01540 if ((global_index+d >= next_break_backward) &&
01541 (global_index+d <= next_break_forward))
01542 return Iterator (*parent, global_index+d, current_block,
01543 index_within_block+d,
01544 next_break_forward, next_break_backward);
01545 else
01546
01547
01548
01549 return Iterator (*parent, global_index+d);
01550 }
01551
01552
01553
01554 template <class BlockVectorType, bool constness>
01555 inline
01556 Iterator<BlockVectorType,constness>
01557 Iterator<BlockVectorType,constness>::
01558 operator - (const difference_type &d) const
01559 {
01560
01561
01562
01563
01564
01565 if ((global_index-d >= next_break_backward) &&
01566 (global_index-d <= next_break_forward))
01567 return Iterator (*parent, global_index-d, current_block,
01568 index_within_block-d,
01569 next_break_forward, next_break_backward);
01570 else
01571
01572
01573
01574 return Iterator (*parent, global_index-d);
01575 }
01576
01577
01578
01579 template <class BlockVectorType, bool constness>
01580 inline
01581 Iterator<BlockVectorType,constness> &
01582 Iterator<BlockVectorType,constness>::
01583 operator += (const difference_type &d)
01584 {
01585
01586
01587
01588
01589
01590 if ((global_index+d >= next_break_backward) &&
01591 (global_index+d <= next_break_forward))
01592 {
01593 global_index += d;
01594 index_within_block += d;
01595 }
01596 else
01597
01598
01599
01600 *this = Iterator (*parent, global_index+d);
01601
01602 return *this;
01603 }
01604
01605
01606
01607 template <class BlockVectorType, bool constness>
01608 inline
01609 Iterator<BlockVectorType,constness> &
01610 Iterator<BlockVectorType,constness>::
01611 operator -= (const difference_type &d)
01612 {
01613
01614
01615
01616
01617
01618 if ((global_index-d >= next_break_backward) &&
01619 (global_index-d <= next_break_forward))
01620 {
01621 global_index -= d;
01622 index_within_block -= d;
01623 }
01624 else
01625
01626
01627
01628 *this = Iterator (*parent, global_index-d);
01629
01630 return *this;
01631 }
01632
01633
01634 template <class BlockVectorType, bool constness>
01635 Iterator<BlockVectorType,constness>::
01636 Iterator (BlockVector &parent,
01637 const unsigned global_index)
01638 :
01639 parent (&parent),
01640 global_index (global_index)
01641 {
01642
01643
01644
01645
01646
01647
01648 if (global_index < parent.size())
01649 {
01650 const std::pair<unsigned int, unsigned int>
01651 indices = parent.block_indices.global_to_local(global_index);
01652 current_block = indices.first;
01653 index_within_block = indices.second;
01654
01655 next_break_backward
01656 = parent.block_indices.local_to_global (current_block, 0);
01657 next_break_forward
01658 = (parent.block_indices.local_to_global (current_block, 0)
01659 +parent.block_indices.block_size(current_block)-1);
01660 }
01661 else
01662
01663
01664 {
01665 this->global_index = parent.size ();
01666 current_block = parent.n_blocks();
01667 index_within_block = 0;
01668 next_break_backward = global_index;
01669 next_break_forward = numbers::invalid_unsigned_int;
01670 };
01671 }
01672
01673
01674
01675 template <class BlockVectorType, bool constness>
01676 void
01677 Iterator<BlockVectorType,constness>::move_forward ()
01678 {
01679 if (global_index != next_break_forward)
01680 ++index_within_block;
01681 else
01682 {
01683
01684
01685 index_within_block = 0;
01686 ++current_block;
01687
01688
01689
01690 next_break_backward = next_break_forward+1;
01691
01692
01693 if (current_block < parent->block_indices.size())
01694 next_break_forward
01695 += parent->block_indices.block_size(current_block);
01696 else
01697
01698
01699
01700
01701 next_break_forward = numbers::invalid_unsigned_int;
01702 };
01703
01704 ++global_index;
01705 }
01706
01707
01708
01709 template <class BlockVectorType, bool constness>
01710 void
01711 Iterator<BlockVectorType,constness>::move_backward ()
01712 {
01713 if (global_index != next_break_backward)
01714 --index_within_block;
01715 else
01716 if (current_block != 0)
01717 {
01718
01719
01720 --current_block;
01721 index_within_block = parent->block_indices.block_size(current_block)-1;
01722
01723
01724
01725 next_break_forward = next_break_backward-1;
01726
01727
01728 next_break_backward
01729 -= parent->block_indices.block_size (current_block);
01730 }
01731 else
01732
01733
01734 {
01735 --current_block;
01736 index_within_block = numbers::invalid_unsigned_int;
01737 next_break_forward = 0;
01738 next_break_backward = 0;
01739 };
01740
01741 --global_index;
01742 }
01743
01744
01745 }
01746
01747 }
01748
01749
01750 template <class VectorType>
01751 inline
01752 BlockVectorBase<VectorType>::BlockVectorBase ()
01753 {}
01754
01755
01756
01757 template <class VectorType>
01758 inline
01759 unsigned int
01760 BlockVectorBase<VectorType>::size () const
01761 {
01762 return block_indices.total_size();
01763 }
01764
01765
01766
01767 template <class VectorType>
01768 inline
01769 unsigned int
01770 BlockVectorBase<VectorType>::n_blocks () const
01771 {
01772 return block_indices.size();
01773 }
01774
01775
01776 template <class VectorType>
01777 inline
01778 typename BlockVectorBase<VectorType>::BlockType &
01779 BlockVectorBase<VectorType>::block (const unsigned int i)
01780 {
01781 Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
01782
01783 return components[i];
01784 }
01785
01786
01787
01788 template <class VectorType>
01789 inline
01790 const typename BlockVectorBase<VectorType>::BlockType &
01791 BlockVectorBase<VectorType>::block (const unsigned int i) const
01792 {
01793 Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
01794
01795 return components[i];
01796 }
01797
01798
01799
01800 template <class VectorType>
01801 inline
01802 const BlockIndices&
01803 BlockVectorBase<VectorType>::get_block_indices () const
01804 {
01805 return block_indices;
01806 }
01807
01808
01809 template <class VectorType>
01810 inline
01811 void
01812 BlockVectorBase<VectorType>::collect_sizes ()
01813 {
01814 std::vector<unsigned int> sizes (n_blocks());
01815
01816 for (unsigned int i=0; i<n_blocks(); ++i)
01817 sizes[i] = block(i).size();
01818
01819 block_indices.reinit(sizes);
01820 }
01821
01822
01823
01824 template <class VectorType>
01825 inline
01826 void
01827 BlockVectorBase<VectorType>::compress ()
01828 {
01829 for (unsigned int i=0; i<n_blocks(); ++i)
01830 block(i).compress ();
01831 }
01832
01833
01834
01835 template <class VectorType>
01836 inline
01837 typename BlockVectorBase<VectorType>::iterator
01838 BlockVectorBase<VectorType>::begin()
01839 {
01840 return iterator(*this, 0U);
01841 }
01842
01843
01844
01845 template <class VectorType>
01846 inline
01847 typename BlockVectorBase<VectorType>::const_iterator
01848 BlockVectorBase<VectorType>::begin() const
01849 {
01850 return const_iterator(*this, 0U);
01851 }
01852
01853
01854 template <class VectorType>
01855 inline
01856 typename BlockVectorBase<VectorType>::iterator
01857 BlockVectorBase<VectorType>::end()
01858 {
01859 return iterator(*this, size());
01860 }
01861
01862
01863
01864 template <class VectorType>
01865 inline
01866 typename BlockVectorBase<VectorType>::const_iterator
01867 BlockVectorBase<VectorType>::end() const
01868 {
01869 return const_iterator(*this, size());
01870 }
01871
01872
01873
01874 template <class VectorType>
01875 bool
01876 BlockVectorBase<VectorType>::all_zero () const
01877 {
01878 for (unsigned int i=0;i<n_blocks();++i)
01879 if (components[i].all_zero() == false)
01880 return false;
01881
01882 return true;
01883 }
01884
01885
01886
01887 template <class VectorType>
01888 bool
01889 BlockVectorBase<VectorType>::is_non_negative () const
01890 {
01891 for (unsigned int i=0;i<n_blocks();++i)
01892 if (components[i].is_non_negative() == false)
01893 return false;
01894
01895 return true;
01896 }
01897
01898
01899
01900 template <class VectorType>
01901 typename BlockVectorBase<VectorType>::value_type
01902 BlockVectorBase<VectorType>::
01903 operator * (const BlockVectorBase<VectorType>& v) const
01904 {
01905 Assert (n_blocks() == v.n_blocks(),
01906 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
01907
01908 value_type sum = 0.;
01909 for (unsigned int i=0;i<n_blocks();++i)
01910 sum += components[i]*v.components[i];
01911
01912 return sum;
01913 }
01914
01915
01916 template <class VectorType>
01917 typename BlockVectorBase<VectorType>::real_type
01918 BlockVectorBase<VectorType>::norm_sqr () const
01919 {
01920 real_type sum = 0.;
01921 for (unsigned int i=0;i<n_blocks();++i)
01922 sum += components[i].norm_sqr();
01923
01924 return sum;
01925 }
01926
01927
01928
01929 template <class VectorType>
01930 typename BlockVectorBase<VectorType>::value_type
01931 BlockVectorBase<VectorType>::mean_value () const
01932 {
01933 value_type sum = 0.;
01934 for (unsigned int i=0;i<n_blocks();++i)
01935 sum += components[i].mean_value() * components[i].size();
01936
01937 return sum/size();
01938 }
01939
01940
01941
01942 template <class VectorType>
01943 typename BlockVectorBase<VectorType>::real_type
01944 BlockVectorBase<VectorType>::l1_norm () const
01945 {
01946 real_type sum = 0.;
01947 for (unsigned int i=0;i<n_blocks();++i)
01948 sum += components[i].l1_norm();
01949
01950 return sum;
01951 }
01952
01953
01954
01955 template <class VectorType>
01956 typename BlockVectorBase<VectorType>::real_type
01957 BlockVectorBase<VectorType>::l2_norm () const
01958 {
01959 return std::sqrt(norm_sqr());
01960 }
01961
01962
01963
01964 template <class VectorType>
01965 typename BlockVectorBase<VectorType>::real_type
01966 BlockVectorBase<VectorType>::linfty_norm () const
01967 {
01968 real_type sum = 0.;
01969 for (unsigned int i=0;i<n_blocks();++i)
01970 {
01971 value_type newval = components[i].linfty_norm();
01972 if (sum<newval)
01973 sum = newval;
01974 }
01975 return sum;
01976 }
01977
01978
01979
01980 template <class VectorType>
01981 BlockVectorBase<VectorType> &
01982 BlockVectorBase<VectorType>::operator += (const BlockVectorBase<VectorType>& v)
01983 {
01984 add (v);
01985 return *this;
01986 }
01987
01988
01989
01990 template <class VectorType>
01991 BlockVectorBase<VectorType> &
01992 BlockVectorBase<VectorType>::operator -= (const BlockVectorBase<VectorType>& v)
01993 {
01994 Assert (n_blocks() == v.n_blocks(),
01995 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
01996
01997 for (unsigned int i=0;i<n_blocks();++i)
01998 {
01999 components[i] -= v.components[i];
02000 }
02001 return *this;
02002 }
02003
02004
02005
02006 template <class VectorType>
02007 template <typename Number>
02008 inline
02009 void
02010 BlockVectorBase<VectorType>::add (const std::vector<unsigned int> &indices,
02011 const std::vector<Number> &values)
02012 {
02013 Assert (indices.size() == values.size(),
02014 ExcDimensionMismatch(indices.size(), values.size()));
02015 add (indices.size(), &indices[0], &values[0]);
02016 }
02017
02018
02019
02020 template <class VectorType>
02021 template <typename Number>
02022 inline
02023 void
02024 BlockVectorBase<VectorType>::add (const std::vector<unsigned int> &indices,
02025 const Vector<Number> &values)
02026 {
02027 Assert (indices.size() == values.size(),
02028 ExcDimensionMismatch(indices.size(), values.size()));
02029 const unsigned int n_indices = indices.size();
02030 for (unsigned int i=0; i<n_indices; ++i)
02031 (*this)(indices[i]) += values(i);
02032 }
02033
02034
02035
02036 template <class VectorType>
02037 template <typename Number>
02038 inline
02039 void
02040 BlockVectorBase<VectorType>::add (const unsigned int n_indices,
02041 const unsigned int *indices,
02042 const Number *values)
02043 {
02044 for (unsigned int i=0; i<n_indices; ++i)
02045 (*this)(indices[i]) += values[i];
02046 }
02047
02048
02049
02050 template <class VectorType>
02051 void BlockVectorBase<VectorType>::add (const value_type a)
02052 {
02053 Assert (numbers::is_finite(a), ExcNumberNotFinite());
02054
02055 for (unsigned int i=0;i<n_blocks();++i)
02056 {
02057 components[i].add(a);
02058 }
02059 }
02060
02061
02062
02063 template <class VectorType>
02064 void BlockVectorBase<VectorType>::add (const BlockVectorBase<VectorType>& v)
02065 {
02066 Assert (n_blocks() == v.n_blocks(),
02067 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02068
02069 for (unsigned int i=0;i<n_blocks();++i)
02070 {
02071 components[i].add(v.components[i]);
02072 }
02073 }
02074
02075
02076
02077 template <class VectorType>
02078 void BlockVectorBase<VectorType>::add (const value_type a,
02079 const BlockVectorBase<VectorType>& v)
02080 {
02081
02082 Assert (numbers::is_finite(a), ExcNumberNotFinite());
02083
02084 Assert (n_blocks() == v.n_blocks(),
02085 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02086
02087 for (unsigned int i=0;i<n_blocks();++i)
02088 {
02089 components[i].add(a, v.components[i]);
02090 }
02091 }
02092
02093
02094
02095 template <class VectorType>
02096 void BlockVectorBase<VectorType>::add (const value_type a,
02097 const BlockVectorBase<VectorType>& v,
02098 const value_type b,
02099 const BlockVectorBase<VectorType>& w)
02100 {
02101
02102 Assert (numbers::is_finite(a), ExcNumberNotFinite());
02103 Assert (numbers::is_finite(b), ExcNumberNotFinite());
02104
02105 Assert (n_blocks() == v.n_blocks(),
02106 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02107 Assert (n_blocks() == w.n_blocks(),
02108 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
02109
02110
02111 for (unsigned int i=0;i<n_blocks();++i)
02112 {
02113 components[i].add(a, v.components[i], b, w.components[i]);
02114 }
02115 }
02116
02117
02118
02119 template <class VectorType>
02120 void BlockVectorBase<VectorType>::sadd (const value_type x,
02121 const BlockVectorBase<VectorType>& v)
02122 {
02123
02124 Assert (numbers::is_finite(x), ExcNumberNotFinite());
02125
02126 Assert (n_blocks() == v.n_blocks(),
02127 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02128
02129 for (unsigned int i=0;i<n_blocks();++i)
02130 {
02131 components[i].sadd(x, v.components[i]);
02132 }
02133 }
02134
02135
02136
02137 template <class VectorType>
02138 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
02139 const BlockVectorBase<VectorType>& v)
02140 {
02141
02142 Assert (numbers::is_finite(x), ExcNumberNotFinite());
02143 Assert (numbers::is_finite(a), ExcNumberNotFinite());
02144
02145 Assert (n_blocks() == v.n_blocks(),
02146 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02147
02148 for (unsigned int i=0;i<n_blocks();++i)
02149 {
02150 components[i].sadd(x, a, v.components[i]);
02151 }
02152 }
02153
02154
02155
02156 template <class VectorType>
02157 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
02158 const BlockVectorBase<VectorType>& v,
02159 const value_type b,
02160 const BlockVectorBase<VectorType>& w)
02161 {
02162
02163 Assert (numbers::is_finite(x), ExcNumberNotFinite());
02164 Assert (numbers::is_finite(a), ExcNumberNotFinite());
02165 Assert (numbers::is_finite(b), ExcNumberNotFinite());
02166
02167 Assert (n_blocks() == v.n_blocks(),
02168 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02169 Assert (n_blocks() == w.n_blocks(),
02170 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
02171
02172 for (unsigned int i=0;i<n_blocks();++i)
02173 {
02174 components[i].sadd(x, a, v.components[i], b, w.components[i]);
02175 }
02176 }
02177
02178
02179
02180 template <class VectorType>
02181 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
02182 const BlockVectorBase<VectorType>& v,
02183 const value_type b,
02184 const BlockVectorBase<VectorType>& w,
02185 const value_type c,
02186 const BlockVectorBase<VectorType>& y)
02187 {
02188
02189 Assert (numbers::is_finite(x), ExcNumberNotFinite());
02190 Assert (numbers::is_finite(a), ExcNumberNotFinite());
02191 Assert (numbers::is_finite(b), ExcNumberNotFinite());
02192 Assert (numbers::is_finite(c), ExcNumberNotFinite());
02193
02194 Assert (n_blocks() == v.n_blocks(),
02195 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02196 Assert (n_blocks() == w.n_blocks(),
02197 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
02198 Assert (n_blocks() == y.n_blocks(),
02199 ExcDimensionMismatch(n_blocks(), y.n_blocks()));
02200
02201 for (unsigned int i=0;i<n_blocks();++i)
02202 {
02203 components[i].sadd(x, a, v.components[i],
02204 b, w.components[i], c, y.components[i]);
02205 }
02206 }
02207
02208
02209
02210 template <class VectorType>
02211 template <class BlockVector2>
02212 void BlockVectorBase<VectorType>::scale (const BlockVector2 &v)
02213 {
02214 Assert (n_blocks() == v.n_blocks(),
02215 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02216 for (unsigned int i=0;i<n_blocks();++i)
02217 components[i].scale(v.block(i));
02218 }
02219
02220
02221
02222 template <class VectorType>
02223 void BlockVectorBase<VectorType>::equ (const value_type a,
02224 const BlockVectorBase<VectorType>& v,
02225 const value_type b,
02226 const BlockVectorBase<VectorType>& w)
02227 {
02228
02229 Assert (numbers::is_finite(a), ExcNumberNotFinite());
02230 Assert (numbers::is_finite(b), ExcNumberNotFinite());
02231
02232 Assert (n_blocks() == v.n_blocks(),
02233 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02234 Assert (n_blocks() == w.n_blocks(),
02235 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
02236
02237 for (unsigned int i=0;i<n_blocks();++i)
02238 {
02239 components[i].equ( a, v.components[i], b, w.components[i]);
02240 }
02241 }
02242
02243
02244
02245 template <class VectorType>
02246 std::size_t
02247 BlockVectorBase<VectorType>::memory_consumption () const
02248 {
02249 std::size_t mem = sizeof(this->n_blocks());
02250 for (unsigned int i=0; i<this->components.size(); ++i)
02251 mem += MemoryConsumption::memory_consumption (this->components[i]);
02252 mem += MemoryConsumption::memory_consumption (this->block_indices);
02253 return mem;
02254 }
02255
02256
02257
02258 template <class VectorType>
02259 template <class BlockVector2>
02260 void BlockVectorBase<VectorType>::equ (const value_type a,
02261 const BlockVector2 &v)
02262 {
02263
02264 Assert (numbers::is_finite(a), ExcNumberNotFinite());
02265
02266 Assert (n_blocks() == v.n_blocks(),
02267 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
02268
02269 for (unsigned int i=0;i<n_blocks();++i)
02270 components[i].equ( a, v.components[i]);
02271 }
02272
02273
02274
02275 template <class VectorType>
02276 void BlockVectorBase<VectorType>::update_ghost_values () const
02277 {
02278 for (unsigned int i=0; i<n_blocks(); ++i)
02279 block(i).update_ghost_values ();
02280 }
02281
02282
02283
02284 template <class VectorType>
02285 BlockVectorBase<VectorType>&
02286 BlockVectorBase<VectorType>::operator = (const value_type s)
02287 {
02288
02289 Assert (numbers::is_finite(s), ExcNumberNotFinite());
02290
02291 for (unsigned int i=0;i<n_blocks();++i)
02292 components[i] = s;
02293
02294 return *this;
02295 }
02296
02297
02298 template <class VectorType>
02299 BlockVectorBase<VectorType>&
02300 BlockVectorBase<VectorType>::operator = (const BlockVectorBase<VectorType>& v)
02301 {
02302 AssertDimension(n_blocks(), v.n_blocks());
02303
02304 for (unsigned int i=0;i<n_blocks();++i)
02305 components[i] = v.components[i];
02306
02307 return *this;
02308 }
02309
02310
02311 template <class VectorType>
02312 template <class VectorType2>
02313 BlockVectorBase<VectorType>&
02314 BlockVectorBase<VectorType>::operator = (const BlockVectorBase<VectorType2> &v)
02315 {
02316 AssertDimension(n_blocks(), v.n_blocks());
02317
02318 for (unsigned int i=0;i<n_blocks();++i)
02319 components[i] = v.components[i];
02320
02321 return *this;
02322 }
02323
02324
02325
02326 template <class VectorType>
02327 BlockVectorBase<VectorType>&
02328 BlockVectorBase<VectorType>::operator = (const VectorType &v)
02329 {
02330 Assert (size() == v.size(),
02331 ExcDimensionMismatch(size(), v.size()));
02332
02333 unsigned int index_v = 0;
02334 for (unsigned int b=0;b<n_blocks();++b)
02335 for (unsigned int i=0; i<block(b).size(); ++i, ++index_v)
02336 block(b)(i) = v(index_v);
02337
02338 return *this;
02339 }
02340
02341
02342
02343 template <class VectorType>
02344 template <class VectorType2>
02345 inline
02346 bool
02347 BlockVectorBase<VectorType>::
02348 operator == (const BlockVectorBase<VectorType2> &v) const
02349 {
02350 Assert (block_indices == v.block_indices, ExcDifferentBlockIndices());
02351
02352 for (unsigned int i=0;i<n_blocks();++i)
02353 if ( ! (components[i] == v.components[i]))
02354 return false;
02355
02356 return true;
02357 }
02358
02359
02360
02361 template <class VectorType>
02362 inline
02363 BlockVectorBase<VectorType> &
02364 BlockVectorBase<VectorType>::operator *= (const value_type factor)
02365 {
02366
02367 Assert (numbers::is_finite(factor), ExcNumberNotFinite());
02368
02369 for (unsigned int i=0;i<n_blocks();++i)
02370 components[i] *= factor;
02371
02372 return *this;
02373 }
02374
02375
02376
02377 template <class VectorType>
02378 inline
02379 BlockVectorBase<VectorType> &
02380 BlockVectorBase<VectorType>::operator /= (const value_type factor)
02381 {
02382
02383 Assert (numbers::is_finite(factor), ExcNumberNotFinite());
02384 Assert (factor > 0., ExcDivideByZero() );
02385
02386 for (unsigned int i=0;i<n_blocks();++i)
02387 components[i] /= factor;
02388
02389 return *this;
02390 }
02391
02392
02393 template <class VectorType>
02394 inline
02395 typename BlockVectorBase<VectorType>::value_type
02396 BlockVectorBase<VectorType>::operator() (const unsigned int i) const
02397 {
02398 const std::pair<unsigned int,unsigned int> local_index
02399 = block_indices.global_to_local (i);
02400 return components[local_index.first](local_index.second);
02401 }
02402
02403
02404
02405 template <class VectorType>
02406 inline
02407 typename BlockVectorBase<VectorType>::reference
02408 BlockVectorBase<VectorType>::operator() (const unsigned int i)
02409 {
02410 const std::pair<unsigned int,unsigned int> local_index
02411 = block_indices.global_to_local (i);
02412 return components[local_index.first](local_index.second);
02413 }
02414
02415
02416
02417 template <class VectorType>
02418 inline
02419 typename BlockVectorBase<VectorType>::value_type
02420 BlockVectorBase<VectorType>::operator[] (const unsigned int i) const
02421 {
02422 return operator()(i);
02423 }
02424
02425
02426
02427 template <class VectorType>
02428 inline
02429 typename BlockVectorBase<VectorType>::reference
02430 BlockVectorBase<VectorType>::operator[] (const unsigned int i)
02431 {
02432 return operator()(i);
02433 }
02434
02435 #endif // DOXYGEN
02436
02437 DEAL_II_NAMESPACE_CLOSE
02438
02439 #endif
02440