FullMatrix::symmetrize

git-svn-id: https://svn.dealii.org/trunk@5030 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html
index a1c704260c838e3d86d41652ffac6b08d7083984..2bae941a21cd285ff0b9f39f27a1da14369fc796 100644 (file)
--- a/deal.II/doc/news/2001/c-3-1.html
+++ b/deal.II/doc/news/2001/c-3-1.html
@@ -356,7 +356,12 @@ documentation, etc</a>.
 
 <ol>
   <li> <p>
+       New: Function <code class="class">FullMatrix::symmetrize()</code>.
+       <br>
+       (WB 2001/09/20)
+       </p>
 
+  <li> <p>
        Improved: the stopping criterion of <code
        class="class">SolverBicgstab</code> without computing the exact
        residual is implemented

diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h
index 40a62d9379de4db3c0711f2d2e30996ea7df1706..de5ac4995bd8714cff32b2c58f144cefc071b135 100644 (file)
--- a/deal.II/lac/include/lac/full_matrix.h
+++ b/deal.II/lac/include/lac/full_matrix.h
@@ -318,6 +318,17 @@ class FullMatrix : public vector2d<number>
     number2 matrix_scalar_product (const Vector<number2> &u,
                                   const Vector<number2> &v) const;
 
+                                    /**
+                                     * Symmetrize the matrix by
+                                     * forming the mean value between
+                                     * the existing matrix and its
+                                     * transpose, $A = \frac 12(A+A^T)$.
+                                     *
+                                     * Obviously the matrix must be
+                                     * square for this operation.
+                                     */
+    void symmetrize ();
+    
                                     /**
                                      * Return the $l_1$-norm of the matrix, i.e.
                                      * $|M|_1=max_{all columns j}\sum_{all 
@@ -374,13 +385,13 @@ class FullMatrix : public vector2d<number>
     
                                     /**
                                      * A=Inverse(A). Inversion of
-                                     * this matrix by
-                                     * Gauss-Jordan algorithm with
-                                     * partial pivoting.  This
-                                     * process is well-behaved for
-                                     * positive definite matrices,
-                                     * but be aware of round-off
-                                     * errors in the indefinite case.
+                                     * this matrix by Gauss-Jordan
+                                     * algorithm with partial
+                                     * pivoting.  This process is
+                                     * well-behaved for positive
+                                     * definite matrices, but be
+                                     * aware of round-off errors in
+                                     * the indefinite case.
                                      *
                                      * The numerical effort to invert
                                      * an @p{n x n} matrix is of the
@@ -389,12 +400,14 @@ class FullMatrix : public vector2d<number>
     void gauss_jordan ();
 
                                     /**
-                                      * Computes the determinant of a matrix.
-                                      * This is only implemented for one two and
-                                      * three dimensions, since for higher
-                                      * dimensions the numerical work explodes.
-                                      * Obviously, the matrix needs to be square
-                                      * for this function.
+                                      * Computes the determinant of a
+                                      * matrix.  This is only
+                                      * implemented for one, two, and
+                                      * three dimensions, since for
+                                      * higher dimensions the
+                                      * numerical work explodes.
+                                      * Obviously, the matrix needs to
+                                      * be square for this function.
                                       */
     double determinant () const;
 
@@ -482,15 +495,17 @@ class FullMatrix : public vector2d<number>
                     const Vector<number3>& b) const;
 
                                     /**
-                                     * Forward elimination of lower triangle.
-                                     * Inverts the lower triangle of a
-                                     * quadratic matrix.
+                                     * Forward elimination of lower
+                                     * triangle.  Inverts the lower
+                                     * triangle of a quadratic
+                                     * matrix.
                                      *
-                                     * If the matrix has more columns than rows,
-                                     * this function only operates on the left
-                                     * square submatrix. If there are more rows,
-                                     * the upper square part of the matrix
-                                     * is considered
+                                     * If the matrix has more columns
+                                     * than rows, this function only
+                                     * operates on the left square
+                                     * submatrix. If there are more
+                                     * rows, the upper square part of
+                                     * the matrix is considered
                                      */
     template<typename number2>
     void forward (Vector<number2>       &dst,
@@ -506,12 +521,15 @@ class FullMatrix : public vector2d<number>
 
                                     /**
                                      * QR-factorization of a matrix.
-                                     * The orthogonal transformation Q is
-                                     * applied to the vector y and this matrix.
+                                     * The orthogonal transformation
+                                     * Q is applied to the vector y
+                                     * and this matrix.
                                      *
-                                     * After execution of householder, the upper
-                                     * triangle contains the resulting matrix R,
-                                     * the lower the incomplete factorization
+                                     * After execution of
+                                     * householder, the upper
+                                     * triangle contains the
+                                     * resulting matrix R, the lower
+                                     * the incomplete factorization
                                      * matrices.
                                      */
     template<typename number2>

diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h
index 0a0494c8f10cd80afc8ca376257e12bb2405412d..0fd628dad934ad07572b6189110610c401d6f3fe 100644 (file)
--- a/deal.II/lac/include/lac/full_matrix.templates.h
+++ b/deal.II/lac/include/lac/full_matrix.templates.h
@@ -292,15 +292,16 @@ double FullMatrix<number>::residual (Vector<number2>& dst,
 }
 
 
+
 template <typename number>
 template <typename number2>
-void FullMatrix<number>::forward (Vector<number2>& dst,
-                                 const Vector<number2>& src) const
+void FullMatrix<number>::forward (Vector<number2>       &dst,
+                                 const Vector<number2> &src) const
 {
   Assert (data() != 0, ExcEmptyMatrix());
   
-  Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
-  Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
+  Assert (dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
+  Assert (src.size() == n(), ExcDimensionMismatch(src.size(), n()));
 
   unsigned int i,j;
   unsigned int nu = ( (m()<n()) ? m() : n());
@@ -314,10 +315,11 @@ void FullMatrix<number>::forward (Vector<number2>& dst,
 }
 
 
+
 template <typename number>
 template <typename number2>
-void FullMatrix<number>::backward (Vector<number2>& dst,
-                                  const Vector<number2>& src) const
+void FullMatrix<number>::backward (Vector<number2>       &dst,
+                                  const Vector<number2> &src) const
 {
   Assert (data() != 0, ExcEmptyMatrix());
   
@@ -333,20 +335,12 @@ void FullMatrix<number>::backward (Vector<number2>& dst,
 }
 
 
-/*  template <typename number> */
-/*  template <typename number2> */
-/*  FullMatrix<number>& */
-/*  FullMatrix<number>::operator = (const FullMatrix<number2>& m)  */
-/*  { */
-  
-/*  } */
-
 
 template <typename number>
 template <typename number2>
-void FullMatrix<number>::fill (const FullMatrix<number2>& src,
-                              const unsigned int i,
-                              const unsigned int j)
+void FullMatrix<number>::fill (const FullMatrix<number2> &src,
+                              const unsigned int         i,
+                              const unsigned int         j)
 {
   Assert (n() >= src.n() + j, ExcInvalidDestination(n(), src.n(), j));
   Assert (m() >= src.m() + i, ExcInvalidDestination(m(), src.m(), i));
@@ -598,6 +592,23 @@ number2 FullMatrix<number>::matrix_scalar_product (const Vector<number2> &u,
 };
 
 
+
+template <typename number>
+void
+FullMatrix<number>::symmetrize ()
+{
+  Assert (m() == n(), ExcNotQuadratic());
+  
+  const unsigned int N = m();
+  for (unsigned int i=0; i<N; ++i)
+    for (unsigned int j=i+1; j<N; ++j)
+      {
+       const number t = (el(i,j) + el(j,i)) / 2;
+       el(i,j) = el(j,i) = t;
+      };
+};
+
+
 template <typename number>
 number FullMatrix<number>::l1_norm () const
 {