From 2174dffb1a7a252077b5930a045bb621cc00f120 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 1 Feb 2000 14:22:06 +0000 Subject: [PATCH] Implement block version of Vanka preconditioner. git-svn-id: https://svn.dealii.org/trunk@2305 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_vanka.h | 215 ++++++++++++++---- .../lac/include/lac/sparse_vanka.templates.h | 132 +++++++++-- 2 files changed, 287 insertions(+), 60 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_vanka.h b/deal.II/lac/include/lac/sparse_vanka.h index 0bbda58ab5..7da845f9e5 100644 --- a/deal.II/lac/include/lac/sparse_vanka.h +++ b/deal.II/lac/include/lac/sparse_vanka.h @@ -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 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 + void apply_preconditioner (Vector &dst, + const Vector &src, + const unsigned int begin, + const unsigned int end) const; private: /** @@ -189,17 +245,17 @@ class SparseVanka */ SmartPointer > matrix; - /** - * Indices of Lagrange - * multipliers. - */ - const vector &selected; - /** * Conserve memory flag. */ const bool conserve_mem; + /** + * Indices of those degrees of + * freedom that we shall work on. + */ + const vector &selected; + /** * Number of threads to be used * when building the @@ -253,49 +309,116 @@ class SparseVanka void compute_inverse (const unsigned int row, map &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 +class SparseBlockVanka : public SparseVanka +{ + 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 - void apply_preconditioner (Vector &dst, - const Vector &src, - const unsigned int begin, - const unsigned int end) const; -}; + SparseBlockVanka (const SparseMatrix &M, + const vector &selected, + const bool conserve_memory = false, + const unsigned int n_threads = 1, + const unsigned int n_blocks = 1); + /** + * Apply the preconditioner. + */ + template + void operator() (Vector &dst, + const Vector &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 > intervals; +}; diff --git a/deal.II/lac/include/lac/sparse_vanka.templates.h b/deal.II/lac/include/lac/sparse_vanka.templates.h index 373beb5d20..598c2092c5 100644 --- a/deal.II/lac/include/lac/sparse_vanka.templates.h +++ b/deal.II/lac/include/lac/sparse_vanka.templates.h @@ -22,14 +22,14 @@ SparseVanka::SparseVanka(const SparseMatrix &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::operator ()(Vector &dst, + template template void @@ -230,7 +231,11 @@ SparseVanka::apply_preconditioner (Vector &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::apply_preconditioner (Vector &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::apply_preconditioner (Vector &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::apply_preconditioner (Vector &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::apply_preconditioner (Vector &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::apply_preconditioner (Vector &dst, }; }; + + + +template +SparseBlockVanka::SparseBlockVanka (const SparseMatrix &M, + const vector &selected, + const bool conserve_memory, + const unsigned int n_threads, + const unsigned int n_blocks) + : + SparseVanka (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 +template +void SparseBlockVanka::operator() (Vector &dst, + const Vector &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 + , Vector&, + const Vector &, unsigned int, unsigned int> MemFunData; + vector mem_fun_data + (n_blocks, + MemFunData (this, + dst, src, 0, 0, + &SparseVanka::template apply_preconditioner)); + + ThreadManager thread_manager; + for (unsigned int block=0; block