/**
* 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
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];
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)));
-// $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
{
* are built up and inverted every time.
*/
void conserve_memory();
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcMatrixNotSquare);
private:
/**
mutable vector<SmartPointer<FullMatrix<float> > > inverses;
};
+
+
+/* ---------------------------- Inline functions -----------------------*/
+
template<typename number>
template<typename number2>
inline
}
+
+
+/*---------------------------- sparse_vanka.h ---------------------------*/
+/* end of #ifndef __sparse_vanka_H */
#endif
+/*---------------------------- sparse_vanka.h ---------------------------*/
#include <map>
+
+
template<typename number>
-SparseVanka<number>::SparseVanka(const SparseMatrix<number>& M,
- const bit_vector& selected)
+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()
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);
}
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);
}
}