]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement block version of Vanka preconditioner.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 1 Feb 2000 14:22:06 +0000 (14:22 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 1 Feb 2000 14:22:06 +0000 (14:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@2305 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 0bbda58ab5eb511d940e0d738cc6f512470db09c..7da845f9e5cf6c78a13b2afa94309db86d644bbb 100644 (file)
@@ -99,7 +99,7 @@
  * whether this might pose some problems in the inversion of the local matrices.
  * Maybe someone would like to check this.
  *
- * @author Guido Kanschat, documentation and extensions by Wolfgang Bangerth; 1999, 2000
+ * @author Guido Kanschat, Wolfgang Bangerth; 1999, 2000
  */
 template<typename number>
 class SparseVanka
@@ -182,6 +182,62 @@ class SparseVanka
                                      * Exception
                                      */
     DeclException0 (ExcMatrixNotSquare);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidRange,
+                   unsigned int, unsigned int,
+                   << "The bounds [" << arg1 << ',' << arg2
+                   << ") do not form a valid range.");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidVectorSize,
+                   unsigned int, unsigned int,
+                   << "The dimensions of vectors and matrices, "
+                   << arg1 << " and " << arg2 << " do not match.");
+
+  protected:
+                                    /**
+                                     * Apply the inverses in the
+                                     * range #[begin,end)# to the
+                                     * #src# vector and move the
+                                     * result into #dst#. Actually,
+                                     * only values of #src# from
+                                     * within the range are taken
+                                     * (all others are set to zero),
+                                     * and only values inside the
+                                     * range are written to #dst#, so
+                                     * the application of this
+                                     * function only does what is
+                                     * announced in the general
+                                     * documentation if the given
+                                     * range is the whole interval.
+                                     *
+                                     * The reason for providing the
+                                     * interval anyway is that in
+                                     * derived classes we may want to
+                                     * apply the preconditioner to
+                                     * blocks of the matrix only, in
+                                     * order to parallelize the
+                                     * application. Then, it is
+                                     * important to only write to
+                                     * some slices of #dst# and only
+                                     * takes values from similar
+                                     * slices of #src#, in order to
+                                     * eliminate the dependencies of
+                                     * threads of each other.
+                                     *
+                                     * The #operator()# of this class
+                                     * of course calls this function
+                                     * with the whole interval
+                                     * #[begin,end)=[0,matrix.m())#.
+                                     */
+    template<typename number2>
+    void apply_preconditioner (Vector<number2>       &dst,
+                              const Vector<number2> &src,
+                              const unsigned int     begin,
+                              const unsigned int     end) const;    
     
   private:
                                     /**
@@ -189,17 +245,17 @@ class SparseVanka
                                      */
     SmartPointer<const SparseMatrix<number> > matrix;
     
-                                    /**
-                                     * Indices of Lagrange
-                                     * multipliers.
-                                     */
-    const vector<bool> &selected;
-    
                                     /**
                                      * Conserve memory flag.
                                      */
     const bool conserve_mem;
 
+                                    /**
+                                     * Indices of those degrees of
+                                     * freedom that we shall work on.
+                                     */
+    const vector<bool> &selected;
+
                                     /**
                                      * Number of threads to be used
                                      * when building the
@@ -253,49 +309,116 @@ class SparseVanka
     void compute_inverse (const unsigned int               row,
                          map<unsigned int, unsigned int> &local_index);
 
+};
+
+
+
+/**
+ * Block version of the sparse Vanka preconditioner. This class
+ * divides the matrix into blocks and works on the diagonal blocks
+ * only, which of course reduces the efficiency as preconditioner, but
+ * is perfectly parallelizable. The constructor takes a parameter into
+ * how many diagonal blocks the matrix shall be subdivided and then
+ * lets the underlying class do the work.
+ *
+ * Division of the matrix is done in a way such that the blocks are
+ * not necessarily of equal size, but such that the number of selected
+ * degrees of freedom for which a local system is to be solved is
+ * equal between blocks. The reason for this strategy to subdivision
+ * is load-balancing for multithreading, but it is necessary to note
+ * that this almost renders the capability as precondition useless if
+ * the degrees of freedom are numbered by component, i.e. all Lagrange
+ * multipliers en bloc.
+ *
+ * This class is probably useless if you don't have a multiprocessor
+ * system, since then the amount of work per preconditioning step is
+ * the same as for the #SparseVanka# class, but preconditioning
+ * properties are worse. On the other hand, if you have a
+ * multiprocessor system, the worse preconditioning quality (leading
+ * to more iterations of the linear solver) usually is well balanced
+ * by the increased speed of application due to the parallelization,
+ * leading to an overall decrease in elapsed wall-time for solving
+ * your linear system. It should be noted that the quality as
+ * preconditioner reduces with growing number of blocks, so there may
+ * be an optimal value (in terms of wall-time per linear solve) for
+ * the number of blocks.
+ *
+ * To facilitate writing portable code, if the number of blocks into
+ * which the matrix is to be subdivided, is set to one, then this
+ * class acts just like the #SparseVanka# class. You may therefore
+ * want to set the number of blocks equal to the number of processors
+ * you have.
+ *
+ * Note that the parallelization is done if #deal.II# was configured
+ * for multithread use and that the number of threads which is spawned
+ * equals the number of blocks. This is reasonable since you will not
+ * want to set the number of blocks unnecessarily large, since, as
+ * mentioned, this reduces the preconditioning properties.
+ *
+ *
+ * \subsection{Typical results}
+ *
+ * As a prototypical test case, we use a nonlinear problem from
+ * optimization, which leads to a series of saddle point problems,
+ * each of which is solved using GMRES with Vanka as
+ * preconditioner. The equation had approx. 850 degrees of
+ * freedom. With the non-blocked version #SparseVanka# (or
+ * #SparseBlockVanka# with #n_blocks==1#), the following numbers of
+ * iterations is needed to solver the linear system in each nonlinear
+ * step: \begin{verbatim} 101 68 64 53 35 21 \end{verbatim} With four
+ * blocks, we need the following numbers of iterations
+ * \begin{verbatim} 124 88 83 66 44 28 \end{verbatim} As can be seen,
+ * more iterations are needed. However, in terms of computing time,
+ * the first version needs 72 seconds wall time (and 79 seconds CPU
+ * time, which is more than wall time since some other parts of the
+ * program were parallelized as well), while the second version needed
+ * 53 second wall time (and 110 seconds CPU time) on a four processor
+ * machine. The total time is in both cases dominated by the linear
+ * solvers. In this case, it is therefore worth while using the
+ * blocked version of the preconditioner if wall time is more
+ * important than CPU time.
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+template<typename number>
+class SparseBlockVanka : public SparseVanka<number>
+{
+  public:
                                     /**
-                                     * Apply the inverses in the
-                                     * range #[begin,end)# to the
-                                     * #src# vector and move the
-                                     * result into #dst#. Actually,
-                                     * only values of #src# from
-                                     * within the range are taken
-                                     * (all others are set to zero),
-                                     * and only values inside the
-                                     * range are written to #dst#, so
-                                     * the application of this
-                                     * function only does what is
-                                     * announced in the general
-                                     * documentation if the given
-                                     * range is the whole interval.
-                                     *
-                                     * The reason for providing the
-                                     * interval anyway is that in
-                                     * derived classes we may want to
-                                     * apply the preconditioner to
-                                     * blocks of the matrix only, in
-                                     * order to parallelize the
-                                     * application. Then, it is
-                                     * important to only write to
-                                     * some slices of #dst# and only
-                                     * takes values from similar
-                                     * slices of #src#, in order to
-                                     * eliminate the dependencies of
-                                     * threads of each other.
-                                     *
-                                     * The #operator()# of this class
-                                     * of course calls this function
-                                     * with the whole interval
-                                     * #[begin,end)=[0,matrix.m())#.
+                                     * Constructor. Pass all
+                                     * arguments except for
+                                     * #n_blocks# to the base class.
                                      */
-    template<typename number2>
-    void apply_preconditioner (Vector<number2>       &dst,
-                              const Vector<number2> &src,
-                              const unsigned int     begin,
-                              const unsigned int     end) const;
-};
+    SparseBlockVanka (const SparseMatrix<number> &M,
+                     const vector<bool>         &selected,
+                     const bool                  conserve_memory = false,
+                     const unsigned int          n_threads       = 1,
+                     const unsigned int          n_blocks        = 1);
 
+                                    /**
+                                     * Apply the preconditioner.
+                                     */
+    template<typename number2>
+    void operator() (Vector<number2>       &dst,
+                    const Vector<number2> &src) const;
+    
+  private:
+                                    /**
+                                     * Store the number of blocks.
+                                     */
+    const unsigned int n_blocks;
 
+                                    /**
+                                     * In this field, we precompute
+                                     * the first and the one after
+                                     * the last index of each
+                                     * block. This computation is
+                                     * done in the constructor, to
+                                     * avoid recomputing each time
+                                     * the preconditioner is called.
+                                     */
+    vector<pair<unsigned int, unsigned int> > intervals;
+};
 
 
 
index 373beb5d2085f15b462e9c229d8c4c3c3da7b710..598c2092c5cc63f50bb788ab5857b127d1d84f51 100644 (file)
@@ -22,14 +22,14 @@ SparseVanka<number>::SparseVanka(const SparseMatrix<number> &M,
                                 const unsigned int          n_threads)
                :
                matrix (&M),
-               selected (selected),
                conserve_mem (conserve_mem),
+               selected (selected),
                n_threads (n_threads),
                inverses (M.m(), 0)
 {
-  Assert (M.m() == M.n(),
-         ExcMatrixNotSquare ());
-
+  Assert (M.m() == M.n(), ExcMatrixNotSquare ());
+  Assert (M.m() == selected.size(), ExcInvalidVectorSize(M.m(), selected.size()));
+  
   if (conserve_mem == false)
     compute_inverses ();
 }
@@ -222,6 +222,7 @@ SparseVanka<number>::operator ()(Vector<number2>       &dst,
 
 
 
+
 template<typename number>
 template<typename number2>
 void
@@ -230,7 +231,11 @@ SparseVanka<number>::apply_preconditioner (Vector<number2>       &dst,
                                           const unsigned int     begin,
                                           const unsigned int     end) const
 {
-  Assert (begin < end, ExcInternalError());
+  Assert (begin < end, ExcInvalidRange(begin, end));
+  Assert (dst.size() == src.size(),
+         ExcInvalidVectorSize(dst.size(), src.size()));
+  Assert (dst.size() == matrix->m(),
+         ExcInvalidVectorSize(dst.size(), src.size()));
   
                                   // first define an alias to the sparsity
                                   // pattern of the matrix, since this
@@ -244,8 +249,7 @@ SparseVanka<number>::apply_preconditioner (Vector<number2>       &dst,
                                   // blocks. this variable is used to
                                   // optimize access to vectors a
                                   // little bit.
-  const bool range_is_restricted = (begin != 0) && (end != matrix->m());
-  
+  const bool range_is_restricted = ((begin != 0) || (end != matrix->m()));
   
                                   // space to be used for local
                                   // systems. allocate as much memory
@@ -262,7 +266,7 @@ SparseVanka<number>::apply_preconditioner (Vector<number2>       &dst,
 
                                   // traverse all rows of the matrix
                                   // which are selected
-  for (unsigned int row=0; row< matrix->m() ; ++row)
+  for (unsigned int row=begin; row<end; ++row)
     if (selected[row] == true)
       {
        const unsigned int row_length = structure.row_length(row);
@@ -280,7 +284,6 @@ SparseVanka<number>::apply_preconditioner (Vector<number2>       &dst,
        
        b.reinit (row_length);
        x.reinit (row_length);
-       
                                         // mapping between:
                                         // 1 column number of all
                                         //   entries in this row, and
@@ -315,11 +318,7 @@ SparseVanka<number>::apply_preconditioner (Vector<number2>       &dst,
            const unsigned int irow_length = structure.row_length(irow);
            
                                             // copy rhs
-           if (!range_is_restricted ||
-               ((begin <= irow) && (irow < end)))
-             b(i) = src(irow);
-           else
-             b(i) = 0;
+           b(i) = src(irow);
            
                                             // for all the DoFs that irow
                                             // couples with
@@ -383,3 +382,108 @@ SparseVanka<number>::apply_preconditioner (Vector<number2>       &dst,
       };
 };
 
+
+
+
+template <typename number>
+SparseBlockVanka<number>::SparseBlockVanka (const SparseMatrix<number> &M,
+                                           const vector<bool>         &selected,
+                                           const bool                  conserve_memory,
+                                           const unsigned int          n_threads,
+                                           const unsigned int          n_blocks)
+               :
+               SparseVanka<number> (M, selected, conserve_memory, n_threads),
+                n_blocks (n_blocks)
+{
+  Assert (n_blocks > 0, ExcInternalError());
+
+                                  // precompute the splitting points
+  intervals.resize (n_blocks);
+  
+  const unsigned int n_inverses = count (selected.begin(),
+                                        selected.end(),
+                                        true);
+  
+  const unsigned int n_inverses_per_block = max(n_inverses / n_blocks,
+                                               1U);
+  
+                                  // set up start and end index for
+                                  // each of the blocks. note that
+                                  // we have to work somewhat to get
+                                  // this appropriate, since the
+                                  // indices for which inverses have
+                                  // to be computed may not be evenly
+                                  // distributed in the vector. as an
+                                  // extreme example consider
+                                  // numbering of DoFs by component,
+                                  // then all indices for which we
+                                  // have to do work will be
+                                  // consecutive, with other
+                                  // consecutive regions where we do
+                                  // not have to do something
+  unsigned int c       = 0;
+  unsigned int block   = 0;
+  intervals[0].first   = 0;
+  
+  for (unsigned int i=0; (i<M.m()) && (block+1<n_blocks); ++i)
+    {
+      if (selected[i] == true)
+       ++c;
+      if (c == n_inverses_per_block)
+       {
+         intervals[block].second  = i;
+         intervals[block+1].first = i;
+         ++block;
+         
+         c = 0;
+       };
+    };
+  intervals[n_blocks-1].second = M.m();
+};
+
+
+
+
+template <typename number>
+template <typename number2>
+void SparseBlockVanka<number>::operator() (Vector<number2>       &dst,
+                                          const Vector<number2> &src) const
+{
+  dst.clear ();
+
+                                  // if no blocking is required, pass
+                                  // down to the underlying class
+  if (n_blocks == 1)
+    apply_preconditioner (dst, src, 0, dst.size());
+  else
+                                    // otherwise: blocking requested
+    {
+#ifdef DEAL_II_USE_MT
+      typedef ThreadManager::Mem_Fun_Data4
+       <const SparseVanka<number>, Vector<number2>&,
+       const Vector<number2> &, unsigned int, unsigned int> MemFunData;
+      vector<MemFunData> mem_fun_data
+       (n_blocks,
+        MemFunData (this,
+                    dst, src, 0, 0,
+                    &SparseVanka<number>::template apply_preconditioner<number2>));
+
+      ThreadManager thread_manager;
+      for (unsigned int block=0; block<n_blocks; ++block)
+       {
+         mem_fun_data[block].arg3 = intervals[block].first;
+         mem_fun_data[block].arg4 = intervals[block].second;
+
+         thread_manager.spawn (&mem_fun_data[block],
+                               THR_SCOPE_SYSTEM | THR_DETACHED);
+       };
+
+      thread_manager.wait ();
+#else
+      for (unsigned int block=0; block<n_blocks; ++block)
+       apply_preconditioner (dst, src,
+                             intervals[block].first,
+                             intervals[block].second);
+#endif
+    };
+};

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.