]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
SparseMatrixEZ const_iterator
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 6 Sep 2002 13:36:49 +0000 (13:36 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 6 Sep 2002 13:36:49 +0000 (13:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@6370 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_matrix_ez.h
tests/lac/sparse_matrices.cc
tests/lac/sparse_matrices.checked

index b74ab7af60aaae482b395c7c0530ea56bacd1195..a04a120a9d4d1b3dbe861cba5b134b1188d7d4d4 100644 (file)
@@ -31,9 +31,50 @@ template<typename number> class FullMatrix;
  * builds the pattern on the fly. Filling the matrix may consume more
  * time as with @p{SparseMatrix}, since large memory movements may be
  * involved. To help optimizing things, an expected row-length may be
- * provided to the constructor, as well as a mininmum size increment
+ * provided to the constructor, as well as an increment size
  * for rows.
  *
+ * The storage structure: like with the usual sparse matrix, it is
+ * attempted to store only non-zero elements. these are stored in a
+ * single data array @p{data}. They are ordered by row and inside each
+ * row by column number. Each row is described by its starting point
+ * in the data array and its length. These are stored in the
+ * @p{row_info} array, together with additional useful information.
+ *
+ * Due to the structure, gaps may occur between rows. Whenever a new
+ * entry must be created, an attempt is made to use the gap in its
+ * row. If the gap is full, the row must be extended and all
+ * subsequent rows must be shifted backwards. This is a very expensive
+ * operation and should be avoided as much as possible.
+ *
+ * This is, where the optimization parameters, provided to the
+ * constructor or to the function @p{reinit} come
+ * in. @p{default_row_length} is the amount of entries that will be
+ * allocated for each row on initialization (the actual length of the
+ * rows is still zero). This means, that @p{default_row_length}
+ * entries can be added to this row without shifting other rows. If
+ * less entries are added, the additional memory will be wasted.
+ *
+ * If the space for a row is not sufficient, then it is enlarged by
+ * @p{default_increment} entries. This way, the subsequent rows are
+ * not shifted by single entries very often.
+ *
+ * Finally, the @p{default_reserve} allocates extra space at the end
+ * of the data array. This space is used whenever a row must be
+ * enlarged. Since @p{std::vector} doubles the capacity everytime it
+ * must increase it, this value should allow for all the growth needed.
+ *
+ * Suggested settings: @p{default_row_length} should be the length of
+ * a typical row, for instance the size of the stencil in regular
+ * parts of the grid. Then, @p{default_increment} may be the expected
+ * amount of entries added to the row by having one hanging node. This
+ * way, a good compromise between memory consumption and speed should
+ * be achieved. @p{default_reserve} should then be an estimate for the
+ * number of hanging nodes times @p{default_increment}.
+ *
+ * If the rows are expected to be filled more or less from first to
+ * last, using a @p{default_row_length} may not be such a bad idea.
+ *
  * The name of this matrix is in reverence to a publication of the
  * Internal Revenue Service of the United States of America. I hope
  * some other aliens will appreciate it. By the way, the suffix makes
@@ -124,9 +165,124 @@ class SparseMatrixEZ : public Subscriptor
        static const unsigned short
        invalid_diagonal = static_cast<unsigned short>(-1);
     };
-    
-    
+
   public:
+                                    /**
+                                     * Accessor class for iterators
+                                     */
+    class Accessor
+    {
+      public:
+                                        /**
+                                         * Constructor. Since we use
+                                         * accessors only for read
+                                         * access, a const matrix
+                                         * pointer is sufficient.
+                                         */
+       Accessor (const SparseMatrixEZ<number>*,
+                 unsigned int row,
+                 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.
+                                         */
+       number value() const;
+       
+       protected:
+                                        /**
+                                         * The matrix accessed.
+                                         */
+       const SparseMatrixEZ<number>* matrix;
+
+                                        /**
+                                         * Current row number.
+                                         */
+       unsigned int a_row;
+
+                                        /**
+                                         * Current index in row.
+                                         */
+       unsigned short a_index;
+      };
+
+                                    /**
+                                     * STL conforming iterator.
+                                     */
+    class const_iterator : private Accessor
+      {
+       public:
+                                          /**
+                                           * Constructor.
+                                           */ 
+       const_iterator(const SparseMatrixEZ<number>*,
+                      unsigned int row,
+                      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 @p{==}.
+                                           */
+       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;
+      };
+    
                                     /**
                                      * Type of matrix entries. In analogy to
                                      * the STL container classes.
@@ -237,31 +393,6 @@ class SparseMatrixEZ : public Subscriptor
                                      */
     unsigned int n () const;
 
-                                    /**
-                                     * Return the number of nonzero
-                                     * elements of this
-                                     * matrix. Actually, it returns
-                                     * the number of entries in the
-                                     * sparsity pattern; if any of
-                                     * the entries should happen to
-                                     * be zero, it is counted anyway.
-                                     */
-//    unsigned int n_nonzero_elements () const;
-
-                                    /**
-                                     * Return the number of actually
-                                     * nonzero elements of this
-                                     * matrix.
-                                     *
-                                     * Note, that this function does
-                                     * (in contrary to the
-                                     * @p{n_nonzero_elements}) NOT
-                                     * count all entries of the
-                                     * sparsity pattern but only the
-                                     * ones that are nonzero.
-                                     */
-//    unsigned int n_actually_nonzero_elements () const;
-    
                                     /**
                                      * Set the element @p{(i,j)} to
                                      * @p{value}. Allocates the entry
@@ -427,6 +558,14 @@ class SparseMatrixEZ : public Subscriptor
     number operator () (const unsigned int i,
                        const unsigned int j) const;
 
+                                    /**
+                                     * Return the value of the entry
+                                     * (i,j). Returns zero for all
+                                     * non-existing entries.
+                                     */
+    number el (const unsigned int i,
+              const unsigned int j) const;
+
                                     /**
                                      * Return the main diagonal element in
                                      * the @p{i}th row. This function throws an
@@ -441,14 +580,14 @@ class SparseMatrixEZ : public Subscriptor
                                      * involve searching for the
                                      * right column number.
                                      */
-    number diag_element (const unsigned int i) const;
+//    number diag_element (const unsigned int i) const;
 
                                     /**
                                      * Same as above, but return a
                                      * writeable reference. You're
                                      * sure you know what you do?
                                      */
-    number & diag_element (const unsigned int i);
+//    number & diag_element (const unsigned int i);
     
                                     /**
                                      * Matrix-vector multiplication:
@@ -592,6 +731,42 @@ class SparseMatrixEZ : public Subscriptor
     void precondition_TSOR (Vector<somenumber>       &dst,
                            const Vector<somenumber> &src,
                            const number              om = 1.) 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 @p{r}.
+                                     */
+    const_iterator begin (unsigned int r) const;
+
+                                    /**
+                                     * Final iterator of row @p{r}.
+                                     */
+    const_iterator end (unsigned int r) const;
+    
+                                    /**
+                                     * Return the number of nonzero
+                                     * elements of this
+                                     * matrix.
+                                     */
+//    unsigned int n_nonzero_elements () const;
+
+                                    /**
+                                     * Return the number of actually
+                                     * nonzero elements of this
+                                     * matrix.
+                                     */
+//    unsigned int n_actually_nonzero_elements () const;
     
                                     /**
                                      * Print the matrix to the given
@@ -667,8 +842,15 @@ class SparseMatrixEZ : public Subscriptor
                                      */
     DeclException0(ExcNoDiagonal);
     
-  private:
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidIndex,
+                   int, int,
+                   << "The entry with index <" << arg1 << ',' << arg2
+                   << "> does not exist.");
 
+  private:
                                     /**
                                      * Find an entry. Return a
                                      * zero-pointer if the entry does
@@ -772,10 +954,9 @@ class SparseMatrixEZ : public Subscriptor
                                      * Increment when a row grows.
                                      */
     unsigned int increment;
-    
-                                    // make all other sparse matrices
-                                    // friends
-    template <typename somenumber> friend class SparseMatrix;
+
+    friend class Accessor;
+    friend class const_iterator;
 };
 
 
@@ -801,31 +982,142 @@ SparseMatrixEZ<number>::Entry::Entry()
 {}
 
 
-// template <typename number>
-// inline
-// bool
-// SparseMatrixEZ<number>::Entry::operator==(const Entry& e) const
-// {
-//   return column == e.column;
-// }
+template <typename number>
+inline
+SparseMatrixEZ<number>::RowInfo::RowInfo(unsigned int start)
+               :
+               start(start), length(0), diagonal(invalid_diagonal)
+{}
+
+
+//----------------------------------------------------------------------//
+template <typename number>
+inline
+SparseMatrixEZ<number>::Accessor::Accessor (
+  const SparseMatrixEZ<number>* matrix,
+  unsigned int r,
+  unsigned short i)
+               :
+               matrix(matrix),
+               a_row(r),
+               a_index(i)
+{}
+
+
+template <typename number>
+inline
+unsigned int
+SparseMatrixEZ<number>::Accessor::row() const
+{
+  return a_row;
+}
+
+
+template <typename number>
+inline
+unsigned int
+SparseMatrixEZ<number>::Accessor::column() const
+{
+  return matrix->data[matrix->row_info[a_row].start+a_index].column;
+}
+
+
+template <typename number>
+inline
+unsigned short
+SparseMatrixEZ<number>::Accessor::index() const
+{
+  return a_index;
+}
 
 
-// template <typename number>
-// inline
-// bool
-// SparseMatrixEZ<number>::Entry::operator<(const Entry& e) const
-// {
-//   return column < e.column;
-// }
 
 template <typename number>
 inline
-SparseMatrixEZ<number>::RowInfo::RowInfo(unsigned int start)
+number
+SparseMatrixEZ<number>::Accessor::value() const
+{
+  return matrix->data[matrix->row_info[a_row].start+a_index].value;
+}
+
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator::const_iterator(
+  const SparseMatrixEZ<number>* matrix,
+  unsigned int r,
+  unsigned short i)
                :
-               start(start), length(0), diagonal(invalid_diagonal)
+               Accessor(matrix, r, i)
 {}
 
 
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator&
+SparseMatrixEZ<number>::const_iterator::operator++ ()
+{
+  Assert (a_row < matrix->m(), ExcIteratorPastEnd());
+  
+  ++a_index;
+  if (a_index >= matrix->row_info[a_row].length)
+    {
+      a_index = 0;
+      a_row++;
+    }
+  return *this;
+}
+
+
+template <typename number>
+inline
+const SparseMatrixEZ<number>::Accessor&
+SparseMatrixEZ<number>::const_iterator::operator* () const
+{
+  return *this;
+}
+
+
+template <typename number>
+inline
+const SparseMatrixEZ<number>::Accessor*
+SparseMatrixEZ<number>::const_iterator::operator-> () const
+{
+  return this;
+}
+
+
+template <typename number>
+inline
+bool
+SparseMatrixEZ<number>::const_iterator::operator == (
+  const const_iterator& other) const
+{
+  return (row() == other->row() && index() == other->index());
+}
+
+
+template <typename number>
+inline
+bool
+SparseMatrixEZ<number>::const_iterator::operator != (
+  const const_iterator& other) const
+{
+  return ! (*this == other);
+}
+
+
+template <typename number>
+inline
+bool
+SparseMatrixEZ<number>::const_iterator::operator < (
+  const const_iterator& other) const
+{
+  return (row() < other->row() ||
+         (row() == other->row() && index() < other->index()));
+}
+
+
 //----------------------------------------------------------------------//
 template <typename number>
 inline
@@ -983,5 +1275,66 @@ void SparseMatrixEZ<number>::add (const unsigned int i,
 };
 
 
+
+template <typename number>
+inline
+number SparseMatrixEZ<number>::el (const unsigned int i,
+                                  const unsigned int j) const
+{
+  Entry* entry = locate(i,j);
+  if (entry)
+    return entry->value;
+  return 0.;
+}
+
+
+
+template <typename number>
+inline
+number SparseMatrixEZ<number>::operator() (const unsigned int i,
+                                          const unsigned int j) const
+{
+  Entry* entry = locate(i,j);
+  if (entry)
+    return entry->value;
+  Assert(false, ExcInvalidEntry(i,j));
+  return 0.;
+}
+
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator
+SparseMatrixEZ<number>::begin () const
+{
+  return const_iterator(this, 0, 0);
+}
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator
+SparseMatrixEZ<number>::end () const
+{
+  return const_iterator(this, m(), 0);
+}
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator
+SparseMatrixEZ<number>::begin (unsigned int r) const
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return const_iterator(this, r, 0);
+}
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator
+SparseMatrixEZ<number>::end (unsigned int r) const
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return const_iterator(this, r+1, 0);
+}
+
 #endif
 /*----------------------------   sparse_matrix.h     ---------------------------*/
index ba74a4d491f451f2b65748b7669968459a6337f6..9cfa6f85c86148046954e300fc4212709d36c2a7 100644 (file)
 
 #include <fstream>
 
+#define PREC_CHECK(solver, method, precond) try \
+ { solver.method (A, u, f, precond);  } catch (...) {} \
+ residuals.push_back(control.last_value())
+
 template<class MATRIX>
 void
-check_vmult_quadratic(const MATRIX& A,
+check_vmult_quadratic(std::vector<double>& residuals,
+                     const MATRIX& A,
                      const char* prefix)
 {
   deallog.push(prefix);
@@ -40,7 +45,7 @@ check_vmult_quadratic(const MATRIX& A,
   SolverControl control(10, 1.e-13, false);
   SolverRichardson<> rich(control, mem, .01);
   SolverRichardson<> prich(control, mem, 1.);
-  PreconditionIdentity prec;
+  PreconditionIdentity identity;
   PreconditionJacobi<MATRIX> jacobi;
   jacobi.initialize(A, .5);
   PreconditionSOR<MATRIX> sor;
@@ -51,20 +56,41 @@ check_vmult_quadratic(const MATRIX& A,
   u = 0.;
   f = 1.;
 
-  try { rich.solve(A, u, f, prec); }    catch (...) {}
-  try { prich.solve(A, u, f, jacobi); } catch (...) {}
-  try { prich.solve(A, u, f, ssor); }   catch (...) {}
-  try { prich.solve(A, u, f, sor); }    catch (...) {}
-
+  PREC_CHECK(rich, solve, identity);
+  PREC_CHECK(prich, solve, jacobi);
+  PREC_CHECK(prich, solve, ssor);
+  PREC_CHECK(prich, solve, sor);
+  
   u = 0.;
   deallog << "Transpose" << std::endl;
-  try { rich.Tsolve(A, u, f, prec); }    catch (...) {}
-  try { prich.Tsolve(A, u, f, jacobi); } catch (...) {}
-  try { prich.Tsolve(A, u, f, ssor); }   catch (...) {}
-  try { prich.Tsolve(A, u, f, sor); }    catch (...) {}
+  PREC_CHECK(rich, Tsolve, identity);
+  PREC_CHECK(prich, Tsolve, jacobi);
+  PREC_CHECK(prich, Tsolve, ssor);
+  PREC_CHECK(prich, Tsolve, sor);
   deallog.pop();
 }
 
+
+template <class MATRIX>
+void
+check_iterator (const MATRIX& A)
+{
+  for (typename MATRIX::const_iterator i = A.begin(); i!= A.end(); ++i)
+    deallog << '\t' << i->row()
+           << '\t' << i->column()
+           << '\t' << i->index()
+           << '\t' << i->value()
+           << std::endl;
+  deallog << "Repeat row 2" << std::endl;
+  for (typename MATRIX::const_iterator i = A.begin(2); i!= A.end(2); ++i)
+    deallog << '\t' << i->row()
+           << '\t' << i->column()
+           << '\t' << i->index()
+           << '\t' << i->value()
+           << std::endl;
+}
+
+
 int main()
 {
   std::ofstream logfile("sparse_matrices.output");
@@ -89,6 +115,9 @@ int main()
   FDMatrix testproblem (size, size);
   unsigned int dim = (size-1)*(size-1);
 
+  std::vector<double> A_res;
+  std::vector<double> E_res;
+  
   deallog << "Structure" << std::endl;
   SparsityPattern structure(dim, dim, 5);
   testproblem.five_point_structure(structure);
@@ -96,13 +125,16 @@ int main()
   SparseMatrix<double>  A(structure);
   deallog << "Assemble" << std::endl;
   testproblem.five_point(A, true);
-  check_vmult_quadratic(A, "5-SparseMatrix<double>");
+  check_vmult_quadratic(A_res, A, "5-SparseMatrix<double>");
 
   SparseMatrixEZ<double> E(dim,dim,row_length,2);
   deallog << "Assemble" << std::endl;
   testproblem.five_point(E, true);
-  check_vmult_quadratic(E, "5-SparseMatrixEZ<double>");
-
+  check_vmult_quadratic(E_res, E, "5-SparseMatrixEZ<double>");
+#ifdef DEBUG
+  check_iterator(E);
+#endif
+  
   A.clear();
   deallog << "Structure" << std::endl;
   structure.reinit(dim, dim, 9);
@@ -111,12 +143,16 @@ int main()
   A.reinit(structure);
   deallog << "Assemble" << std::endl;
   testproblem.nine_point(A);
-  check_vmult_quadratic(A, "9-SparseMatrix<double>");
+  check_vmult_quadratic(A_res, A, "9-SparseMatrix<double>");
 
   E.clear();
   E.reinit(dim,dim,row_length,2);
   deallog << "Assemble" << std::endl;
   testproblem.nine_point(E);
-  check_vmult_quadratic(E, "9-SparseMatrixEZ<double>");
-  
+  check_vmult_quadratic(E_res, E, "9-SparseMatrixEZ<double>");
+
+  for (unsigned int i=0;i<A_res.size();++i)
+    if (std::fabs(A_res[i] - E_res[i]) > 1.e-14)
+      deallog << "SparseMatrix and SparseMatrixEZ differ!!!"
+             << std::endl;
 }
index 15304de6ce44b5b777f61c48a85a0f0e7d0e4dea..7b9387f15b615ba34aca3729a44a03eb94ecb11b 100644 (file)
@@ -40,6 +40,75 @@ DEAL:5-SparseMatrixEZ<double>:RichardsonT::Starting
 DEAL:5-SparseMatrixEZ<double>:RichardsonT::Failure step 10 
 DEAL::GrowingVectorMemory:Overall allocated vectors: 16
 DEAL::GrowingVectorMemory:Maximum allocated vectors: 2
+DEAL:: 0       0       0       5.00
+DEAL:: 0       1       1       -2.00
+DEAL:: 0       4       2       -1.00
+DEAL:: 1       0       0       -1.00
+DEAL:: 1       1       1       5.00
+DEAL:: 1       2       2       -2.00
+DEAL:: 1       5       3       -1.00
+DEAL:: 2       1       0       -1.00
+DEAL:: 2       2       1       5.00
+DEAL:: 2       3       2       -2.00
+DEAL:: 2       6       3       -1.00
+DEAL:: 3       2       0       -1.00
+DEAL:: 3       3       1       5.00
+DEAL:: 3       7       2       -1.00
+DEAL:: 4       0       0       -1.00
+DEAL:: 4       4       1       5.00
+DEAL:: 4       5       2       -2.00
+DEAL:: 4       8       3       -1.00
+DEAL:: 5       1       0       -1.00
+DEAL:: 5       4       1       -1.00
+DEAL:: 5       5       2       5.00
+DEAL:: 5       6       3       -2.00
+DEAL:: 5       9       4       -1.00
+DEAL:: 6       2       0       -1.00
+DEAL:: 6       5       1       -1.00
+DEAL:: 6       6       2       5.00
+DEAL:: 6       7       3       -2.00
+DEAL:: 6       10      4       -1.00
+DEAL:: 7       3       0       -1.00
+DEAL:: 7       6       1       -1.00
+DEAL:: 7       7       2       5.00
+DEAL:: 7       11      3       -1.00
+DEAL:: 8       4       0       -1.00
+DEAL:: 8       8       1       5.00
+DEAL:: 8       9       2       -2.00
+DEAL:: 8       12      3       -1.00
+DEAL:: 9       5       0       -1.00
+DEAL:: 9       8       1       -1.00
+DEAL:: 9       9       2       5.00
+DEAL:: 9       10      3       -2.00
+DEAL:: 9       13      4       -1.00
+DEAL:: 10      6       0       -1.00
+DEAL:: 10      9       1       -1.00
+DEAL:: 10      10      2       5.00
+DEAL:: 10      11      3       -2.00
+DEAL:: 10      14      4       -1.00
+DEAL:: 11      7       0       -1.00
+DEAL:: 11      10      1       -1.00
+DEAL:: 11      11      2       5.00
+DEAL:: 11      15      3       -1.00
+DEAL:: 12      8       0       -1.00
+DEAL:: 12      12      1       5.00
+DEAL:: 12      13      2       -2.00
+DEAL:: 13      9       0       -1.00
+DEAL:: 13      12      1       -1.00
+DEAL:: 13      13      2       5.00
+DEAL:: 13      14      3       -2.00
+DEAL:: 14      10      0       -1.00
+DEAL:: 14      13      1       -1.00
+DEAL:: 14      14      2       5.00
+DEAL:: 14      15      3       -2.00
+DEAL:: 15      11      0       -1.00
+DEAL:: 15      14      1       -1.00
+DEAL:: 15      15      2       5.00
+DEAL::Repeat row 2
+DEAL:: 2       1       0       -1.00
+DEAL:: 2       2       1       5.00
+DEAL:: 2       3       2       -2.00
+DEAL:: 2       6       3       -1.00
 DEAL::Structure
 DEAL::Assemble
 DEAL:9-SparseMatrix<double>:Richardson::Starting 

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.