]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add functions symmetrize to SparsityPattern and SparseMatrix.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 2 Dec 2000 21:30:13 +0000 (21:30 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 2 Dec 2000 21:30:13 +0000 (21:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@3520 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/2000/c-3-0.html
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix.templates.h
deal.II/lac/include/lac/sparsity_pattern.h
deal.II/lac/source/sparsity_pattern.cc

index 218cac5a1072671fa4c7f81d48ceb51e826ba9c2..b7885f62c8be51d55ee21005f509a266217b81e2 100644 (file)
@@ -250,6 +250,15 @@ documentation, etc</a>.
 <h3>lac</h3>
 
 <ol>
+  <li> <p>
+       New: There are now functions 
+       <code class="member">SparsityPattern::symmetrize</code> and
+       <code class="member">SparseMatrix::symmetrize</code> that
+       generate a symmetric matrix from a non-symmetric one.
+       <br>
+       (WB 2000/12/02)
+       </p>
+
   <li> <p>
        New: almost all classes that store data now have a function
        <code class="member">memory_consumption</code> that returns an
index ffd0d5f5cad3e094703bd8c097eb109b6541ee96..09d9ad660aa6ae92a05cd5011a3659646df1a95e 100644 (file)
@@ -254,6 +254,34 @@ class SparseMatrix : public Subscriptor
     void add (const unsigned int i, const unsigned int j,
              const number value);
 
+                                    /**
+                                     * Symmetrize the matrix by
+                                     * forming the mean value between
+                                     * the existing matrix and its
+                                     * transpose, $A = \frac 12(A+A^T)$.
+                                     *
+                                     * This operation assumes that
+                                     * the underlying sparsity
+                                     * pattern represents a symmetric
+                                     * object. If this is not the
+                                     * case, then the result of this
+                                     * operation will not be a
+                                     * symmetric matrix, since it
+                                     * only explicitely symmetrizes
+                                     * by looping over the lower left
+                                     * triangular part for efficiency
+                                     * reasons; if there are entries
+                                     * in the upper right triangle,
+                                     * then these elements are missed
+                                     * in the
+                                     * symmetrization. Symmetrization
+                                     * of the sparsity pattern can be
+                                     * obtain by the
+                                     * @ref{SparsityPattern}@p{::symmetrize}
+                                     * function.
+                                     */
+    void symmetrize ();
+    
                                     /**
                                      * Copy the given matrix to this
                                      * one.  The operation throws an
index c1db80dc5515ac84bb1d37f8250c623225ef932b..4bf1714ddeb327e3500c598c805025fba1cdfaf9 100644 (file)
@@ -40,6 +40,7 @@ SparseMatrix<number>::SparseMatrix () :
 {};
 
 
+
 template <typename number>
 SparseMatrix<number>::SparseMatrix (const SparseMatrix &m) :
                Subscriptor (m),
@@ -53,6 +54,7 @@ SparseMatrix<number>::SparseMatrix (const SparseMatrix &m) :
 };
 
 
+
 template <typename number>
 SparseMatrix<number>&
 SparseMatrix<number>::operator = (const SparseMatrix<number> &m)
@@ -65,6 +67,7 @@ SparseMatrix<number>::operator = (const SparseMatrix<number> &m)
 };
 
 
+
 template <typename number>
 SparseMatrix<number>::SparseMatrix (const SparsityPattern &c) :
                cols(&c),
@@ -75,6 +78,7 @@ SparseMatrix<number>::SparseMatrix (const SparsityPattern &c) :
 };
 
 
+
 template <typename number>
 SparseMatrix<number>::~SparseMatrix ()
 {
@@ -85,6 +89,7 @@ SparseMatrix<number>::~SparseMatrix ()
 };
 
 
+
 template <typename number>
 void
 SparseMatrix<number>::reinit ()
@@ -113,6 +118,7 @@ SparseMatrix<number>::reinit ()
 }
 
 
+
 template <typename number>
 void
 SparseMatrix<number>::reinit (const SparsityPattern &sparsity)
@@ -122,6 +128,7 @@ SparseMatrix<number>::reinit (const SparsityPattern &sparsity)
 };
 
 
+
 template <typename number>
 void
 SparseMatrix<number>::clear ()
@@ -133,6 +140,7 @@ SparseMatrix<number>::clear ()
 };
 
 
+
 template <typename number>
 bool
 SparseMatrix<number>::empty () const
@@ -144,6 +152,7 @@ SparseMatrix<number>::empty () const
 };
 
 
+
 template <typename number>
 unsigned int
 SparseMatrix<number>::n_nonzero_elements () const
@@ -153,6 +162,7 @@ SparseMatrix<number>::n_nonzero_elements () const
 };
 
 
+
 template <typename number>
 unsigned int
 SparseMatrix<number>::n_actually_nonzero_elements () const
@@ -163,6 +173,44 @@ SparseMatrix<number>::n_actually_nonzero_elements () const
 };
 
 
+
+template <typename number>
+void
+SparseMatrix<number>::symmetrize ()
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (cols->rows == cols->cols, ExcMatrixNotSquare());
+  
+  const unsigned int n_rows = m();
+  for (unsigned int row=0; row<n_rows; ++row)
+    {
+                                      // first skip diagonal entry
+      number             *val_ptr = &val[cols->rowstart[row]+1];
+      const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[row]+1];      
+      const number       *const val_end_of_row = &val[cols->rowstart[row+1]];
+
+                                      // treat lower left triangle
+      while ((val_ptr != val_end_of_row) && (*colnum_ptr<row))
+       {
+                                          // compute the mean of this
+                                          // and the transpose value
+         const number mean_value = (*val_ptr +
+                                    val[(*cols)(*colnum_ptr,row)]) / 2.0;
+                                          // set this value and the
+                                          // transpose one to the
+                                          // mean
+         *val_ptr = mean_value;
+         set (*colnum_ptr, row, mean_value);
+
+                                          // advance pointers
+         ++val_ptr;
+         ++colnum_ptr;
+       };
+    };
+};
+
+
+
 template <typename number>
 template <typename somenumber>
 SparseMatrix<number> &
index fbf44987ce2cf0a0140cc006ed8ed0b220afeb4e..c29fe1110bbde6500d9828fc2629b3608838c5e6 100644 (file)
@@ -351,7 +351,8 @@ class SparsityPattern : public Subscriptor
                                      * because the column numbers are
                                      * sorted.
                                      */
-    unsigned int operator() (const unsigned int i, const unsigned int j) const;
+    unsigned int operator() (const unsigned int i, 
+                            const unsigned int j) const;
 
                                     /**
                                      * Add a nonzero entry to the matrix.
@@ -361,7 +362,21 @@ class SparsityPattern : public Subscriptor
                                      * If the entry already exists, nothing
                                      * bad happens.
                                      */
-    void add (const unsigned int i, const unsigned int j);
+    void add (const unsigned int i, 
+             const unsigned int j);
+    
+                                    /**
+                                     * Make the sparsity pattern
+                                     * symmetric by adding the
+                                     * sparsity pattern of the
+                                     * transpose object.
+                                     *
+                                     * This function throws an
+                                     * exception if the sparsity
+                                     * pattern does not represent a
+                                     * square matrix.
+                                     */
+    void symmetrize ();
     
                                     /**
                                      * Print the sparsity of the matrix
index 81bb2999ee0172e9c2ee3823c5fcf306ad183f13..d8752837686bc8cc68b70848f97cbabf8a99d8eb 100644 (file)
@@ -567,7 +567,42 @@ SparsityPattern::add (const unsigned int i,
                                   // wrong: there was not enough space
                                   // in this line
   Assert (false, ExcNotEnoughSpace(i, rowstart[i+1]-rowstart[i]));
-}
+};
+
+
+
+void
+SparsityPattern::symmetrize () 
+{
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
+  Assert (compressed==false, ExcMatrixIsCompressed());
+  Assert (rows==cols, ExcNotSquare());
+
+                                  // loop over all elements presently
+                                  // in the sparsity pattern and add
+                                  // the transpose element. note:
+                                  //
+                                  // 1. that the sparsity pattern
+                                  // changes which we work on
+                                  //
+                                  // 2. that the @p{add} function can
+                                  // be called on elements that
+                                  // already exist without any harm
+  for (unsigned int row=0; row<rows; ++row)
+    for (unsigned int k=rowstart[row]; k<rowstart[row+1]; k++)
+      {
+                                        // check whether we are at
+                                        // the end of the entries of
+                                        // this row. if so, go to
+                                        // next row
+       if (colnums[k] == invalid_entry)
+         break;
+
+                                        // otherwise add the
+                                        // transpose entry
+       add (colnums[k], row);
+      };
+};
 
 
 

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.