]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Mostly doc updates.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 17 Jul 1999 21:57:49 +0000 (21:57 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 17 Jul 1999 21:57:49 +0000 (21:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@1587 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_vanka.h
deal.II/lac/include/lac/sparse_vanka.templates.h

index 9056e40c8ad086b629aac056e52718461acfa502..ea46c58104fe4344f4d4d946522a2d11fd7788e5 100644 (file)
@@ -290,14 +290,15 @@ class SparseMatrixStruct : public Subscriptor
                                     /**
                                      * Number of entries in a specific row.
                                      */
-    unsigned int row_length(unsigned int row) const;
+    unsigned int row_length (const unsigned int row) const;
 
                                     /**
                                      * Access to column nuber field.
                                      * Return the column number of
                                      * the #index#th entry in #row#.
                                      */
-    unsigned int column_number(unsigned int row, unsigned int index) const;
+    unsigned int column_number (const unsigned int row,
+                               const unsigned int index) const;
 
                                     /**
                                      * Compute the bandwidth of the matrix
@@ -1202,7 +1203,7 @@ SparseMatrixStruct::get_column_numbers () const
 
 inline
 unsigned int
-SparseMatrixStruct::row_length(unsigned int row) const
+SparseMatrixStruct::row_length (const unsigned int row) const
 {
   Assert(row<rows, ExcIndexRange(row,0,rows));
   return rowstart[row+1]-rowstart[row];
@@ -1210,7 +1211,8 @@ SparseMatrixStruct::row_length(unsigned int row) const
 
 inline
 unsigned int
-SparseMatrixStruct::column_number(unsigned int row, unsigned int index) const
+SparseMatrixStruct::column_number (const unsigned int row,
+                                  const unsigned int index) const
 {
   Assert(row<rows, ExcIndexRange(row,0,rows));
   Assert(index<row_length(row), ExcIndexRange(index,0,row_length(row)));
index 3a75a7ef92ce721c02a1162e61e53174616c20a9..8e638aeafdb87029816675b4bee1726521024000 100644 (file)
@@ -1,8 +1,11 @@
-// $Id$
-// Copyright Guido Kanschat, 1999
+/*----------------------------   sparse_vanka.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __sparse_vanka_H
+#define __sparse_vanka_H
+/* Copyright Guido Kanschat, 1999 */
+/*----------------------------   sparse_vanka.h     ---------------------------*/
+
 
-#ifndef __lac_sparse_vanka_H
-#define __lac_sparse_vanka_H
 
 #include <base/smartpointer.h>
 #include <lac/forward-declarations.h>
 /**
  * Point-wise Vanka preconditioning.
  * This class does Vanka preconditioning  on a point-wise base.
- * Vanka preconditioners are used for saddle point problems. There the
- * application of Jacobi or Gauß-Seidel methods is impossible, because
- * the diagonal elements are zero in the rows of the Lagrange multiplier.
+ * Vanka preconditioners are used for saddle point problems like Stoke's
+ * problem or problems arising in optimization where Lagrange multiplier
+ * occur and let Netwon's matrix have a zero block. With these matrices the
+ * application of Jacobi or Gauss-Seidel methods is impossible, because
+ * some diagonal elements are zero in the rows of the Lagrange multiplier.
+ * The approach of Vanka is to solve a small (usually indefinite) system
+ * of equations for each Langrange multiplie variable (we will also call
+ * the pressure in Stoke's equation a Langrange multiplier since it
+ * can be interpreted as such).
  *
- * It is constructed initializing a vector of indices to the degrees of
- * freedom of the Lagrange multiplier.
+ * Objects of this class are constructed by passing a vector of indices
+ * of the degrees of freedom of the Lagrange multiplier. In the actual
+ * preconditioning method, these rows are traversed in the order in which
+ * the appear in the matrix. Since this is a Gauß-Seidel like procedure,
+ * remember to have a good ordering in advance (for transport dominated
+ * problems, Cuthill-McKee algorithms are a good means for this, if points
+ * on the inflow boundary are chosen as starting points for the renumbering).
  *
- * In the actual preconditioning method, these rows are traversed in
- * original order. Since this is a Gauß-Seidel like procedure,
- * remember to have a good ordering in advance.
+ * For each selected degree of freedom, a local system of equations is built
+ * by the degree of freedom itself and all other values coupling immediately,
+ * i.e. the set of degrees of freedom considered for the local system of
+ * equations for degree of freedom #i# is #i# itself and all #j# such that
+ * the element #(i,j)# is a nonzero entry in the sparse matrix under 
+ * consideration. The elements #(j,i)# are not considered. We now pick all
+ * matrix entries from rows and columns out of the set of degrees of freedom
+ * just described out of the global matrix and put it into a local matrix,
+ * which is subsequently inverted. This system may be of different size for
+ * each degree of freedom, depending for example on the local neighborhood of
+ * the respective node on a computational grid.
  *
- * For each row, a local system of equations is built by the degree of
- * freedom itself and all other values coupling immediately. The right
- * hand side is augmented by all further couplings.
+ * The right hand side is built up in the same way, i.e. by copying all entries
+ * that coupled with the one under present consideration, but it is augmented
+ * by all degrees of freedom coupling with the degrees from the set described
+ * above (i.e. the DoFs coupling second order to the present one).
  *
  * This local system is solved and the values are updated into the
  * destination vector.
  *
  * Remark: the Vanka method is a non-symmetric preconditioning method.
- * @author Guido Kanschat */
+ *
+ *
+ * \subsection{Example of Use}
+ * This little example is taken from a program doing parameter optimization.
+ * The Lagrange multiplier is the third component of the finite element
+ * used. The system is solved by the GMRES method.
+ * \begin{verbatim}
+ *                        // tag the Lagrange multiplier variable
+ *    vector<bool> signature(3);
+ *    signature[0] = signature[1] = false;
+ *    signature[2] = true;
+ *
+ *                        // tag all dofs belonging to the
+ *                        // Lagrange multiplier
+ *    vector<bool> selected_dofs (dof.n_dofs(), false);
+ *    DoFTools::extract_dofs(dof, signature, p_select);
+ *                        // create the Vanka object
+ *    SparseVanka<double> vanka (global_matrix, selected_dofs);
+ *
+ *                        // create the solver
+ *    SolverGMRES<PreconditionedSparseMatrix<double>,
+ *                Vector<double> >    gmres(control,memory,504);
+ *    
+ *                       // solve
+ *    gmres.solve (global_matrix, solution, right_hand_side,
+ *                vanka);
+ * \end{verbatim}
+ *
+ * @author Guido Kanschat, 1999
+ */
 template<typename number>
 class SparseVanka
 {
@@ -97,6 +149,11 @@ class SparseVanka
                                      * are built up and inverted every time.
                                      */
     void conserve_memory();
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcMatrixNotSquare);
     
   private:
                                     /**
@@ -119,6 +176,10 @@ class SparseVanka
     mutable vector<SmartPointer<FullMatrix<float> > > inverses;
 };
 
+
+
+/* ---------------------------- Inline functions -----------------------*/
+
 template<typename number>
 template<typename number2>
 inline
@@ -131,4 +192,9 @@ SparseVanka<number>::operator() (Vector<number2>& dst,
 }
 
 
+
+
+/*----------------------------   sparse_vanka.h     ---------------------------*/
+/* end of #ifndef __sparse_vanka_H */
 #endif
+/*----------------------------   sparse_vanka.h     ---------------------------*/
index 3f100bcd729dfe40c7c85cbc0ac3c26c205ae25f..dd24436823b4d41035697b984a11433f90ade130 100644 (file)
@@ -6,12 +6,21 @@
 
 #include <map>
 
+
+
 template<typename number>
-SparseVanka<number>::SparseVanka(const SparseMatrix<number>M,
-                                const bit_vectorselected)
+SparseVanka<number>::SparseVanka(const SparseMatrix<number> &M,
+                                const bit_vector           &selected)
                :
-               matrix(&M), selected(selected), conserve_mem(false), inverses(M.m(),0)
-{}
+               matrix(&M),
+               selected(selected),
+               conserve_mem(false),
+               inverses(M.m(),0)
+{
+  Assert (M.m() == M.n(),
+         ExcMatrixNotSquare ());
+}
+
 
 template<typename number>
 SparseVanka<number>::~SparseVanka()
@@ -33,68 +42,112 @@ void
 SparseVanka<number>::forward(Vector<number2>& dst,
                           const Vector<number2>& src) const
 {
+                                  // first define an alias to the sparsity
+                                  // pattern of the matrix, since this
+                                  // will be used quite often
+  const SparseMatrixStruct &structure
+    = matrix->get_sparsity_pattern();
+                                  // space to be used for local systems
   FullMatrix<float> local_matrix;
+
+                                  // traverse all rows of the matrix
   for (unsigned int row=0; row< matrix->m() ; ++row)
     {
-      bool build_matrix = true;
-      
+                                      // but skip those that were not selected
       if (!selected[row])
        continue;
 
-             
-      const SparseMatrixStruct& structure = matrix->get_sparsity_pattern();
-      unsigned int n = structure.row_length(row);
+                                      // shall we reconstruct the system
+                                      // of equations for this DoF?
+      bool build_matrix = true;
+      const unsigned int row_length = structure.row_length(row);
       
       if (!conserve_mem)
        {
          if (inverses[row] == 0)
+                                            // inverse not yet built
            {
-             FullMatrix<float>p = new FullMatrix<float>;
+             FullMatrix<float> *p = new FullMatrix<float>;
              inverses[row] = p;
            } else {
+                                              // inverse already built
              build_matrix = false;
            }
        }
 
-      FullMatrix<float>& A = (conserve_mem) ? local_matrix : (*inverses[row]);
+                                      // alias to the matrix which is
+                                      // to be used
+      FullMatrix<float> &A = (conserve_mem ?
+                             local_matrix :
+                             (*inverses[row]));
 
-      if (build_matrix) A.reinit(n);
-      
-      Vector<float> b(n);
-      Vector<float> x(n);
+      if (build_matrix)
+       A.reinit(row_length);
       
-      map<unsigned int, unsigned int> local_index;
+      Vector<float> b(row_length);
+      Vector<float> x(row_length);
 
-                                      // Build local index
-
-      for (unsigned int i=0;i<n;++i)
+                                      // mapping between:
+                                      // 1 column number of all
+                                      //   entries in this row, and
+                                      // 2 the position within this
+                                      //   row (as stored in the
+                                      //   sparsematrixstruct object
+                                      //
+                                      // since we do not explicitely
+                                      // consider nonsysmmetric sparsity
+                                      // patterns, the first element
+                                      // of each entry simply denotes
+                                      // all degrees of freedom that
+                                      // couple with #row#.
+      map<unsigned int, unsigned int> local_index;
+      for (unsigned int i=0; i<row_length; ++i)
        local_index.insert(pair<unsigned int, unsigned int>
                           (structure.column_number(row, i), i));
-
-//       for (map<unsigned int, unsigned int>::iterator is=local_index.begin();
-//        is!=local_index.end();++is)
-//     cerr << "map " << is->first << '\t' << is->second << endl;
       
-                                      // Build local matrix
-
-      for (map<unsigned int, unsigned int>::iterator is=local_index.begin();
-          is!=local_index.end();++is)
+                                      // Build local matrix and rhs
+      for (map<unsigned int, unsigned int>::const_iterator is=local_index.begin();
+          is!=local_index.end(); ++is)
        {
-         unsigned int irow = is->first;
-         unsigned int i = is->second;
-         unsigned int n = structure.row_length(irow);
-         
+                                          // irow loops over all DoFs that
+                                          // couple with the present DoF
+         const unsigned int irow = is->first;
+                                          // index of DoF irow in the matrix
+                                          // row corresponding to DoF #row#.
+                                          // runs between 0 and row_length
+         const unsigned int i = is->second;
+                                          // number of DoFs coupling to
+                                          // irow (including irow itself)
+         const unsigned int irow_length = structure.row_length(irow);
+
+                                          // copy rhs
          b(i) = src(irow);
-         
-         for (unsigned int j=0;j<n;++j)
+
+                                          // for all the DoFs that irow
+                                          // couples with
+         for (unsigned int j=0; j<irow_length; ++j)
            {
-             unsigned int col = structure.column_number(irow, j);
-             map<unsigned int, unsigned int>::iterator js
+                                              // col is the number of
+                                              // this dof
+             const unsigned int col = structure.column_number(irow, j);
+                                              // find out whether this DoF
+                                              // (that couples with #irow#,
+                                              // which itself couples with
+                                              // #row#) also couples with
+                                              // #row#.
+             const map<unsigned int, unsigned int>::const_iterator js
                = local_index.find(col);
+                                              // if not, then still use
+                                              // this dof to modify the rhs
+                                              //
+                                              // note that if so, we already
+                                              // have copied the entry above
              if (js == local_index.end())
                {
                  b(i) -= matrix->raw_entry(irow,j) * dst(col);
                } else {
+                                                  // if so, then build the
+                                                  // matrix out of it
                  if (build_matrix)
                    A(i,js->second) = matrix->raw_entry(irow,j);
                }
@@ -106,11 +159,11 @@ SparseVanka<number>::forward(Vector<number2>& dst,
       A.vmult(x,b);
       
                                       // Distribute new values
-      for (map<unsigned int, unsigned int>::iterator is=local_index.begin();
-          is!=local_index.end();++is)
+      for (map<unsigned int, unsigned int>::const_iterator is=local_index.begin();
+          is!=local_index.end(); ++is)
        {
-         unsigned int irow = is->first;
-         unsigned int i = is->second;
+         const unsigned int irow = is->first;
+         const unsigned int i = is->second;
          dst(irow) = x(i);
        }
     }

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.