From f7cd843259797f354d071264d614af620473d75b Mon Sep 17 00:00:00 2001 From: kanschat Date: Tue, 25 Sep 2007 20:21:19 +0000 Subject: [PATCH] new FilteredMatrix git-svn-id: https://svn.dealii.org/trunk@15238 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/grid/grid_generator.cc | 11 +- deal.II/doc/news/changes.h | 10 + deal.II/lac/include/lac/filtered_matrix.h | 595 ++++++++++++------ .../include/lac/filtered_matrix.templates.h | 272 -------- deal.II/lac/source/filtered_matrix.cc | 37 +- 5 files changed, 412 insertions(+), 513 deletions(-) delete mode 100644 deal.II/lac/include/lac/filtered_matrix.templates.h diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index aaf94e5372..70adfd55c7 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -2287,17 +2287,20 @@ void GridGenerator::laplace_solve (const SparseMatrix &S, Vector &u) { const unsigned int n_dofs=S.n(); - FilteredMatrix > SF (S); + FilteredMatrix > SF (S); + PreconditionJacobi > prec; + prec.initialize(S, 1.2); + FilteredMatrix > PF (prec); + SolverControl control (1000, 1.e-10, false, false); PrimitiveVectorMemory > mem; SolverCG > solver (control, mem); - PreconditionJacobi > > prec; + Vector f(n_dofs); SF.add_constraints(m); - prec.initialize (SF); SF.apply_constraints (f, true); - solver.solve(SF, u, f, prec); + solver.solve(SF, u, f, PF); } diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index e08bb5d0b3..7ebb011822 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -30,6 +30,16 @@ inconvenience this causes.

    + +
  1. Changed: FilteredMatrix now can be applied to any matrix having the standard +set of vmult functions. In order to achieve this, its interface had to be overhauled. +Only the VECTOR template argument remains. Furthermore, instead of +PreconditionJacobi being applied to FilteredMatrix, FilteredMatrix +can now be applied to any preconditioner. +
    +(GK 2007/09/25) +

    +
  2. Changed: The deprecated typedefs internal::Triangulation::Line, internal::Triangulation::Quad, and diff --git a/deal.II/lac/include/lac/filtered_matrix.h b/deal.II/lac/include/lac/filtered_matrix.h index d9c3394f70..cbefe2b29f 100644 --- a/deal.II/lac/include/lac/filtered_matrix.h +++ b/deal.II/lac/include/lac/filtered_matrix.h @@ -18,6 +18,9 @@ #include #include #include +#include +#include +#include #include #include @@ -35,7 +38,82 @@ template class Vector; * This class is a wrapper for linear systems of equations with simple * equality constraints fixing individual degrees of freedom to a * certain value such as when using Dirichlet boundary - * values. Mathematically speaking, it is used to represent a system + * values. + * + * In order to accomplish this, the vmult(), Tvmult(), vmult_add() and + * Tvmult_add functions modify the same function of the original + * matrix such as if all constrained entries of the source vector were + * zero. Additionally, all constrained entries of the destination + * vector are set to zero. + * + *

    Usage

    + * + * Usage is simple: create an object of this type, point it to a + * matrix that shall be used for $A$ above (either through the + * constructor, the copy constructor, or the + * set_referenced_matrix() function), specify the list of boundary + * values or other constraints (through the add_constraints() + * function), and then for each required solution modify the right + * hand side vector (through apply_constraints()) and use this + * object as matrix object in a linear solver. As linear solvers + * should only use vmult() and residual() functions of a + * matrix class, this class should be as a good a matrix as any other + * for that purpose. + * + * Furthermore, also the precondition_Jacobi() function is + * provided (since the computation of diagonal elements of the + * filtered matrix $A_X$ is simple), so you can use this as a + * preconditioner. Some other function useful for matrices are also + * available. + * + * A typical code snippet showing the above steps is as follows: + * @verbatim + * ... // set up sparse matrix A and right hand side b somehow + * + * // initialize filtered matrix with + * // matrix and boundary value constraints + * FilteredMatrix > filtered_A (A); + * filtered_A.add_constraints (boundary_values); + * + * // set up a linear solver + * SolverControl control (1000, 1.e-10, false, false); + * PrimitiveVectorMemory > mem; + * SolverCG > solver (control, mem); + * + * // set up a preconditioner object + * PreconditionJacobi > > prec; + * prec.initialize (filtered_A, 1.2); + * + * // compute modification of right hand side + * filtered_A.apply_constraints (b, true); + * + * // solve for solution vector x + * solver.solve (filtered_A, x, b, prec); + * @endverbatim + * + *

    Connection to other classes

    + * + * The function MatrixTools::apply_boundary_values() does exactly + * the same that this class does, except for the fact that that + * function actually modifies the matrix. Due to this, it is only + * possible to solve with a matrix onto which + * MatrixTools::apply_boundary_values() was applied for one right + * hand side and one set of boundary values since the modification of + * the right hand side depends on the original matrix. + * + * While this is fine (and the recommended way) in cases where only + * one solution of the linear system is required, for example in + * solving linear stationary systems, one would often like to have the + * ability to solve multiply with the same matrix in nonlinear + * problems (where one often does not want to update the Hessian + * between Newton steps, despite having different right hand sides in + * subsequent steps) or time dependent problems, without having to + * re-assemble the matrix or copy it to temporary matrices with which + * one then can work. For these cases, this class is meant. + * + * + *

    Some background

    + * Mathematically speaking, it is used to represent a system * of linear equations $Ax=b$ with the constraint that $B_D x = g_D$, * where $B_D$ is a rectangular matrix with exactly one $1$ in each * row, and these $1$s in those columns representing constrained @@ -100,72 +178,6 @@ template class Vector; * hand side, through the apply_constraints() function. * * - *

    Connection to other classes

    - * - * The function MatrixTools::apply_boundary_values() does exactly - * the same that this class does, except for the fact that that - * function actually modifies the matrix. Due to this, it is only - * possible to solve with a matrix onto which - * MatrixTools::apply_boundary_values() was applied for one right - * hand side and one set of boundary values since the modification of - * the right hand side depends on the original matrix. - * - * While this is fine (and the recommended way) in cases where only - * one solution of the linear system is required, for example in - * solving linear stationary systems, one would often like to have the - * ability to solve multiply with the same matrix in nonlinear - * problems (where one often does not want to update the Hessian - * between Newton steps, despite having different right hand sides in - * subsequent steps) or time dependent problems, without having to - * re-assemble the matrix or copy it to temporary matrices with which - * one then can work. For these cases, this class is meant. - * - * - *

    Usage

    - * - * Usage is simple: create an object of this type, point it to a - * matrix that shall be used for $A$ above (either through the - * constructor, the copy constructor, or the - * set_referenced_matrix() function), specify the list of boundary - * values or other constraints (through the add_constraints() - * function), and then for each required solution modify the right - * hand side vector (through apply_constraints()) and use this - * object as matrix object in a linear solver. As linear solvers - * should only use vmult() and residual() functions of a - * matrix class, this class should be as a good a matrix as any other - * for that purpose. - * - * Furthermore, also the precondition_Jacobi() function is - * provided (since the computation of diagonal elements of the - * filtered matrix $A_X$ is simple), so you can use this as a - * preconditioner. Some other function useful for matrices are also - * available. - * - * A typical code snippet showing the above steps is as follows: - * @verbatim - * ... // set up sparse matrix A and right hand side b somehow - * - * // initialize filtered matrix with - * // matrix and boundary value constraints - * FilteredMatrix > filtered_A (A); - * filtered_A.add_constraints (boundary_values); - * - * // set up a linear solver - * SolverControl control (1000, 1.e-10, false, false); - * PrimitiveVectorMemory > mem; - * SolverCG > solver (control, mem); - * - * // set up a preconditioner object - * PreconditionJacobi > > prec; - * prec.initialize (filtered_A, 1.2); - * - * // compute modification of right hand side - * filtered_A.apply_constraints (b, true); - * - * // solve for solution vector x - * solver.solve (filtered_A, x, b, prec); - * @endverbatim - * * *

    Template arguments

    * @@ -191,26 +203,19 @@ template class Vector; * bottleneck. If you don't want this serialization of operations, you * have to use several objects of this type. * - * @author Wolfgang Bangerth 2001, Luca Heltai 2006 + * @author Wolfgang Bangerth 2001, Luca Heltai 2006, Guido Kanschat 2007 */ -template > +template class FilteredMatrix : public Subscriptor { public: - /** - * Type of matrix entries. In - * analogy to the STL container - * classes. - */ - typedef typename MATRIX::value_type value_type; - /** * Typedef defining a type that * represents a pair of degree of * freedom index and the value it * shall have. */ - typedef std::pair IndexValuePair; + typedef std::pair IndexValuePair; /** * Default constructor. You will @@ -228,12 +233,21 @@ class FilteredMatrix : public Subscriptor */ FilteredMatrix (const FilteredMatrix &fm); - /** + /** * Constructor. Use the given * matrix for future operations. + * + * @arg @p m: The matrix being used in multiplications. + * + * @arg @p + * expect_constrained_source: See + * documentation of + * #expect_constrained_source. */ - FilteredMatrix (const MATRIX &matrix); - + template + FilteredMatrix (const MATRIX &matrix, + bool expect_constrained_source = false); + /** * Copy operator. Take over * matrix and constraints from @@ -248,16 +262,18 @@ class FilteredMatrix : public Subscriptor * clear_constraints() * function if constraits were * previously added. + * + * @arg @p m: The matrix being used in multiplications. + * + * @arg @p + * expect_constrained_source: See + * documentation of + * #expect_constrained_source. */ - void set_referenced_matrix (const MATRIX &m); + template + void initialize (const MATRIX &m, + bool expect_constrained_source = false); - /** - * Return a reference to the - * matrix that is used by this - * object. - */ - const MATRIX & get_referenced_matrix () const; - /** * Add a list of constraints to * the ones already managed by @@ -275,6 +291,13 @@ class FilteredMatrix : public Subscriptor * but also a * std::map. * + * The second component of these + * pairs will only be used in + * apply_constraints(). The first + * is used to set values to zero + * in matrix vector + * multiplications. + * * It is an error if the argument * contains an entry for a degree * of freedom that has already @@ -305,22 +328,6 @@ class FilteredMatrix : public Subscriptor void apply_constraints (VECTOR &v, const bool matrix_is_symmetric) const; - /** - * Return the dimension of the - * image space. To remember: the - * matrix is of dimension - * $m \times n$. - */ - unsigned int m () const; - - /** - * Return the dimension of the - * range space. To remember: the - * matrix is of dimension - * $m \times n$. - */ - unsigned int n () const; - /** * Matrix-vector multiplication: * let $dst = M*src$ with $M$ @@ -352,81 +359,42 @@ class FilteredMatrix : public Subscriptor const VECTOR &src) const; /** - * Return the square of the norm - * of the vector $v$ with respect - * to the norm induced by this - * matrix, - * i.e. $\left(v,Mv\right)$. This - * is useful, e.g. in the finite - * element context, where the - * $L_2$ norm of a function - * equals the matrix norm with - * respect to the mass matrix of - * the vector representing the - * nodal values of the finite - * element function. - * - * Obviously, the matrix needs to - * be square for this operation. - * - * Note that in many cases, you - * will not want to compute the - * norm with respect to the - * filtered matrix, but with - * respect to the original - * one. For example, if you want - * to compute the $L^2$ norm of a - * vector by forming the matrix - * norm with the mass matrix, - * then you want to use the - * original mass matrix, not the - * filtered one where you might - * have eliminated Dirichlet - * boundary values. - */ - value_type matrix_norm_square (const VECTOR &v) const; - - /** - * Compute the residual of an - * equation Mx=b, where the - * residual is defined to be - * r=b-Mx with @p x - * typically being an approximate - * of the true solution of the - * equation. Write the residual - * into @p dst. The l2 norm of - * the residual vector is - * returned. + * Adding matrix-vector multiplication. * - * Note that it is assumed that - * @p b is a vector that has been - * treated by the - * modify_rhs() function, - * since we can then assume that - * the components of the residual - * which correspond to - * constrained degrees of freedom - * do not contribute to the - * residual at all. + * @note The result vector of + * this multiplication will have + * the constraint entries set to + * zero, independent of the + * previous value of + * dst. We excpect that + * in most cases this is the + * required behavior. */ - value_type residual (VECTOR &dst, - const VECTOR &x, - const VECTOR &b) const; + void vmult_add (VECTOR &dst, + const VECTOR &src) const; /** - * Apply the Jacobi - * preconditioner, which - * multiplies every element of - * the @p src vector by the - * inverse of the respective - * diagonal element and - * multiplies the result with the - * damping factor @p omega. + * Adding transpose matrix-vector multiplication: + * + * Because we need to use a + * temporary variable and since + * we only allocate that each + * time the matrix changed, this + * function only works for + * quadratic matrices. + * + * @note The result vector of + * this multiplication will have + * the constraint entries set to + * zero, independent of the + * previous value of + * dst. We excpect that + * in most cases this is the + * required behavior. */ - void precondition_Jacobi (VECTOR &dst, - const VECTOR &src, - const value_type omega = 1.) const; - + void Tvmult_add (VECTOR &dst, + const VECTOR &src) const; + /** * Determine an estimate for the * memory consumption (in bytes) @@ -438,6 +406,27 @@ class FilteredMatrix : public Subscriptor unsigned int memory_consumption () const; private: + /** + * Determine, whether + * multiplications can expect + * that the source vector has all + * constrained entries set to + * zero. + * + * If so, the auxiliary vector + * can be avoided and memory as + * well as time can be saved. + * + * We expect this for instance in + * Newton's method, where the + * residual already should be + * zero on constrained + * nodes. This is, because there + * is no testfunction in these + * nodes. + */ + bool expect_constrained_source; + /** * Declare an abbreviation for an * iterator into the array @@ -474,7 +463,7 @@ class FilteredMatrix : public Subscriptor * it using the SmartPointer * class. */ - SmartPointer matrix; + boost::shared_ptr > matrix; /** * Sorted list of pairs denoting @@ -525,32 +514,16 @@ class FilteredMatrix : public Subscriptor */ void post_filter (const VECTOR &in, VECTOR &out) const; - - /** - * Based on the size of the - * matrix and type of the matrix - * and vector, allocate a - * temporary vector. This - * function has to be overloaded - * for the various template - * parameter choices. Since the - * allocated vector will be - * filled by the site that calls - * this function, no - * initialization is necessary. - */ - void allocate_tmp_vector (); - }; /*@}*/ /*---------------------- Inline functions -----------------------------------*/ -template +template inline bool -FilteredMatrix::PairComparison:: +FilteredMatrix::PairComparison:: operator () (const IndexValuePair &i1, const IndexValuePair &i2) const { @@ -559,10 +532,66 @@ operator () (const IndexValuePair &i1, -template +template +template +inline +void +FilteredMatrix::initialize (const MATRIX &m, bool ecs) +{ + matrix = boost::shared_ptr > ( + new_pointer_matrix_base(m, VECTOR())); + + expect_constrained_source = ecs; +} + + + +template +FilteredMatrix::FilteredMatrix () +{} + + + +template +FilteredMatrix:: +FilteredMatrix (const FilteredMatrix &fm) + : + Subscriptor(), + constraints (fm.constraints) +{ + initialize (*fm.matrix, fm.expect_constrained_source); +} + + + +template +template +inline +FilteredMatrix:: +FilteredMatrix (const MATRIX &m, bool ecs) +{ + initialize (m, ecs); +} + + + +template +inline +FilteredMatrix & +FilteredMatrix::operator = (const FilteredMatrix &fm) +{ + matrix = fm.matrix; + expect_constrained_source = fm.expect_constrained_source; + constraints = fm.constraints; + return *this; +} + + + +template template void -FilteredMatrix:: +FilteredMatrix:: add_constraints (const ConstraintList &new_constraints) { // add new constraints to end @@ -581,34 +610,190 @@ add_constraints (const ConstraintList &new_constraints) -template +template +inline +void +FilteredMatrix::clear_constraints () +{ + // swap vectors to release memory + std::vector empty; + constraints.swap (empty); +} + + + +template +inline +void +FilteredMatrix:: +apply_constraints (VECTOR &v, + const bool /* matrix_is_symmetric */) const +{ + tmp_vector.reinit(v); + const_index_value_iterator i = constraints.begin(); + const const_index_value_iterator e = constraints.end(); + for (; i!=e; ++i) + tmp_vector(i->first) = -i->second; + + // This vmult is without bc, to get the rhs correction in a correct way. + matrix->vmult_add(v, tmp_vector); + + // finally set constrained entries themselves + for (i=constraints.begin(); i!=e; ++i) + v(i->first) = i->second; +} + + + +template inline -const MATRIX & -FilteredMatrix::get_referenced_matrix () const +void +FilteredMatrix::pre_filter (VECTOR &v) const { - return *matrix; + // iterate over all constraints and + // zero out value + const_index_value_iterator i = constraints.begin(); + const const_index_value_iterator e = constraints.end(); + for (; i!=e; ++i) + v(i->first) = 0; } -template +template inline -unsigned int FilteredMatrix::m () const +void +FilteredMatrix::post_filter (const VECTOR &in, + VECTOR &out) const { - return matrix->m(); + // iterate over all constraints and + // set value correctly + const_index_value_iterator i = constraints.begin(); + const const_index_value_iterator e = constraints.end(); + for (; i!=e; ++i) + out(i->first) = in(i->first); } -template +template inline -unsigned int FilteredMatrix::n () const +void +FilteredMatrix::vmult (VECTOR& dst, const VECTOR& src) const { - return matrix->n(); + if (!expect_constrained_source) + { + tmp_mutex.acquire (); + // first copy over src vector and + // pre-filter + tmp_vector.reinit(src, true); + tmp_vector = src; + pre_filter (tmp_vector); + // then let matrix do its work + matrix->vmult (dst, tmp_vector); + // tmp_vector now no more needed + tmp_mutex.release (); + } + else + matrix->vmult (dst, src); + // finally do post-filtering + post_filter (src, dst); } +template +inline +void +FilteredMatrix::Tvmult (VECTOR& dst, const VECTOR& src) const +{ + if (!expect_constrained_source) + { + tmp_mutex.acquire (); + // first copy over src vector and + // pre-filter + tmp_vector.reinit(src, true); + tmp_vector = src; + pre_filter (tmp_vector); + // then let matrix do its work + matrix->Tvmult (dst, tmp_vector); + // tmp_vector now no more needed + tmp_mutex.release (); + } + else + matrix->Tvmult (dst, src); + // finally do post-filtering + post_filter (src, dst); +} + + + +template +inline +void +FilteredMatrix::vmult_add (VECTOR& dst, const VECTOR& src) const +{ + if (!expect_constrained_source) + { + tmp_mutex.acquire (); + // first copy over src vector and + // pre-filter + tmp_vector.reinit(src, true); + tmp_vector = src; + pre_filter (tmp_vector); + // then let matrix do its work + matrix->vmult_add (dst, tmp_vector); + // tmp_vector now no more needed + tmp_mutex.release (); + } + else + matrix->vmult_add (dst, src); + // finally do post-filtering + post_filter (src, dst); +} + + + +template +inline +void +FilteredMatrix::Tvmult_add (VECTOR& dst, const VECTOR& src) const +{ + if (!expect_constrained_source) + { + tmp_mutex.acquire (); + // first copy over src vector and + // pre-filter + tmp_vector.reinit(src, true); + tmp_vector = src; + pre_filter (tmp_vector); + // then let matrix do its work + matrix->Tvmult_add (dst, tmp_vector); + // tmp_vector now no more needed + tmp_mutex.release (); + } + else + matrix->Tvmult_add (dst, src); + // finally do post-filtering + post_filter (src, dst); +} + + + +template +inline +unsigned int +FilteredMatrix::memory_consumption () const +{ + return (MemoryConsumption::memory_consumption (matrix) + + MemoryConsumption::memory_consumption (constraints) + + MemoryConsumption::memory_consumption (tmp_vector)); +} + + + + + /*---------------------------- filtered_matrix.h ---------------------------*/ DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/lac/include/lac/filtered_matrix.templates.h b/deal.II/lac/include/lac/filtered_matrix.templates.h deleted file mode 100644 index 6b9f20ebdc..0000000000 --- a/deal.II/lac/include/lac/filtered_matrix.templates.h +++ /dev/null @@ -1,272 +0,0 @@ -//--------------------------------------------------------------------------- -// $Id$ -// Version: $Name$ -// -// Copyright (C) 2001, 2002, 2003, 2005, 2006 by the deal.II authors -// -// This file is subject to QPL and may not be distributed -// without copyright and license information. Please refer -// to the file deal.II/doc/license.html for the text and -// further information on this license. -// -//--------------------------------------------------------------------------- -#ifndef __deal2__filtered_matrix_templates_h -#define __deal2__filtered_matrix_templates_h - - -#include -#include -#include -#include -#include -#include -#include - -DEAL_II_NAMESPACE_OPEN - - -template -FilteredMatrix::FilteredMatrix () -{} - -template -FilteredMatrix:: -FilteredMatrix (const FilteredMatrix &fm) - : - Subscriptor (), - constraints (fm.constraints) -{ - set_referenced_matrix (*fm.matrix); -} - - - -template -FilteredMatrix:: -FilteredMatrix (const MATRIX &m) -{ - set_referenced_matrix (m); -} - - - -template -FilteredMatrix & -FilteredMatrix::operator = (const FilteredMatrix &fm) -{ - set_referenced_matrix (*fm.matrix); - constraints = fm.constraints; - return *this; -} - - - -template -void -FilteredMatrix:: -set_referenced_matrix (const MATRIX &m) -{ - matrix = &m; - allocate_tmp_vector(); -} - - - -template -void -FilteredMatrix::clear_constraints () -{ - // swap vectors to release memory - std::vector empty; - constraints.swap (empty); -} - - - -template -void -FilteredMatrix:: -apply_constraints (VECTOR &v, - const bool /* matrix_is_symmetric */) const -{ - tmp_vector = 0; - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) - tmp_vector(i->first) = -i->second; - - // This vmult is without bc, to get the rhs correction in a correct way. - matrix->vmult_add(v, tmp_vector); - - // finally set constrained entries themselves - for (i=constraints.begin(); i!=e; ++i) - v(i->first) = i->second; -} - - - -template -void -FilteredMatrix::pre_filter (VECTOR &v) const -{ - // iterate over all constraints and - // zero out value - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) - v(i->first) = 0; -} - - - -template -void -FilteredMatrix::post_filter (const VECTOR &in, - VECTOR &out) const -{ - // iterate over all constraints and - // set value correctly - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) - out(i->first) = in(i->first); -} - - - -template -void -FilteredMatrix::vmult (VECTOR &dst, - const VECTOR &src) const -{ - tmp_mutex.acquire (); - // first copy over src vector and - // pre-filter - tmp_vector = src; - pre_filter (tmp_vector); - // then let matrix do its work - matrix->vmult (dst, tmp_vector); - // tmp_vector now no more needed - tmp_mutex.release (); - // finally do post-filtering - post_filter (src, dst); -} - - - -template -typename FilteredMatrix::value_type -FilteredMatrix::residual (VECTOR &dst, - const VECTOR &x, - const VECTOR &b) const -{ - tmp_mutex.acquire (); - // first copy over x vector and - // pre-filter - tmp_vector = x; - pre_filter (tmp_vector); - // then let matrix do its work - value_type res = matrix->residual (dst, tmp_vector, b); - value_type res2 = res*res; - // tmp_vector now no more needed - tmp_mutex.release (); - // finally do post-filtering. here, - // we set constrained indices to - // zero, but have to subtract their - // contributions to the residual - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) - { - const value_type v = dst(i->first); - res2 -= v*v; - dst(i->first) = 0; - }; - - Assert (res2>=0, ExcInternalError()); - return std::sqrt (res2); -} - - - -template -void -FilteredMatrix::Tvmult (VECTOR &dst, - const VECTOR &src) const -{ - tmp_mutex.acquire (); - // first copy over src vector and - // pre-filter - tmp_vector = src; - pre_filter (tmp_vector); - // then let matrix do its work - matrix->Tvmult (dst, tmp_vector); - // tmp_vector now no more needed - tmp_mutex.release (); - // finally do post-filtering - post_filter (src, dst); -} - - - -template -typename FilteredMatrix::value_type -FilteredMatrix::matrix_norm_square (const VECTOR &v) const -{ - tmp_mutex.acquire (); - tmp_vector = v; - - // zero out constrained entries and - // form matrix norm with original - // matrix. this is equivalent to - // forming the matrix norm of the - // original vector with the matrix - // where we have zeroed out rows - // and columns - pre_filter (tmp_vector); - const value_type ret = matrix->matrix_norm_square (tmp_vector); - tmp_mutex.release (); - return ret; -} - - - -template -void -FilteredMatrix:: -precondition_Jacobi (VECTOR &dst, - const VECTOR &src, - const value_type omega) const -{ - // first precondition as usual, - // using the fast algorithms of the - // matrix class - matrix->precondition_Jacobi (dst, src, omega); - - // then modify the constrained - // degree of freedom. as the - // diagonal entries of the filtered - // matrix would be 1.0, simply copy - // over old and new values - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) - dst(i->first) = src(i->first); -} - - - -template -unsigned int -FilteredMatrix::memory_consumption () const -{ - return (MemoryConsumption::memory_consumption (matrix) + - MemoryConsumption::memory_consumption (constraints) + - MemoryConsumption::memory_consumption (tmp_vector)); -} - - - -DEAL_II_NAMESPACE_CLOSE - -#endif diff --git a/deal.II/lac/source/filtered_matrix.cc b/deal.II/lac/source/filtered_matrix.cc index e97391c828..a216440c2f 100644 --- a/deal.II/lac/source/filtered_matrix.cc +++ b/deal.II/lac/source/filtered_matrix.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2006 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -16,36 +16,9 @@ DEAL_II_NAMESPACE_OPEN -#define FILT(MM,VV) \ -template <> \ -void FilteredMatrix::allocate_tmp_vector() \ -{\ - Threads::ThreadMutex::ScopedLock lock (tmp_mutex); \ - tmp_vector.reinit (matrix->n(), true); \ -} - -#define BFILT(MM,VV) \ -template <> \ -void \ -FilteredMatrix:: \ -allocate_tmp_vector () \ -{ \ - std::vector block_sizes (matrix->n_block_rows()); \ - for (unsigned int i=0; iblock(i,i).n(); \ - \ - Threads::ThreadMutex::ScopedLock lock (tmp_mutex); \ - tmp_vector.reinit (block_sizes, true); \ -} - -FILT(SparseMatrix, Vector) -BFILT(BlockSparseMatrix, BlockVector) -template class FilteredMatrix,Vector >; -template class FilteredMatrix,BlockVector >; - -FILT(SparseMatrix, Vector) -BFILT(BlockSparseMatrix, BlockVector) -template class FilteredMatrix,Vector >; -template class FilteredMatrix,BlockVector >; +template class FilteredMatrix >; +template class FilteredMatrix >; +template class FilteredMatrix >; +template class FilteredMatrix >; DEAL_II_NAMESPACE_CLOSE -- 2.39.5