include/deal.II/lac/block_matrix_base.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: block_matrix_base.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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 /* ------------------------- Template functions ---------------------- */
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                                      // check if this is a regular row or
01304                                      // the end of the matrix
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                                          // find the first block that does
01311                                          // have an entry in this row
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                                          // hm, there is no block that has
01326                                          // an entry in this column. we need
01327                                          // to take the next entry then,
01328                                          // which may be the first entry of
01329                                          // the next row, or recursively the
01330                                          // next row, or so on
01331         *this = Accessor (matrix, row+1, 0);
01332       }
01333     else
01334       {
01335                                          // we were asked to create the end
01336                                          // iterator for this matrix
01337         this->row_block = numbers::invalid_unsigned_int;
01338         this->col_block = numbers::invalid_unsigned_int;
01339       }
01340   }
01341 
01342 
01343 //   template <class BlockMatrix>
01344 //   inline
01345 //   Accessor<BlockMatrix, true>::Accessor (const Accessor<BlockMatrix, true>& other)
01346 //                :
01347 //                matrix(other.matrix),
01348 //                base_iterator(other.base_iterator)
01349 //   {
01350 //     this->row_block = other.row_block;
01351 //     this->col_block = other.col_block;
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                                      // Remember current row inside block
01419     unsigned int local_row = base_iterator->row();
01420 
01421                                      // Advance one element inside the
01422                                      // current block
01423     ++base_iterator;
01424 
01425                                      // while we hit the end of the row of a
01426                                      // block (which may happen multiple
01427                                      // times if rows inside a block are
01428                                      // empty), we have to jump to the next
01429                                      // block and take the
01430     while (base_iterator ==
01431            matrix->block(this->row_block, this->col_block).end(local_row))
01432       {
01433                                          // jump to next block in this block
01434                                          // row, if possible, otherwise go
01435                                          // to next row
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                                              // jump back to next row in
01445                                              // first block column
01446             this->col_block = 0;
01447             ++local_row;
01448 
01449                                              // see if this has brought us
01450                                              // past the number of rows in
01451                                              // this block. if so see
01452                                              // whether we've just fallen
01453                                              // off the end of the whole
01454                                              // matrix
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                                        // end iterators do not necessarily
01485                                        // have to have the same
01486                                        // base_iterator representation, but
01487                                        // valid iterators have to
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                                      // check if this is a regular row or
01512                                      // the end of the matrix
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                                          // find the first block that does
01519                                          // have an entry in this row
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                                          // hm, there is no block that has
01534                                          // an entry in this column. we need
01535                                          // to take the next entry then,
01536                                          // which may be the first entry of
01537                                          // the next row, or recursively the
01538                                          // next row, or so on
01539         *this = Accessor (matrix, row+1, 0);
01540       }
01541     else
01542       {
01543                                          // we were asked to create the end
01544                                          // iterator for this matrix
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                                      // Remember current row inside block
01618     unsigned int local_row = base_iterator->row();
01619 
01620                                      // Advance one element inside the
01621                                      // current block
01622     ++base_iterator;
01623 
01624                                      // while we hit the end of the row of a
01625                                      // block (which may happen multiple
01626                                      // times if rows inside a block are
01627                                      // empty), we have to jump to the next
01628                                      // block and take the
01629     while (base_iterator ==
01630            matrix->block(this->row_block, this->col_block).end(local_row))
01631       {
01632                                          // jump to next block in this block
01633                                          // row, if possible, otherwise go
01634                                          // to next row
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                                              // jump back to next row in
01644                                              // first block column
01645             this->col_block = 0;
01646             ++local_row;
01647 
01648                                              // see if this has brought us
01649                                              // past the number of rows in
01650                                              // this block. if so see
01651                                              // whether we've just fallen
01652                                              // off the end of the whole
01653                                              // matrix
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                                        // end iterators do not necessarily
01685                                        // have to have the same
01686                                        // base_iterator representation, but
01687                                        // valid iterators have to
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                                    // reset block indices to empty
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                                  // Write the single set manually,
01840                                  // since the other function has a lot
01841                                  // of overhead in that case.
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                                  // This is a very messy function, since
01922                                  // we need to calculate to each position
01923                                  // the location in the global array.
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                                    // Resize scratch arrays
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                                    // Resize sub-arrays to n_cols. This
01945                                    // is a bit wasteful, but we resize
01946                                    // only a few times (then the maximum
01947                                    // row length won't increase that
01948                                    // much any more). At least we know
01949                                    // that all arrays are going to be of
01950                                    // the same size, so we can check
01951                                    // whether the size of one is large
01952                                    // enough before actually going
01953                                    // through all of them.
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                                    // Reset the number of added elements
01964                                    // in each block to zero.
01965   for (unsigned int i=0; i<this->n_block_cols(); ++i)
01966     counter_within_block[i] = 0;
01967 
01968                                    // Go through the column indices to
01969                                    // find out which portions of the
01970                                    // values should be set in which
01971                                    // block of the matrix. We need to
01972                                    // touch all the data, since we can't
01973                                    // be sure that the data of one block
01974                                    // is stored contiguously (in fact,
01975                                    // indices will be intermixed when it
01976                                    // comes from an element matrix).
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                                    // If in debug mode, do a check whether
01995                                    // the right length has been obtained.
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                                    // Now we found out about where the
02003                                    // individual columns should start and
02004                                    // where we should start reading out
02005                                    // data. Now let's write the data into
02006                                    // the individual blocks!
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                                    // save some cycles for zero additions, but
02038                                    // only if it is safe for the matrix we are
02039                                    // working with
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                                  // This is a very messy function, since
02115                                  // we need to calculate to each position
02116                                  // the location in the global array.
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                                    // TODO: Look over this to find out
02131                                    // whether we can do that more
02132                                    // efficiently.
02133   if (col_indices_are_sorted == true)
02134     {
02135 #ifdef DEBUG
02136                                        // check whether indices really are
02137                                        // sorted.
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                                    // Resize scratch arrays
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                                    // Resize sub-arrays to n_cols. This
02190                                    // is a bit wasteful, but we resize
02191                                    // only a few times (then the maximum
02192                                    // row length won't increase that
02193                                    // much any more). At least we know
02194                                    // that all arrays are going to be of
02195                                    // the same size, so we can check
02196                                    // whether the size of one is large
02197                                    // enough before actually going
02198                                    // through all of them.
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                                    // Reset the number of added elements
02209                                    // in each block to zero.
02210   for (unsigned int i=0; i<this->n_block_cols(); ++i)
02211     counter_within_block[i] = 0;
02212 
02213                                    // Go through the column indices to
02214                                    // find out which portions of the
02215                                    // values should be written into
02216                                    // which block of the matrix. We need
02217                                    // to touch all the data, since we
02218                                    // can't be sure that the data of one
02219                                    // block is stored contiguously (in
02220                                    // fact, data will be intermixed when
02221                                    // it comes from an element matrix).
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                                    // If in debug mode, do a check whether
02240                                    // the right length has been obtained.
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                                    // Now we found out about where the
02248                                    // individual columns should start and
02249                                    // where we should start reading out
02250                                    // data. Now let's write the data into
02251                                    // the individual blocks!
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                                    // in block notation, the residual is
02651                                    // r_i = b_i - \sum_j A_ij x_j.
02652                                    // this can be written as
02653                                    // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
02654                                    //
02655                                    // for the first two terms, we can
02656                                    // call the residual function of
02657                                    // A_i0. for the other terms, we
02658                                    // use vmult_add. however, we want
02659                                    // to subtract, so in order to
02660                                    // avoid a temporary vector, we
02661                                    // perform a sign change of the
02662                                    // first two term before, and after
02663                                    // adding up
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                                    // first find out the row sizes
02781                                    // from the first block column
02782   for (unsigned int r=0; r<this->n_block_rows(); ++r)
02783     row_sizes[r] = sub_objects[r][0]->m();
02784                                    // then check that the following
02785                                    // block columns have the same
02786                                    // sizes
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                                    // finally initialize the row
02793                                    // indices with this array
02794   this->row_block_indices.reinit (row_sizes);
02795 
02796 
02797                                    // then do the same with the columns
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                                    // finally initialize the row
02806                                    // indices with this array
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Thu May 17 2012 20:04:39 by doxygen 1.7.3