]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
preconditioners for SparseMatrixEZ
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 5 Sep 2002 17:43:24 +0000 (17:43 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 5 Sep 2002 17:43:24 +0000 (17:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@6363 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_matrix_ez.h
deal.II/lac/include/lac/sparse_matrix_ez.templates.h
tests/lac/sparse_matrices.cc
tests/lac/testmatrix.cc
tests/lac/testmatrix.h

index a3eb58fa1c79d991ca3d2871dd65cbdfcf4b8b5f..b74ab7af60aaae482b395c7c0530ea56bacd1195 100644 (file)
@@ -55,7 +55,7 @@ class SparseMatrixEZ : public Subscriptor
                                         /**
                                          * Standard constructor. Sets
                                          * @p{column} to
-                                         * @p{invalid_entry}.
+                                         * @p{invalid}.
                                          */
        Entry();
 
@@ -77,18 +77,55 @@ class SparseMatrixEZ : public Subscriptor
                                         /**
                                          * Comparison operator for finding.
                                          */
-       bool operator==(const Entry&) const;
+//     bool operator==(const Entry&) const;
 
                                         /**
                                          * Less than operator for sorting.
                                          */
-       bool operator < (const Entry&) const;
+//     bool operator < (const Entry&) const;
                                         /**
                                          * Non-existent column number.
                                          */
-       static const unsigned int invalid_entry = static_cast<unsigned int>(-1);
+       static const unsigned int invalid = static_cast<unsigned int>(-1);
     };
 
+                                    /**
+                                     * Structure for storing
+                                     * information on a matrix
+                                     * row. One object for each row
+                                     * will be stored in the matrix.
+                                     */
+    struct RowInfo
+    {
+                                        /**
+                                         * Constructor.
+                                         */
+       RowInfo (unsigned int start = Entry::invalid);
+       
+                                        /**
+                                         * Index of first entry of
+                                         * the row in the data field.
+                                         */
+       unsigned int start;
+                                        /**
+                                         * Number of entries in this
+                                         * row.
+                                         */
+       unsigned short length;
+                                        /**
+                                         * Position of the diagonal
+                                         * element relative tor the
+                                         * start index.
+                                         */
+       unsigned short diagonal;
+                                        /**
+                                         * Value for non-existing diagonal.
+                                         */
+       static const unsigned short
+       invalid_diagonal = static_cast<unsigned short>(-1);
+    };
+    
+    
   public:
                                     /**
                                      * Type of matrix entries. In analogy to
@@ -120,12 +157,20 @@ class SparseMatrixEZ : public Subscriptor
                                     /**
                                      * Constructor. Generates a
                                      * matrix of the given size,
-                                     * ready to be filled. The optional parameters 
+                                     * ready to be filled. The
+                                     * optional parameters
+                                     * @p{default_row_length} and
+                                     * @p{default_increment} allow
+                                     * for preallocating
+                                     * memory. Providing these
+                                     * properly is essential for an
+                                     * efficient assembling of the
+                                     * matrix.
                                      */
     explicit SparseMatrixEZ (unsigned int n_rows,
                             unsigned int n_columns = n_rows,
-                            unsigned int default_row_length = Entry::invalid_entry,
-                            unsigned int default_increment = Entry::invalid_entry);
+                            unsigned int default_row_length = Entry::invalid,
+                            unsigned int default_increment = Entry::invalid);
     
                                     /**
                                      * Destructor. Free all memory, but do not
@@ -144,12 +189,20 @@ class SparseMatrixEZ : public Subscriptor
                                      * Reinitialize the sparse matrix
                                      * to the dimensions provided.
                                      * The matrix will have no
-                                     * entries at this point.
+                                     * entries at this point. The
+                                     * optional parameters
+                                     * @p{default_row_length} and
+                                     * @p{default_increment} allow
+                                     * for preallocating
+                                     * memory. Providing these
+                                     * properly is essential for an
+                                     * efficient assembling of the
+                                     * matrix.
                                      */
     virtual void reinit (unsigned int n_rows,
                         unsigned int n_columns = n_rows,
-                        unsigned int default_row_length = Entry::invalid_entry,
-                        unsigned int default_increment = Entry::invalid_entry);
+                        unsigned int default_row_length = Entry::invalid,
+                        unsigned int default_increment = Entry::invalid);
 
                                     /**
                                      * Release all memory and return
@@ -540,67 +593,6 @@ class SparseMatrixEZ : public Subscriptor
                            const Vector<somenumber> &src,
                            const number              om = 1.) const;
     
-                                    /**
-                                     * Perform SSOR preconditioning
-                                     * in-place.  Apply the
-                                     * preconditioner matrix without
-                                     * copying to a second vector.
-                                     * @p{omega} is the relaxation
-                                     * parameter.
-                                     */
-    template <typename somenumber>
-    void SSOR (Vector<somenumber> &v,
-              const number        omega = 1.) const;
-
-                                    /**
-                                     * Perform an SOR preconditioning in-place.
-                                     * The result is $v = (\omega D - L)^{-1} v$.
-                                     * @p{omega} is the damping parameter.
-                                     */
-    template <typename somenumber>
-    void SOR (Vector<somenumber> &v,
-             const number        om = 1.) const;
-
-                                    /**
-                                     * Perform a transpose SOR preconditioning in-place.
-                                     * The result is $v = (\omega D - L)^{-1} v$.
-                                     * @p{omega} is the damping parameter.
-                                     */
-    template <typename somenumber>
-    void TSOR (Vector<somenumber> &v,
-             const number        om = 1.) const;
-
-                                    /**
-                                     * Do one SOR step on @p{v}.
-                                     * Performs a direct SOR step
-                                     * with right hand side @p{b}.
-                                     */
-    template <typename somenumber>
-    void SOR_step (Vector<somenumber> &v,
-                  const Vector<somenumber> &b,
-                  const number        om = 1.) const;
-
-                                    /**
-                                     * Do one adjoint SOR step on
-                                     * @p{v}.  Performs a direct TSOR
-                                     * step with right hand side @p{b}.
-                                     */
-    template <typename somenumber>
-    void TSOR_step (Vector<somenumber> &v,
-                   const Vector<somenumber> &b,
-                   const number        om = 1.) const;
-
-                                    /**
-                                     * Do one adjoint SSOR step on
-                                     * @p{v}.  Performs a direct SSOR
-                                     * step with right hand side @p{b}
-                                     * by performing TSOR after SOR.
-                                     */
-    template <typename somenumber>
-    void SSOR_step (Vector<somenumber> &v,
-                   const Vector<somenumber> &b,
-                   const number        om = 1.) const;
-
                                     /**
                                      * Print the matrix to the given
                                      * stream, using the format
@@ -662,29 +654,19 @@ class SparseMatrixEZ : public Subscriptor
                                      * of this object.
                                      */
     unsigned int memory_consumption () const;
-    
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcInvalidIndex,
-                   int, int,
-                   << "The entry with index <" << arg1 << ',' << arg2
-                   << "> does not exist.");
 
                                     /**
-                                     * Exception
+                                     * Exception for applying
+                                     * inverse-type operators to
+                                     * rectangular matrices.
                                      */
-    DeclException0 (ExcMatrixNotSquare);
-
+    DeclException0(ExcNoSquareMatrix);
+    
                                     /**
-                                     * Exception
+                                     * Exception for missing diagonal entry.
                                      */
-    DeclException2 (ExcIteratorRange,
-                   int, int,
-                   << "The iterators denote a range of " << arg1
-                   << " elements, but the given number of rows was " << arg2);
+    DeclException0(ExcNoDiagonal);
     
-
   private:
 
                                     /**
@@ -694,6 +676,7 @@ class SparseMatrixEZ : public Subscriptor
                                      */
     const Entry* locate (unsigned int row,
                         unsigned int col) const;
+
                                     /**
                                      * Find an entry or generate it.
                                      */
@@ -776,10 +759,9 @@ class SparseMatrixEZ : public Subscriptor
     unsigned int n_columns;
 
                                     /**
-                                     * Start of indices rows. Points
-                                     * into the data field.
+                                     * Info structure for each row.
                                      */
-    std::vector<unsigned int> row_start;
+    std::vector<RowInfo> row_info;
     
                                     /**
                                      * Data storage.
@@ -814,28 +796,34 @@ template <typename number>
 inline
 SparseMatrixEZ<number>::Entry::Entry()
                :
-               column(invalid_entry),
+               column(invalid),
   value(0)
 {}
 
 
-template <typename number>
-inline
-bool
-SparseMatrixEZ<number>::Entry::operator==(const Entry& e) const
-{
-  return column == e.column;
-}
+// template <typename number>
+// inline
+// bool
+// SparseMatrixEZ<number>::Entry::operator==(const Entry& e) const
+// {
+//   return column == e.column;
+// }
 
 
+// template <typename number>
+// inline
+// bool
+// SparseMatrixEZ<number>::Entry::operator<(const Entry& e) const
+// {
+//   return column < e.column;
+// }
+
 template <typename number>
 inline
-bool
-SparseMatrixEZ<number>::Entry::operator<(const Entry& e) const
-{
-  return column < e.column;
-}
-
+SparseMatrixEZ<number>::RowInfo::RowInfo(unsigned int start)
+               :
+               start(start), length(0), diagonal(invalid_diagonal)
+{}
 
 
 //----------------------------------------------------------------------//
@@ -843,7 +831,7 @@ template <typename number>
 inline
 unsigned int SparseMatrixEZ<number>::m () const
 {
-  return row_start.size() - 1;
+  return row_info.size();
 };
 
 
@@ -864,13 +852,14 @@ const SparseMatrixEZ<number>::Entry* SparseMatrixEZ<number>::locate (
   Assert (row<m(), ExcIndexRange(row,0,m()));
   Assert (col<n(), ExcIndexRange(col,0,n()));
 
-  const unsigned int end = row_start[row+1];
-  for (unsigned int i=row_start[row];i<end;++i)
+  const RowInfo& r = row_info[row];
+  const unsigned int end = r.start + r.length;
+  for (unsigned int i=r.start;i<end;++i)
     {
       const Entry * const entry = &data[i];
       if (entry->column == col)
        return entry;
-      if (entry->column == Entry::invalid_entry)
+      if (entry->column == Entry::invalid)
        return 0;
     }
   return 0;
@@ -887,63 +876,82 @@ SparseMatrixEZ<number>::Entry* SparseMatrixEZ<number>::allocate (
   Assert (row<m(), ExcIndexRange(row,0,m()));
   Assert (col<n(), ExcIndexRange(col,0,n()));
 
-  const unsigned int end = row_start[row+1];
-  for (unsigned int i=row_start[row];i<end;++i)
+  RowInfo& r = row_info[row];
+  const unsigned int end = r.start + r.length;
+
+  unsigned int i = r.start;
+  if (r.diagonal != RowInfo::invalid_diagonal && col >= row)
+    i += r.diagonal;
+                                  // Find position of entry
+  while (i<end && data[i].column < col) ++i;
+  
+  Entry* entry = &data[i];
+                                  // entry found
+  if (entry->column == col)
+    return entry;
+
+                                  // Now, we must insert the new
+                                  // entry and move all successive
+                                  // entries back.
+  
+                                  // If no more space is available
+                                  // for this row, insert new
+                                  // elements into the vector.
+  if (row != row_info.size()-1)
     {
-      Entry* entry = &data[i];
-                                      // entry found
-      if (entry->column == col)
-       return entry;
-                                      // entry does not exist,
-                                      // create it
-      if (entry->column > col)
+      if (end >= row_info[row+1].start)
+       {
+                                          // Insert new entries
+         data.insert(data.begin()+end, increment, Entry());
+         entry = &data[i];
+                                          // Update starts of
+                                          // following rows
+         for (unsigned int rn=row+1;rn<row_info.size();++rn)
+           row_info[rn].start += increment;
+       }
+    } else {
+      if (end >= data.size())
        {
-                                          // Save original entry
-         Entry temp = data[i];
-                                          // Insert new entry here to
-                                          // make sure all entries
-                                          // are ordered by column
-                                          // index
-         entry->column = col;
-         entry->value = 0;
-
-                                          // Move all entries in this
-                                          // row up by one
-         for (unsigned int j = i+1;j<end;++j)
-           {
-             Entry temp2 = data[j];
-             data[j] = temp;
-             temp = temp2;
-             if (temp.column == Entry::invalid_entry)
-               break;
-           }
-                                          // Extend row if there is
-                                          // still a valid entry left
-         if (temp.column != Entry::invalid_entry)
-           {
-                                              // Insert new entries
-             data.insert(data.begin()+end, increment, Entry());
-                                              // Update starts of
-                                              // following rows
-             for (unsigned int r=row+1;r<row_start.size();++r)
-               row_start[r] += increment;
-                                              // Insert missing entry
-             data[end] = temp;
-           }
-         return entry;
+                                          // Here, appending a block
+                                          // does not increase
+                                          // performance.
+         data.push_back(Entry());
+         entry = &data[i];
        }
     }
-                                  // Row is full, but entry not
-                                  // found. Therefore, we must
-                                  // extend again
-  data.insert(data.begin()+end, increment, Entry());
-                                  // Update starts of
-                                  // following rows
-  for (unsigned int r=row+1;r<row_start.size();++r)
-    row_start[r] += increment;
-                                  // Insert missing entry
-  data[end].column = col;
-  return &data[end];
+                                  // Save original entry
+  Entry temp = *entry;
+                                  // Insert new entry here to
+                                  // make sure all entries
+                                  // are ordered by column
+                                  // index
+  entry->column = col;
+  entry->value = 0;
+                                  // Update row_info
+  ++r.length;
+  if (col == row)
+    r.diagonal = i - r.start;
+  else if (col<row && r.diagonal!= RowInfo::invalid_diagonal)
+    ++r.diagonal;
+
+  if (i == end)
+      return entry;
+  
+                                  // Move all entries in this
+                                  // row up by one
+  for (unsigned int j = i+1; j < end; ++j)
+    {
+                                      // There should be no invalid
+                                      // entry below end
+      Assert (data[j].column != Entry::invalid, ExcInternalError());
+      Entry temp2 = data[j];
+      data[j] = temp;
+      temp = temp2;
+    }
+  Assert (data[end].column == Entry::invalid, ExcInternalError());
+  data[end] = temp;
+
+  return entry;
 }
 
 
index f3c5d5694a5c9ceae7fbc762fbaf4357bb66faad..f54b66ec1da8cfce5fa3f71e745d00732638bd90 100644 (file)
@@ -49,21 +49,22 @@ SparseMatrixEZ<number>::reinit(unsigned int n_rows,
                               unsigned int default_row_length,
                               unsigned int default_increment)
 {
-  if (default_row_length == Entry::invalid_entry)
+  clear();
+  if (default_row_length == Entry::invalid)
     default_row_length = 5;
-  if (default_increment == Entry::invalid_entry)
+  if (default_increment == Entry::invalid)
     default_increment = 4;
   if (default_increment == 0)
     default_increment = 4;
   increment = default_increment;
   
   n_columns = n_cols;
-  row_start.resize(n_rows+1);
+  row_info.resize(n_rows);
   data.reserve(default_row_length * n_rows + n_rows * increment);
   data.resize(default_row_length * n_rows);
 
-  for (unsigned int i=0;i<=n_rows;++i)
-    row_start[i] = i * default_row_length;
+  for (unsigned int i=0;i<n_rows;++i)
+    row_info[i].start = i * default_row_length;
 }
 
 
@@ -72,7 +73,7 @@ void
 SparseMatrixEZ<number>::clear()
 {
   n_columns = 0;
-  row_start.resize(0);
+  row_info.resize(0);
   data.resize(0);
 }
 
@@ -94,18 +95,20 @@ SparseMatrixEZ<number>::vmult (Vector<somenumber>& dst,
   Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
 
-  const unsigned int end_row = row_start.size() - 1;
-  for (unsigned int i = 0; i < end_row;++i)
+  const unsigned int end_row = row_info.size();
+  for (unsigned int row = 0; row < end_row; ++row)
     {
-      unsigned int index = row_start[i];
-      unsigned int end = row_start[i+1];
+      const RowInfo& ri = row_info[row];
+      typename std::vector<Entry>::const_iterator
+       entry = data.begin() + ri.start;
       double s = 0.;
-      for (;index != end && data[index].column != Entry::invalid_entry;++index)
+      for (unsigned short i=0;i<ri.length;++i,++entry)
        {
-         const Entry& entry = data[index];
-         s += entry.value * src(entry.column);
+         Assert (entry->column != Entry::invalid,
+                 ExcInternalError());
+         s += entry->value * src(entry->column);
        }
-      dst(i) = s;
+      dst(row) = s;
     }
 }
 
@@ -130,18 +133,20 @@ SparseMatrixEZ<number>::vmult_add (Vector<somenumber>& dst,
   Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
 
-  const unsigned int end_row = row_start.size() - 1;
-  for (unsigned int i = 0; i < end_row;++i)
+  const unsigned int end_row = row_info.size();
+  for (unsigned int row = 0; row < end_row; ++row)
     {
-      unsigned int index = row_start[i];
-      unsigned int end = row_start[i+1];
+      const RowInfo& ri = row_info[i];
+      typename std::vector<Entry>::const_iterator
+       entry = data.begin() + ri.start;
       double s = 0.;
-      for (;index != end && data[index].column != Entry::invalid_entry;++index)
+      for (unsigned short i=0;i<ri.length;++i,++entry)
        {
-         const Entry& entry = data[index];
-         s += entry.value * src(entry.column);
+         Assert (entry->column != Entry::invalid,
+                 ExcInternalError());
+         s += entry->value * src(entry->column);
        }
-      dst(i) += s;
+      dst(row) += s;
     }
 }
 
@@ -154,19 +159,157 @@ SparseMatrixEZ<number>::Tvmult_add (Vector<somenumber>& dst,
   Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size()));
   Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size()));
 
-  const unsigned int end_row = row_start.size() - 1;
-  for (unsigned int i = 0; i < end_row;++i)
+  const unsigned int end_row = row_info.size();
+  for (unsigned int row = 0; row < end_row; ++row)
     {
-      unsigned int index = row_start[i];
-      unsigned int end = row_start[i+1];
-      for (;index != end && data[index].column != Entry::invalid_entry;++index)
+      const RowInfo& ri = row_info[row];
+      typename std::vector<Entry>::const_iterator
+       entry = data.begin() + ri.start;
+      for (unsigned short i=0;i<ri.length;++i,++entry)
        {
-         const Entry& entry = data[index];
-         dst(entry.column) += entry.value * src(i);
+         Assert (entry->column != Entry::invalid,
+                 ExcInternalError());
+         dst(entry->column) += entry->value * src(row);
        }
     }
 }
 
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrixEZ<number>::precondition_Jacobi (Vector<somenumber>       &dst,
+                                            const Vector<somenumber> &src,
+                                            const number              om) const
+{
+  Assert (m() == n(), ExcNoSquareMatrix());
+  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+  somenumber              *dst_ptr = dst.begin();
+  const somenumber        *src_ptr = src.begin();
+  typename std::vector<RowInfo>::const_iterator ri = row_info.begin();
+  const typename std::vector<RowInfo>::const_iterator end = row_info.end();
+  
+  for (; ri != end; ++dst_ptr, ++src_ptr, ++ri)
+    {
+      Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
+      *dst_ptr = om * *src_ptr / data[ri->start + ri->diagonal].value;
+    }
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrixEZ<number>::precondition_SOR (Vector<somenumber>       &dst,
+                                         const Vector<somenumber> &src,
+                                         const number              om) const
+{
+  Assert (m() == n(), ExcNoSquareMatrix());
+  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+  somenumber       *dst_ptr = dst.begin();
+  const somenumber *src_ptr = src.begin();
+  typename std::vector<RowInfo>::const_iterator ri = row_info.begin();
+  const typename std::vector<RowInfo>::const_iterator end = row_info.end();
+  
+  for (; ri != end; ++dst_ptr, ++src_ptr, ++ri)
+    {
+      Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
+      number s = *src_ptr;
+      const unsigned int end_row = ri->start + ri->diagonal;
+      for (unsigned int i=ri->start;i<end_row;++i)
+       s -= data[i].value * dst(data[i].column);
+      
+      *dst_ptr = om * s / data[ri->start + ri->diagonal].value;
+    }
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrixEZ<number>::precondition_TSOR (Vector<somenumber>       &dst,
+                                         const Vector<somenumber> &src,
+                                         const number              om) const
+{
+  Assert (m() == n(), ExcNoSquareMatrix());
+  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+  somenumber       *dst_ptr = dst.begin()+dst.size()-1;
+  const somenumber *src_ptr = src.begin()+src.size()-1;
+  typename std::vector<RowInfo>::const_reverse_iterator
+    ri = row_info.rbegin();
+  const typename std::vector<RowInfo>::const_reverse_iterator
+    end = row_info.rend();
+  
+  for (; ri != end; --dst_ptr, --src_ptr, ++ri)
+    {
+      Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
+      number s = *src_ptr;
+      const unsigned int end_row = ri->start + ri->length;
+      for (unsigned int i=ri->start+ri->diagonal+1;i<end_row;++i)
+       s -= data[i].value * dst(data[i].column);
+      
+      *dst_ptr = om * s / data[ri->start + ri->diagonal].value;
+    }
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrixEZ<number>::precondition_SSOR (Vector<somenumber>       &dst,
+                                          const Vector<somenumber> &src,
+                                          const number              om) const
+{
+  Assert (m() == n(), ExcNoSquareMatrix());
+  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+  somenumber       *dst_ptr = dst.begin();
+  const somenumber *src_ptr = src.begin();
+  typename std::vector<RowInfo>::const_iterator ri;
+  const typename std::vector<RowInfo>::const_iterator end = row_info.end();
+
+                                  // Forward
+  for (ri = row_info.begin(); ri != end; ++dst_ptr, ++src_ptr, ++ri)
+    {
+      Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
+      number s = *src_ptr;
+      const unsigned int end_row = ri->start + ri->diagonal;
+      for (unsigned int i=ri->start;i<end_row;++i)
+       s -= om * data[i].value * dst(data[i].column);
+      
+      *dst_ptr = s / data[ri->start + ri->diagonal].value;
+    }
+                                  // Diagonal
+  dst_ptr = dst.begin();
+  for (ri = row_info.begin(); ri != end; ++dst_ptr, ++ri)
+    *dst_ptr *= (2.-om) * data[ri->start + ri->diagonal].value;
+
+                                  // Backward
+  typename std::vector<RowInfo>::const_reverse_iterator rri;
+   const typename std::vector<RowInfo>::const_reverse_iterator
+     rend = row_info.rend();
+   dst_ptr = dst.begin()+dst.size()-1;
+   for (rri = row_info.rbegin(); rri != rend; --dst_ptr, ++rri)
+     {
+      const unsigned int end_row = rri->start + rri->length;
+      for (unsigned int i=rri->start+rri->diagonal+1;i<end_row;++i)
+       *dst_ptr -= om * data[i].value * dst(data[i].column);
+      *dst_ptr /= data[rri->start + rri->diagonal].value;
+     }
+};
+
+
+
 template <typename number>
 unsigned int
 SparseMatrixEZ<number>::memory_consumption() const
index df6e0b4b0860823a5e4d01f035d829eea79998e3..0775afe138c5eea44e012fa9552ef49f87aa5978 100644 (file)
@@ -37,28 +37,31 @@ check_vmult_quadratic(const MATRIX& A,
   Vector<double> f(A.m());
   GrowingVectorMemory<> mem;
 
-  SolverControl control(20, 1.e-3, false);
+  SolverControl control(10, 1.e-13, false);
   SolverRichardson<> rich(control, mem, .01);
+  SolverRichardson<> prich(control, mem, 1.);
   PreconditionIdentity prec;
-
+  PreconditionJacobi<MATRIX> jacobi;
+  jacobi.initialize(A, .5);
+  PreconditionSOR<MATRIX> sor;
+  sor.initialize(A, 1.2);
+  PreconditionSSOR<MATRIX> ssor;
+  ssor.initialize(A, 1.2);
+  
   u = 0.;
   f = 1.;
 
-  try
-    {
-      rich.solve(A, u, f, prec);
-    }
-  catch (...)
-    {
-    }
+  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 (...) {}
+
+  u = 0.;
   deallog << "Transpose" << std::endl;
-  try
-    {
-      rich.Tsolve(A, u, f, prec);
-    }
-  catch (...)
-    {
-    }
+  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 (...) {}
   deallog.pop();
 }
 
@@ -72,13 +75,15 @@ int main()
                                   // Switch between regression test
                                   // and benchmark
 #ifdef DEBUG  
-  deallog.depth_console(0);
-  const unsigned int size = 10;
+  deallog.depth_console(3);
+  const unsigned int size = 5;
+  const unsigned int row_length = 3;
 #else
   deallog.depth_console(1000);
   deallog.log_execution_time(true);
   deallog.log_time_differences(true);
   const unsigned int size = 500;
+  const unsigned int row_length = 9;
 #endif
   
   FDMatrix testproblem (size, size);
@@ -90,12 +95,12 @@ int main()
   structure.compress();
   SparseMatrix<double>  A(structure);
   deallog << "Assemble" << std::endl;
-  testproblem.five_point(A);
+  testproblem.five_point(A, true);
   check_vmult_quadratic(A, "5-SparseMatrix<double>");
 
-  SparseMatrixEZ<double> E(dim,dim,5,1);
+  SparseMatrixEZ<double> E(dim,dim,row_length,2);
   deallog << "Assemble" << std::endl;
-  testproblem.five_point(E);
+  testproblem.five_point(E, true);
   check_vmult_quadratic(E, "5-SparseMatrixEZ<double>");
 
   A.clear();
@@ -109,7 +114,7 @@ int main()
   check_vmult_quadratic(A, "9-SparseMatrix<double>");
 
   E.clear();
-  E.reinit(dim,dim,9,2);
+  E.reinit(dim,dim,row_length,2);
   deallog << "Assemble" << std::endl;
   testproblem.nine_point(E);
   check_vmult_quadratic(E, "9-SparseMatrixEZ<double>");
index 5b1336216369de96257a7c0f829c45a063441ce2..7ca6bdb4c308400ef86db9570c155d54fde9df08 100644 (file)
@@ -114,7 +114,7 @@ FDMatrix::nine_point_structure(SparsityPattern& structure) const
 
 template<typename MATRIX>
 void
-FDMatrix::nine_point(MATRIX& A) const
+FDMatrix::nine_point(MATRIX& A, bool nonsymmetric) const
 {
    for(unsigned int i=0;i<=ny-2;i++)
     {
@@ -170,7 +170,7 @@ FDMatrix::nine_point(MATRIX& A) const
 
 template<typename MATRIX>
 void
-FDMatrix::five_point(MATRIX& A) const
+FDMatrix::five_point(MATRIX& A, bool nonsymmetric) const
 {
    for(unsigned int i=0;i<=ny-2;i++)
     {
@@ -178,10 +178,16 @@ FDMatrix::five_point(MATRIX& A) const
        {
                                           // Number of the row to be entered
          unsigned int row = j+(nx-1)*i;
-         A.set(row, row, 4.);
+         if (nonsymmetric)
+           A.set(row, row, 5.);
+         else
+           A.set(row, row, 4.);
          if (j>0)
            {
-             A.set(row-1, row, -1.);
+             if (nonsymmetric)
+               A.set(row-1, row, -2.);
+             else
+               A.set(row-1, row, -1.);
              A.set(row, row-1, -1.);
            }
          if (j<nx-2)
index d671d9c4dd0705d45dbba1481ae27ee98eef466b..bbff9c7dcf8ab6dd8b6990e78d087d8c1e85d7e2 100644 (file)
@@ -30,13 +30,13 @@ class FDMatrix
                                      * Fill the matrix with values.
                                      */
     template <typename MATRIX>
-    void five_point(MATRIX&) const;
+    void five_point(MATRIX&, bool nonsymmetric = false) const;
 
                                     /**
                                      * Fill the matrix with values.
                                      */
     template <typename MATRIX>
-    void nine_point(MATRIX&) const;
+    void nine_point(MATRIX&, bool nonsymmetric = false) const;
 
     template <typename number>
     void gnuplot_print(std::ostream&, const Vector<number>&) const;
@@ -52,4 +52,3 @@ class FDMatrix
                                      */
     unsigned int ny;
 };
-

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.