]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement petsc matrix iterators.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 7 Mar 2004 22:31:45 +0000 (22:31 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 7 Mar 2004 22:31:45 +0000 (22:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@8670 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/petsc_matrix_base.h
deal.II/lac/include/lac/petsc_sparse_matrix.h
deal.II/lac/source/petsc_matrix_base.cc

index 223398e078430606e90807f834afb3c89cead61d..4878b4eb70b62509225a6fcdb76127f1b1ccdb04 100644 (file)
 #ifdef DEAL_II_USE_PETSC
 
 #include <petscmat.h>
+#include <boost/shared_ptr.hpp>
+#include <vector>
+
 
 
 namespace PETScWrappers
 {
-                                   // forward declaration
+                                   // forward declarations
   class VectorBase;
+  class MatrixBase;
+
+  namespace MatrixIterators
+  {
+/**
+ * STL conforming iterator. This class acts as an iterator walking over the
+ * elements of PETSc matrices. Since PETSc offers a uniform interface for all
+ * types of matrices, this iterator can be used to access both sparse and full
+ * matrices.
+ *
+ * Note that PETSc does not give any guarantees as to the order of elements
+ * within each row. Note also that accessing the elements of a full matrix
+ * surprisingly only shows the nonzero elements of the matrix, not all
+ * elements.
+ *
+ * @author Guido Kanschat, Roy Stogner, Wolfgang Bangerth, 2004
+ */
+    class const_iterator
+    {
+      private:
+                                         /**
+                                          * Accessor class for iterators
+                                          */
+        class Accessor
+        {
+          public:
+                                             /**
+                                              * Constructor. Since we use
+                                              * accessors only for read
+                                              * access, a const matrix
+                                              * pointer is sufficient.
+                                              */
+            Accessor (const MatrixBase    *matrix,
+                      const unsigned int   row,
+                      const unsigned short index);
+
+                                             /**
+                                              * Row number of the element
+                                              * represented by this
+                                              * object.
+                                              */
+            unsigned int row() const;
+
+                                             /**
+                                              * Index in row of the element
+                                              * represented by this
+                                              * object.
+                                              */
+            unsigned short index() const;
+
+                                             /**
+                                              * Column number of the
+                                              * element represented by
+                                              * this object.
+                                              */
+            unsigned int column() const;
+
+                                             /**
+                                              * Value of this matrix entry.
+                                              */
+            PetscScalar value() const;
+
+                                             /**
+                                              * Exception
+                                              */
+            DeclException0 (ExcBeyondEndOfMatrix);
+            
+          private:
+                                             /**
+                                              * The matrix accessed.
+                                              */
+            const MatrixBase *matrix;
+
+                                             /**
+                                              * Current row number.
+                                              */
+            unsigned int a_row;
+
+                                             /**
+                                              * Current index in row.
+                                              */
+            unsigned int a_index;
+
+                                             /**
+                                              * Cache where we store the
+                                              * column indices of the present
+                                              * row. This is necessary, since
+                                              * PETSc makes access to the
+                                              * elements of its matrices
+                                              * rather hard, and it is much
+                                              * more efficient to copy all
+                                              * column entries of a row once
+                                              * when we enter it than
+                                              * repeatedly asking PETSc for
+                                              * individual ones. This also
+                                              * makes some sense since it is
+                                              * likely that we will access
+                                              * them sequentially anyway.
+                                              *
+                                              * In order to make copying of
+                                              * iterators/accessor of
+                                              * acceptable performance, we
+                                              * keep a shared pointer to these
+                                              * entries so that more than one
+                                              * accessor can access this data
+                                              * if necessary.
+                                              */
+            boost::shared_ptr<const std::vector<unsigned int> > colnum_cache;
+
+                                             /**
+                                              * Similar cache for the values
+                                              * of this row.
+                                              */
+            boost::shared_ptr<const std::vector<PetscScalar> > value_cache;
+            
+                                             /**
+                                              * Discard the old row caches
+                                              * (they may still be used by
+                                              * other accessors) and generate
+                                              * new ones for the row pointed
+                                              * to presently by this accessor.
+                                              */
+            void visit_present_row ();
+
+                                             /**
+                                              * Make enclosing class a
+                                              * friend.
+                                              */
+            friend class const_iterator;
+        };
+        
+      public:
+          
+                                         /**
+                                          * Constructor. Create an iterator
+                                          * into the matrix @p matrix for the
+                                          * given row and the index within it.
+                                          */ 
+        const_iterator (const MatrixBase    *matrix,
+                        const unsigned int   row,
+                        const unsigned short index);
+          
+                                         /**
+                                          * Prefix increment.
+                                          */
+        const_iterator& operator++ ();
+
+                                         /**
+                                          * Postfix increment.
+                                          */
+        const_iterator operator++ (int);
+
+                                         /**
+                                          * Dereferencing operator.
+                                          */
+        const Accessor& operator* () const;
+
+                                         /**
+                                          * Dereferencing operator.
+                                          */
+        const Accessor* operator-> () const;
+
+                                         /**
+                                          * Comparison. True, if
+                                          * both iterators point to
+                                          * the same matrix
+                                          * position.
+                                          */
+        bool operator == (const const_iterator&) const;
+                                         /**
+                                          * Inverse of <tt>==</tt>.
+                                          */
+        bool operator != (const const_iterator&) 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 const_iterator&) const;
+
+                                         /**
+                                          * Exception
+                                          */
+        DeclException2 (ExcInvalidIndexWithinRow,
+                        int, int,
+                        << "Attempt to access element " << arg2
+                        << " of row " << arg1
+                        << " which doesn't have that many elements.");
+        
+      private:
+                                         /**
+                                          * Store an object of the
+                                          * accessor class.
+                                          */
+        Accessor accessor;
+    };
+    
+  }
   
   
 /**
@@ -61,6 +267,8 @@ namespace PETScWrappers
   class MatrixBase 
   {
     public:
+      typedef MatrixIterators::const_iterator const_iterator;
+      
                                        /**
                                         * Default constructor.
                                         */
@@ -337,6 +545,29 @@ namespace PETScWrappers
                             const VectorBase &x,
                             const VectorBase &b) const;
 
+                                       /**
+                                        * STL-like iterator with the
+                                        * first entry.
+                                        */
+      const_iterator begin () const;
+
+                                       /**
+                                        * Final iterator.
+                                        */
+      const_iterator end () const;
+
+                                       /**
+                                        * STL-like iterator with the
+                                        * first entry of row <tt>r</tt>.
+                                        */
+      const_iterator begin (const unsigned int r) const;
+
+                                       /**
+                                        * Final iterator of row
+                                        * <tt>r</tt>.
+                                        */
+      const_iterator end (const unsigned int r) const;
+      
                                        /**
                                         * Conversion operator to gain access
                                         * to the underlying PETSc type. If you
@@ -408,6 +639,157 @@ namespace PETScWrappers
 
 // -------------------------- inline and template functions ----------------------
 
+
+  namespace MatrixIterators
+  {
+
+    inline
+    const_iterator::Accessor::
+    Accessor (const MatrixBase    *matrix,
+              const unsigned int   row,
+              const unsigned short index)
+                    :
+                    matrix(matrix),
+                    a_row(row),
+                    a_index(index)
+    {
+      visit_present_row ();
+      Assert (index < colnum_cache->size(),
+              ExcInvalidIndexWithinRow (index, row));
+    }
+
+
+    inline
+    unsigned int
+    const_iterator::Accessor::row() const
+    {
+      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
+      return a_row;
+    }
+
+
+    inline
+    unsigned int
+    const_iterator::Accessor::column() const
+    {
+      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
+      return (*colnum_cache)[a_index];
+    }
+
+
+    inline
+    unsigned short
+    const_iterator::Accessor::index() const
+    {
+      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
+      return a_index;
+    }
+
+
+    inline
+    PetscScalar
+    const_iterator::Accessor::value() const
+    {
+      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
+      return (*value_cache)[a_index];
+    }
+
+
+    inline
+    const_iterator::
+    const_iterator(const MatrixBase    *matrix,
+                   const unsigned int   row,
+                   const unsigned short index)
+                    :
+                    accessor(matrix, row, index)
+    {}
+
+
+
+    inline
+    const_iterator &
+    const_iterator::operator++ ()
+    {
+      Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
+
+      ++accessor.a_index;
+      if (accessor.a_index >= accessor.colnum_cache->size())
+        {
+          accessor.a_index = 0;
+          accessor.a_row++;
+          accessor.visit_present_row();
+        }
+      return *this;
+    }
+
+
+    inline
+    const_iterator
+    const_iterator::operator++ (int)
+    {
+      const const_iterator old_state = *this;
+    
+      Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
+
+      ++accessor.a_index;
+      if (accessor.a_index >= accessor.colnum_cache->size())
+        {
+          accessor.a_index = 0;
+          accessor.a_row++;
+          accessor.visit_present_row();
+        }
+      return old_state;
+    }
+
+
+    inline
+    const const_iterator::Accessor &
+    const_iterator::operator* () const
+    {
+      return accessor;
+    }
+
+
+    inline
+    const const_iterator::Accessor *
+    const_iterator::operator-> () const
+    {
+      return &accessor;
+    }
+
+
+    inline
+    bool
+    const_iterator::
+    operator == (const const_iterator& other) const
+    {
+      return (accessor.row() == other.accessor.row() &&
+              accessor.index() == other.accessor.index());
+    }
+
+
+    inline
+    bool
+    const_iterator::
+    operator != (const const_iterator& other) const
+    {
+      return ! (*this == other);
+    }
+
+
+    inline
+    bool
+    const_iterator::
+    operator < (const const_iterator& other) const
+    {
+      return (accessor.row() < other.accessor.row() ||
+              (accessor.row() == other.accessor.row() &&
+               accessor.index() < other.accessor.index()));
+    }
+    
+  }
+  
+  
   inline
   PetscScalar
   MatrixBase::operator() (const unsigned int i,
@@ -415,8 +797,42 @@ namespace PETScWrappers
   {
     return el(i,j);
   }
+
   
-  
+
+  inline
+  MatrixBase::const_iterator
+  MatrixBase::begin() const
+  {
+    return const_iterator(this, 0, 0);
+  }
+
+
+  inline
+  MatrixBase::const_iterator
+  MatrixBase::end() const
+  {
+    return const_iterator(this, m(), 0);
+  }
+
+
+  inline
+  MatrixBase::const_iterator
+  MatrixBase::begin(const unsigned int r) const
+  {
+    Assert (r < m(), ExcIndexRange(r, 0, m()));
+    return const_iterator(this, r, 0);
+  }
+
+
+  inline
+  MatrixBase::const_iterator
+  MatrixBase::end(const unsigned int r) const
+  {
+    Assert (r < m(), ExcIndexRange(r, 0, m()));
+    return const_iterator(this, r+1, 0);
+  }
+      
 }
 
 #endif // DEAL_II_USE_PETSC
index a7abf350612f5016fe0d2b19b57a9fb63f667e02..0c4a797eebb8de201b5bc28a5a9ac32a071bf992 100644 (file)
@@ -24,7 +24,6 @@
 #include <vector>
 
 
-
 namespace PETScWrappers
 {
 /**
@@ -99,7 +98,6 @@ namespace PETScWrappers
                     const unsigned int               n,
                     const std::vector<unsigned int> &row_lengths,
                     const bool                       is_symmetric = false);
-      
   };
 }
 
index 819c7a7e6a32145165548c4456224d9c95f11e6b..bff066f87b6734aa7f52b83827b7ef9e77556d0c 100644 (file)
 
 namespace PETScWrappers
 {
+  namespace MatrixIterators
+  {
+    void
+    MatrixBase::const_iterator::Accessor::
+    visit_present_row ()
+    {
+      const int a_row_int = this->a_row;  // need signed int for PETSc?
+
+                                       // get a representation of the present
+                                       // row
+      int          ncols;
+      int         *colnums;
+      PetscScalar *values;
+      int ierr;
+      ierr = MatGetRow(*matrix, a_row_int, &ncols, &colnums, &values);
+      AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr));
+
+                                       // copy it into our caches
+      colnum_cache =
+        boost::shared_ptr<std::vector<unsigned int> >
+        ( new std::vector<unsigned int> (colnums, colnums+ncols));
+      value_cache =
+        boost::shared_ptr<std::vector<PetscScalar> >
+        ( new std::vector<PetscScalar> (values, values+ncols));
+
+                                       // and finally restore the matrix
+      ierr = MatRestoreRow(*matrix, a_row_int, &ncols, &colnums, &values);
+      AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr));
+    }
+  }
+
+
+  
   MatrixBase::MatrixBase ()
                   :
                   last_action (LastAction::none)

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.