]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add row_length(). Make sure that matrix iterators also work if there are empty rows...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 3 May 2005 03:02:23 +0000 (03:02 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 3 May 2005 03:02:23 +0000 (03:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@10603 0785d39b-7218-0410-832d-ea1e28bc413d

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

index df58d602c7bd3e0ece11e9a646da3d13a0a3f751..ace5bf7effe66f897c8c8885ef85f9c8d768d864 100644 (file)
@@ -491,6 +491,12 @@ namespace PETScWrappers
                                         */
       unsigned int n_nonzero_elements () const;
       
+
+                                       /**
+                                        * Number of entries in a specific row.
+                                        */
+      unsigned int row_length (const unsigned int row) const;
+      
                                        /**
                                         * Return the l1-norm of the matrix, that is
                                         * $|M|_1=max_{all columns j}\sum_{all 
@@ -669,12 +675,27 @@ namespace PETScWrappers
                                        /**
                                         * STL-like iterator with the
                                         * first entry of row @p r.
+                                        *
+                                        * Note that if the given row is empty,
+                                        * i.e. does not contain any nonzero
+                                        * entries, then the iterator returned by
+                                        * this function equals
+                                        * <tt>end(r)</tt>. Note also that the
+                                        * iterator may not be dereferencable in
+                                        * that case.
                                         */
       const_iterator begin (const unsigned int r) const;
 
                                        /**
-                                        * Final iterator of row
-                                        * @p r.
+                                        * Final iterator of row <tt>r</tt>. It
+                                        * points to the first element past the
+                                        * end of line @p r, or past the end of
+                                        * the entire sparsity pattern.
+                                        *
+                                        * Note that the end iterator is not
+                                        * necessarily dereferencable. This is in
+                                        * particular the case if it is the end
+                                        * iterator for the last row of a matrix.
                                         */
       const_iterator end (const unsigned int r) const;
       
@@ -764,9 +785,6 @@ namespace PETScWrappers
                     a_index(index)
     {
       visit_present_row ();
-      Assert ((row == matrix->m()) ||
-              (index < colnum_cache->size()),
-              ExcInvalidIndexWithinRow (index, row));
     }
 
 
@@ -824,10 +842,25 @@ namespace PETScWrappers
       Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
 
       ++accessor.a_index;
+
+                                       // if at end of line: do one step, then
+                                       // cycle until we find a row with a
+                                       // nonzero number of entries
       if (accessor.a_index >= accessor.colnum_cache->size())
         {
           accessor.a_index = 0;
-          accessor.a_row++;
+          ++accessor.a_row;
+      
+          while (accessor.a_index >= accessor.matrix->row_length(accessor.a_row))
+            {
+              ++accessor.a_row;
+
+                                           // if we happened to find the end
+                                           // of the matrix, then stop here
+              if (accessor.a_row == accessor.matrix->m())
+                break;
+            }
+
           accessor.visit_present_row();
         }
       return *this;
@@ -839,16 +872,7 @@ namespace PETScWrappers
     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();
-        }
+      ++(*this);
       return old_state;
     }
 
@@ -932,7 +956,10 @@ namespace PETScWrappers
   MatrixBase::begin(const unsigned int r) const
   {
     Assert (r < m(), ExcIndexRange(r, 0, m()));
-    return const_iterator(this, r, 0);
+    if (row_length(r) > 0)
+      return const_iterator(this, r, 0);
+    else
+      return end (r);
   }
 
 
@@ -941,7 +968,17 @@ namespace PETScWrappers
   MatrixBase::end(const unsigned int r) const
   {
     Assert (r < m(), ExcIndexRange(r, 0, m()));
-    return const_iterator(this, r+1, 0);
+
+                                     // place the iterator on the first entry
+                                     // past this line, or at the end of the
+                                     // matrix
+    for (unsigned int i=r+1; i<m(); ++i)
+      if (row_length(i) > 0)
+        return const_iterator(this, i, 0);
+    
+                                     // if there is no such line, then take the
+                                     // end iterator of the matrix
+    return end();
   }
 
 
index 7e5f0c9c03fb59150b502175f28e270ee53bcfae..e09f59fd2262da9c3655655f367b966f7335e17d 100644 (file)
@@ -63,18 +63,15 @@ namespace PETScWrappers
       AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr));
 
                                        // copy it into our caches if the line
-                                       // isn't empty
-      if (ncols != 0)
-        {
-          colnum_cache.reset (new std::vector<unsigned int> (colnums,
-                                                             colnums+ncols));
-          value_cache.reset (new std::vector<PetscScalar> (values, values+ncols));
-        }
-      else
-        {
-          colnum_cache.reset ();
-          value_cache.reset ();
-        }
+                                       // isn't empty. if it is, then we've
+                                       // done something wrong, since we
+                                       // shouldn't have initialized an
+                                       // iterator for an empty line (what
+                                       // would it point to?)
+      Assert (ncols != 0, ExcInternalError());
+      colnum_cache.reset (new std::vector<unsigned int> (colnums,
+                                                         colnums+ncols));
+      value_cache.reset (new std::vector<PetscScalar> (values, values+ncols));
       
                                        // and finally restore the matrix
       ierr = MatRestoreRow(*matrix, this->a_row, &ncols, &colnums, &values);
@@ -309,6 +306,48 @@ namespace PETScWrappers
 
 
   
+  unsigned int
+  MatrixBase::
+  row_length (const unsigned int row) const
+  {
+//TODO: this function will probably only work if compress() was called on the
+//matrix previously. however, we can't do this here, since it would impose
+//global communication and one would have to make sure that this function is
+//called the same number of times from all processors, something that is
+//unreasonable. there should simply be a way in PETSc to query the number of
+//entries in a row bypassing the call to compress(), but I can't find one
+    Assert (row < m(), ExcInternalError());
+
+                                     // get a representation of the present
+                                     // row
+    int ncols;
+
+#if (PETSC_VERSION_MAJOR <= 2) && \
+    ((PETSC_VERSION_MINOR < 2) ||  \
+     ((PETSC_VERSION_MINOR == 2) && (PETSC_VERSION_SUBMINOR == 0)))
+    int         *colnums;
+    PetscScalar *values;
+#else
+    const int         *colnums;
+    const PetscScalar *values;
+#endif
+
+//TODO: this is probably horribly inefficient; we should lobby for a way to
+//query this information from PETSc
+    int ierr;
+    ierr = MatGetRow(*this, row, &ncols, &colnums, &values);
+    Assert (ierr == 0, MatrixBase::ExcPETScError(ierr));
+
+                                     // then restore the matrix and return the
+                                     // number of columns in this row as
+                                     // queried previously
+    ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values);
+    AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr));
+
+    return ncols;
+  }
+
+
   PetscReal
   MatrixBase::l1_norm () const
   {

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.