From db9c30bab5e63fc463b92007d92c6e2afac656ea Mon Sep 17 00:00:00 2001 From: wolf Date: Sat, 17 Jul 1999 21:57:49 +0000 Subject: [PATCH] Mostly doc updates. git-svn-id: https://svn.dealii.org/trunk@1587 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_matrix.h | 10 +- deal.II/lac/include/lac/sparse_vanka.h | 98 ++++++++++--- .../lac/include/lac/sparse_vanka.templates.h | 131 ++++++++++++------ 3 files changed, 180 insertions(+), 59 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index 9056e40c8a..ea46c58104 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -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 #include @@ -13,26 +16,75 @@ /** * 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 signature(3); + * signature[0] = signature[1] = false; + * signature[2] = true; + * + * // tag all dofs belonging to the + * // Lagrange multiplier + * vector selected_dofs (dof.n_dofs(), false); + * DoFTools::extract_dofs(dof, signature, p_select); + * // create the Vanka object + * SparseVanka vanka (global_matrix, selected_dofs); + * + * // create the solver + * SolverGMRES, + * Vector > gmres(control,memory,504); + * + * // solve + * gmres.solve (global_matrix, solution, right_hand_side, + * vanka); + * \end{verbatim} + * + * @author Guido Kanschat, 1999 + */ template 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 > > inverses; }; + + +/* ---------------------------- Inline functions -----------------------*/ + template template inline @@ -131,4 +192,9 @@ SparseVanka::operator() (Vector& dst, } + + +/*---------------------------- sparse_vanka.h ---------------------------*/ +/* end of #ifndef __sparse_vanka_H */ #endif +/*---------------------------- sparse_vanka.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/sparse_vanka.templates.h b/deal.II/lac/include/lac/sparse_vanka.templates.h index 3f100bcd72..dd24436823 100644 --- a/deal.II/lac/include/lac/sparse_vanka.templates.h +++ b/deal.II/lac/include/lac/sparse_vanka.templates.h @@ -6,12 +6,21 @@ #include + + template -SparseVanka::SparseVanka(const SparseMatrix& M, - const bit_vector& selected) +SparseVanka::SparseVanka(const SparseMatrix &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 SparseVanka::~SparseVanka() @@ -33,68 +42,112 @@ void SparseVanka::forward(Vector& dst, const Vector& 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 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* p = new FullMatrix; + FullMatrix *p = new FullMatrix; inverses[row] = p; } else { + // inverse already built build_matrix = false; } } - FullMatrix& A = (conserve_mem) ? local_matrix : (*inverses[row]); + // alias to the matrix which is + // to be used + FullMatrix &A = (conserve_mem ? + local_matrix : + (*inverses[row])); - if (build_matrix) A.reinit(n); - - Vector b(n); - Vector x(n); + if (build_matrix) + A.reinit(row_length); - map local_index; + Vector b(row_length); + Vector x(row_length); - // Build local index - - for (unsigned int i=0;i local_index; + for (unsigned int i=0; i (structure.column_number(row, i), i)); - -// for (map::iterator is=local_index.begin(); -// is!=local_index.end();++is) -// cerr << "map " << is->first << '\t' << is->second << endl; - // Build local matrix - - for (map::iterator is=local_index.begin(); - is!=local_index.end();++is) + // Build local matrix and rhs + for (map::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::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::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::forward(Vector& dst, A.vmult(x,b); // Distribute new values - for (map::iterator is=local_index.begin(); - is!=local_index.end();++is) + for (map::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); } } -- 2.39.5