00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__block_matrix_base_h
00013 #define __deal2__block_matrix_base_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/table.h>
00018 #include <deal.II/base/utilities.h>
00019 #include <deal.II/base/smartpointer.h>
00020 #include <deal.II/base/memory_consumption.h>
00021 #include <deal.II/lac/block_indices.h>
00022 #include <deal.II/lac/exceptions.h>
00023 #include <deal.II/lac/full_matrix.h>
00024 #include <deal.II/lac/matrix_iterator.h>
00025
00026 #include <cmath>
00027
00028 DEAL_II_NAMESPACE_OPEN
00029
00030
00031 template <typename> class MatrixIterator;
00032
00033
00034
00044 namespace BlockMatrixIterators
00045 {
00051 template <class BlockMatrix>
00052 class AccessorBase
00053 {
00054 public:
00059 typedef typename BlockMatrix::value_type value_type;
00060
00065 AccessorBase ();
00066
00072 unsigned int block_row() const;
00073
00079 unsigned int block_column() const;
00080
00081 protected:
00086 unsigned int row_block;
00087
00092 unsigned int col_block;
00093
00098 template <typename>
00099 friend class MatrixIterator;
00100 };
00105 template <class BlockMatrix, bool ConstNess>
00106 class Accessor;
00107
00112 template <class BlockMatrix>
00113 class Accessor<BlockMatrix, false>
00114 :
00115 public AccessorBase<BlockMatrix>
00116 {
00117 public:
00122 typedef BlockMatrix MatrixType;
00123
00128 typedef typename BlockMatrix::value_type value_type;
00129
00142 Accessor (BlockMatrix* m,
00143 const unsigned int row,
00144 const unsigned int col);
00145
00151 unsigned int row() const;
00152
00158 unsigned int column() const;
00159
00164 value_type value() const;
00165
00169 void set_value(value_type newval) const;
00170
00171 protected:
00175 BlockMatrix * matrix;
00176
00181 typename BlockMatrix::BlockType::iterator base_iterator;
00182
00186 void advance ();
00187
00192 bool operator == (const Accessor &a) const;
00193
00194 template <typename> friend class MatrixIterator;
00195 friend class Accessor<BlockMatrix, true>;
00196 };
00197
00203 template <class BlockMatrix>
00204 class Accessor<BlockMatrix, true>
00205 :
00206 public AccessorBase<BlockMatrix>
00207 {
00208 public:
00213 typedef const BlockMatrix MatrixType;
00214
00219 typedef typename BlockMatrix::value_type value_type;
00220
00233 Accessor (const BlockMatrix* m,
00234 const unsigned int row,
00235 const unsigned int col);
00236
00241 Accessor(const Accessor<BlockMatrix, false>&);
00242
00248 unsigned int row() const;
00249
00255 unsigned int column() const;
00256
00261 value_type value() const;
00262 protected:
00266 const BlockMatrix * matrix;
00267
00272 typename BlockMatrix::BlockType::const_iterator base_iterator;
00273
00277 void advance ();
00278
00283 bool operator == (const Accessor &a) const;
00284
00289 template <typename>
00290 friend class ::MatrixIterator;
00291 };
00292 }
00293
00294
00295
00355 template <typename MatrixType>
00356 class BlockMatrixBase : public Subscriptor
00357 {
00358 public:
00363 typedef MatrixType BlockType;
00364
00369 typedef typename BlockType::value_type value_type;
00370 typedef value_type *pointer;
00371 typedef const value_type *const_pointer;
00372 typedef value_type &reference;
00373 typedef const value_type &const_reference;
00374 typedef std::size_t size_type;
00375
00376 typedef
00377 MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, false> >
00378 iterator;
00379
00380 typedef
00381 MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, true> >
00382 const_iterator;
00383
00384
00388 BlockMatrixBase ();
00389
00417 template <class BlockMatrixType>
00418 BlockMatrixBase &
00419 copy_from (const BlockMatrixType &source);
00420
00425 BlockType &
00426 block (const unsigned int row,
00427 const unsigned int column);
00428
00429
00435 const BlockType &
00436 block (const unsigned int row,
00437 const unsigned int column) const;
00438
00445 unsigned int m () const;
00446
00453 unsigned int n () const;
00454
00455
00462 unsigned int n_block_rows () const;
00463
00470 unsigned int n_block_cols () const;
00471
00481 void set (const unsigned int i,
00482 const unsigned int j,
00483 const value_type value);
00484
00512 template <typename number>
00513 void set (const std::vector<unsigned int> &indices,
00514 const FullMatrix<number> &full_matrix,
00515 const bool elide_zero_values = false);
00516
00524 template <typename number>
00525 void set (const std::vector<unsigned int> &row_indices,
00526 const std::vector<unsigned int> &col_indices,
00527 const FullMatrix<number> &full_matrix,
00528 const bool elide_zero_values = false);
00529
00548 template <typename number>
00549 void set (const unsigned int row,
00550 const std::vector<unsigned int> &col_indices,
00551 const std::vector<number> &values,
00552 const bool elide_zero_values = false);
00553
00570 template <typename number>
00571 void set (const unsigned int row,
00572 const unsigned int n_cols,
00573 const unsigned int *col_indices,
00574 const number *values,
00575 const bool elide_zero_values = false);
00576
00586 void add (const unsigned int i,
00587 const unsigned int j,
00588 const value_type value);
00589
00617 template <typename number>
00618 void add (const std::vector<unsigned int> &indices,
00619 const FullMatrix<number> &full_matrix,
00620 const bool elide_zero_values = true);
00621
00629 template <typename number>
00630 void add (const std::vector<unsigned int> &row_indices,
00631 const std::vector<unsigned int> &col_indices,
00632 const FullMatrix<number> &full_matrix,
00633 const bool elide_zero_values = true);
00634
00652 template <typename number>
00653 void add (const unsigned int row,
00654 const std::vector<unsigned int> &col_indices,
00655 const std::vector<number> &values,
00656 const bool elide_zero_values = true);
00657
00675 template <typename number>
00676 void add (const unsigned int row,
00677 const unsigned int n_cols,
00678 const unsigned int *col_indices,
00679 const number *values,
00680 const bool elide_zero_values = true,
00681 const bool col_indices_are_sorted = false);
00682
00694 value_type operator () (const unsigned int i,
00695 const unsigned int j) const;
00696
00714 value_type el (const unsigned int i,
00715 const unsigned int j) const;
00716
00734 value_type diag_element (const unsigned int i) const;
00735
00740 void compress ();
00741
00746 BlockMatrixBase & operator *= (const value_type factor);
00747
00752 BlockMatrixBase & operator /= (const value_type factor);
00753
00771 template <class BlockMatrixType>
00772 void add (const value_type factor,
00773 const BlockMatrixType &matrix);
00774
00775
00782 template <class BlockVectorType>
00783 void vmult_add (BlockVectorType &dst,
00784 const BlockVectorType &src) const;
00785
00796 template <class BlockVectorType>
00797 void Tvmult_add (BlockVectorType &dst,
00798 const BlockVectorType &src) const;
00799
00824 template <class BlockVectorType>
00825 value_type
00826 matrix_norm_square (const BlockVectorType &v) const;
00827
00832 template <class BlockVectorType>
00833 value_type
00834 matrix_scalar_product (const BlockVectorType &u,
00835 const BlockVectorType &v) const;
00836
00842 template <class BlockVectorType>
00843 value_type residual (BlockVectorType &dst,
00844 const BlockVectorType &x,
00845 const BlockVectorType &b) const;
00846
00851 iterator begin ();
00852
00856 iterator end ();
00857
00862 iterator begin (const unsigned int r);
00863
00867 iterator end (const unsigned int r);
00872 const_iterator begin () const;
00873
00877 const_iterator end () const;
00878
00883 const_iterator begin (const unsigned int r) const;
00884
00888 const_iterator end (const unsigned int r) const;
00889
00894 const BlockIndices & get_row_indices () const;
00895
00900 const BlockIndices & get_column_indices () const;
00901
00910 std::size_t memory_consumption () const;
00911
00918 DeclException4 (ExcIncompatibleRowNumbers,
00919 int, int, int, int,
00920 << "The blocks [" << arg1 << ',' << arg2 << "] and ["
00921 << arg3 << ',' << arg4 << "] have differing row numbers.");
00925 DeclException4 (ExcIncompatibleColNumbers,
00926 int, int, int, int,
00927 << "The blocks [" << arg1 << ',' << arg2 << "] and ["
00928 << arg3 << ',' << arg4 << "] have differing column numbers.");
00930 protected:
00950 void clear ();
00951
00955 BlockIndices row_block_indices;
00956 BlockIndices column_block_indices;
00957
00961 Table<2,SmartPointer<BlockType, BlockMatrixBase<MatrixType> > > sub_objects;
00962
00998 void collect_sizes ();
00999
01016 template <class BlockVectorType>
01017 void vmult_block_block (BlockVectorType &dst,
01018 const BlockVectorType &src) const;
01019
01038 template <class BlockVectorType,
01039 class VectorType>
01040 void vmult_block_nonblock (BlockVectorType &dst,
01041 const VectorType &src) const;
01042
01061 template <class BlockVectorType,
01062 class VectorType>
01063 void vmult_nonblock_block (VectorType &dst,
01064 const BlockVectorType &src) const;
01065
01084 template <class VectorType>
01085 void vmult_nonblock_nonblock (VectorType &dst,
01086 const VectorType &src) const;
01087
01107 template <class BlockVectorType>
01108 void Tvmult_block_block (BlockVectorType &dst,
01109 const BlockVectorType &src) const;
01110
01129 template <class BlockVectorType,
01130 class VectorType>
01131 void Tvmult_block_nonblock (BlockVectorType &dst,
01132 const VectorType &src) const;
01133
01152 template <class BlockVectorType,
01153 class VectorType>
01154 void Tvmult_nonblock_block (VectorType &dst,
01155 const BlockVectorType &src) const;
01156
01175 template <class VectorType>
01176 void Tvmult_nonblock_nonblock (VectorType &dst,
01177 const VectorType &src) const;
01178
01179
01180 protected:
01181
01192 void prepare_add_operation();
01193
01199 void prepare_set_operation();
01200
01201
01202 private:
01209 std::vector<unsigned int> counter_within_block;
01210
01217 std::vector<std::vector<unsigned int> > column_indices;
01218
01225 std::vector<std::vector<double> > column_values;
01226
01227
01233 #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG
01234 template <typename, bool>
01235 friend class BlockMatrixIterators::Accessor;
01236
01237 template <typename>
01238 friend class MatrixIterator;
01239 #else
01240 typedef BlockMatrixIterators::Accessor<BlockMatrixBase, true> ConstAccessor;
01241 typedef BlockMatrixIterators::Accessor<BlockMatrixBase, false> Accessor;
01242 friend class ConstAccessor;
01243
01244 friend class const_iterator;
01245 #endif
01246
01247 };
01248
01249
01252 #ifndef DOXYGEN
01253
01254
01255
01256 namespace BlockMatrixIterators
01257 {
01258 template <class BlockMatrix>
01259 inline
01260 AccessorBase<BlockMatrix>::AccessorBase()
01261 :
01262 row_block(0),
01263 col_block(0)
01264 {}
01265
01266
01267 template <class BlockMatrix>
01268 inline
01269 unsigned int
01270 AccessorBase<BlockMatrix>::block_row() const
01271 {
01272 Assert (row_block != numbers::invalid_unsigned_int,
01273 ExcIteratorPastEnd());
01274
01275 return row_block;
01276 }
01277
01278
01279 template <class BlockMatrix>
01280 inline
01281 unsigned int
01282 AccessorBase<BlockMatrix>::block_column() const
01283 {
01284 Assert (col_block != numbers::invalid_unsigned_int,
01285 ExcIteratorPastEnd());
01286
01287 return col_block;
01288 }
01289
01290
01291 template <class BlockMatrix>
01292 inline
01293 Accessor<BlockMatrix, true>::Accessor (
01294 const BlockMatrix *matrix,
01295 const unsigned int row,
01296 const unsigned int col)
01297 :
01298 matrix(matrix),
01299 base_iterator(matrix->block(0,0).begin())
01300 {
01301 Assert(col==0, ExcNotImplemented());
01302
01303
01304
01305 if (row < matrix->m())
01306 {
01307 const std::pair<unsigned int,unsigned int> indices
01308 = matrix->row_block_indices.global_to_local(row);
01309
01310
01311
01312 for (unsigned int bc=0; bc<matrix->n_block_cols(); ++bc)
01313 {
01314 base_iterator
01315 = matrix->block(indices.first, bc).begin(indices.second);
01316 if (base_iterator !=
01317 matrix->block(indices.first, bc).end(indices.second))
01318 {
01319 this->row_block = indices.first;
01320 this->col_block = bc;
01321 return;
01322 }
01323 }
01324
01325
01326
01327
01328
01329
01330
01331 *this = Accessor (matrix, row+1, 0);
01332 }
01333 else
01334 {
01335
01336
01337 this->row_block = numbers::invalid_unsigned_int;
01338 this->col_block = numbers::invalid_unsigned_int;
01339 }
01340 }
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355 template <class BlockMatrix>
01356 inline
01357 Accessor<BlockMatrix, true>::Accessor (const Accessor<BlockMatrix, false>& other)
01358 :
01359 matrix(other.matrix),
01360 base_iterator(other.base_iterator)
01361 {
01362 this->row_block = other.row_block;
01363 this->col_block = other.col_block;
01364 }
01365
01366
01367 template <class BlockMatrix>
01368 inline
01369 unsigned int
01370 Accessor<BlockMatrix, true>::row() const
01371 {
01372 Assert (this->row_block != numbers::invalid_unsigned_int,
01373 ExcIteratorPastEnd());
01374
01375 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
01376 base_iterator->row());
01377 }
01378
01379
01380 template <class BlockMatrix>
01381 inline
01382 unsigned int
01383 Accessor<BlockMatrix, true>::column() const
01384 {
01385 Assert (this->col_block != numbers::invalid_unsigned_int,
01386 ExcIteratorPastEnd());
01387
01388 return (matrix->column_block_indices.local_to_global(this->col_block,0) +
01389 base_iterator->column());
01390 }
01391
01392
01393 template <class BlockMatrix>
01394 inline
01395 typename Accessor<BlockMatrix, true>::value_type
01396 Accessor<BlockMatrix, true>::value () const
01397 {
01398 Assert (this->row_block != numbers::invalid_unsigned_int,
01399 ExcIteratorPastEnd());
01400 Assert (this->col_block != numbers::invalid_unsigned_int,
01401 ExcIteratorPastEnd());
01402
01403 return base_iterator->value();
01404 }
01405
01406
01407
01408 template <class BlockMatrix>
01409 inline
01410 void
01411 Accessor<BlockMatrix, true>::advance ()
01412 {
01413 Assert (this->row_block != numbers::invalid_unsigned_int,
01414 ExcIteratorPastEnd());
01415 Assert (this->col_block != numbers::invalid_unsigned_int,
01416 ExcIteratorPastEnd());
01417
01418
01419 unsigned int local_row = base_iterator->row();
01420
01421
01422
01423 ++base_iterator;
01424
01425
01426
01427
01428
01429
01430 while (base_iterator ==
01431 matrix->block(this->row_block, this->col_block).end(local_row))
01432 {
01433
01434
01435
01436 if (this->col_block < matrix->n_block_cols()-1)
01437 {
01438 ++this->col_block;
01439 base_iterator
01440 = matrix->block(this->row_block, this->col_block).begin(local_row);
01441 }
01442 else
01443 {
01444
01445
01446 this->col_block = 0;
01447 ++local_row;
01448
01449
01450
01451
01452
01453
01454
01455 if (local_row == matrix->block(this->row_block, this->col_block).m())
01456 {
01457 local_row = 0;
01458 ++this->row_block;
01459 if (this->row_block == matrix->n_block_rows())
01460 {
01461 this->row_block = numbers::invalid_unsigned_int;
01462 this->col_block = numbers::invalid_unsigned_int;
01463 return;
01464 }
01465 }
01466
01467 base_iterator
01468 = matrix->block(this->row_block, this->col_block).begin(local_row);
01469 }
01470 }
01471 }
01472
01473
01474 template <class BlockMatrix>
01475 inline
01476 bool
01477 Accessor<BlockMatrix, true>::operator == (const Accessor &a) const
01478 {
01479 if (matrix != a.matrix)
01480 return false;
01481
01482 if (this->row_block == a.row_block
01483 && this->col_block == a.col_block)
01484
01485
01486
01487
01488 return (((this->row_block == numbers::invalid_unsigned_int)
01489 &&
01490 (this->col_block == numbers::invalid_unsigned_int))
01491 ||
01492 (base_iterator == a.base_iterator));
01493
01494 return false;
01495 }
01496
01497
01498
01499
01500 template <class BlockMatrix>
01501 inline
01502 Accessor<BlockMatrix, false>::Accessor (
01503 BlockMatrix *matrix,
01504 const unsigned int row,
01505 const unsigned int col)
01506 :
01507 matrix(matrix),
01508 base_iterator(matrix->block(0,0).begin())
01509 {
01510 Assert(col==0, ExcNotImplemented());
01511
01512
01513 if (row < matrix->m())
01514 {
01515 const std::pair<unsigned int,unsigned int> indices
01516 = matrix->row_block_indices.global_to_local(row);
01517
01518
01519
01520 for (unsigned int bc=0; bc<matrix->n_block_cols(); ++bc)
01521 {
01522 base_iterator
01523 = matrix->block(indices.first, bc).begin(indices.second);
01524 if (base_iterator !=
01525 matrix->block(indices.first, bc).end(indices.second))
01526 {
01527 this->row_block = indices.first;
01528 this->col_block = bc;
01529 return;
01530 }
01531 }
01532
01533
01534
01535
01536
01537
01538
01539 *this = Accessor (matrix, row+1, 0);
01540 }
01541 else
01542 {
01543
01544
01545 this->row_block = numbers::invalid_unsigned_int;
01546 this->col_block = numbers::invalid_unsigned_int;
01547 }
01548 }
01549
01550
01551 template <class BlockMatrix>
01552 inline
01553 unsigned int
01554 Accessor<BlockMatrix, false>::row() const
01555 {
01556 Assert (this->row_block != numbers::invalid_unsigned_int,
01557 ExcIteratorPastEnd());
01558
01559 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
01560 base_iterator->row());
01561 }
01562
01563
01564 template <class BlockMatrix>
01565 inline
01566 unsigned int
01567 Accessor<BlockMatrix, false>::column() const
01568 {
01569 Assert (this->col_block != numbers::invalid_unsigned_int,
01570 ExcIteratorPastEnd());
01571
01572 return (matrix->column_block_indices.local_to_global(this->col_block,0) +
01573 base_iterator->column());
01574 }
01575
01576
01577 template <class BlockMatrix>
01578 inline
01579 typename Accessor<BlockMatrix, false>::value_type
01580 Accessor<BlockMatrix, false>::value () const
01581 {
01582 Assert (this->row_block != numbers::invalid_unsigned_int,
01583 ExcIteratorPastEnd());
01584 Assert (this->col_block != numbers::invalid_unsigned_int,
01585 ExcIteratorPastEnd());
01586
01587 return base_iterator->value();
01588 }
01589
01590
01591
01592 template <class BlockMatrix>
01593 inline
01594 void
01595 Accessor<BlockMatrix, false>::set_value (typename Accessor<BlockMatrix, false>::value_type newval) const
01596 {
01597 Assert (this->row_block != numbers::invalid_unsigned_int,
01598 ExcIteratorPastEnd());
01599 Assert (this->col_block != numbers::invalid_unsigned_int,
01600 ExcIteratorPastEnd());
01601
01602 base_iterator->value() = newval;
01603 }
01604
01605
01606
01607 template <class BlockMatrix>
01608 inline
01609 void
01610 Accessor<BlockMatrix, false>::advance ()
01611 {
01612 Assert (this->row_block != numbers::invalid_unsigned_int,
01613 ExcIteratorPastEnd());
01614 Assert (this->col_block != numbers::invalid_unsigned_int,
01615 ExcIteratorPastEnd());
01616
01617
01618 unsigned int local_row = base_iterator->row();
01619
01620
01621
01622 ++base_iterator;
01623
01624
01625
01626
01627
01628
01629 while (base_iterator ==
01630 matrix->block(this->row_block, this->col_block).end(local_row))
01631 {
01632
01633
01634
01635 if (this->col_block < matrix->n_block_cols()-1)
01636 {
01637 ++this->col_block;
01638 base_iterator
01639 = matrix->block(this->row_block, this->col_block).begin(local_row);
01640 }
01641 else
01642 {
01643
01644
01645 this->col_block = 0;
01646 ++local_row;
01647
01648
01649
01650
01651
01652
01653
01654 if (local_row == matrix->block(this->row_block, this->col_block).m())
01655 {
01656 local_row = 0;
01657 ++this->row_block;
01658 if (this->row_block == matrix->n_block_rows())
01659 {
01660 this->row_block = numbers::invalid_unsigned_int;
01661 this->col_block = numbers::invalid_unsigned_int;
01662 return;
01663 }
01664 }
01665
01666 base_iterator
01667 = matrix->block(this->row_block, this->col_block).begin(local_row);
01668 }
01669 }
01670 }
01671
01672
01673
01674 template <class BlockMatrix>
01675 inline
01676 bool
01677 Accessor<BlockMatrix, false>::operator == (const Accessor &a) const
01678 {
01679 if (matrix != a.matrix)
01680 return false;
01681
01682 if (this->row_block == a.row_block
01683 && this->col_block == a.col_block)
01684
01685
01686
01687
01688 return (((this->row_block == numbers::invalid_unsigned_int)
01689 &&
01690 (this->col_block == numbers::invalid_unsigned_int))
01691 ||
01692 (base_iterator == a.base_iterator));
01693
01694 return false;
01695 }
01696 }
01697
01698
01699
01700
01701
01702 template <typename MatrixType>
01703 inline
01704 BlockMatrixBase<MatrixType>::BlockMatrixBase ()
01705 {}
01706
01707
01708
01709 template <class MatrixType>
01710 template <class BlockMatrixType>
01711 inline
01712 BlockMatrixBase<MatrixType> &
01713 BlockMatrixBase<MatrixType>::
01714 copy_from (const BlockMatrixType &source)
01715 {
01716 for (unsigned int r=0; r<n_block_rows(); ++r)
01717 for (unsigned int c=0; c<n_block_cols(); ++c)
01718 block(r,c).copy_from (source.block(r,c));
01719
01720 return *this;
01721 }
01722
01723
01724 template <class MatrixType>
01725 std::size_t
01726 BlockMatrixBase<MatrixType>::memory_consumption () const
01727 {
01728 std::size_t mem =
01729 MemoryConsumption::memory_consumption(row_block_indices)+
01730 MemoryConsumption::memory_consumption(column_block_indices)+
01731 MemoryConsumption::memory_consumption(sub_objects)+
01732 MemoryConsumption::memory_consumption(counter_within_block)+
01733 MemoryConsumption::memory_consumption(column_indices)+
01734 MemoryConsumption::memory_consumption(column_values);
01735
01736 for (unsigned int r=0; r<n_block_rows(); ++r)
01737 for (unsigned int c=0; c<n_block_cols(); ++c)
01738 {
01739 MatrixType *p = this->sub_objects[r][c];
01740 mem += MemoryConsumption::memory_consumption(*p);
01741 }
01742
01743 return mem;
01744 }
01745
01746
01747
01748 template <class MatrixType>
01749 inline
01750 void
01751 BlockMatrixBase<MatrixType>::clear ()
01752 {
01753 for (unsigned int r=0; r<n_block_rows(); ++r)
01754 for (unsigned int c=0; c<n_block_cols(); ++c)
01755 {
01756 MatrixType *p = this->sub_objects[r][c];
01757 this->sub_objects[r][c] = 0;
01758 delete p;
01759 }
01760 sub_objects.reinit (0,0);
01761
01762
01763 row_block_indices = column_block_indices = BlockIndices ();
01764 }
01765
01766
01767
01768 template <class MatrixType>
01769 inline
01770 typename BlockMatrixBase<MatrixType>::BlockType &
01771 BlockMatrixBase<MatrixType>::block (const unsigned int row,
01772 const unsigned int column)
01773 {
01774 Assert (row<n_block_rows(),
01775 ExcIndexRange (row, 0, n_block_rows()));
01776 Assert (column<n_block_cols(),
01777 ExcIndexRange (column, 0, n_block_cols()));
01778
01779 return *sub_objects[row][column];
01780 }
01781
01782
01783
01784 template <class MatrixType>
01785 inline
01786 const typename BlockMatrixBase<MatrixType>::BlockType &
01787 BlockMatrixBase<MatrixType>::block (const unsigned int row,
01788 const unsigned int column) const
01789 {
01790 Assert (row<n_block_rows(),
01791 ExcIndexRange (row, 0, n_block_rows()));
01792 Assert (column<n_block_cols(),
01793 ExcIndexRange (column, 0, n_block_cols()));
01794
01795 return *sub_objects[row][column];
01796 }
01797
01798
01799 template <class MatrixType>
01800 inline
01801 unsigned int
01802 BlockMatrixBase<MatrixType>::m () const
01803 {
01804 return row_block_indices.total_size();
01805 }
01806
01807
01808
01809 template <class MatrixType>
01810 inline
01811 unsigned int
01812 BlockMatrixBase<MatrixType>::n () const
01813 {
01814 return column_block_indices.total_size();
01815 }
01816
01817
01818
01819 template <class MatrixType>
01820 inline
01821 unsigned int
01822 BlockMatrixBase<MatrixType>::n_block_cols () const
01823 {
01824 return column_block_indices.size();
01825 }
01826
01827
01828
01829 template <class MatrixType>
01830 inline
01831 unsigned int
01832 BlockMatrixBase<MatrixType>::n_block_rows () const
01833 {
01834 return row_block_indices.size();
01835 }
01836
01837
01838
01839
01840
01841
01842 template <class MatrixType>
01843 inline
01844 void
01845 BlockMatrixBase<MatrixType>::set (const unsigned int i,
01846 const unsigned int j,
01847 const value_type value)
01848 {
01849 prepare_set_operation();
01850
01851 Assert (numbers::is_finite(value), ExcNumberNotFinite());
01852
01853 const std::pair<unsigned int,unsigned int>
01854 row_index = row_block_indices.global_to_local (i),
01855 col_index = column_block_indices.global_to_local (j);
01856 block(row_index.first,col_index.first).set (row_index.second,
01857 col_index.second,
01858 value);
01859 }
01860
01861
01862
01863 template <class MatrixType>
01864 template <typename number>
01865 inline
01866 void
01867 BlockMatrixBase<MatrixType>::set (const std::vector<unsigned int> &row_indices,
01868 const std::vector<unsigned int> &col_indices,
01869 const FullMatrix<number> &values,
01870 const bool elide_zero_values)
01871 {
01872 Assert (row_indices.size() == values.m(),
01873 ExcDimensionMismatch(row_indices.size(), values.m()));
01874 Assert (col_indices.size() == values.n(),
01875 ExcDimensionMismatch(col_indices.size(), values.n()));
01876
01877 for (unsigned int i=0; i<row_indices.size(); ++i)
01878 set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
01879 elide_zero_values);
01880 }
01881
01882
01883
01884 template <class MatrixType>
01885 template <typename number>
01886 inline
01887 void
01888 BlockMatrixBase<MatrixType>::set (const std::vector<unsigned int> &indices,
01889 const FullMatrix<number> &values,
01890 const bool elide_zero_values)
01891 {
01892 Assert (indices.size() == values.m(),
01893 ExcDimensionMismatch(indices.size(), values.m()));
01894 Assert (values.n() == values.m(), ExcNotQuadratic());
01895
01896 for (unsigned int i=0; i<indices.size(); ++i)
01897 set (indices[i], indices.size(), &indices[0], &values(i,0),
01898 elide_zero_values);
01899 }
01900
01901
01902
01903 template <class MatrixType>
01904 template <typename number>
01905 inline
01906 void
01907 BlockMatrixBase<MatrixType>::set (const unsigned int row,
01908 const std::vector<unsigned int> &col_indices,
01909 const std::vector<number> &values,
01910 const bool elide_zero_values)
01911 {
01912 Assert (col_indices.size() == values.size(),
01913 ExcDimensionMismatch(col_indices.size(), values.size()));
01914
01915 set (row, col_indices.size(), &col_indices[0], &values[0],
01916 elide_zero_values);
01917 }
01918
01919
01920
01921
01922
01923
01924 template <class MatrixType>
01925 template <typename number>
01926 inline
01927 void
01928 BlockMatrixBase<MatrixType>::set (const unsigned int row,
01929 const unsigned int n_cols,
01930 const unsigned int *col_indices,
01931 const number *values,
01932 const bool elide_zero_values)
01933 {
01934 prepare_set_operation();
01935
01936
01937 if (column_indices.size() < this->n_block_cols())
01938 {
01939 column_indices.resize (this->n_block_cols());
01940 column_values.resize (this->n_block_cols());
01941 counter_within_block.resize (this->n_block_cols());
01942 }
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954 if (column_indices[0].size() < n_cols)
01955 {
01956 for (unsigned int i=0; i<this->n_block_cols(); ++i)
01957 {
01958 column_indices[i].resize(n_cols);
01959 column_values[i].resize(n_cols);
01960 }
01961 }
01962
01963
01964
01965 for (unsigned int i=0; i<this->n_block_cols(); ++i)
01966 counter_within_block[i] = 0;
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977 for (unsigned int j=0; j<n_cols; ++j)
01978 {
01979 double value = values[j];
01980
01981 if (value == 0 && elide_zero_values == true)
01982 continue;
01983
01984 const std::pair<unsigned int, unsigned int>
01985 col_index = this->column_block_indices.global_to_local(col_indices[j]);
01986
01987 const unsigned int local_index = counter_within_block[col_index.first]++;
01988
01989 column_indices[col_index.first][local_index] = col_index.second;
01990 column_values[col_index.first][local_index] = value;
01991 }
01992
01993 #ifdef DEBUG
01994
01995
01996 unsigned int length = 0;
01997 for (unsigned int i=0; i<this->n_block_cols(); ++i)
01998 length += counter_within_block[i];
01999 Assert (length <= n_cols, ExcInternalError());
02000 #endif
02001
02002
02003
02004
02005
02006
02007 const std::pair<unsigned int,unsigned int>
02008 row_index = this->row_block_indices.global_to_local (row);
02009 for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
02010 {
02011 if (counter_within_block[block_col] == 0)
02012 continue;
02013
02014 block(row_index.first, block_col).set
02015 (row_index.second,
02016 counter_within_block[block_col],
02017 &column_indices[block_col][0],
02018 &column_values[block_col][0],
02019 false);
02020 }
02021 }
02022
02023
02024
02025 template <class MatrixType>
02026 inline
02027 void
02028 BlockMatrixBase<MatrixType>::add (const unsigned int i,
02029 const unsigned int j,
02030 const value_type value)
02031 {
02032
02033 Assert (numbers::is_finite(value), ExcNumberNotFinite());
02034
02035 prepare_add_operation();
02036
02037
02038
02039
02040 typedef typename MatrixType::Traits MatrixTraits;
02041 if ((MatrixTraits::zero_addition_can_be_elided == true)
02042 &&
02043 (value == value_type()))
02044 return;
02045
02046 const std::pair<unsigned int,unsigned int>
02047 row_index = row_block_indices.global_to_local (i),
02048 col_index = column_block_indices.global_to_local (j);
02049 block(row_index.first,col_index.first).add (row_index.second,
02050 col_index.second,
02051 value);
02052 }
02053
02054
02055
02056 template <class MatrixType>
02057 template <typename number>
02058 inline
02059 void
02060 BlockMatrixBase<MatrixType>::add (const std::vector<unsigned int> &row_indices,
02061 const std::vector<unsigned int> &col_indices,
02062 const FullMatrix<number> &values,
02063 const bool elide_zero_values)
02064 {
02065 Assert (row_indices.size() == values.m(),
02066 ExcDimensionMismatch(row_indices.size(), values.m()));
02067 Assert (col_indices.size() == values.n(),
02068 ExcDimensionMismatch(col_indices.size(), values.n()));
02069
02070 for (unsigned int i=0; i<row_indices.size(); ++i)
02071 add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
02072 elide_zero_values);
02073 }
02074
02075
02076
02077 template <class MatrixType>
02078 template <typename number>
02079 inline
02080 void
02081 BlockMatrixBase<MatrixType>::add (const std::vector<unsigned int> &indices,
02082 const FullMatrix<number> &values,
02083 const bool elide_zero_values)
02084 {
02085 Assert (indices.size() == values.m(),
02086 ExcDimensionMismatch(indices.size(), values.m()));
02087 Assert (values.n() == values.m(), ExcNotQuadratic());
02088
02089 for (unsigned int i=0; i<indices.size(); ++i)
02090 add (indices[i], indices.size(), &indices[0], &values(i,0),
02091 elide_zero_values);
02092 }
02093
02094
02095
02096 template <class MatrixType>
02097 template <typename number>
02098 inline
02099 void
02100 BlockMatrixBase<MatrixType>::add (const unsigned int row,
02101 const std::vector<unsigned int> &col_indices,
02102 const std::vector<number> &values,
02103 const bool elide_zero_values)
02104 {
02105 Assert (col_indices.size() == values.size(),
02106 ExcDimensionMismatch(col_indices.size(), values.size()));
02107
02108 add (row, col_indices.size(), &col_indices[0], &values[0],
02109 elide_zero_values);
02110 }
02111
02112
02113
02114
02115
02116
02117 template <class MatrixType>
02118 template <typename number>
02119 inline
02120 void
02121 BlockMatrixBase<MatrixType>::add (const unsigned int row,
02122 const unsigned int n_cols,
02123 const unsigned int *col_indices,
02124 const number *values,
02125 const bool elide_zero_values,
02126 const bool col_indices_are_sorted)
02127 {
02128 prepare_add_operation();
02129
02130
02131
02132
02133 if (col_indices_are_sorted == true)
02134 {
02135 #ifdef DEBUG
02136
02137
02138 unsigned int before = col_indices[0];
02139 for (unsigned int i=1; i<n_cols; ++i)
02140 if (col_indices[i] <= before)
02141 Assert (false, ExcMessage ("Flag col_indices_are_sorted is set, but "
02142 "indices appear to not be sorted."))
02143 else
02144 before = col_indices[i];
02145 #endif
02146 const std::pair<unsigned int,unsigned int>
02147 row_index = this->row_block_indices.global_to_local (row);
02148
02149 if (this->n_block_cols() > 1)
02150 {
02151 const unsigned int * first_block = Utilities::lower_bound (col_indices,
02152 col_indices+n_cols,
02153 this->column_block_indices.block_start(1));
02154
02155 const unsigned int n_zero_block_indices = first_block - col_indices;
02156 block(row_index.first, 0).add (row_index.second,
02157 n_zero_block_indices,
02158 col_indices,
02159 values,
02160 elide_zero_values,
02161 col_indices_are_sorted);
02162
02163 if (n_zero_block_indices < n_cols)
02164 this->add(row, n_cols - n_zero_block_indices, first_block,
02165 values + n_zero_block_indices, elide_zero_values,
02166 false);
02167 }
02168 else
02169 {
02170 block(row_index.first, 0). add (row_index.second,
02171 n_cols,
02172 col_indices,
02173 values,
02174 elide_zero_values,
02175 col_indices_are_sorted);
02176 }
02177
02178 return;
02179 }
02180
02181
02182 if (column_indices.size() < this->n_block_cols())
02183 {
02184 column_indices.resize (this->n_block_cols());
02185 column_values.resize (this->n_block_cols());
02186 counter_within_block.resize (this->n_block_cols());
02187 }
02188
02189
02190
02191
02192
02193
02194
02195
02196
02197
02198
02199 if (column_indices[0].size() < n_cols)
02200 {
02201 for (unsigned int i=0; i<this->n_block_cols(); ++i)
02202 {
02203 column_indices[i].resize(n_cols);
02204 column_values[i].resize(n_cols);
02205 }
02206 }
02207
02208
02209
02210 for (unsigned int i=0; i<this->n_block_cols(); ++i)
02211 counter_within_block[i] = 0;
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221
02222 for (unsigned int j=0; j<n_cols; ++j)
02223 {
02224 double value = values[j];
02225
02226 if (value == 0 && elide_zero_values == true)
02227 continue;
02228
02229 const std::pair<unsigned int, unsigned int>
02230 col_index = this->column_block_indices.global_to_local(col_indices[j]);
02231
02232 const unsigned int local_index = counter_within_block[col_index.first]++;
02233
02234 column_indices[col_index.first][local_index] = col_index.second;
02235 column_values[col_index.first][local_index] = value;
02236 }
02237
02238 #ifdef DEBUG
02239
02240
02241 unsigned int length = 0;
02242 for (unsigned int i=0; i<this->n_block_cols(); ++i)
02243 length += counter_within_block[i];
02244 Assert (length <= n_cols, ExcInternalError());
02245 #endif
02246
02247
02248
02249
02250
02251
02252 const std::pair<unsigned int,unsigned int>
02253 row_index = this->row_block_indices.global_to_local (row);
02254 for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
02255 {
02256 if (counter_within_block[block_col] == 0)
02257 continue;
02258
02259 block(row_index.first, block_col).add
02260 (row_index.second,
02261 counter_within_block[block_col],
02262 &column_indices[block_col][0],
02263 &column_values[block_col][0],
02264 false,
02265 col_indices_are_sorted);
02266 }
02267 }
02268
02269
02270
02271 template <class MatrixType>
02272 inline
02273 typename BlockMatrixBase<MatrixType>::value_type
02274 BlockMatrixBase<MatrixType>::operator () (const unsigned int i,
02275 const unsigned int j) const
02276 {
02277 const std::pair<unsigned int,unsigned int>
02278 row_index = row_block_indices.global_to_local (i),
02279 col_index = column_block_indices.global_to_local (j);
02280 return block(row_index.first,col_index.first) (row_index.second,
02281 col_index.second);
02282 }
02283
02284
02285
02286 template <class MatrixType>
02287 inline
02288 typename BlockMatrixBase<MatrixType>::value_type
02289 BlockMatrixBase<MatrixType>::el (const unsigned int i,
02290 const unsigned int j) const
02291 {
02292 const std::pair<unsigned int,unsigned int>
02293 row_index = row_block_indices.global_to_local (i),
02294 col_index = column_block_indices.global_to_local (j);
02295 return block(row_index.first,col_index.first).el (row_index.second,
02296 col_index.second);
02297 }
02298
02299
02300
02301 template <class MatrixType>
02302 inline
02303 typename BlockMatrixBase<MatrixType>::value_type
02304 BlockMatrixBase<MatrixType>::diag_element (const unsigned int i) const
02305 {
02306 Assert (n_block_rows() == n_block_cols(),
02307 ExcNotQuadratic());
02308
02309 const std::pair<unsigned int,unsigned int>
02310 index = row_block_indices.global_to_local (i);
02311 return block(index.first,index.first).diag_element(index.second);
02312 }
02313
02314
02315
02316 template <class MatrixType>
02317 inline
02318 void
02319 BlockMatrixBase<MatrixType>::compress ()
02320 {
02321 for (unsigned int r=0; r<n_block_rows(); ++r)
02322 for (unsigned int c=0; c<n_block_cols(); ++c)
02323 block(r,c).compress ();
02324 }
02325
02326
02327
02328 template <class MatrixType>
02329 inline
02330 BlockMatrixBase<MatrixType> &
02331 BlockMatrixBase<MatrixType>::operator *= (const value_type factor)
02332 {
02333 Assert (n_block_cols() != 0, ExcNotInitialized());
02334 Assert (n_block_rows() != 0, ExcNotInitialized());
02335
02336 for (unsigned int r=0; r<n_block_rows(); ++r)
02337 for (unsigned int c=0; c<n_block_cols(); ++c)
02338 block(r,c) *= factor;
02339
02340 return *this;
02341 }
02342
02343
02344
02345 template <class MatrixType>
02346 inline
02347 BlockMatrixBase<MatrixType> &
02348 BlockMatrixBase<MatrixType>::operator /= (const value_type factor)
02349 {
02350 Assert (n_block_cols() != 0, ExcNotInitialized());
02351 Assert (n_block_rows() != 0, ExcNotInitialized());
02352 Assert (factor !=0, ExcDivideByZero());
02353
02354 const value_type factor_inv = 1. / factor;
02355
02356 for (unsigned int r=0; r<n_block_rows(); ++r)
02357 for (unsigned int c=0; c<n_block_cols(); ++c)
02358 block(r,c) *= factor_inv;
02359
02360 return *this;
02361 }
02362
02363
02364
02365 template <class MatrixType>
02366 const BlockIndices &
02367 BlockMatrixBase<MatrixType>::get_row_indices () const
02368 {
02369 return this->row_block_indices;
02370 }
02371
02372
02373
02374 template <class MatrixType>
02375 const BlockIndices &
02376 BlockMatrixBase<MatrixType>::get_column_indices () const
02377 {
02378 return this->column_block_indices;
02379 }
02380
02381
02382
02383 template <class MatrixType>
02384 template <class BlockVectorType>
02385 void
02386 BlockMatrixBase<MatrixType>::
02387 vmult_block_block (BlockVectorType &dst,
02388 const BlockVectorType &src) const
02389 {
02390 Assert (dst.n_blocks() == n_block_rows(),
02391 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
02392 Assert (src.n_blocks() == n_block_cols(),
02393 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
02394
02395 for (unsigned int row=0; row<n_block_rows(); ++row)
02396 {
02397 block(row,0).vmult (dst.block(row),
02398 src.block(0));
02399 for (unsigned int col=1; col<n_block_cols(); ++col)
02400 block(row,col).vmult_add (dst.block(row),
02401 src.block(col));
02402 };
02403 }
02404
02405
02406
02407 template <class MatrixType>
02408 template <class BlockVectorType,
02409 class VectorType>
02410 void
02411 BlockMatrixBase<MatrixType>::
02412 vmult_nonblock_block (VectorType &dst,
02413 const BlockVectorType &src) const
02414 {
02415 Assert (n_block_rows() == 1,
02416 ExcDimensionMismatch(1, n_block_rows()));
02417 Assert (src.n_blocks() == n_block_cols(),
02418 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
02419
02420 block(0,0).vmult (dst, src.block(0));
02421 for (unsigned int col=1; col<n_block_cols(); ++col)
02422 block(0,col).vmult_add (dst, src.block(col));
02423 }
02424
02425
02426
02427 template <class MatrixType>
02428 template <class BlockVectorType,
02429 class VectorType>
02430 void
02431 BlockMatrixBase<MatrixType>::
02432 vmult_block_nonblock (BlockVectorType &dst,
02433 const VectorType &src) const
02434 {
02435 Assert (dst.n_blocks() == n_block_rows(),
02436 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
02437 Assert (1 == n_block_cols(),
02438 ExcDimensionMismatch(1, n_block_cols()));
02439
02440 for (unsigned int row=0; row<n_block_rows(); ++row)
02441 block(row,0).vmult (dst.block(row),
02442 src);
02443 }
02444
02445
02446
02447 template <class MatrixType>
02448 template <class VectorType>
02449 void
02450 BlockMatrixBase<MatrixType>::
02451 vmult_nonblock_nonblock (VectorType &dst,
02452 const VectorType &src) const
02453 {
02454 Assert (1 == n_block_rows(),
02455 ExcDimensionMismatch(1, n_block_rows()));
02456 Assert (1 == n_block_cols(),
02457 ExcDimensionMismatch(1, n_block_cols()));
02458
02459 block(0,0).vmult (dst, src);
02460 }
02461
02462
02463
02464 template <class MatrixType>
02465 template <class BlockVectorType>
02466 void
02467 BlockMatrixBase<MatrixType>::vmult_add (BlockVectorType &dst,
02468 const BlockVectorType &src) const
02469 {
02470 Assert (dst.n_blocks() == n_block_rows(),
02471 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
02472 Assert (src.n_blocks() == n_block_cols(),
02473 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
02474
02475 for (unsigned int row=0; row<n_block_rows(); ++row)
02476 {
02477 block(row,0).vmult_add (dst.block(row),
02478 src.block(0));
02479 for (unsigned int col=1; col<n_block_cols(); ++col)
02480 block(row,col).vmult_add (dst.block(row),
02481 src.block(col));
02482 };
02483 }
02484
02485
02486
02487
02488 template <class MatrixType>
02489 template <class BlockVectorType>
02490 void
02491 BlockMatrixBase<MatrixType>::
02492 Tvmult_block_block (BlockVectorType &dst,
02493 const BlockVectorType &src) const
02494 {
02495 Assert (dst.n_blocks() == n_block_cols(),
02496 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
02497 Assert (src.n_blocks() == n_block_rows(),
02498 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
02499
02500 dst = 0.;
02501
02502 for (unsigned int row=0; row<n_block_rows(); ++row)
02503 {
02504 for (unsigned int col=0; col<n_block_cols(); ++col)
02505 block(row,col).Tvmult_add (dst.block(col),
02506 src.block(row));
02507 };
02508 }
02509
02510
02511
02512 template <class MatrixType>
02513 template <class BlockVectorType,
02514 class VectorType>
02515 void
02516 BlockMatrixBase<MatrixType>::
02517 Tvmult_block_nonblock (BlockVectorType &dst,
02518 const VectorType &src) const
02519 {
02520 Assert (dst.n_blocks() == n_block_cols(),
02521 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
02522 Assert (1 == n_block_rows(),
02523 ExcDimensionMismatch(1, n_block_rows()));
02524
02525 dst = 0.;
02526
02527 for (unsigned int col=0; col<n_block_cols(); ++col)
02528 block(0,col).Tvmult_add (dst.block(col), src);
02529 }
02530
02531
02532
02533 template <class MatrixType>
02534 template <class BlockVectorType,
02535 class VectorType>
02536 void
02537 BlockMatrixBase<MatrixType>::
02538 Tvmult_nonblock_block (VectorType &dst,
02539 const BlockVectorType &src) const
02540 {
02541 Assert (1 == n_block_cols(),
02542 ExcDimensionMismatch(1, n_block_cols()));
02543 Assert (src.n_blocks() == n_block_rows(),
02544 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
02545
02546 block(0,0).Tvmult (dst, src.block(0));
02547
02548 for (unsigned int row=1; row<n_block_rows(); ++row)
02549 block(row,0).Tvmult_add (dst, src.block(row));
02550 }
02551
02552
02553
02554 template <class MatrixType>
02555 template <class VectorType>
02556 void
02557 BlockMatrixBase<MatrixType>::
02558 Tvmult_nonblock_nonblock (VectorType &dst,
02559 const VectorType &src) const
02560 {
02561 Assert (1 == n_block_cols(),
02562 ExcDimensionMismatch(1, n_block_cols()));
02563 Assert (1 == n_block_rows(),
02564 ExcDimensionMismatch(1, n_block_rows()));
02565
02566 block(0,0).Tvmult (dst, src);
02567 }
02568
02569
02570
02571 template <class MatrixType>
02572 template <class BlockVectorType>
02573 void
02574 BlockMatrixBase<MatrixType>::Tvmult_add (BlockVectorType& dst,
02575 const BlockVectorType& src) const
02576 {
02577 Assert (dst.n_blocks() == n_block_cols(),
02578 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
02579 Assert (src.n_blocks() == n_block_rows(),
02580 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
02581
02582 for (unsigned int row=0; row<n_block_rows(); ++row)
02583 {
02584 for (unsigned int col=0; col<n_block_cols(); ++col)
02585 block(row,col).Tvmult_add (dst.block(col),
02586 src.block(row));
02587 };
02588 }
02589
02590
02591
02592 template <class MatrixType>
02593 template <class BlockVectorType>
02594 typename BlockMatrixBase<MatrixType>::value_type
02595 BlockMatrixBase<MatrixType>::matrix_norm_square (const BlockVectorType &v) const
02596 {
02597 Assert (n_block_rows() == n_block_cols(), ExcNotQuadratic());
02598 Assert (v.n_blocks() == n_block_rows(),
02599 ExcDimensionMismatch(v.n_blocks(), n_block_rows()));
02600
02601 value_type norm_sqr = 0;
02602 for (unsigned int row=0; row<n_block_rows(); ++row)
02603 for (unsigned int col=0; col<n_block_cols(); ++col)
02604 if (row==col)
02605 norm_sqr += block(row,col).matrix_norm_square (v.block(row));
02606 else
02607 norm_sqr += block(row,col).matrix_scalar_product (v.block(row),
02608 v.block(col));
02609 return norm_sqr;
02610 }
02611
02612
02613
02614 template <class MatrixType>
02615 template <class BlockVectorType>
02616 typename BlockMatrixBase<MatrixType>::value_type
02617 BlockMatrixBase<MatrixType>::
02618 matrix_scalar_product (const BlockVectorType &u,
02619 const BlockVectorType &v) const
02620 {
02621 Assert (u.n_blocks() == n_block_rows(),
02622 ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
02623 Assert (v.n_blocks() == n_block_cols(),
02624 ExcDimensionMismatch(v.n_blocks(), n_block_cols()));
02625
02626 value_type result = 0;
02627 for (unsigned int row=0; row<n_block_rows(); ++row)
02628 for (unsigned int col=0; col<n_block_cols(); ++col)
02629 result += block(row,col).matrix_scalar_product (u.block(row),
02630 v.block(col));
02631 return result;
02632 }
02633
02634
02635
02636 template <class MatrixType>
02637 template <class BlockVectorType>
02638 typename BlockMatrixBase<MatrixType>::value_type
02639 BlockMatrixBase<MatrixType>::
02640 residual (BlockVectorType &dst,
02641 const BlockVectorType &x,
02642 const BlockVectorType &b) const
02643 {
02644 Assert (dst.n_blocks() == n_block_rows(),
02645 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
02646 Assert (b.n_blocks() == n_block_rows(),
02647 ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
02648 Assert (x.n_blocks() == n_block_cols(),
02649 ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
02650
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664 for (unsigned int row=0; row<n_block_rows(); ++row)
02665 {
02666 block(row,0).residual (dst.block(row),
02667 x.block(0),
02668 b.block(row));
02669
02670 for (unsigned int i=0; i<dst.block(row).size(); ++i)
02671 dst.block(row)(i) = -dst.block(row)(i);
02672
02673 for (unsigned int col=1; col<n_block_cols(); ++col)
02674 block(row,col).vmult_add (dst.block(row),
02675 x.block(col));
02676
02677 for (unsigned int i=0; i<dst.block(row).size(); ++i)
02678 dst.block(row)(i) = -dst.block(row)(i);
02679 };
02680
02681 value_type res = 0;
02682 for (unsigned int row=0; row<n_block_rows(); ++row)
02683 res += dst.block(row).norm_sqr ();
02684 return std::sqrt(res);
02685 }
02686
02687
02688
02689 template <class MatrixType>
02690 inline
02691 typename BlockMatrixBase<MatrixType>::const_iterator
02692 BlockMatrixBase<MatrixType>::begin () const
02693 {
02694 return const_iterator(this, 0);
02695 }
02696
02697
02698
02699 template <class MatrixType>
02700 inline
02701 typename BlockMatrixBase<MatrixType>::const_iterator
02702 BlockMatrixBase<MatrixType>::end () const
02703 {
02704 return const_iterator(this, m());
02705 }
02706
02707
02708
02709 template <class MatrixType>
02710 inline
02711 typename BlockMatrixBase<MatrixType>::const_iterator
02712 BlockMatrixBase<MatrixType>::begin (const unsigned int r) const
02713 {
02714 Assert (r<m(), ExcIndexRange(r,0,m()));
02715 return const_iterator(this, r);
02716 }
02717
02718
02719
02720 template <class MatrixType>
02721 inline
02722 typename BlockMatrixBase<MatrixType>::const_iterator
02723 BlockMatrixBase<MatrixType>::end (const unsigned int r) const
02724 {
02725 Assert (r<m(), ExcIndexRange(r,0,m()));
02726 return const_iterator(this, r+1);
02727 }
02728
02729
02730
02731 template <class MatrixType>
02732 inline
02733 typename BlockMatrixBase<MatrixType>::iterator
02734 BlockMatrixBase<MatrixType>::begin ()
02735 {
02736 return iterator(this, 0);
02737 }
02738
02739
02740
02741 template <class MatrixType>
02742 inline
02743 typename BlockMatrixBase<MatrixType>::iterator
02744 BlockMatrixBase<MatrixType>::end ()
02745 {
02746 return iterator(this, m());
02747 }
02748
02749
02750
02751 template <class MatrixType>
02752 inline
02753 typename BlockMatrixBase<MatrixType>::iterator
02754 BlockMatrixBase<MatrixType>::begin (const unsigned int r)
02755 {
02756 Assert (r<m(), ExcIndexRange(r,0,m()));
02757 return iterator(this, r);
02758 }
02759
02760
02761
02762 template <class MatrixType>
02763 inline
02764 typename BlockMatrixBase<MatrixType>::iterator
02765 BlockMatrixBase<MatrixType>::end (const unsigned int r)
02766 {
02767 Assert (r<m(), ExcIndexRange(r,0,m()));
02768 return iterator(this, r+1);
02769 }
02770
02771
02772
02773 template <class MatrixType>
02774 void
02775 BlockMatrixBase<MatrixType>::collect_sizes ()
02776 {
02777 std::vector<unsigned int> row_sizes (this->n_block_rows());
02778 std::vector<unsigned int> col_sizes (this->n_block_cols());
02779
02780
02781
02782 for (unsigned int r=0; r<this->n_block_rows(); ++r)
02783 row_sizes[r] = sub_objects[r][0]->m();
02784
02785
02786
02787 for (unsigned int c=1; c<this->n_block_cols(); ++c)
02788 for (unsigned int r=0; r<this->n_block_rows(); ++r)
02789 Assert (row_sizes[r] == sub_objects[r][c]->m(),
02790 ExcIncompatibleRowNumbers (r,0,r,c));
02791
02792
02793
02794 this->row_block_indices.reinit (row_sizes);
02795
02796
02797
02798 for (unsigned int c=0; c<this->n_block_cols(); ++c)
02799 col_sizes[c] = sub_objects[0][c]->n();
02800 for (unsigned int r=1; r<this->n_block_rows(); ++r)
02801 for (unsigned int c=0; c<this->n_block_cols(); ++c)
02802 Assert (col_sizes[c] == sub_objects[r][c]->n(),
02803 ExcIncompatibleRowNumbers (0,c,r,c));
02804
02805
02806
02807 this->column_block_indices.reinit (col_sizes);
02808 }
02809
02810
02811
02812 template <class MatrixType>
02813 void
02814 BlockMatrixBase<MatrixType>::prepare_add_operation ()
02815 {
02816 for (unsigned int row=0; row<n_block_rows(); ++row)
02817 for (unsigned int col=0; col<n_block_cols(); ++col)
02818 block(row, col).prepare_add();
02819 }
02820
02821
02822
02823 template <class MatrixType>
02824 void
02825 BlockMatrixBase<MatrixType>::prepare_set_operation ()
02826 {
02827 for (unsigned int row=0; row<n_block_rows(); ++row)
02828 for (unsigned int col=0; col<n_block_cols(); ++col)
02829 block(row, col).prepare_set();
02830 }
02831
02832 #endif // DOXYGEN
02833
02834
02835 DEAL_II_NAMESPACE_CLOSE
02836
02837 #endif // __deal2__block_matrix_base_h
02838