]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new generalized iterator for block matrices
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 20 Feb 2006 17:15:59 +0000 (17:15 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 20 Feb 2006 17:15:59 +0000 (17:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@12429 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/block_matrix_base.h
deal.II/lac/include/lac/matrix_iterator.h [new file with mode: 0644]

index 2beb4499b35b10bb7aa8d56ab3b33302210fbc73..a99a311a8f4d4684f4e4b8d30d81befb8071d1b6 100644 (file)
@@ -19,6 +19,7 @@
 #include <base/smartpointer.h>
 #include <lac/block_indices.h>
 #include <lac/exceptions.h>
+#include <lac/matrix_iterator.h>
 
 #include <cmath>
 
  *@{
  */
 
-                                /**
-                                 * @cond internal
-                                 */
-namespace internal
-{
-
 /**
  * Namespace in which iterators in block matrices are implemented.
  *
  * @author Wolfgang Bangerth, 2004
  */
-  namespace BlockMatrixIterators
+namespace BlockMatrixIterators
+{
+                                  /**
+                                   * Base class for block matrix
+                                   * accessors, implementing the
+                                   * stepping through a matrix.
+                                   */
+  template <class BlockMatrix>
+  class AccessorBase
   {
-    template <typename> class ConstIterator;
-    
-                                    /**
-                                     * Accessor class for iterators
-                                     */
-    template <class BlockMatrix>
-    class Accessor
-    {
-      public:
-                                         /**
-                                          * Typedef the value type of the
-                                          * matrix we point into.
-                                          */
-        typedef typename BlockMatrix::value_type value_type;
-        
-                                        /**
-                                         * Constructor. Since we use
-                                         * accessors only for read
-                                         * access, a const matrix
-                                         * pointer is sufficient.
-                                         *
-                                         * Place the iterator at the
-                                         * beginning of the given row of the
-                                         * matrix, or create the end pointer
-                                         * if @p row equals the total number
-                                         * of rows in the matrix.
-                                         */
-       Accessor (const BlockMatrix  *m,
-                 const unsigned int  row);
-       
-                                        /**
-                                         * Row number of the element
-                                         * represented by this
-                                         * object.
-                                         */
-       unsigned int row() const;
-       
-                                        /**
-                                         * Column number of the
-                                         * element represented by
-                                         * this object.
-                                         */
-       unsigned int column() const;
-       
-                                        /**
-                                         * Value of this matrix entry.
-                                         */
-       value_type value() const;
-
-                                        /**
-                                         * Block row of the
-                                         * element represented by
-                                         * this object.
-                                         */
-       unsigned int block_row() const;
-       
-                                        /**
-                                         * Block column of the
-                                         * element represented by
-                                         * this object.
-                                         */
-       unsigned int block_column() const;
-        
-      protected:
-                                        /**
-                                         * The matrix accessed.
-                                         */
-       const BlockMatrix * matrix;
-       
-                                        /**
-                                         * Iterator of the underlying matrix
-                                         * class.
-                                         */
-       typename BlockMatrix::BlockType::const_iterator base_iterator;
-       
-                                        /**
-                                         * Block row into which we presently
-                                         * point.
-                                         */
-       unsigned int row_block;
-       
-                                        /**
-                                         * Block column into which we
-                                         * presently point.
-                                         */
-       unsigned int col_block;
-
-                                         /**
-                                          * Move ahead one element.
-                                          */
-        void advance ();
-
-                                         /**
-                                          * Compare this accessor with another
-                                          * one for equality.
-                                          */
-        bool operator == (const Accessor &a) const;
-        
-                                         /**
-                                          * Let the iterator class be a
-                                          * friend.
-                                          */
-        template <typename>
-        friend class ConstIterator;
-    };
-    
-                                    /**
-                                     * STL conforming iterator.
-                                     */
-    template <class BlockMatrix>
-    class ConstIterator
-    {
-      public:
-                                         /**
-                                          * Constructor.
-                                         *
-                                         * Place the iterator at the
-                                         * beginning of the given row of the
-                                         * matrix, or create the end pointer
-                                         * if @p row equals the total number
-                                         * of rows in the matrix.
-                                          */ 
-       ConstIterator(const BlockMatrix  *matrix,
-                      const unsigned int  row);
-         
-                                         /**
-                                          * Prefix increment.
-                                          */
-       ConstIterator& operator++ ();
-
-                                         /**
-                                          * Postfix increment.
-                                          */
-       ConstIterator& operator++ (int);
-
-                                         /**
-                                          * Dereferencing operator.
-                                          */
-       const Accessor<BlockMatrix> & operator* () const;
-
-                                         /**
-                                          * Dereferencing operator.
-                                          */
-       const Accessor<BlockMatrix> * operator-> () const;
-
-                                         /**
-                                          * Comparison. True, if
-                                          * both iterators point to
-                                          * the same matrix
-                                          * position.
-                                          */
-       bool operator == (const ConstIterator&) const;
-                                         /**
-                                          * Inverse of operator==().
-                                          */
-       bool operator != (const ConstIterator&) const;
-
-                                         /**
-                                          * Comparison
-                                          * operator. Result is true
-                                          * if either the first row
-                                          * number is smaller or if
-                                          * the row numbers are
-                                          * equal and the first
-                                          * index is smaller.
-                                          */
-       bool operator < (const ConstIterator&) const;
-
-      private:
-                                         /**
-                                          * The accessor with which we work.
-                                          */
-        Accessor<BlockMatrix> accessor;
-    };
-  }
+    public:
+                                      /**
+                                       * Typedef the value type of the
+                                       * matrix we point into.
+                                       */
+      typedef typename BlockMatrix::value_type value_type;
+
+                                      /**
+                                       * Initialize data fields to
+                                       * default values.
+                                       */
+      AccessorBase ();
+      
+                                      /**
+                                       * Block row of the
+                                       * element represented by
+                                       * this object.
+                                       */
+      unsigned int block_row() const;
+      
+                                      /**
+                                       * Block column of the
+                                       * element represented by
+                                       * this object.
+                                       */
+      unsigned int block_column() const;
+      
+    protected:
+                                      /**
+                                       * Block row into which we presently
+                                       * point.
+                                       */
+      unsigned int row_block;
+      
+                                      /**
+                                       * Block column into which we
+                                       * presently point.
+                                       */
+      unsigned int col_block;
+      
+                                      /**
+                                       * Let the iterator class be a
+                                       * friend.
+                                       */
+      template <typename>
+      friend class MatrixIterator;
+  };
+                                  /**
+                                   * Accessor classes in
+                                   * block matrices.
+                                   */
+  template <class BlockMatrix, bool ConstNess>
+  class Accessor;
+  
+                                  /**
+                                   * Block matrix accessor for non
+                                   * const matrices.
+                                   */
+  template <class BlockMatrix>
+  class Accessor<BlockMatrix, false>
+    :
+      public AccessorBase<BlockMatrix>
+  {
+    public:
+                                      /**
+                                       * Type of the matrix used in
+                                       * this accessor.
+                                       */
+      typedef BlockMatrix MatrixType;
+      
+                                      /**
+                                       * Typedef the value type of the
+                                       * matrix we point into.
+                                       */
+      typedef typename BlockMatrix::value_type value_type;
+      
+                                      /**
+                                       * Constructor. Since we use
+                                       * accessors only for read
+                                       * access, a const matrix
+                                       * pointer is sufficient.
+                                       *
+                                       * Place the iterator at the
+                                       * beginning of the given row of the
+                                       * matrix, or create the end pointer
+                                       * if @p row equals the total number
+                                       * of rows in the matrix.
+                                       */
+      Accessor (BlockMatrix* m,
+               const unsigned int row,
+               const unsigned int col);
+      
+                                      /**
+                                       * Row number of the element
+                                       * represented by this
+                                       * object.
+                                       */
+      unsigned int row() const;
+      
+                                      /**
+                                       * Column number of the
+                                       * element represented by
+                                       * this object.
+                                       */
+      unsigned int column() const;
+      
+                                      /**
+                                       * Value of the entry at the
+                                       * current position.
+                                       */
+      value_type value() const;
+      
+                                      /**
+                                       * Set new value.
+                                       */
+      void set_value(value_type newval) const;
+      
+    protected:
+                                      /**
+                                       * The matrix accessed.
+                                       */
+      BlockMatrix * matrix;
+      
+                                      /**
+                                       * Iterator of the underlying matrix
+                                       * class.
+                                       */
+      typename BlockMatrix::BlockType::iterator base_iterator;
+      
+                                      /**
+                                       * Move ahead one element.
+                                       */
+      void advance ();
+      
+                                      /**
+                                       * Compare this accessor with another
+                                       * one for equality.
+                                       */
+      bool operator == (const Accessor &a) const;
+      
+      template <typename> friend class MatrixIterator;
+      friend class Accessor<BlockMatrix, true>;
+  };
+  
+                                  /**
+                                   * Block matrix accessor for
+                                   * constant matrices, implementing
+                                   * the stepping through a matrix.
+                                   */
+  template <class BlockMatrix>
+  class Accessor<BlockMatrix, true>
+    :
+      public AccessorBase<BlockMatrix>
+  {
+    public:
+                                      /**
+                                       * Thype of the matrix used in
+                                       * this accessor.
+                                       */
+      typedef const BlockMatrix MatrixType;
+      
+                                      /**
+                                       * Typedef the value type of the
+                                       * matrix we point into.
+                                       */
+      typedef typename BlockMatrix::value_type value_type;
+      
+                                      /**
+                                       * Constructor. Since we use
+                                       * accessors only for read
+                                       * access, a const matrix
+                                       * pointer is sufficient.
+                                       *
+                                       * Place the iterator at the
+                                       * beginning of the given row of the
+                                       * matrix, or create the end pointer
+                                       * if @p row equals the total number
+                                       * of rows in the matrix.
+                                       */
+      Accessor (const BlockMatrix* m,
+               const unsigned int row,
+               const unsigned int col);
+      
+                                      /**
+                                       * Initalize const accessor
+                                       * from non const accessor.
+                                       */
+      Accessor(const Accessor<BlockMatrix, false>&);
+      
+                                      /**
+                                       * Row number of the element
+                                       * represented by this
+                                       * object.
+                                       */
+      unsigned int row() const;
+      
+                                      /**
+                                       * Column number of the
+                                       * element represented by
+                                       * this object.
+                                       */
+      unsigned int column() const;
+      
+                                      /**
+                                       * Value of the entry at the
+                                       * current position.
+                                       */
+      value_type value() const;
+    protected:
+                                      /**
+                                       * The matrix accessed.
+                                       */
+      const BlockMatrix * matrix;
+      
+                                      /**
+                                       * Iterator of the underlying matrix
+                                       * class.
+                                       */
+      typename BlockMatrix::BlockType::const_iterator base_iterator;
+      
+                                      /**
+                                       * Move ahead one element.
+                                       */
+      void advance ();
+      
+                                      /**
+                                       * Compare this accessor with another
+                                       * one for equality.
+                                       */
+      bool operator == (const Accessor &a) const;
+      
+                                      /**
+                                       * Let the iterator class be a
+                                       * friend.
+                                       */
+      template <typename>
+      friend class MatrixIterator;
+  };
 }
-                                   /**
-                                    * @endcond
-                                    */
+
 
 
 /**
@@ -302,11 +366,11 @@ class BlockMatrixBase : public Subscriptor
     typedef size_t                  size_type;
 
     typedef
-    internal::BlockMatrixIterators::ConstIterator<BlockMatrixBase>
+    MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, false> >
     iterator;
 
     typedef
-    internal::BlockMatrixIterators::ConstIterator<BlockMatrixBase>
+    MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, true> >
     const_iterator;
 
 
@@ -604,6 +668,27 @@ class BlockMatrixBase : public Subscriptor
                         const BlockVectorType &x,
                         const BlockVectorType &b) const;
 
+                                    /**
+                                     * STL-like iterator with the
+                                     * first entry.
+                                     */
+    iterator begin ();
+
+                                    /**
+                                     * Final iterator.
+                                     */
+    iterator end ();
+    
+                                    /**
+                                     * STL-like iterator with the
+                                     * first entry of row <tt>r</tt>.
+                                     */
+    iterator begin (const unsigned int r);
+
+                                    /**
+                                     * Final iterator of row <tt>r</tt>.
+                                     */
+    iterator end (const unsigned int r);
                                     /**
                                      * STL-like iterator with the
                                      * first entry.
@@ -900,14 +985,15 @@ class BlockMatrixBase : public Subscriptor
                                      * a compiler bug here again.
                                      */
 #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG
-    template <typename>
-    friend class internal::BlockMatrixIterators::Accessor;
+    template <typename, bool>
+    friend class BlockMatrixIterators::Accessor;
 
     template <typename>
-    friend class internal::BlockMatrixIterators::ConstIterator;
+    friend class MatrixIterator;
 #else
-    typedef internal::BlockMatrixIterators::Accessor<BlockMatrixBase> Accessor;
-    friend class Accessor;
+    typedef BlockMatrixIterators::Accessor<BlockMatrixBase, true> ConstAccessor;
+    typedef BlockMatrixIterators::Accessor<BlockMatrixBase, false> Accessor;
+    friend class ConstAccessor;
     
     friend class const_iterator;
 #endif    
@@ -920,316 +1006,449 @@ class BlockMatrixBase : public Subscriptor
 /* ------------------------- Template functions ---------------------- */
 
 
-namespace internal
+namespace BlockMatrixIterators
 {
-  namespace BlockMatrixIterators
-  {
-    template <class BlockMatrix>
-    inline
-    Accessor<BlockMatrix>::
-    Accessor (const BlockMatrix  *matrix,
-              const unsigned int  row)
-                    :
-                    matrix(matrix),
-                    base_iterator(matrix->block(0,0).begin()),
-                    row_block(0),
-                    col_block(0)
-    {
-                                       // check if this is a regular row or
-                                       // the end of the matrix
-      if (row < matrix->m())
-        {
-          const std::pair<unsigned int,unsigned int> indices
-            = matrix->row_block_indices.global_to_local(row);
-
-                                           // find the first block that does
-                                           // have an entry in this row
-          for (unsigned int bc=0; bc<matrix->n_block_cols(); ++bc)
-            {
-              base_iterator
-                = matrix->block(indices.first, bc).begin(indices.second);
-              if (base_iterator !=
-                  matrix->block(indices.first, bc).end(indices.second))
-                {
-                  row_block = indices.first;
-                  col_block = bc;
-                  return;
-                }
-            }
-
-                                           // hm, there is no block that has
-                                           // an entry in this column. we need
-                                           // to take the next entry then,
-                                           // which may be the first entry of
-                                           // the next row, or recursively the
-                                           // next row, or so on
-          *this = Accessor (matrix, row+1);
-        }
-      else
-        {
-                                           // we were asked to create the end
-                                           // iterator for this matrix
-          row_block = deal_II_numbers::invalid_unsigned_int;
-          col_block = deal_II_numbers::invalid_unsigned_int;
-        }
-    }
-
-
-    template <class BlockMatrix>
-    inline
-    unsigned int
-    Accessor<BlockMatrix>::row() const
-    {
-      Assert (row_block != deal_II_numbers::invalid_unsigned_int,
-              ExcIteratorPastEnd());
-      
-      return (matrix->row_block_indices.local_to_global(row_block, 0) +
-              base_iterator->row());
-    }
-
-
-    template <class BlockMatrix>
-    inline
-    unsigned int
-    Accessor<BlockMatrix>::column() const
-    {
-      Assert (col_block != deal_II_numbers::invalid_unsigned_int,
-              ExcIteratorPastEnd());
-
-      return (matrix->column_block_indices.local_to_global(col_block,0) +
-              base_iterator->column());
-    }
-
-
-    template <class BlockMatrix>
-    inline
-    unsigned int
-    Accessor<BlockMatrix>::block_row() const
-    {
-      Assert (row_block != deal_II_numbers::invalid_unsigned_int,
-              ExcIteratorPastEnd());
-
-      return row_block;
-    }
-
-
-    template <class BlockMatrix>
-    inline
-    unsigned int
-    Accessor<BlockMatrix>::block_column() const
-    {
-      Assert (col_block != deal_II_numbers::invalid_unsigned_int,
-              ExcIteratorPastEnd());
-
-      return col_block;
-    }
-
-
-    template <class BlockMatrix>
-    inline
-    typename Accessor<BlockMatrix>::value_type
-    Accessor<BlockMatrix>::value () const
-    {
-      Assert (row_block != deal_II_numbers::invalid_unsigned_int,
-              ExcIteratorPastEnd());
-      Assert (col_block != deal_II_numbers::invalid_unsigned_int,
-              ExcIteratorPastEnd());
-
-      return base_iterator->value();
-    }
-
+  template <class BlockMatrix>
+  inline
+  AccessorBase<BlockMatrix>::AccessorBase()
+                 :
+                 row_block(0),
+                 col_block(0)
+  {}
+  
 
+  template <class BlockMatrix>
+  inline
+  unsigned int
+  AccessorBase<BlockMatrix>::block_row() const
+  {
+    Assert (row_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+    return row_block;
+  }
+  
+  
+  template <class BlockMatrix>
+  inline
+  unsigned int
+  AccessorBase<BlockMatrix>::block_column() const
+  {
+    Assert (col_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+    return col_block;
+  }
+  
+  
+  template <class BlockMatrix>
+  inline
+  Accessor<BlockMatrix, true>::Accessor (
+    const BlockMatrix  *matrix,
+    const unsigned int  row,
+    const unsigned int  col)
+                 :
+                 matrix(matrix),
+                 base_iterator(matrix->block(0,0).begin())
+  {
+    Assert(col==0, ExcNotImplemented());
+    
+                                    // check if this is a regular row or
+                                    // the end of the matrix
+    if (row < matrix->m())
+      {
+       const std::pair<unsigned int,unsigned int> indices
+         = matrix->row_block_indices.global_to_local(row);
+       
+                                        // find the first block that does
+                                        // have an entry in this row
+       for (unsigned int bc=0; bc<matrix->n_block_cols(); ++bc)
+         {
+           base_iterator
+             = matrix->block(indices.first, bc).begin(indices.second);
+           if (base_iterator !=
+               matrix->block(indices.first, bc).end(indices.second))
+             {
+               this->row_block = indices.first;
+               this->col_block = bc;
+               return;
+             }
+         }
+       
+                                        // hm, there is no block that has
+                                        // an entry in this column. we need
+                                        // to take the next entry then,
+                                        // which may be the first entry of
+                                        // the next row, or recursively the
+                                        // next row, or so on
+       *this = Accessor (matrix, row+1, 0);
+      }
+    else
+      {
+                                        // we were asked to create the end
+                                        // iterator for this matrix
+       this->row_block = deal_II_numbers::invalid_unsigned_int;
+       this->col_block = deal_II_numbers::invalid_unsigned_int;
+      }
+  }
+  
 
-    template <class BlockMatrix>
-    inline
-    void
-    Accessor<BlockMatrix>::advance ()
-    {      
-      Assert (row_block != deal_II_numbers::invalid_unsigned_int,
-              ExcIteratorPastEnd());
-      Assert (col_block != deal_II_numbers::invalid_unsigned_int,
-              ExcIteratorPastEnd());
+//   template <class BlockMatrix>
+//   inline
+//   Accessor<BlockMatrix, true>::Accessor (const Accessor<BlockMatrix, true>& other)
+//               :
+//               matrix(other.matrix),
+//               base_iterator(other.base_iterator)
+//   {
+//     this->row_block = other.row_block;
+//     this->col_block = other.col_block;
+//   }
+  
 
-                                       // Remember current row inside block
-      unsigned int local_row = base_iterator->row();
+  template <class BlockMatrix>
+  inline
+  Accessor<BlockMatrix, true>::Accessor (const Accessor<BlockMatrix, false>& other)
+                 :
+                 matrix(other.matrix),
+                 base_iterator(other.base_iterator)
+  {
+     this->row_block = other.row_block;
+     this->col_block = other.col_block;    
+  }
+  
 
-                                       // Advance one element inside the
-                                       // current block
-      ++base_iterator;
+  template <class BlockMatrix>
+  inline
+  unsigned int
+  Accessor<BlockMatrix, true>::row() const
+  {
+    Assert (this->row_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+    return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
+           base_iterator->row());
+  }
   
-                                       // while we hit the end of the row of a
-                                       // block (which may happen multiple
-                                       // times if rows inside a block are
-                                       // empty), we have to jump to the next
-                                       // block and take the
-      while (base_iterator ==
-             matrix->block(row_block,
-                           col_block).end(local_row))
-        {
-                                           // jump to next block in this block
-                                           // row, if possible, otherwise go
-                                           // to next row
-          if (col_block < matrix->n_block_cols()-1)
-            {
-              ++col_block;
-              base_iterator
-                = matrix->block(row_block, col_block).begin(local_row);
-            }
-          else
-            {
-                                               // jump back to next row in
-                                               // first block column
-              col_block = 0;
-              ++local_row;
-
-                                               // see if this has brought us
-                                               // past the number of rows in
-                                               // this block. if so see
-                                               // whether we've just fallen
-                                               // off the end of the whole
-                                               // matrix
-              if (local_row == matrix->block(row_block, col_block).m())
-                {
-                  local_row = 0;
-                  ++row_block;
-                  if (row_block == matrix->n_block_rows())
-                    {
-                      row_block = deal_II_numbers::invalid_unsigned_int;
-                      col_block = deal_II_numbers::invalid_unsigned_int;
-                      return;
-                    }
-                }
+  
+  template <class BlockMatrix>
+  inline
+  unsigned int
+  Accessor<BlockMatrix, true>::column() const
+  {
+    Assert (this->col_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+    return (matrix->column_block_indices.local_to_global(this->col_block,0) +
+           base_iterator->column());
+  }
+  
+  
+  template <class BlockMatrix>
+  inline
+  typename Accessor<BlockMatrix, true>::value_type
+  Accessor<BlockMatrix, true>::value () const
+  {
+    Assert (this->row_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    Assert (this->col_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+    return base_iterator->value();
+  }
+  
+  
+  
+  template <class BlockMatrix>
+  inline
+  void
+  Accessor<BlockMatrix, true>::advance ()
+  {      
+    Assert (this->row_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    Assert (this->col_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+                                    // Remember current row inside block
+    unsigned int local_row = base_iterator->row();
+    
+                                    // Advance one element inside the
+                                    // current block
+    ++base_iterator;
+    
+                                    // while we hit the end of the row of a
+                                    // block (which may happen multiple
+                                    // times if rows inside a block are
+                                    // empty), we have to jump to the next
+                                    // block and take the
+    while (base_iterator ==
+          matrix->block(this->row_block, this->col_block).end(local_row))
+      {
+                                        // jump to next block in this block
+                                        // row, if possible, otherwise go
+                                        // to next row
+       if (this->col_block < matrix->n_block_cols()-1)
+         {
+           ++this->col_block;
+           base_iterator
+             = matrix->block(this->row_block, this->col_block).begin(local_row);
+         }
+       else
+         {
+                                            // jump back to next row in
+                                            // first block column
+           this->col_block = 0;
+           ++local_row;
+
+                                            // see if this has brought us
+                                            // past the number of rows in
+                                            // this block. if so see
+                                            // whether we've just fallen
+                                            // off the end of the whole
+                                            // matrix
+           if (local_row == matrix->block(this->row_block, this->col_block).m())
+             {
+               local_row = 0;
+               ++this->row_block;
+               if (this->row_block == matrix->n_block_rows())
+                 {
+                   this->row_block = deal_II_numbers::invalid_unsigned_int;
+                   this->col_block = deal_II_numbers::invalid_unsigned_int;
+                   return;
+                 }
+             }
               
-              base_iterator
-                = matrix->block(row_block, col_block).begin(local_row);
-            }
-        }
-    }
+           base_iterator
+             = matrix->block(this->row_block, this->col_block).begin(local_row);
+         }
+      }
+  }
     
-      
-
-    template <class BlockMatrix>
-    inline
-    bool
-    Accessor<BlockMatrix>::operator == (const Accessor &a) const
-    {      
-      if (matrix != a.matrix)
-        return false;
-
-      if (row_block == a.row_block
-          && col_block == a.col_block)
-                                         // end iterators do not necessarily
-                                         // have to have the same
-                                         // base_iterator representation, but
-                                         // valid iterators have to
-        return (((row_block == deal_II_numbers::invalid_unsigned_int)
-                 &&
-                 (col_block == deal_II_numbers::invalid_unsigned_int))
-                ||
-                (base_iterator == a.base_iterator));
 
+  template <class BlockMatrix>
+  inline
+  bool
+  Accessor<BlockMatrix, true>::operator == (const Accessor &a) const
+  {      
+    if (matrix != a.matrix)
       return false;
-    }
-    
-
-//---------------------------------------------------------------------------
-
-
-    template <class BlockMatrix>
-    inline
-    ConstIterator<BlockMatrix>::
-    ConstIterator (const BlockMatrix *m,
-                   const unsigned int r)
-                    :
-                    accessor (m, r)
-    {}
-
-
-
-    template <class BlockMatrix>
-    inline
-    ConstIterator<BlockMatrix> &
-    ConstIterator<BlockMatrix>::operator++ ()
-    {
-      accessor.advance ();
-      return *this;
-    }
-
-
-
-//  template <class BlockMatrix>
-//  inline
-//  const_iterator&
-//  ConstIterator::operator++ (int)
-//  {
-//    Assert (false, ExcNotImplemented());
-//  }
-
-
-
-
-    template <class BlockMatrix>
-    inline
-    const Accessor<BlockMatrix> &
-    ConstIterator<BlockMatrix>::operator* () const
-    {
-      return accessor;
-    }
-
-
-    template <class BlockMatrix>
-    inline
-    const Accessor<BlockMatrix> *
-    ConstIterator<BlockMatrix>::operator-> () const
-    {
-      return &accessor;
-    }
-
-
-
-    template <class BlockMatrix>
-    inline
-    bool
-    ConstIterator<BlockMatrix>::
-    operator == (const ConstIterator& i) const
-    {
-      return (accessor == i.accessor);
-    }
-
-
-
-    template <class BlockMatrix>
-    inline
-    bool
-    ConstIterator<BlockMatrix>::
-    operator != (const ConstIterator& i) const
-    {
-      return !(*this == i);
-    }
 
+    if (this->row_block == a.row_block
+       && this->col_block == a.col_block)
+                                      // end iterators do not necessarily
+                                      // have to have the same
+                                      // base_iterator representation, but
+                                      // valid iterators have to
+      return (((this->row_block == deal_II_numbers::invalid_unsigned_int)
+              &&
+              (this->col_block == deal_II_numbers::invalid_unsigned_int))
+             ||
+             (base_iterator == a.base_iterator));
+
+    return false;
+  }
+    
+//----------------------------------------------------------------------//      
 
+  
+  template <class BlockMatrix>
+  inline
+  Accessor<BlockMatrix, false>::Accessor (
+    BlockMatrix  *matrix,
+    const unsigned int  row,
+    const unsigned int  col)
+                 :
+                 matrix(matrix),
+                 base_iterator(matrix->block(0,0).begin())
+  {
+    Assert(col==0, ExcNotImplemented());
+                                    // check if this is a regular row or
+                                    // the end of the matrix
+    if (row < matrix->m())
+      {
+       const std::pair<unsigned int,unsigned int> indices
+         = matrix->row_block_indices.global_to_local(row);
+       
+                                        // find the first block that does
+                                        // have an entry in this row
+       for (unsigned int bc=0; bc<matrix->n_block_cols(); ++bc)
+         {
+           base_iterator
+             = matrix->block(indices.first, bc).begin(indices.second);
+           if (base_iterator !=
+               matrix->block(indices.first, bc).end(indices.second))
+             {
+               this->row_block = indices.first;
+               this->col_block = bc;
+               return;
+             }
+         }
+       
+                                        // hm, there is no block that has
+                                        // an entry in this column. we need
+                                        // to take the next entry then,
+                                        // which may be the first entry of
+                                        // the next row, or recursively the
+                                        // next row, or so on
+       *this = Accessor (matrix, row+1, 0);
+      }
+    else
+      {
+                                        // we were asked to create the end
+                                        // iterator for this matrix
+       this->row_block = deal_II_numbers::invalid_unsigned_int;
+       this->col_block = deal_II_numbers::invalid_unsigned_int;
+      }
+  }
+  
+  
+  template <class BlockMatrix>
+  inline
+  unsigned int
+  Accessor<BlockMatrix, false>::row() const
+  {
+    Assert (this->row_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+    return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
+           base_iterator->row());
+  }
+  
+  
+  template <class BlockMatrix>
+  inline
+  unsigned int
+  Accessor<BlockMatrix, false>::column() const
+  {
+    Assert (this->col_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+    return (matrix->column_block_indices.local_to_global(this->col_block,0) +
+           base_iterator->column());
+  }
+  
+  
+  template <class BlockMatrix>
+  inline
+  typename Accessor<BlockMatrix, false>::value_type
+  Accessor<BlockMatrix, false>::value () const
+  {
+    Assert (this->row_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    Assert (this->col_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+    return base_iterator->value();
+  }
+  
+  
+  
+  template <class BlockMatrix>
+  inline
+  void
+  Accessor<BlockMatrix, false>::set_value (typename Accessor<BlockMatrix, false>::value_type newval) const
+  {
+    Assert (this->row_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    Assert (this->col_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+    base_iterator->value() = newval;
+  }
+  
+  
+  
+  template <class BlockMatrix>
+  inline
+  void
+  Accessor<BlockMatrix, false>::advance ()
+  {      
+    Assert (this->row_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    Assert (this->col_block != deal_II_numbers::invalid_unsigned_int,
+           ExcIteratorPastEnd());
+    
+                                    // Remember current row inside block
+    unsigned int local_row = base_iterator->row();
+    
+                                    // Advance one element inside the
+                                    // current block
+    ++base_iterator;
+    
+                                    // while we hit the end of the row of a
+                                    // block (which may happen multiple
+                                    // times if rows inside a block are
+                                    // empty), we have to jump to the next
+                                    // block and take the
+    while (base_iterator ==
+          matrix->block(this->row_block, this->col_block).end(local_row))
+      {
+                                        // jump to next block in this block
+                                        // row, if possible, otherwise go
+                                        // to next row
+       if (this->col_block < matrix->n_block_cols()-1)
+         {
+           ++this->col_block;
+           base_iterator
+             = matrix->block(this->row_block, this->col_block).begin(local_row);
+         }
+       else
+         {
+                                            // jump back to next row in
+                                            // first block column
+           this->col_block = 0;
+           ++local_row;
+
+                                            // see if this has brought us
+                                            // past the number of rows in
+                                            // this block. if so see
+                                            // whether we've just fallen
+                                            // off the end of the whole
+                                            // matrix
+           if (local_row == matrix->block(this->row_block, this->col_block).m())
+             {
+               local_row = 0;
+               ++this->row_block;
+               if (this->row_block == matrix->n_block_rows())
+                 {
+                   this->row_block = deal_II_numbers::invalid_unsigned_int;
+                   this->col_block = deal_II_numbers::invalid_unsigned_int;
+                   return;
+                 }
+             }
+              
+           base_iterator
+             = matrix->block(this->row_block, this->col_block).begin(local_row);
+         }
+      }
+  }
+    
+      
 
-    template <class BlockMatrix>
-    inline
-    bool
-    ConstIterator<BlockMatrix>::
-    operator < (const ConstIterator& i) const
-    {
-      if (accessor.row_block < i.accessor.row_block)
-        return true;
-      if (accessor.row_block == i.accessor.row_block)
-        {
-          if (accessor.base_iterator->row() < i.accessor.base_iterator->row())
-            return true;
-          if (accessor.base_iterator->row() == i.accessor.base_iterator->row())
-            return (accessor.base_iterator < i.accessor.base_iterator);
-        }
+  template <class BlockMatrix>
+  inline
+  bool
+  Accessor<BlockMatrix, false>::operator == (const Accessor &a) const
+  {      
+    if (matrix != a.matrix)
       return false;
-    }
 
+    if (this->row_block == a.row_block
+       && this->col_block == a.col_block)
+                                      // end iterators do not necessarily
+                                      // have to have the same
+                                      // base_iterator representation, but
+                                      // valid iterators have to
+      return (((this->row_block == deal_II_numbers::invalid_unsigned_int)
+              &&
+              (this->col_block == deal_II_numbers::invalid_unsigned_int))
+             ||
+             (base_iterator == a.base_iterator));
+
+    return false;
   }
 }
 
+
 //---------------------------------------------------------------------------
 
 
@@ -1843,6 +2062,48 @@ BlockMatrixBase<MatrixType>::end (const unsigned int r) const
 
 
 
+template <class MatrixType>
+inline
+typename BlockMatrixBase<MatrixType>::iterator
+BlockMatrixBase<MatrixType>::begin ()
+{
+  return iterator(this, 0);
+}
+
+
+
+template <class MatrixType>
+inline
+typename BlockMatrixBase<MatrixType>::iterator
+BlockMatrixBase<MatrixType>::end ()
+{
+  return iterator(this, m());
+}
+
+
+
+template <class MatrixType>
+inline
+typename BlockMatrixBase<MatrixType>::iterator
+BlockMatrixBase<MatrixType>::begin (const unsigned int r)
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return iterator(this, r);
+}
+
+
+
+template <class MatrixType>
+inline
+typename BlockMatrixBase<MatrixType>::iterator
+BlockMatrixBase<MatrixType>::end (const unsigned int r)
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return iterator(this, r+1);
+}
+
+
+
 template <class MatrixType>
 void
 BlockMatrixBase<MatrixType>::collect_sizes ()
diff --git a/deal.II/lac/include/lac/matrix_iterator.h b/deal.II/lac/include/lac/matrix_iterator.h
new file mode 100644 (file)
index 0000000..55a43e1
--- /dev/null
@@ -0,0 +1,236 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__matrix_iterator_h
+#define __deal2__matrix_iterator_h
+
+
+#include <base/config.h>
+
+                                /**
+                                 * STL conforming iterator for constant
+                                 * and non-constant matrices.
+                                 *
+                                 * This iterator is abstracted from
+                                 * the actual matrix type and can be
+                                 * used for any matrix having the
+                                 * required ACCESSOR type.
+                                 *
+                                 * @author Guido Kanschat, 2006, based on previous a implementation
+                                 */
+template <class ACCESSOR>
+class MatrixIterator
+{
+  public:
+                                    /**
+                                     * Typedef for the matrix type
+                                     * (including constness) we are to
+                                     * operate on.
+                                     */
+    typedef typename ACCESSOR::MatrixType MatrixType;
+    
+                                    /**
+                                     * Constructor. Create an
+                                     * iterator into the matrix
+                                     * <tt>matrix</tt> for the given
+                                     * <tt>row</tt> and the
+                                     * <tt>index</tt> within it.
+                                     */ 
+    MatrixIterator (MatrixType        *matrix,
+                   const unsigned int row = 0,
+                   const unsigned int index = 0);
+
+                                    /**
+                                     * Copy from another matrix
+                                     * iterator. Mostly implemented
+                                     * to allow initialization of a
+                                     * constant iterator from a non
+                                     * constant, this function only
+                                     * requires that a conversion
+                                     * from the other iterator's
+                                     * accessor to this accessor
+                                     * object is possible.
+                                     */
+    template <class OtherAccessor>
+    MatrixIterator(const MatrixIterator<OtherAccessor>& other);
+    
+                                    /**
+                                     * Prefix increment.
+                                     */
+    MatrixIterator & operator++ ();
+
+                                    /**
+                                     * Postfix increment.
+                                     */
+    MatrixIterator operator++ (int);
+
+                                    /**
+                                     * Dereferencing operator.
+                                     */
+    const ACCESSOR & operator* () const;
+
+                                    /**
+                                     * Dereferencing operator.
+                                     */
+    const ACCESSOR * operator-> () const;
+
+                                    /**
+                                     * Comparison. True, if
+                                     * both accessors are equal.
+                                     */
+    bool operator == (const MatrixIterator &) const;
+
+                                    /**
+                                     * Inverse of <tt>==</tt>.
+                                     */
+    bool operator != (const MatrixIterator &) const;
+
+                                    /**
+                                     * Comparison operator. Result is
+                                     * true if either the first row
+                                     * number is smaller or if the row
+                                     * numbers are equal and the first
+                                     * index is smaller.
+                                     *
+                                     * This function is only valid if
+                                     * both iterators point into the same
+                                     * matrix.
+                                     */
+    bool operator < (const MatrixIterator &) const;
+
+                                    /**
+                                     * Comparison operator. Works in the
+                                     * same way as above operator, just
+                                     * the other way round.
+                                     */
+    bool operator > (const MatrixIterator &) const;
+        
+  private:
+                                    /**
+                                     * Store an object of the
+                                     * accessor class.
+                                     */
+    ACCESSOR accessor;
+
+                                    /**
+                                     * Allow other iterators access
+                                     * to private data.
+                                     */
+    template <class OtherAccessor> friend class MatrixIterator;
+};
+
+
+//----------------------------------------------------------------------//
+
+template <class ACCESSOR>
+inline
+MatrixIterator<ACCESSOR>::
+MatrixIterator (MatrixType        *matrix,
+               const unsigned int r,
+               const unsigned int i)
+               :
+               accessor(matrix, r, i)
+{}
+
+
+template <class ACCESSOR>
+template <class OtherAccessor>
+inline
+MatrixIterator<ACCESSOR>::
+MatrixIterator (const MatrixIterator<OtherAccessor>& other)
+               :
+               accessor(other.accessor)
+{}
+
+
+template <class ACCESSOR>
+inline
+MatrixIterator<ACCESSOR> &
+MatrixIterator<ACCESSOR>::operator++ ()
+{
+  accessor.advance ();
+  return *this;
+}
+
+
+template <class ACCESSOR>
+inline
+MatrixIterator<ACCESSOR>
+MatrixIterator<ACCESSOR>::operator++ (int)
+{
+  const MatrixIterator iter = *this;
+  accessor.advance ();
+  return iter;
+}
+
+
+template <class ACCESSOR>
+inline
+const ACCESSOR &
+MatrixIterator<ACCESSOR>::operator* () const
+{
+  return accessor;
+}
+
+
+template <class ACCESSOR>
+inline
+const ACCESSOR *
+MatrixIterator<ACCESSOR>::operator-> () const
+{
+  return &accessor;
+}
+
+
+template <class ACCESSOR>
+inline
+bool
+MatrixIterator<ACCESSOR>::
+operator == (const MatrixIterator& other) const
+{
+  return (accessor == other.accessor);
+}
+
+
+template <class ACCESSOR>
+inline
+bool
+MatrixIterator<ACCESSOR>::
+operator != (const MatrixIterator& other) const
+{
+  return ! (*this == other);
+}
+
+
+template <class ACCESSOR>
+inline
+bool
+MatrixIterator<ACCESSOR>::
+operator < (const MatrixIterator& other) const
+{
+  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
+         ExcInternalError());
+      
+  return (accessor < other.accessor);
+}
+
+
+template <class ACCESSOR>
+inline
+bool
+MatrixIterator<ACCESSOR>::
+operator > (const MatrixIterator& other) const
+{
+  return (other < *this);
+}
+
+#endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.